UNIT
Numerical Methods – II 2
2.1 NUMERICAL SOLUTION OF SIMULTANEOUS
FIRST ORDER ORDINARY DIFFERENTIAL
EQUATIONS
2.1.1 Runge-Kutta Method of Fourth Order
dy dz
1. Let = f (x, y, z), = g(x, y, z) be the simultaneous first order ordi-
dx dx
nary differential equations with the initial values y(x0 ) = y0 , z(x0 ) = z0 .
The modified form of expression for computation are as follows.
k1 = h f (x0 , y0 , z0 ); l1 = hg(x0 , y0 , z0 )
h k1 l1 h k1 l1
k2 = h f x0 + , y0 + , z0 + ; l2 = hg x0 + , y0 + , z0 +
2 2 2 2 2 2
h k2 l2 h k2 l2
k3 = h f x0 + , y0 + , z0 + ; l3 = hg x0 + , y0 + , z0 +
2 2 2 2 2 2
k4 = h f x0 + h, y0 + k3 , z0 + l3 ; l4 = hg x0 + h, y0 + k3 , z0 + l3
Then according to second order Runge-Kutta method
y(x0 + h) = y0 + k2
z(x0 + h) = z0 + l2
According to fourth order Runge-Kutta method
1
y(x0 + h) = y0 + k1 + 2k2 + 2k3 + k4 and
6
1
z(x0 + h) = z0 + l1 + 2l2 + 2l3 + l4
6
46 Engineering Mathematics-IV
= .1 [1.11 – 0.1] = 0.101
l4 = hg t0 + h, x0 + k3 , y0 + l3
= 0.1 [1.1005 + 0.1] = 0.12005
1
x(t0 + h) = x0 + k1 + 2k2 + 2k3 + k4
∴
6
1
=1+ [0.1 + 2 × 0.1 + 2 × 0.1005 + 0.101]
6
= 1.100333
and
1
y(t0 + h) = y0 + l1 + 2l2 + 2l3 + l4
6
1
=1+ [0.1 + 2(0.11) + 2(0.11) + 0.12005]
6
= 1.1100083
Thus,
x(0.1) = 1.100333 1.1
y(0.1) = 1.1100083 1.11.
Example 2: Use fourth order Runge-Kutta method to solve the system of
equations
dy dz
= 1 + xz; = – xy given y(0) = 0, z(0) = 1 at x = 0.3.
dx dx
Solution: Given simultaneous equations are
dy dz
= 1 + xz, = – xy and y(0) = 0, z(0) = 1
dx dx
or, f (x, y, z) = 1 + xz, g(x, y, z) = – xy
Here x0 = 0, y0 = 0, z0 = 1 and given h = 0.3.
Numerical Methods–II 47
Now, we shall compute the following.
k1 = h f (x0 , y0 , z0 ) = 0.3 [1 + 1 × 0] = 0.3
l1 = hg(x0 , y0 , z0 ) = 0.3 [ – 0 × 0] = 0
h k l
k2 = h f x0 + , y0 + 1 , z0 + 1
2 2 2
= 0.3 [1 + (1)(0.15)] = 0.345
h k l
l2 = hg x0 + , y0 + 1 , z0 + 1
2 2 2
= 0.3 [ – (0.15)(0.15)] = – 0.00675
h k l
k3 = h f x0 + , y0 + 2 , z0 + 2
2 2 2
= 0.3 [(1 + 0.996625)(0.15)] = 0.34485
h k l
l3 = h f x0 + , y0 + 2 , z0 + 2
2 2 2
= 0.3 [ – (0.15)(0.1725)] = – 0.0077625
k4 = h f x0 + h, y0 + k3 , z0 + l3
= 0.3 [1 + (0.99224)(0.3)] = 0.3893
l4 = h f x0 + h, y0 + k3 , z0 + l3
= 0.3 [ – 0.3(.34485)] = – 0.03104
1
y(x0 + h) = y0 + k1 + 2k2 + 2k3 + k4
∴
6
1
∴ y(0.3) = 0 + [.3 + 2(.345) + 2(0.34485) + 0.3893]
6
= 0.34483
48 Engineering Mathematics-IV
and
1
z(x0 + h) = z0 + l1 + 2l2 + 2l3 + l4
6
1
=1+ [0 + 2( – 0.00675) + 2( – 0.0077625) + ( – 0.03104)]
6
∴ z(0.3) = 0.98999.
2.1.2 Picard’s Method
Let the simultaneous differential equations be
dy dz
= f (x, y, z) and = g(x, y, z)
dx dx
with initial conditions y(x0 ) = y0 and z(x0 ) = z0 . Then
y1 = y0 + f (x, y0 , z0 ) dx, z1 = z 0 + g(x, y0 , z0 ) dx
y2 = y0 + f (x, y1 , z1 ) dx, z2 = z 0 + g(x, y1 , z1 ) dx
y3 = y0 + f (x, y2 , z2 ) dx, z3 = z 0 + g(x, y2 , z2 ) dx
and so on.
Continuing this process, a sequence of functions of x, i.e., y1 , y2 , y3 . . . and
z1 , z2 , z3 . . . are obtained each giving a better approximations of the desired
solution.
Example 1: Using Picard’s method find approximate values of y and z corre-
dy dz
sponding to x = 0.1 given that y(0) = 2, z(0) = 1 and = x + z, = x – y2 .
dx dx
Solution: Given equations are
dy dz 2
= f (x, y, z) = x + z and = g(x, y, z) = x – y
dx dx
x x
∴ y = y0 + f (x, y, z) dx and z = z0 + g(x, y, z) dx
x0 x0
Numerical Methods–II 49
First approximation
x x
x2
y1 = y0 + f (x, y0 , z0 ) dx = 2 + (x + 1) dx = 2 + x +
2
x0 0
x x
1
z1 = z0 + g(x, y0 , z0 ) dx = 1 + (x – 4) dx = 1 – 4x + x2
2
x0 0
Second approximation
x x
1
y2 = y0 + f (x, y1 , z1 )dx = 2 + x + 1 – 4x + x2 dx
2
x0 0
3 2 x3
=2+x– x +
2 6
x x
1 2 2
z2 = z0 + g(x, y1 , z1 )dx = 1 + x– 2+x+ x dx
2
x0 0
3 2 3 x4 x5
= 1 + 4x + x –x – –
2 4 20
Third approximation
x
y3 = y0 + f (x, y2 , z2 )dx
x0
3 2 1 3 1 4 1 5 1 6
=2+x– x – x – x – x – x
2 2 4 20 120
x
z3 = z0 + g(x, y2 , z2 )dx
0
3 2 5 3 7 4 31 5 1 6 1 7
= 1 – 4x – x – x + x – x + x – x
2 3 12 60 12 252
and so on.
50 Engineering Mathematics-IV
when x = 0.1
y1 = 2.105, y2 = 2.08517 y3 = 2.08447
z1 = 0.605, z2 = 0.58397 z3 = 0.58672
Thus, y(0.1) = 2.0845, z(0.1) = 0.5867 correct to four decimal places.
2.1.3 Taylor’s Series Method
Let the simultaneous differential equations be
dy
= f (x, y, z) (1)
dx
dz
and = g(x, y, z) (2)
dx
with initial conditions y(x0 ) = y0 and z(x0 ) = z0 .
If h be the step-size, y1 = y(x0 + h) and z1 = z(x0 + h)
Then Taylor’s algorithm for (1) and (2) gives
h2 h3
y1 = y0 + hy0 + y + y ... (3)
2! 0 3! 0
h2 h3
z1 = z0 + hz0 + z + z ... (4)
2! 0 3! 0
Differentiating (1) and (2) successively we get y , z , etc.
So the values y0 , y0 , y0 . . . and z0 , z0 , z0 . . . are known. Substituting these
in (3) and (4), we obtain y1 , z1 for the next step.
Similarly,
h2 h3
y2 = y1 + hy1 + y + y ... (5)
2! 1 3! 1
h2 h3
z2 = z1 + hz1 + z + z ... (6)
2! 1 3! 1
Since y1 and z1 are known, we calculate y1 , y1 . . . and z1 , z1 . . . substituting
these in (5) and (6) we get y2 and z2 .
Proceeding further, we can calculate the other values of y and z step by
step.
Numerical Methods–II 51
dy dz
Example 1: Given = z and = – xz – y with y(0) = 1, z(0) = 0, obtain
dx dx
y and z for x = 0.1, 0.2, 0.3 by Taylor’s series method.
Hint:
We have
y = z and z = – xz – y
We use Taylor’s series method to find y and z.
2.2 NUMERICAL SOLUTION OF SECOND ORDER
ORDINARY DIFFERENTIAL EQUATIONS
2.2.1 Picard’s Method
Let the second order differential equation be
d2y dy
= f x, y, (1)
dx2 dx
with y(x0 ) = y0 and y (x0 ) = y0 .
dy d 2 y dz
Now, let = z then = .
dx dx2 dx
Substituting in (1)
dz
= f (x, y, z) with y(x0 ) = y0 and z(x0 ) = z0
dx
The problem reduces to solving the simultaneous equations
dy
= z = f 1 (x, y, z)
dx
dz
= f (x, y, z) subject to y(x0 ) = y0 and z(x0 ) = z0 .
dx 2
By Picard’s method we have
x x
y1 = y0 + f1 (x, y0 , z0 ) dx, z1 = z0 + f2 (x, y0 , z0 ) dx
x0 x0
x x
y2 = y0 + f1 (x, y1 , z1 ) dx, z2 = z0 + f2 (x, y1 , z1 ) dx
x0 x0
Numerical Methods–II 53
x x
x2
.2 2
= 0.5 + .1 – 0.2 – 0.5x dx, z2 = 0.1 + – 2x 0.1 – x – .5x
2 2
0 0
– (0.5 + 0.1x) dx
x3 x2 x2 x4 x3
= 0.5 + 0.1x – 0.2 – 0.5 , = 0.1 – .2 + 0.2 +1
6 2 2 4 3
x2
– 0.5x – 0.1
2
and so on.
when x = 0.1
y1 = 0.51 y2 = 0.5074
z1 = y1 = 0.049 y2 = z2 = 0.0488
∴ y(0.1) = 0.5074
y (0.1) = 0.0488
Example 2: Using Picard’s method, obtain the second approximation to the
solution of
d 2 y 3 dy 3 1
=x +x y so that y(0) = 1, y (0) = .
dx2 dx 2
Solution: Given the second order differential equation is
d 2 y 3 dy 3 1
=x +x y with y(0) = 1, y (0) =
dx2 dx 2
dy d 2 y dz
Let =z so that =
dx dx2 dx
dy dz 3 3 1
∴ =z and =x z+x y with y(0) = 1, z(0) = .
dx dx 2
= f (x, y, z) = f 1 (x, y, z)
By Picard’s methods,
54 Engineering Mathematics-IV
First approximation
x x
y1 = y0 + f (x y0 z0 )dx z1 = z0 + f2 (x, y0 z0 ) dx
x0 x0
x x
1 1 1 3 1 3 3
=1+ dx = + x + 1x3 dx = + x dx
2 2 2 2 2
0 0
x 1 3x4
=1+ = +
2 2 8
Second approximation
x x
y2 = y0 + f (x, y1 , z1 ) dx z 2 = z0 + f (x, y1 , z1 ) dx
x0
x0
x
1 3x4 4
x
1 1 3x x
=1+ + dx = + x3 +
3
+x 1+ dx
2 8 2 0 2 8 2
0
x
x 3x5 x3 3 7 3 x4
1
=1+ + = + + x +x + dx
2 x0 2 2 8 2
0
1 x4 3x8 x4 x5
= + + + +
2 8 64 4 10
1 3x4 x5 3x8
+ + + .=
2 8 10 64
Example 3: Using Picard’s method find the solution of the equation
y + xy + y = 0, with y(0) = 1 y (0) = 0 for x = 0.1 and 0.3.
Solution: Given equation is
y + xy + y = 0 with y(0) = 1, y (0) = 0.
dy d2y dz
Let = y = z so that
2
= y =
dx dx dx
dy dz
∴ =z and = – xz + y with y(0) = 1, z(0) = 0
dx dx
= f (x, y, z), = f 1 (x, y, z).
56 Engineering Mathematics-IV
when x = 0.1
y1 (0.1) = 1 y2 (0.1) = 0.955, y3 (0.1) = 0.99501
z1 (0.1) = – 0.1 z2 (0.1) = – 0.099, z3 (0.1) = – 0.0995
∴ y(0.1) = 0.99501
z(0.1) = y (0.1) = – 0.0995
when x = 0.3
y1 (0.3) = 1 y2 (0.3) = 0.955 y3 (0.3) = 0.9556
z1 (0.3) = – 0.3 z2 (0.3) = – 0.291, z3 (0.3) = – 0.2911
∴ y(0.3) = 0.9556
z(0.3) = y (0.3) = – 0.2911.
2.2.2 Runge-Kutta Method
Let the second order differential equation be
d2y dy
= f x, y, . . . (1) with y(x0 ) = y0 , y (x0 ) = y0 .
dx2 dx
dy d 2 y dz
Now, let = z so that =
dx dx2 dx
dy dz
∴ =z and = f (x, y, z)
dx dx
= f 1 (x, y, z), = f 2 (x, y, z) with y(x0 ) = y0 , z(x0 ) = z0 .
The problem reduces to solving the simultaneous equations
dy dz
= z = f 1 (x, y, z) and = f (x, y, z) with y(x0 ) = y0 and z(x0 ) = z0 .
dx dx 2
Starting at x0 , y0 , z0 and taking the step-sizes for x, y, z to be h, k, l respec-
tively, the Runge-Kutta method gives
Numerical Methods–II 57
K1 = h f1 (x0 , y0 , z0 ) l1 = h f2 (x0 , y0 , z0 )
1 k1 l1 h k1 l1
K2 = h f1 x0 + h, y0 + , z0 + l 2 = h f 2 x 0 + , y0 + , z 0 +
2 2 2 2 2 2
h k2 l2 h k2 l2
k3 = h f1 x0 + , y0 + , z0 + l 3 = h f 2 x 0 + , y0 + , z 0 +
2 2 2 2 2 2
k4 = h f1 x0 + h, y0 + k3 , z0 + l3 l4 = h f2 x0 + h, y0 + h3 , z0 + l3
1 h
y1 = y0 + k1 + 2k2 + 2k3 + k4 z 1 = z0 + l1 + 2l2 + 2l3 + l4
∴
6 6
Example 1: Using Runge-Kutta method solve
y = xy 2 – y2 for x = 0.2 correct to 4 decimal places initial conditions are
x = 0, y = 1, y = 0.
Solution: Given equation of second order is
y = xy 2 – y2 with x = 0, y = 1, y = 0
dy
Let = z = f 1 (x, y, z), then
dx
dz 2 2
= xz – y = f 2 (x, y, z) with x0 = 0, y0 = 1, z0 = 0
dx
we compute k1 , k2 , k3 , k4 for f1 (x, y, z) and l1 , l2 , l3 , l4 for f2 (x, y, z).
∴ by Runge-Kutta formulae, we have
k1 = h f1 (x0 , y0 , z0 ) l1 = h f2 (x0 , y0 , z0 )
= 0.2(0) = 0 = 0.2( – 1) = – 0.2
h k l h k1 l1
k2 = h f 1 x0 + , y 0 + 1 , z 0 + 1 l 2 = h f 2 x0 + , y 0 + , z 0 +
2 2 2 2 2 2
= 0.2( – 0.1) = – 0.02 = 0.2( – 0.999) = – 0.1998
h k l h k2 l2
k3 = h f1 x0 + , y 0 + 2 , z 0 + 2 l 3 = h f 2 x0 + , y 0 + , z 0 +
2 2 2 2 2 2
= 0.2( – 0.0999) = – 0.02 = 0.2( – 0.9791) = – 0.1958
Numerical Methods–II 59
dy d 2 y dz
Let = z then = .
dx dx2 dx
Substituting in (1)
dz
= f (x, y, z) with y(x0 ) = y0 and z(x0 ) = z0
dx
The problem reduces to solving the simultaneous equations.
dy
= z = f 1 (x, y, z)
dx
dz
= f (x, y, z) subject to y(x0 ) = y0 , and z(x0 ) = z0 .
dx 2
The predictor and corrector formulae can be written as
4h
y(P)
n+1
= yn – 3 + 2y n–2
–y
n–1
+ 2yn
3
4h
z(P)
n+1
= zn – 3 + 2z n–2
–z
n–1
+ 2zn
3
and
h
y(c) = yn – 1 + y + 4yn + y
n+1 3 n–1 n+1
h
zn(c)+ 1 = zn – 1 + zn – 1 + 4zn + zn + 1
3
In particular when n = 3
(P) 4h
y4 = y0 + 2y1 – y2 + 2y3
3
(2)
(P) 4h
z 4 = z0 + 2z1 – z2 + 2z3
3
and
(c) h
y4 = y2 + y + 4y3 + y4
3 2
(3)
(c) h
z 4 = z2 + z2 + 4z3 + z4
3
(2) is called Milne’s predictor formula and (3) is called Milne’s corrector for-
mula.
60 Engineering Mathematics-IV
d2y dy
Example 1: Given 2
+ x + y = 0, y(0) = 1, y (0) = 0 obtain y for
dx dx
x = 0, 0.1, 0.2, 0.3 by any method. Further compute y(0.4) by Milne’s method.
Solution: The given equation is
y + xy + y = 0 . . . (1) with y(0) = 1, y(0)
=1
Let y = z then y = z
∴ z + xz + y = 0 and y = z
dy
= y = z = f 1 (x, y, z) with y(0) = 1
∴
dx
dz
= z = f 2 (x, y, z) = – xz – y with z(0) = 0
dx
We use Taylor’s series method to find y. Differentiating (1) n times we get
yn + 2 + xyn + 1 + nyn + yn = 0
At x = 0, yn + 2 = – (n + 1)yn
0
∴ y(0) = 1 gives y2 (0) = – 1, y(0)
4
= 32, y
(0)
6
= – 5×3
and y1 (0) = 0 gives y3 (0) = y5 (0) = . . . = 0.
Expand y(x) by Taylor’s series
x2 x3
y(x) = y(0) + xy1 (0) + y2 (0) + y3 (0) + . . .
2! 3!
x2 3 4 5 × 3 6
∴ y(x) = 1 –
+ x – x + ... (2)
2! 4! 6!
1 1
and z(x) = y (x) = – x + x3 – x5 . . . = – xy (3)
2 8
From (2) we have
(0.1)2 1 4
y(0.1) = 1 – + (.1) + . . . = 0.995
2 8
(0.2)2 (0.2)4
y(0.2) = 1 – + + . . . = 0.9802
2 8
(0.3)2 (0.3)4 (0.3)6
y(0.3) = 1 – + – . . . = 0.956.
2 8 48
62 Engineering Mathematics-IV
and
h
y(0.4) = y(0.2) + y (0.2) + 4y (0.3) + y (0.4)
3
0.1
= – 0.9802 + [ – 0.196 – 1.1452 – 0.3692]
3
= 0.9232
Hence, y(0.4) = 0.9232 and y(0.4) = – 0.3692.
d2y dy
Example 2: Given 2
+ 2x + y = 0 and y(0) = 0.5, y (0) = 0.1 obtain y
dx dx
for x = 0, 0.1, 0.2, 0.3 by any suitable method. Hence, apply Milne’s method
to compute y(0.4).
2
Example 3: Given y = xy – y2 and y(0) = 1, y (0) = 0 obtain y for
x = 0, 0.1, 0.2, 0.3 by Taylor’s series method. Hence, apply Milne’s method
to compute y(0.4).
Ans: y(0.4) = 0.9232
2.2.4 Taylor’s Series Method
Example 1: Given y + xy + y = 0, y(0) = 1, y (0) = 0. Obtain y for
x = 0, 0.1, 0.2, 0.3 by Taylor’s series method.
Solution: Given second order differential equation is
y + xy + y = 0 (1)
Let y = z then y = z
∴ z + xz + y = 0; y = z (2)
We use Taylor’s series method to find y.
Differentiating the given equation n times we get
yn + 2 + xyn + 1 + nyn + yn = 0
At x = 0,
(yn + 2 )0 = – (n + 1)yn
∴ y(0) = 1 gives y2 (0) = – 1, y4 (0) = 3, y6 (0) = – 5 × 3 and y1 (0) gives
y3 (0) = y5 (0) = . . . = 0
