Nwnerical Solution of Ordinary DUJerenllal & 111atlons 1209

,,.,,,--
.. l + O. l + 2 (0.000 3333) + 0.01 + 0.000001 333 + 0.00001666
: 1.110683993

io,!S EULER'S METHOD

conSl·der a first order differential equation
dy
= f(x, y)
dx ... (11)

1lieri ,rbere Y (xJ = Yo ... {1 2J

rot us take points Xo, XI' Xi, ...... xn, ... .. such that X = X t+l + h.
J.i', • I
r,et the actual solution of the differential equation be den oted by
the curve given in Fig. 10. 1 With point (Xa = ye) lying on the curve.
We will calculate the value of y of the cmve at x = x 1• The equation
of tangent at (Xa, ye) to the cmve is given by
Y - Yo = Y (Xo, Ye) [Xa]
= f (Xo, Ye) [x-Xo]
Y = Yo + f (Xo, Yo) [x - Xo]
This gives the value of y on the tangent corresponding to value
)I x. In the interval (Xa, yI), the cmve is approximated by the tangent.
Therefore, the value of y on the cmve is approximately equal to
the value of y on the tangent at (Xa, ye) corresponding to Xo, x I.
Therefore, YI = Yo + f (Xo, Y0 ) [x1 - Xa] = Yo + hf (Xo, Yol

--+----------Y

Fig. 10.1

Next, we approximate the curve by the line through (x1 • Y 1) for

Wh.tch slope is J (xI, y ). We get
1
Y2 = Y1 + f (x1, YI) [~ - X1]

= Y 1 hf(x1, Y1l
 210 N11mc rfcnl M et11ods 11.1fth Program ln MA 11.AU

Proceeding in the s am e way, W<' get the gener

as al itf'te1
Y,H 1 == Yn + 1!1 (xn, Y,) for rt ::: 0
This fonnula is called ~ulcr's AlgorJthrn. ' 12 J

In this method, we approximate the actua1 <'IJ

of short straight lines. As the length of Int rve by
straight lines deViates much from the actual
error becomes quite significant.
•:al hit.
rve arid
The error at x == x 1 is given by

Error - (x1 - Xo )2 "( h2

- 2I y X1,Yi) ::: 2 if'fx1,y I
Therefore, it is of order h 2 •

10.6 MODIFIED EULER'S METHOD

In Euler's method, we approxunated the actual cunre of the
in the interval (Xa, x1) by the tangent taking
Y1 == Y1 + hf (XrJ, Yo)
Then we obtained the slope of the curve of the solute:
y 1 ) and hence the tangent at (XrJ, y ). •
0
In modified Euler's method, we find a better apprOXin:z:t:
of by y I taking the slope of the curve as the mean of the s::__
the tangents at Xa, y 0 ) and (x1 , y ), i.e.
1

