30125 CE Computational Methods: Assignment 5
December 8, 2010
Problem 1
A fourth order accurate central approximation to the second derivative (p=2, n=4) requires p + n - 1 =
2 + 4 - 1 = 6 - 1 = 5 nodes. The second central derivative can therefore be approximate to O(h4 ) as:
(2)
fi
E =
1 fi2 + 2 fi1 + 3 fi + 4 fi+1 + 5 fi+2
h2
T.S. expansions about xi are:
(1)
fi2 = fi 2hfi
(2)
+ 2h2 fi
(1)
+
(1)
+
fi1 = fi hfi
fi+1 = fi + hfi
(1)
fi+2 = fi + 2hfi
4h3 (3) 2h4 (4) 4h5 (5) 4h6 (6)
f+
f
f+
f + O(h7 )
3i
3i
15 i
45 i
h2 (2) h3 (3) h4 (4)
h5 (5)
h6 (6)
fi fi + fi
fi +
f + O(h7 )
2
6
24
120
720 i
fi = fi
h2 (2) h3 (3) h4 (4)
h5 (5)
h6 (6)
fi + fi + fi +
fi +
f + O(h7 )
2
6
24
120
720 i
(2)
+ 2h2 fi
+
4h3 (3) 2h4 (4) 4h5 (5) 4h6 (6)
f+
f+
f+
f + O(h7 )
3i
3i
15 i
45 i
Substituting into our assumed form of fi (2) and re-arranging:
(2)
fi
E =
1 fi2 + 2 fi1 3 fi + 4 fi+1 + 5 fi+2
h2
(1 + 2 + 3 + 4 + 5 )
(21 2 + 4 + 25 ) (1)
2
4
(2)
fi +
fi + (21 +
+
+ 25 )fi
2
h
h
2
2
41 2 4 45
21 2 4 25 2 (4)
41
2
4
45 3 (5)
(3)
+(
+
+
)hfi + (
+
+
+
)h fi + (
+
+
)h fi
3
6
6
3
3
24
24
3
15
120 120
15
41
2
4
45 4 (6)
+(
+
+
+
) h fi
45
720 720
45
=
We want fi (2) and 4th order accuracy, so we note that the coecients must satisfy the following
equalities:
1 + 2 + 3 + 4 + 5 = 0
21 2 + 4 + 25 = 0
2
4
21 +
+
+ 25 = 1
2
2
41
2
4
45
+
+
=0
3
6
6
3
21
2
4
25
+
+
+
=0
3
24
24
3
1
We can then write this system of equations in matrix form AX = B
0
1
1
1
1
1
1
2
1
0
1
2 2 0
2
0.5
0
0.5
2 3 = 1
1.33333 0.16667 0 0.16667 1.33333 4 0
0
5
0.66667
0.041667 0 0.041667 0.66667
1
Solving this equation in Matlab (see asgn5 prob1.m) yields 1 = - 12 , 2 = 4 , 3 = - 5 , 4 = 4 , and
3
2
3
1
5 = - 12 .
The equation for our fourth order accurate central approximation to the second derivative now becomes
1 (6)
(1)
(2)
(3)
(4)
(5)
)f + O(h5 )
= (0)fi + (0)fi + (1)fi + (0)fi + (0)fi + (0)fi + (
90 i
The fourth order accurate central approximation to the second derivative is as follows:
fi2 + 16fi1 30fi + 16fi+1 fi+2
+E
12h2
where
1 4 (6)
E
h fi
=
90
2
Problem 2
For the third derivative p = 3 and we must use our N+1 = 4 nodes to establish an interpolating
polynomial of degree N 3. We let the four points be x0 = 0, x1 = h, x2 = 2h, and x3 = 3h. Using
Newton forward interpolation, we then establish the following interpolating polynomial:
g ( x) = f0 +
3 f0
2 f0
f0
(x x0 )(x x1 ) +
(x x0 )(x x1 )(x x2 )
(x x0 ) +
2
h
2h
6h3
with the following associated error term:
e(x) =
f (3) ( )
4 f0
(x x0 )(x x1 )(x x2 )(x x3 )
(x x0 )(x x1 )(x x2 )(x x3 )
=
24
24h4
where
f0 = f1 f0
2
f0 = f2 2f1 + f0
3 f0 = f3 3f2 + 3f1 f0
Since we are looking for the third derivatives of the interpolating function and the error term, we
note that these terms dierentiate to zero, and simply focus in on the third and fourth order terms.
Dierentiating these terms yields:
g (1) (x) =
3 f0
3 f0
3 f 0
(x x1 )(x x2 ) +
(x x0 )(x x2 ) +
(x x0 )(x x1 )
3
3
6h
6h
6h3
with the following associated error term:
e(1) (x) =
f (3) ( )
f (3) ( )
(x x1 )(x x2 )(x x3 ) +
(x x0 )(x x2 )(x x3 )
24
24
f (3) ( )
f (3) ( )
(x x0 )(x x1 )(x x3 ) +
(x x0 )(x x1 )(x x2 )
24
24
Dierenting a second time yields:
+
g (2) (x) =
3 f 0
3 f0
3 f 0
(x x1 ) +
( x x2 ) +
(x x0 )
6h3
6h3
6h3
3 f0
3 f 0
3 f0
( x x2 ) +
(x x0 ) +
( x x1 )
3
3
6h
6h
6h3
with the following associated error term:
+
e(2) (x) =
f (3) ( )
f (3) ( )
f (3) ( )
(x x2 )(x x3 ) +
(x x1 )(x x3 ) +
(x x1 )(x x2 )
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
(x x2 )(x x3 ) +
(x x0 )(x x3 ) +
(x x0 )(x x2 )
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
(x x1 )(x x3 ) +
(x x0 )(x x3 ) +
(x x0 )(x x1 )
24
24
24
f (3) ( )
f (3) ( )
f (3) ( )
(x x1 )(x x2 ) +
(x x0 )(x x2 ) +
(x x0 )(x x1 )
24
24
24
Finally, dierentiating a third time and evaluating at x = x3 gives us the backward approximations
to the third derivative and the error associated with the formula.
+
g (3) (x) = 6
3 f0
f3 3f2 + 3f1 f0
3 f0
=
=
3
3
6h
h
h3
3
Dierentiating the error term one more time leads to
e(3) (x) =
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x3 ) +
(x x2 ) +
( x x3 ) +
(x x1 )
24
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x2 ) +
(x x1 ) +
( x x3 ) +
(x x2 )
24
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x0 ) +
(x x3 ) +
( x x2 ) +
(x x0 )
24
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x3 ) +
(x x1 ) +
( x x3 ) +
(x x0 )
24
24
24
24
+
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x1 ) +
(x x0 ) +
( x x2 ) +
(x x1 )
24
24
24
24
f (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x2 ) +
(x x0 ) +
( x x1 ) +
(x x0 )
24
24
24
24
The error term simplies to the following:
+
e(3) =
f (x) (3) ( )
f (3) ( )
f (3) ( )
f (3) ( )
(x x0 ) +
(x x1 ) +
( x x2 ) +
(x x3 )
4
4
4
4
f (3) ( )
[(x x3 ) + (x x2 ) + (x x1 ) + (x x0 )]
4
We then evaluate this at x3 in order to get the error associated with the backward approximation:
=
e(3) (x3 ) =
f (3) ( )
f (3) ( )
3f (3) ( )
[(0) + (h) + (2h) + (3h)] =
[6h] =
h
4
4
2
Finally, we re-index so that the node j corresponds with x3
d3 f
dx3
=
xj
fj 3fj 1 + 3fj 2 fj 3
h3
3 (3)
E fj h
=
2
4
Problem 3
For this problem, we evaluate the approximations for f(4) and veried them using a Python script (see
prob 3.py ).
(a)
3 (fi )
3 (fi+1 fi )
2 (fi+1 fi )
2 (fi+2 2fi+1 + fi ))
(fi+2 2fi+1 + fi ))
4 f i
=
=
=
=
=
4
4
4
4
h
h
h
h
h4
h4
=
(fi+3 fi+2 2fi+2 + 2fi+1 + fi+1 fi )
(fi+3 3fi+2 + 3fi+1 fi )
fi+3 3fi+2 + 3fi+1 fi
=
=
h4
h4
h4
fi+4 fi+3 3fi+3 + 3fi+2 + 3fi+2 3fi+1 fi+1 + fi
fi+4 4fi+3 + 6fi+2 4fi+1 + fi
=
=
h4
h4
(b)
3 f i
3 (fi fi1 )
2 (fi fi1 )
2 (fi+1 2fi + fi1 )
(fi+1 2fi + fi1 )
=
=
=
=
4
4
4
4
h
h
h
h
h4
(fi+2 3fi+1 + 3fi fi1 )
(fi+2 3fi+1 + 3fi fi1 )
=
4
h
h4
fi+3 fi+2 3fi+2 + 3fi+1 + 3fi+1 3fi fi + fi1
fi+3 4fi+2 + 6fi+1 4fi + fi1
=
=
h4
h4
=
(c)
2
2 fi
=
h4
2
(fi+1 fi )
=
h4
2
(fi+1 fi )
h4
( fi+2 2 fi+1 + fi )
(fi+2 2fi+1 + fi )
=
4
h
h4
fi+2 4fi+1 + 6fi 4fi1 + fi2
(fi+2 3fi+1 + 3fi fi1 )
=
4
h
h4
=
=
2
5
Problem 4
For this problem we are working with the following function:
f (x) = 3x5 4x3 + x 6
which has the following derivatives:
f (1) (x) = 15x4 12x2 + 1
f (2) (x) = 60x3 24x
f (3) (x) = 180x2 24
(a) The matlab script computes all the values asked for in the homework assignment. The results are
shown in Table 1.
h
1.0
0.5
0.1
0.05
0.01
f2
60.000000
-6.000000
-5.831040
-5.903970
-5.980032
f1
-6.000000
-5.906250
-5.903970
-5.950499
-5.990004
f0
-6.000000
-6.000000
-6.000000
-6.000000
-6.000000
f0 (1) approx
-33.000000
0.375000
1.075800
1.019739
1.000800
Eactual
34.000000
0.625000
-0.075800
-0.019738
-0.000800
E1
-8.000000
-2.000000
-0.080000
-0.020000
-0.00080
Table 1: Functional values, approximations, and errors for problem 4
(b)
(i)
As can be seen in Figure 1, the estimated error curve has a slope of 2 on the log-log plot, and the
actual error curve generally has a slope of 2. The slope of the estimated error makes sense since taking
logs of both sides and then rearranging yields the following:
1
1 (3)
1 (3)
(3)
log(E1 ) = log( h2 f0 ) = log( f0 ) + log(h2 ) = 2 log(h) + log( f0 )
3
3
3
1 (3)
log(E1 ) = 2 log(h) + log( f0 )
3
Since log(E1 ) is on the y-axis and log(h) is on the x-axis in the plot, this equation is in the form y = mx
+ b. The slope (m) is 2 which is what we see in the plot. Simply put, we expect a second order accurate
formula to result in a straight line with a slope of 2 on the log-log plot.
(ii)
When we have larger h, we see that the two curves may dier. This is because the leading order
truncated error term may not dominate the rest of the Taylor series in this case since the higher order
terms may be non-negligible. We would have to inspect the actual functional values, and their derivatives
in relation to grid spacing, to understand why the actual error is worse than expected as h approaches
1, but is actually better than expected for some range. Regardless, we see that as h becomes suciently
small, the actual error approaches the estimated error since the error term eventually dominates the rest
of the series.
6
Figure 1: Error Curves.
7
Figure 2: Molecule for problem 5 part (a).
Problem 5
(a) For this problem, we use the fourth-order accurate central dierence approximation to the rst
derivative:
1 4 (6)
fi+2 + 16fi+1 30fi + 16fi1 fi2
(2)
+ E, E
h fi
fi =
=
2
12h
90
With this, we have the following partial derivatives with respect to x and y:
fi+2,j + 16fi+1,j 30fi,j + 16fi1,j fi2,j
2f
=
x2
12x2
2f
fi,j +2 + 16fi,j +1 30fi,j + 16fi,j 1 fi,j 2
=
y 2
12y 2
We then consider Laplaces equation and x = y = h
2
f=
2f
2f
+2
x2
y
1
(fi+2,j + 16fi+1,j + 16fi1,j fi2,j fi,j +2 + 16fi,j +1 + 16fi,j 1 fi,j 2 60fi,j )
12h2
This equation can then be represented by the molecule shown in Figure 2.
(b)
By simply adjusting the i and j values in the above equation, we can derive the following two approximations:
2
2
fi,j =
fi1,j =
1
(fi+1,j +16fi,j +16fi2,j fi3,j fi1,j +2 +16fi1,j +1 +16fi1,j 1 fi1,j 2 60fi1,j )
12h2
1
(fi+4,j +1 +16fi+3,j +1 +16fi+1,j +1 fi,j +1 fi+2,j +3 +16fi+2,j +2 +16fi+2,j fi+2,j 1 60fi+2,j +1 )
12h2
These approximations can be represented by the molecules shown in Figure 3. With this, we can see
that no nodal values would be shared in both approximations.
2
fi+2,j +1 =
8
Figure 3: Molecule for problem 5 part (b). Red molecule is for
9
2
fi1,j . Blue molecule is for
2
fi+2,j +1 .