APM3711/201/2/2014
Tutorial Letter 201/2/2014
Numerical Methods II
APM3711
Semester 2
Department of Mathematical Sciences
IMPORTANT INFORMATION:
This tutorial letter contains solutions
to assignment 01.
BAR CODE
university
Learn without limits. of south africa
ONLY FOR SEMESTER 2 STUDENTS
ASSIGNMENT 01
FIXED CLOSING DATE: 15 AUGUST 2014
Unique Number: 897324
Question 1
Solve the differential equation
dy
= 3x + 2y + xy, y(0) = −1
dx
by means of the Taylor–series expansion to get the value of y at x = 0.1. Use terms up to x6 .
Solution:
We approximate y (0.1) by the truncated Taylor series
0 (0.1)2 00 (0.1)6 (6)
y (0.1) ' y (0) + (0.1) y (0) + y (0) + · · · + y (0) (1)
2 6!
The derivatives at x = 0 are calculated in the way illustrated on p. 398 [354] of Gerald :
y 0 (x) = 3x + (2 + x) y, y 0 (0) = 2 (−1) = −2
y 00 (x) = 3 + y + (2 + x) y 0 , y 00 (0) = 3 − 1 + 2 (−2) = −2
y 000 (x) = 2y 0 + (2 + x) y 00 , y 000 (0) = 2 (−2) + 2 (−2) = −8
y (4) (x) = 3y 00 + (2 + x) y 000 , y (4) (0) = 3 (−2) + 2 (−8) = −22
y (5) (x) = 4y 000 + (2 + x) y (4) , y (5) (0) = 4 (−8) + 2 (−22) = −76
y (6) (x) = 5y (4) + (2 + x) y (5) , y (6) (0) = −5 · 22 − 2 · 76 = −262
Substituting these values in (1) we get
y (0.1) ' −1.211431697.
The exact value, rounded to 10 significant numbers, can be calculated by adding more terms of the
Taylor series to (1) until it is clear that further addition of terms will not alter the tenth digit. You
may verify that this gives
y (0.1) = −1.211431726.
2
APM3711/201
Question 2
Consider the system of coupled second-order differential equations
2
u00 − (t + 1) (u0 ) + 2uv − u3 = cos t
2v 00 + (sin t) u0 v 0 − 6u = 2t + 3
with initial conditions
u(0) = 1, u0 (0) = 2, v(0) = 3, v 0 (0) = 4.
Use the second-order Runge-Kutta method with h = 0.2 and a = 2/3, b = 1/3, α = β = 3/2, to find
u, u0 , v and v 0 at t = 0.2.
Solution:
Given
2
u00 − (t + 1) (u0 ) + 2uv − u3 = cos t
2v 00 + (sin t) u0 v 0 − 6u = 2t + 3
u (0) = 1, u0 (0) = 2, v (0) = 3, v 0 (0) = 4
First, we convert the system of second–order differential equations to a system of first–order differential
equations, by introducing two new variables which equal the derivatives of the original variables. Let
u0 = w,
v 0 = x,
then the original differential equations can be written as
w0 = (t + 1) w2 − 2uv + u3 + cos t
1 3
x0 = 3u − (sin t) wx + t +
2 2
Therefore, the corresponding system of first–order differential equations is
0
u = w u (0) = 1
0
v = x v (0) = 3
0 (1)
w = (t + 1) w2 − 2uv + u3 + cos t w (0) = 2
0
x = 3u − 12 (sin t) wx + t + 32 x (0) = 4
To see how the algorithms can be adapted to systems such as (1), note first that the system (1) can
be written as 0
u = U (t, u, v, w, x) u (0) = 1
0
v = V (t, u, v, w, x) v (0) = 3
0 (2)
w = W (t, u, v, w, x) w (0) = 2
0
x = X (t, u, v, w, x) x (0) = 4
3
where
U (t, u, v, w, x) = w,
V (t, u, v, w, x) = x,
W (t, u, v, w, x) = (t + 1) w2 − 2uv + u3 + cos t,
X (t, u, v, w, x) = 3u − 12 (sin t) wx + t + 23 .
The basic idea is to adapt the algorithm for solving problems of the form
y 0 = f (t, y) , y (t0 ) = y0 (3)
to solving problems of the form (2), by replacing y by u, v, w, and x, and replacing f by U, V, W
and X. Note that the algorithm is not carried out separately for each unknown, but simultaneously
for all four.
2 1 3
In this case of the second–order Runge–Kutta method with a = , b = , α = β = , the algorithm
3 3 2
for (3) is
2 1
yn+1 = yn + k1 + k2 ,
3 3
k1 = hf (tn , yn ) ,
3 3
k2 = hf tn + h, yn + k1 .
2 2
Adapted to the system (2), this gives us the algorithm
2 1
un+1 = un + ku1 + ku2
3 3
2 1
vn+1 = vn + kv1 + kv2
3 3
2 1
wn+1 = wn + kw1 + kw2
3 3
2 1
xn+1 = xn + kx1 + kx2
3 3
ku1 = h U (tn , un , vn , wn , xn )
kv1 = hV (tn , un , vn , wn , xn )
kw1 = hW (tn , un , vn , wn , xn )
kx1 = hX (tn , un , vn , wn , xn )
4
APM3711/201
3 3 3 3 3
ku2 = hU tn + h, un + kv1 , vn + kv1 , wn + kw1 , xn + kx1
2 2 2 2 2
3 3 3 3 3
kv2 = hV tn + h, un + kv1 , vn + kv1 , wn + kw1 , xn + kx1
2 2 2 2 2
3 3 3 3 3
kw2 = hW tn + h, un + kv1 , vn + kv1 , wn + kw1 , xn + kx1
2 2 2 2 2
3 3 3 3 3
kx2 = hX tn + h, un + kv1 , vn + kv1 , wn + kw1 , xn + kx1 .
2 2 2 2 2
Hence, with h = 0.2, we calculate u (0.2) , v (0.2) , u0 (0.2) and v 0 (0.2) from the given initial values as
follows:
ku1 = (0.2) U (0, 1, 3, 2, 4) = (0.2) (2) = 0.4
kv1 = (0.2) V (0, 1, 3, 2, 4) = (0.2) (4) = 0.8
kw1 = (0.2) W (0, 1, 3, 2, 4)
= (0.2) (0 + 1) 22 − 2 · 1 · 3 + 13 + 1 = 0
kx1 = (0.2) X (0, 1, 3, 2, 4)
= (0.2) (3 · 1 − (0.5) · 0 · 2 · 4 + 0 + 1.5) = 0.9
3 3 3 3 3
ku2 = (0.2) U 0 + (0.2) , 1 + (0.4) , 3 + (0.8) , 2 + (0) , 4 + (0.9)
2 2 2 2 2
= (0.2) U (0.3, 1.6, 4.2, 2.0, 5.35)
= (0.2) · 2 = 0.4
kv2 = (0.2) V (0.3, 1.6, 4.2, 2.0, 5.35)
= (0.2) · 5.35 = 1.07
kw2 = (0.2) W (0.3, 1.6, 4.2, 2.0, 5.35)
(0.2) (0.3 + 1) 22 − 2(1.6)(4.2) + (1.6)3 + cos (0.3)
=
= −0.6377
kx2 = (0.2) X (0.3, 1.6, 4.2, 2.0, 5.35)
= (0.2) [3(1.6) − (0.5) · sin(0.3) · 2 · (5.35) + 0.3 + 1.5]
= 1.0038
2 1
u (0.2) = 1 + (0.4) + (0.4) ≈ 1.4
3 3
2 1
v (0.2) = 3 + (0.8) + (1.07) ≈ 3.89
3 3
5
2 1
u0 (0.2) = w (0.2) = 2 + (0) + (−0.6377) ≈ 1.79
3 3
2 1
v 0 (0.2) = x (0.2) = 4 + (0.9) + (1.0038) 9 ≈ 4.93
3 3
Question 3
The differential equation
y 0 = y − x2 y(0) = 1
and starting values
y(0.2) = 1.2186, y(0.4) = 1.4682, y(0.6) = 1.7379.
Use the fourth–order Adams–Moulton predictor–corrector method with h = 0.2 to solve the
equation up to x = 1.2.
Compare with the analytical solution,
y = x2 + 2x + 2 − ex .
Solution:
Given
y 0 = y − x2
y(0) = 1; y(0.2) = 1.2186;
y(0.4) = 1.4682; y(0.6) = 1.7379
Analytical solution
dy
− y = −x2
dx
This is a first order linear equation and can be put into the form
y 0 + P (x)y = Q(x)
P (x) = −1
and
Q(x) = −x2
The integrating factor is
R R
−dx
e P (x)dx
=e = e−x
∴ y 0 e−x − ye−x = −x2 e−x
6
APM3711/201
Z
−x
ye = − x2 e−x dx
= x2 e−x + 2xe−x + 2e−x + c
(using integration by parts)
∴ y = x2 + 2x + 2 + cex
x = 0 : 1 = 2 + ce0
∴ c = −1
∴ y = x2 + 2x + 2 − ex
Numerical solution
The Adams–Moulton method is a multistep method, i.e. more than one past value is used to
approximate the new value. In Gerald four multistep methods are discussed.
(1) The 3rd order Adams method (see (5.13) [(5.10)]) has a local error of order 4 and a global
error of order 3:
h
yi+1 = yi + (23fi − 16fi−1 + 5fi−2 )
12
where fi means f (xi , yi ).
(2) The 4th order Adams method (also called the Adams–Bashforth method; see (5.15) [p. 363])
is
h
yi+1 = yi + (55fi − 59fi−1 + 37fi−2 − 9fi−3 )
24
(3) Milne’s predictor–corrector method (see p. 413 [364]) is
4h
yp = (2fi − fi−1 + 2fi−2 ) + yi−3
3
h
yi+1 = yi−1 + (f (xi + h, yp ) + 4fi + fi−1 )
3
(4) The Adams–Moulton method (also called the Adams predictor–corrector method; see (5.24)
– (5.25) [(5.22)]) is
h
yp = yi + (55fi − 59fi−1 + 37fi−2 − 9fi−3 )
24
h
yi+1 = yi + (9f (xi + h, yp ) + 19fi − 5fi−1 + fi−2 )
24
The last three methods all have a local error of order h5 and a global error of order h4 , but Milne’s
method is sometimes unstable, and the error constant for the Adams–Moulton method is smaller than
that of the 4th order Adams method (19/720 acs opposed to 251/720).
See pages 409, 414, 417–418 [361, 365, 369–370] for a discussion of the relative merits of these multistep
methods and the Runge–Kutta methods.
7
The program and results follows. You should make sure that you can reproduce the table of results
by using a calculator.
PROGRAM A99 1 3(output);
CONST
size = 10;
VAR
i : integer;
h, x, y, yp, fp, exact : real;
fn : array[0..size] of real;
z : text;
FUNCTION f(x, y : real) : real;
BEGIN
f := y - x*x;
END; {f}
FUNCTION yexact(x : real) : real;
BEGIN
yexact := x*x + 2*x + 2 - exp(x);
END; {yexact}
BEGIN
assign(z,’AS99 1 3.DAT’);
rewrite(z);
writeln(z);
writeln(z,’ ***** ASSIGNMENT 1, ’,
’QUESTION 3(b) *****’);
writeln(z);
writeln(z,’ *********** ADAMS-MOULTON’,
’ METHOD **********’);
h := 0.2;
writeln(z,’ h = ’,h:4:2);
writeln(z,’ x exact y y error’,
’ f(x,y) yp f(x+h,yp)’);
fn[0] := f(0.0, 1.0);
fn[1] := f(0.2, 1.2186);
fn[2] := f(0.4, 1.4682);
x := 0.6;
y := 1.7379;
write(z,x:6:2,’ ’,y:7:4,’ ’,y:7:4,’ ’,
(0.0):7:4);
FOR i := 3 to 5 DO
BEGIN
fn[i] := f(x,y);
yp := y + h*(55*fn[i] - 59*fn[i-1]
8
APM3711/201
+ 37*fn[i-2] - 9*fn[i-3])/24;
fp := f(x + h,yp);
y := y + h*(9*fp + 19*fn[i] - 5*fn[i-1]
+ fn[i-2])/24;
x := x + h;
writeln(z,’ ’,fn[i]:7:4,’ ’,yp:7:4,’ ’,
fp:7:4);
exact := yexact(x);
write(z,x:6:2,’ ’,exact:7:4,’ ’,y:7:4,’ ’,
(exact - y):7:4);
END; {for i}
writeln(z);
close(z);
END.
***** ASSIGNMENT 1, QUESTION 3(b) *****
*********** ADAMS-MOULTON METHOD **********
h = 0.20
x exact y y error f(x,y) yp f(x+h,yp)
0.60 1.7379 1.7379 0.0000 1.3779 2.0146 1.3746
0.80 2.0145 2.0145 -0.0000 1.3745 2.2819 1.2819
1.00 2.2817 2.2817 -0.0000 1.2817 2.5202 1.0802
1.20 2.5199 2.5199 -0.0000
9
Question 4
(a) The first three Chebyshev polynomials are
T0 = 1
T1 = x
T2 = 2x2 − 1
T3 (x) = 4x3 − 3x
T4 (x) = 8x4 − 8x2 + 1.
(i) Economize the truncated power series
p(x) = 1 + 2x − x3 + 3x4 .
(ii) Give an upper limit for the absolute value of the difference between the original truncated
power series and the economized one.
(b) Find the Padé approximation R4 (x), with denominator of degree 2, to the function
f (x) = 4 + 4x2 − 2x6 .
Solution:
(a) (i)
p (x) = 1 + 2x − x3 + 3x4
Subtract: 43 T4 (x) = 34 − 3x4 − 3x2 + 38 from p (x) , in order to remove the highest order
term as to get the economized series, g (x) = −x3 + 3x2 + 2x + 85
(ii) For all Chebysshev polynomials Ti we have |Ti | ≤ 1 for all x ∈ [−1, 1] , and therefore
|p (x) − g (x)| = 83 T4 (x) ≤ 38 = 0.375
(b) The Padé approximation is
a0 + a1 x + a2 x 2
R4 (x) =
1 + b1 x + b2 x 2
This gives
a0 + a1 x + a2 x 2
R4 (x) − f (x) = − 4 − 4x2 + 2x6
1 + b1 x + b2 x 2
which gives:
4 + 4b1 x + 4 (1 + b2 ) x2 + (4b1 − a3 ) x3 + 4b2 x4 − a0 − a1 x − a2 x2 .
Equate coefficients of like powers of x :
10
APM3711/201
x0 : 4 − a0 = 0
x1 : 4b1 − a1 = 0
x2 : 4 (1 + b2 ) − a2 = 0
x3 : 4b1 = 0
x4 : 4b2 = 0
Therefore,
b1 = 0, a0 = 4, a1 = 0, b2 = 0, a2 = 4
Thus,
R4 (x) = 4 + 4x2
11