y(l\ = Yo+~[J(xo,Yo)+f(x1,Yil]
As the slope of the tangent at (x1 , y 1) is not krio\\11, we ~
the value of y 1 obtained from Eq. (14) in f(x • Y) and then
1 1
it in Eq. (15) to find the modified value of Yi·
To find next modified value of Y,, we put the ,-a1u,.,r •
f( • JI) ) to obtain the modified value of the slope ~t [. 1,'
XI' !:F I
then substitute . f (x- • y ) by J(xl'.•
it in Eq. (15), rep Iacmg 1 1

. h[ ·(
.1Jll
y(2\ = Yo+ 2 J(xo,Yo) + J x,,!t ~

We repeat this step, tI·11 two consecutive '--·deredof~ i

. valut'S

ble The last obtained val ue Of y 1(,) is cons1

j inten·~11 1, ·(
agreea
is .
then taken as the starting point for the nex
next value Y2 is qbtained as
 Nunierlcal Solutto11 of Or<lf,,w I) ,
I
~--~ -- !I IJJ, ,, rttlr,t Eqr.taJl.on.s [ 21
Y2 - Y0 r IU (x
1
, y J
1
,i1e. b,,tter
,. approxi1natlons are obtain er as
1
1
Y2° = Yi+~ ff(x 1 ,!Ji) + Jfx:1.,YzJ]

~ + ~ [ J(x1 ,y1J+ J(x,,y~')]

2
Y2 [ ) = Y1
~
id so on. We continue this process till Y2 becomes stauo
~ we proceed
~1en to obtain y3 as above and so on • naiy.

This ts modified Euler's met~od Which provides great

iJJJprovement in accuracy over Euler s method.

txample 10.6 Using Euler's method, Solve dy = x + y where

dx
iOJ = I and obtain ~e values of y(l,0) taking step length h = o. .
Check your answer with the exact solution. 2
Solution: We have j(x, y) = x + y, Xo = 0, y = I
0
Since h x 1 = 0.2, ~
= 0.2, = 0.4, X":3 = 0.6, x
4
= 0.8, X's = 1.0
By Euler's algortthm
Y1 = Yo+ hf (Xo, Yo) = l + (0.2) [0 +l] = 1.2
Y2 = Y1 + hf (xl, Y1 ) = 1.2 + (0.2) [0.2 + 1.2]
= 1.48
Y3 = Y2 + hf(~. Y2 ) = 1.48 + (0.2) [0.4 + 1.48]
= 1.856
Y4 = Y3 + hf (X:3, y 3 ) = 1.856 + (0.2) [0.6 + 1.856]
= 2.3472
y 5 = y4 + hf (x4 , y 4 ) = 2.3472 + (0.2) (0.8 + 2.3472]
= 2.94664
Now the exact solution of the equation is
y(x) = 2eK - x - l
y1 = y(0.2) = 1.2428
y2 = y(0.4) = 1.5836
y 3 = y(0.6) = 2.0442
y = y(l.O) = 3.4366
Note that the 5values of y deviates irom
r, the exact values as x
creases.
 Numerical Solution of Orrltnrmi f)i'{f' Ii L ,.•,. [
,,.---------- ----.:--..:.:
. <'rr•n a a.Ai1Wl.lons 21
.:=:~==~!! 5 _---
J0,7 aUNGE-KUTTA METHODS

111ese methods .refer ~o ~

class of one-step methods used for the
,. n1erical solution of initial value problems. The previous methods
I, :: solve the differential equations have the disadvantage of slow
~(\ n\•ergence. Runge- Kutta methods deals with this issue by not
' ~~ploying higher order derivatives and involving the function values
at some selected points. Runge-Kutta methods are known by their
order. In ,th order Runge-Kutta method, we consider Taylor's series
solution up to the terms in hr .These methods are all based on
the general form of the extrapolation equation
Yi+l = Yi + Slope x Inteival size = yi +
kh ... (16)
where k represents the slope that is the weighted averages •of
the slopes at various points in the inteival h.
First order Runge-Kutta method
From Euler's method, we have
Y1 = Yo + hf (Xo, Yo)
By Taylor's expansion

Y1 = Y (Xo+h) = Yo +hy'(xo)+ ~ Yo"(xo)+ ...... .

i~ Note that the Euler's method agrees with the Taylor's series
solution up to term in h.
Therefore, Euler's method is the Runge-Kutta method of order
one.
Second order Runge-Kutta method
From Modified Euler's method, we have

y 1 = Yo+; [J(xo,Yo) + f(x1,Y1)] ... (1 7)

Substituting the value of y 1 = Yo + hf (x0 • YJ in RHS of Eq. (171.

we get

Y1 = Yo+;[f(xo,Yo)+j(x1,Y o +hH-\o,tJo))] ... (18)

EXJ)andmg
. .flx y 0 + hif (Xo, Yo)) in Eq. (l8) using Taylor's selies
1, •
for a function of two variables, we get
J( Xi,Yo + hif(Xo,Yo )) = f(x 0 + h, Yo+ /1fo)

•
 ri al Methods wit. h Pr<ogram in MA7LAB
216 I
-------__ ,
Nume C r ' -------

= f(x0 ,Yol+ h(t~l + lifo ( ~l + ...... "•fl )

9
Substituting Eq. (19) in Eq. (18), we get

Yi = Yo+:[ f(x0 ,y0 )+ f(xo,Yo)+h(t~l +hfo(t1 + ... ]

[(!~1 +(~)J~h'J
2

= Yo +lif(xo,Yo)+ ~

= Yo + hfo + h 2 fo + O(h 3)
2

... (20)

Also expanding y 1 = y (Xo + h) on LHS of Eq. (18), we get

2
Y1 = y(xo + h) =Yo + hy'(xo) + : Yo" (xo) + .....

... (21)
Comparing Eq. (20) and (21), it can be seen that the modified
Euler's
term in method
h2 • agrees With the Taylor's senes solution up to the

For Runge-Kutta method of order two the expression in Eq. [16)

can be written as

where
k1 = hf (Xo, Yo) and
Js = hf (Xo + h. Yo + k )
Third order ~UDge-Kutta method
1

The algont1nn for Runge-Kutta method of order three is written as

l
Y1 == Yo+ (k + 4k + k )
6 1 2 3
 Numerical Solution of
0
rdinary Di[[ere11t(al Equations 217
,, ,,
'•
\,.11ere k 1 = hf(~. Yo)

/G, = hf ( Xo + ~ ,Yo + i)
., j
'•,
IC:3 = hf(~ + h. Yo + 21<:i - k )
rourth order Runge-Kutta method 1

This method is mostly used in b

Jllentioned. It is often referred a p1; lems unless otherwise
algorithm for Runge-Kutta method ~ udngefi-Ku~ta method. The
0
or er our 1s written as
1
Y1 = Yo + (k1 +2k2 +2k +k )
6 3 4
where k 1 = hf (Xo, Yo)

k, = hf ( Xo + ~ ' Yo + i)
1s = hf ( Xo + ; . Yo + i)
k4 = hf(Xo + h. Yo + Js)
An important advantage of these methods is that the operation
ls identical whether the differential equation is linear or non-linear.
An important advantage of these methods is that the operation
..[21) ls identical whether he differential equation is linear or non-linear.
ed
Uie Example IO.IO Solve dy = -y where y(O) = lusing Runge-Kutta
dx
method of; (i) Second order; (ii) Third order and (iii) Fourth order
and obtain the values of y(O. l) and y(0.2).
Solution: Here, J(x, y) = - y, Xo = 0, Yo = l. x 1 = 0.1, X:i= 0.2
Q) Second Order Runge-Kutta method
k 1 = hf {x0 , ye) =(0.1) (- Yo) =- (0. l)
k;_ = hf (Xo + h. Yo + k1 )
= _ h (Yo+ k 1) = _ (0.1) (1-0.1) =- 0.09
1
y(O. l) = Y1 = Yo + 2 (ki + k 2 )

= l + .!..(-0.1 - 0.09) = 0.905

2
 . I I [).,(VJr(lTII l\ft\'fl.AIJ
Numcn·,·al Met I rexJ.;;
i11
118 • tl'lf ,, '~

Again starting from (0 . l. 0.905) 1·cplacing (x0 • Y,,J by (x,

I
y, • ~
1
get le1 = hf (xi' y 1/ = It (- Y1l (0. l) (- 0,90-;J
= _ 0.0905
l<:i = hf (x1 + h, Y l + Jcl)
= -h (yl + k1J= - (0.1) (0.905- 0 .09051
=- 0.01845
1
y (0.2) = Y2 = Yi + 2 (ki + k 2 )
1
= 0.905 + 2 (-0.0905- 0.08145) = 0,819025
(ii) Third order Runge-Kutta method

k 1 = lif(xo,Yo)= h(-yo) = -0.l

Ts = !if ( Xo + ; 'Yo + :i)

= -h(Yo+';'.)= -0. 1(1- 0;1) =-0.095
Js = h_r (Xo + h, Yo + 2Js. - kl)
=-0.l(y0 +2k2 -ki)=-0.1(1-0.19+0.l)=-0.091
1
y(0.1) = Yi =yo +-(k1 +4k2 +k3)
6
1
= 1 + 6 (-0.1- 0.38 - 0.091) =0.904

getAgain starting from (0.I, 0.904) replacing (x,,. yJ by (x • yiJ. ~•

1

k1 = hf (x1 , Y1 ) = h (- y 1 ) =- 0.0904
ki = lif( X1 ~
+ •Y1 + :
1
) =-h(Yi+ i)
0
== -0.1(0.904- · 0904 ) = -0.08588
2
's == kh (xl + h. Y1 + kl) 2Js-
== - h (yl + kl) 2Js -
== - O. l (0.94- 0. 1717+ (0.0904) = - 0.0B2?"
-1
\ Numerical Solution 01 0
:1
d.
,-, mary Differenltal. Equations 219
\

y (0.2) == Y2 ::: Y1 + ~ (k1 + 4k2 + k3)

==
0 904
· + ~ (-0.0904 - 0.34352 - 0.08227)
== 0.8179
(iii) Fourth order Runge-Kutta method
kl == hf (JCo, Yo) = - 0.1

Js. == lif(Xo +~,Yo+ i) = -0.095

Is = lif( xo +~,Yo+'; )=-(o.1i(1- 0-~95 )
== - 0.09525
k4 == hf (Xo + h. Yo + Js) = - (0.1) (1- 0.09525)
== - 0.090475
l
y(O.l) == Y1=Yo+-(k1+2k2+2k 3+k4)
6
l
= l +-(-O.l-0.181-0.1905 -0.090475)
6
= 0.90633
Again starting from (0.1, 0.90633) replacing (JCo, yo) by (x , y 1),
1
we get
k 1 = 1if(x1 ,y1 )=-0.090633

Is_ = hf ( X1 + ~ ·Y1 + i)
= -0.1 ( 0.90633- 0.090633) = -0.086101
2

Js = hf (xi + ~ ,Y1 + i)
0.086101)
= -(0.1) ( 0.90633- 2 = -0.08632

k = hf(X1 + h,Y1 + k3)

4
= -(O. l)(0.90633- 0.08632) = -0.082001
 diJ
-
. flit n.vv,ram in MA TI.AB
Numerlcnl M<!lltCK1s w ',~,

y (Xo + h) = Yi , Y (Xo + 211) = Y2, Y (x 0 + 3 h ~

by Picard's method or Taylor's series m e thod. Y:j
Next we calculate
fo = f(Xo, Yo), fi = f (Xo + h, Y 1) ,
fz = f(Xo + 2 h. Y2) , /2 = f (Xo + 3h, Y:J
By Newton's forward interpolation formula, we h ave
2
L\ Yo a
y(x0 + nh) = Yn =Yo+ nL\yo + n(n - 1 ) - - + n(n - l)(n - 2) .1 Yo
2I ~+
3! ...
Changing y toy', the formula can be written as
A2 I
, , ~A , ( ) 0 Yo L\ 3 ,
Y =Yo + 1 LL.1Yo +n n-1 --+n(n-I)(n-2)_1!..Q_
21 3! + ......
Integrating both sides from Xo to x 4 ,we get

Jy'dx = Jf(x,y) dx = j
X4

x0
X4

xo
X

xo
[
Yo+
I
nuyo
A /

+n(n-l)(n-2J_}lQ+
,::\ 2
+-n(n-I)--iF-
,13 ,
'j dx
· 31 .....

j4 y'dx =
7[f(xo,Yo) + n4f(xo,~o) + n(n-1) .12f(xo,Yo)
XQ X
0
3 2f dx
+n(n-I)(n -2) .1 f(xo,Yo)
31 + .....
Using x = Xo + nh, dx - hd
- n, we get

= hj fo + n4fo + n(n- I) .12fo

2
0 3
+n(n - I)(n _ 2) .1 fo ' dx
31 + ..... .

2
nfo + -n A~o
2
f' + l
( n3
- - - -
n 2) 11,
2f'
4

l (n4 2 3)· 2 JO
+-
6 -4- n 3 +n 2 113 fo + ..... .
0
 Numc1ical Solution if
o OrcU11nrr1 lJtfJ ~·, ~
,,,,,.-- . . u , nlial l!.</1wttons L223

Y,1 == Yo +h
[
4J + 8
o
·
3
20
l.Vo +--t/· .t0 +~L'\::i
3 JO r j
+.!il\''
45 Jro-....
we neglect fourth and higher d
,-,tlues of differences or er differences and also use the
l1fo == (E - 1) fr == r _ r
1F O J1 Jo
!1:;o == (E - 1)2 fr == f.: _ 2 r r
3 0 2 U1 + Jo
l11o == (E - 1)3 fr - F 3tr
o - J3 - 'J2 + 3fi1 - ro
We get, J,

i
Y4 = Yo + h [ f1 - : f2 + f3] i
4
Y4 = Yo+ h[2f1 - f2 +2J3 ] ... (24)
3
Equation (24) is called the predictor formula which provides the
value of y 4 . To find a better approximation of y 4 , we find the first
a~proximation to /4 = f (Xo + 4h, yJ. Using this value of L in
Sunpson's rule we obtain a better approximation to y4 as

... (25)

Equation (25) is the corrector formula.

Then again an improved value of Ji is calculated and then again
corrector is applied to obtain further appropriate approximation to
Y4, The process continues till we get two consecutive values of Y.1
that SUitably agree.
After getting the values of y4 and Ji to the desired degree of
accuracy, we calculate Ys from predictor Eq. (24) as
 -1 Sh. !}5), Using
and then tl11d •f,,5 = f(.x-v U f corrector Eq. (
251
\IJ
obtain a better approxima on o Us as
h
Ys = Y3 +
3
(J3 + 4 f4 + f s)

\Ve repeat the process till y 5 becomes stationary and th

proceed to find fs and so on. The predictor and corrector fi en W
can be written in general form as orrnu1a
4
Yn+l = Yn-3 + 3 h[2fn-2 - fn-1 +2fn]
(Predictor
h
Yn+l = Yn-1 + 3 (f_n-1 + 4fn + fn+d (Corrector
Example 10.12 Find y (2) using Milne's method if y (.xj is
d I the
solution of ,}; = 2 {x + y) given y (0) = 2, y (0.5)= 2. 636. y (I) ,
3.595, y (1.5), = 4.968.
Solution: Here
Xo = 0, X1 = 0.5, Xi = l.0, ~ ::: 1.5, x
4
:: 2.0
Yo = 2, Y1 =2.636, !h, =3.595, y3 = 4.968, h == 0.5
Since dy I
dx = f(x, y) = y' = (x + y)
2
Therefore,
fa = f(Xo, YJ = l, Ji =f (x1, Y1) = 1.5680.
/2 = f (Xi, Y2) = 2.2975, fs = f(Xs , Y3) = 3.2340
Using ~ e Predictor formula
4
Yn+l = Yn-3 + h [2fn-2 - fn-1 + 2fn]
3

We get 4
Y4 = Yo+ h[2f1 - f2 + 2J3]
3
4 5
= 2+ (0. ) (2(1.5680)- 2.2975 ;+- 2(3, 2340)1
3
\\
= 6.8710 . l )
Taking X4 = =
2, y 4 6.8710, we get 1t = f(.\p y~
= 4.4355
-, Numerical Solution of Ordin
ary DUferenttol Eq11r11irms
. 0 Milne Corrector formula 225
tJ~tll~

h
Yn+l = Yn- 1 + (fn- i + 4f. + f. )
3 n 11+1

,re get
= 3 595 o. 5
. + 3 (2.2975 + 4(3.2340) + 4.4355)
= 6.8732
V ing corrected value of y 4 , we get ii = f (x , _
rj s al f F4 • 4 4, Y,v - 4.4336
·•~iJlg this v ue o J 1n corrector formula, we get a bette;
•· .,;niation of y 4 as
appro~•·
h
Y4 = Y2 + 3 (f2 + 4f3 + f4)= 3.595 + 0. 5 (2 2975
3 .

+4(3.2340) + 4.4366) = 6.8733

Since two consecutive values of y4 agree to three decimal places,
we consider 6.8733 as the solution.

Example 10.13 Using Milne's method find y(4.4) if y(x) is the

dy 2-y2
solution of dx = Sx given y (4) = 1, y (4.1) = 1.0049, y (4.2),
-~=4.3 = 1.0143.
Solution: Here
,Xo = 4, XI = 4.1, X;i = 4.2, ~ = 4.3, X4 = 4.4
Yo =l,Y1 = 1.0049, Y2 = 1.0097, Y3 = 1.0143, h = 0.1
dy I 2-y2
We have dx = f(x,y) = Y = 5x
Therefore,
_ l = 0.05, J;1 =f (X1, Y1) = 0.0493,
fo - f(Xo, YO' ) - o 0452
}; = f (X;i, Y2) = 0.0467, fs = j{~, Y3 - .
Using Milne Predictor formula
4
Yn+l = Yn-3
+ 3 h [2fn-2 - fn-1 + 2 fn]
We get
 1111111

226 '''""' ncol M< /hod u lth l>rnr1m1n In MA TlAB

HO. I)[
- 1 f-
3 2(0 0 l') J) 0 0-167 ;;,

- 1.01 897
Taking
x4 - 4 .4, Y 4 = 1.01897. we~

2- tJ~ 2 (1.0 l 7
5x4 5(-1.41
- 0 .04371
Us ing Milne Corrector formula

h
Y n+l - Yn - 1 + 3 (fn -J + 4fn + Jr, 41 )

We get h
Y4 - Y2 + (f2 + 4f 3 + 14 )
3

- 1.0097 + O~ l (0.0467 + 4(0.045 +

21
- 1.01874
Using corrected value of y , we get
4

= 0.04373.
Using this value of -4 in corrector for mula, we get a
approximation of y 4 as

- 1.0097 + O. l (0.0 467 + 4(0.0452),. 0. "tJ

3
- 1.01874
Since two consecutive values of y 4 agree to five decimJI P
we consider 1.01874 as the solution.