Numerical Solution of ODE

MD. MEHEDI HASAN

Lecturer(Mathematics)
Department of General Education,
BDU

1
 Numerical Solution of ODE
 A first order initial value problem of ODE may be written in
the form
y '(t)  f ( y,t), y(a)  y0 , a  t  b
 Example:
y'(t)  3y  5, y(0)  1
y'(t)  ty 1, y(0)  0
 Numerical methods for ordinary differential equations
calculate solution on the points,tn  tn1  h where h is the steps
size
How to write Ordinary Differential Equation
How does one write a first order differential equation in
the form of
dy
 f x, y 
dx

Example
 2 y  1.3e x , y 0   5
dy is rewritten as
dx

 1.3e x  2 y, y 0   5
dy
dx
In this case

f x, y   1.3e x  2 y
 Numerical Methods for ODE
Euler Methods
 Euler Methods
 Modified Euler Method
Higher order Taylor’s Methods
 Second Order
 Fourth Order
Runge-Kutta Methods
 Second Order
 Fourth Order

4
 Euler’s Method
y

 f x, y , y 0   y 0
dy
dx
True value

Rise

y1,
Slope x0,y0 Φ
Predicted
Run value

y1  y0
 Step size, h

x1  x0 x

 f x0 , y 0 
Figure 1 Graphical interpretation of the first step of Euler’s method

y1  y 0  f x0 , y 0 x1  x0 
 y0  f x0 , y0 h
6
 Euler’s Method

yi1  yi  f xi , yi h True Value

h  xi1  xi yi+1, Predicted value

Φ
yi
h
Step size
x
xi xi+1

Figure 2. General graphical interpretation of Euler’s method

 Euler Method
 Consider the forward difference approximation for first
derivative yn1  y n
yn ' , t(i) t(i-1)+(i-1)h i= 1,2,3…n
h
 Rewriting the above equation we have
yn1  yn  hyn ', yn '  f ( yn ,t n )
 So, yn is recursively calculated as 1
y1  y0  hy0 '  y0  h f ( y0 ,t0 )
y2  y1  h f ( y1 ,t1 )
⁝
yn  yn1  h f ( yn1 ,t n1 )
 Example: solve
y'  ty 1, y0  y(0)  1, 0  t  1, h  0.25

Solution:
for t0  0, y0  y(0)  1
for t1  0.25, y1  y0  hy0 '
 y0  h(t0 y0 1)
1 0.25(0*11)  1.25
for t2  0.5, y2  y1  hy1 '
 y1  h(t1 y1 1)
1.25  0.25(0.25*1.25 1)  1.5781
etc
Graphthesolution
 Modified Euler’s Method


yn1  yn   f (xn , yn )  f (xn1 , yn1 )
h *
(3)
2

where *
yn1  yn  hf (xn , yn ) (4)

is commonly known as the Improved Euler’s

method or Modified Euler Method.
In general, the improved Euler’s method is an
example of predictor-corrector method.
Modified Euler method: Heun’s method
Modified Euler’s Method
 Why a modification? C
yn1
f n  f n1
P

2
error fPn
modify
yn1
P
f n1

tn tn1 tn tn1

Predictor
P
yn1  yn  hf n
yn'  y n1
'
f (t n , y n )  f (t n1 , y n1
P
)
Average slope y 
'

2 2
h
Corrector C
yn1  y n  hy '  y n  [ f (t n , y n )  f (t n1 , yn1 )]
2
 Example 2 Improved

Use the improved Euler’s method to obtain the

approximate value y(1.5) for the solution of y' 2xy, y(1)  1.
Compare the results for h = 0.1 and h = 0.05.
Solution:
With x0 = 1, y0 = 1, f(xn, yn) = 2xnyn, n = 0, and h = 0.1
y1* = y0 + (0.1)(2x0y0) = 1.2
Using (3) with x1 = 1 + h = 1.1
2x0 y0  2x1 y1* 2(1)(1)  2(1.1)(1.2)
y1  y0  (0.1)  1 (0.1)  1.232
2 2

The results are given in Table 6.1.3 and 6.1.4.

Consider the initial value problem.

y  y  x , y(0)  0.5

Find the value of y(0.8).

Again, let h  0.1

f ( x0 , y0 )  y0  x0  0.5  0  0.5

y  0.5  0.1 0.5  0.55

*
1

f ( x , y )  y  x  0.55  0.1  0.45

* *
1 1 1 1

* 
1
 
y1  y0  h  f x0 , y0  f x1 , y1  
2 
 * 
1
  
y1  y0  h  f x0 , y0  f x1 , y1 
2 
 
y1  0.5  0.1 0.5  0.45  0.5475
1
2 
f ( x1 , y1 )  y1  x1  0.5475  0.1
 0.4475

y  0.5475  0.1 0.4475

*
2

 0.59225
f ( x , y )  y  x  0.59225  0.2
* *
2 2 2 2

 0.39225
 * 
1
  
y2  y1  h  f x1 , y1  f x2 , y2 
2 

y2  0.5475  0.1 0.4475  0.39225
1 
2 
 0.589487
 Find y(0.1) and y(0.2), if y\ = y – x , y(0)= 2, take h = 0.1, three
decimal places are required using Modified Euler Formula

To get y1 = y(0.1)
y10  y 0  hf(x 0 , y 0 )  2  0.1(2  0)  2.2
h 0.1
y11  y 0  [f(x 0 , y 0 )  f(x 1 , y10 )]  2  [(2  0)  (2.2  0.1)]  2.205
2 2
h 0.1
y1  y 0  [f(x 0 , y 0 )  f(x 1 , y1 )]  2 
2 1
[(2  0)  (2.205  0.1)]  2.20525
2 2
Then y(0.1) = 2.205
To get y2 = y(0.2)
y12  y1  hf(x 1 , y1 )  2.205  0.1(2.205  0.1)  2.416
h 0.1
y12  y1  [f(x 1 , y1 )  f(x 2 , y 02 )]  2.205  [(2.205  0.1)  (2.416  0.2)]  2.42105
2 2
h 0.1
y 22  y1  [f(x 1 , y1 )  f(x 2 , y12 )]  2.205  [(2.205  0.1)  (2.421 0.2)]  2.4213
2 2
Then y(0.2) = 2.241
Graphthesolution
 Example
A ball at 1200K is allowed to cool down in air at an ambient
temperature of 300K. Assuming heat is lost only due to radiation,
the differential equation for the temperature of the ball is given by
d
dt
 
 2.2067 10 12  4  8110 8 ,  0   1200K
Find the temperature at t  480 seconds using Euler’s method
Assume a step size of h  240 seconds.
Solution

Step 1:
d
dt

 2.2067 10 12  4  8110 8 

f t,   2.2067 10 12  4  8110 8 
i 1  i  f ti , i h
1   0  f t0 ,  0 h
 1200  f 0,1200 240
 1200   2.2067 10 12 1200 4  81108 240
 1200   4.5579 240
 106.09K
1 is the approximate temperature at t  t1  t0  h  0  240  240
 240 1  106.09K
 Solution Cont

Step 2: For i  1, t1  240, 1  106.09

 2  1  f t1 , 1 h
 106.09  f 240,106.09 240
 
 106.09   2.2067 10 12 106.09 4  8110 8 240 
 106.09  0.017595 240
 110.32K

 2 is the approximate temperature at t  t2  t1  h  240  240  480

 480  2  110.32K

34
 Solution Cont
The exact solution of the ordinary differential equation is
given by the solution of a non-linear equation as

0.92593ln
  300 1.8519 tan1 0.00333   0.22067 103 t  2.9282
  300

The solution to this nonlinear equation at t=480 seconds is

 (480)  647.57K
Comparison of Exact and Numerical Solutions

1400

1200
Temperature, θ(K
)

1000
E xa c t S ol u tio n
800

600

400
h= 240
200

0
0 100 200 300 400 500
Time, t(sec)

Figure 3. Comparing exact and Euler’s method

 Effect of step size

Table 1. Temperature at 480 seconds as a function of step size, h

Step, h (480) Et |єt|%

480 −987.81 1635.4 252.54
240 110.32 537.26 82.964
120 546.77 100.80 15.566
60 614.97 32.607 5.0352
30 632.77 14.806 2.2864

 (480)  647.57K (exact)

 Comparison with exact results

1500

1000 Exact solution

500
Temperature,

h=120
h=240
0
θ(K)

0 100 200 300 400 500

-500
Time, t (sec) h=480
-1000

-1500

Figure 4. Comparison of Euler’s method with exact solution for different step sizes
Effects of step size on Euler’s Method

800

400
Temperature,θ(K)

0
0 100 200 300 400 500
-400
Step size, h (s)
-800

-1200

Figure 5. Effect of step size in Euler’s method.

39
 xtt3|%
|h
E

step size,

3 362.5 373.22 3483.0

1.5 720.31 709.60 6622.2
0.75 284.65 273.93 2556.5
0.375 10.718 0.0024912 0.023249
0.1875 10.714 0.0010803 0.010082
 Algorithm for Euler’s Method
Given
f(x,y) (* expression for x and y *)
x0 (* initial value for x *)
y0 (* initial value for y *)
xn (* terminal value for x *)
n (* number of partitions of the interval *)
deltax = (xn-x0)/n
xi = x0
xprev = xi
yi = y0
for(i=1, in, i++,
xi = xi + deltax
yi = yi + f(xprev,yi)*deltax)
xprev=xi
Return yi
 It is a numerical methods for solving ordinary differential
equations (O.D.Es) with a given conditions.
y\ = f(x, y), y(x0) = y0.
Taylor series expansion of the function f(x) around a point of expansion x
has the following form: (x x 0 ) 2
(x x 0 ) 
y(x)  y(x 0 )  y(x 0 ) y (x 0 )  .......
1! 2!
Assume h = x – x0
h h2 h3
Hint:  y(x 0  h)  y(x 0 )  y(x 0 )  y (x 0 )  y(x 0 )  .......
1! 2! 3!
1- When the step h becomes small, the numerical solution becomes nearly
the same as the exact solution.
2- y = f(x).
3- computing approximate values of the solution f(x) at the points:
x1 = x0 +h, x2 = x0 +2h, x3 = x0 + 3h,…………….
y 1  y(x 0  h), y 2  y(x1  h), y 3  y(x 2  h),..........
 Solve the following boundary value problem y\ = - y at x = 0.2, 0.4 given
that y(0) = 1

Taylor expansion given

2
by: 3
h h h
y(x 0  h)  y(x 0 )  y(x 0 )  y (x 0 )  y(x 0 )  .......
1! 2! 3!
Consider the point of expansion zero
 h  x  x0  0.2  0  0.2
Evaluate the derivatives at the point of expansion zero
y(x)  y  y(x 0 )  y(x 0 )  1,
y (x)   y(x)  y (x 0 )   y(x 0 )  1,
y(x)   y (x)  y(x 0 )   y (x 0 )  1,
Then 0.2 (0.2)2 (0.2)3
y(0.2)  y(0  0.2)  1 (1)  (1)  (1)  .......  0.81867
1! 2! 3!
0.2 (0.2)2 (0.2)3
y(0.4)  y(0.2  0.2)  0.81867  (0.81867)  (0.81867)  (0.81867)  .......  0.67022
1! 2! 3!
 Runge-Kutta Method
 The Taylor methods outlined in the previous section have the
desirable property of high-order local truncation error,

 but the disadvantage of requiring the computation and

evaluation of the derivatives of f (t, y).

 This is a complicated and time-consuming procedure for most

problems, so the Taylor methods are seldom used in practice.

 Runge-Kutta methods have the high-order local truncation

error of the Taylor methods but eliminate the need to compute
and evaluate the derivatives of f (t, y).
 Second Order Runge-Kutta Method
 Consider the IVP

dy(x)
 f (x, y), y(x0 )  y0 (1)
dx
 Suppose the increment of y as the weighted average of two
estimates of the increment, say K1 and K2.
i.e. yi1  yi  y  yi  w1K1  w2 K2 (2)
K1  hf (xi , yi )
K2  hf (xi   h, yi  K1 )
where  ,  , w1, w2 are to be determined
 Second Order Runge-Kutta Method

Second Order Runge-Kutta method becomes:

yi1  yi  K1  K 2 
1
2
where
K1  hf (xi , yi )
K 2  hf (xi  h, yi  K1 )
 2nd Order Runge-Kutta
RK2

Typical value of   1, Know as RK2

Equivalent to Heun's method with a single corrector
k1  f ( xi , yi )
k2  f ( xi  h, yi  k1 h)

yi1  yi  k1  k2 
h
2
Local error is O(h 3 ) and global error is O(h 2 )
 Second order Runge-Kutta Method
Example
Solve the following system to find x(1.02) using RK2
x(t)  1  x 2  t 3 , x(1)  4, h  0.01,   1

STEP1:
K1  h f (t0  1, x0  4)  0.01(1  x02  t03 )  0.18
K 2  h f (t0  h, x0  K1 )
 0.01(1  (x0  0.18)2  (t0  .01)3 )  0.1662
x(1  0.01)  x(1)  K1  K 2 / 2
 4  (0.18  0.1662) / 2  3.8269
 Second order Runge-Kutta Method
Example
STEP 2

K1  hf (t 1  1.01, x1  3.8269)  0.01(1 x1

2
 t1 )  0.1668
3

K 2  hf (t1  h, x1  K1 )
 0.01(1 (x1  0.1668)2  (t1  .01)3 )  0.1546

x(1.01 0.01)  x(1.01) 

1
K1  K2 
2
1
 3.8269  (0.1668 0.1546)  3.6662
2
 x(t) 1  x (t)  t , x(1)  4,
2 3

Solution for t [1,2]

Using RK2,   1

59
CISE301D_Tr.opiLitcan8L4Ku&5mar Saha, Department of Applied Mathematics
 Higher-Order Runge-Kutta

Higher order Runge-Kutta methods are available.

Derived similar to second-order Runge-Kutta.

Higher order methods are more accurate but

require more calculations.

60
 4th Order R unge-K utta Method
 Consider the IVP

dy(x)
 f (x, y), y(x0 )  y0 (1)
dx
 The Runge-Kutta 4th order method is based on the following
yi1  yi  a1k1  a2 k 2  a3 k3  a4 k 4 h
 where knowing the value of y  yi at xi , we can find the value
of y  yi1 at xi1 , and h  xi1  xi

 Equation (1) is equated to the first five terms of Taylor series

 4th Order R unge-K utta Method
1 ' 1 1
y i1  y i  f x i , y i h f xi , y i  h  f
2 ''
xi , y i  h  f
3 '''
xi , y i  h4
2! 3! 4!
....

 Deriving in similar fashion of the Runge-Kutta method of order

2, we will get eleven equations and thirteen unknowns. One of
the popular solutions used for Runge-Kutta method of order 4 is
 4th Order R unge-K utta Method

 y i  k1  2k 2  2k 3  k 4 h
1
y i1
6
 where

k1  f xi , y i 

k2  f  xi  h, yi  k1h 
1 1
 2 2 
k 3  f  x i  h, y i  k 2 h 
1 1
 2 2 
k 4  f xi  h, y i  k 3 h 
 Example
A ball at 1200K is allowed to cool down in air at an ambient
temperature of 300K. Assuming heat is lost only due to radiation,
the differential equation for the temperature of the ball is given by
d
dt
 
 2.2067 10 12  4  8110 8 ,  0   1200K
Find the temperature at t  480 seconds using R-K method
Assume a step size of h  240 seconds.

d
dt

 2.2067 10 12  4  8110 8 

f t,  2.206710 12  4  8110 8 
 Solution

  i  k1  2k 2  2k3  k4 h
1
i 1
6
Step 1: i  0, t0  0, 0   (0)  1200

 
k1  f t0 , o   f 0,1200   2.2067 10 12 1200 4  81108  4.5579

k 2  f  t0  h,  0  k1h   f  0  240 ,1200   4.5579 240 

1 1 1 1
 2 2   2 2 
 
 f 120,653.05 2.2067 10 12 653.05 4  81108  0.38347

k3  f  t 0  h,  0  k 2 h   f  0  240 ,1200   0.38347 240 

1 1 1 1
 2 2   2 2 
 f 120,1154.0   2.2067 10 12 1154.0 4  81108  3.8954

k 4  f t0  h, 0  k3h   f 0  240 ,1200   3.984 240

 f 240,265.10  2.206710 12 265.10 4  81108  0.0069750
 Solution Cont

  
1 0
1
k1  2k2  2k 3  k 4 h
6
 1200 
1
 4.5579  2 0.38347 23.8954 0.069750240
6
 1200   2.1848240
1
6
 675.65K

1 is the approximate temperature at

t  t1  t0  h  0  240  240
 Solution Cont

Step 2: i  1,t1  240,1  675.65K

 
k1  f t1 , 1   f 240,675.65  2.2067 10 12 675.65 4  81108  0.44199

k 2  f  t1  h,1  k1h   f  240  240,675.65  0.44199240

1 1 1 1
 2 2   2 2 
 
 f 360,622.61 2.2067 10 12 622.614  81108  0.31372

k3  f t1  h,1  k 2h   f  240  240,675.65  0.31372240

1 1 1 1
 2 2   2 2 
 f 360, 638.00   2.2067 10 12 638.00 4  81108  0.34775
k 4  f t1  h, 1  k3 h   f 240  240 ,675.65   0.34775240 
 f 480,592.19   2.2067 10 12 592.19 4  81108  0.25351
 Solution Cont

1
 2  1  k1  2k2  2k3  k4 h
6
 675.65  0.44199  2 0.31372 2 0.34775 0.25351240
1
6
 675.65   2.0184 240
1
6
 594.91K
 is the approximate temperature at

t2  t1  h  240  240  480

 480   2  594.91K
 Comparison with exact results
Temperature, θ(K) 1600

1200

h= 120
800 E xa c t
h= 240
400

h=480
0
0 200 400 600
-400
Time,t(sec)

Figure 1. Comparison of Runge-Kutta 4th order method with exact solution

 Effect of step size

Table 1. Temperature at 480 seconds as a function of step size, h

Step size, h  (480) Et |єt|%

480 −90.278 737.85 113.94

240 594.91 52.660 8.1319
120 646.16 1.4122 0.21807
60 647.54 0.033626 0.0051926
30 647.57 0.00086900 0.00013419
 (480)  647.57K (exact)
 Comparison of Runge-Kutta Methods

1400
θ(K)

1200
4th order
1000
Temperature,

800
Exact
600
Heun
400

200 Euler

0
0 100 200 300 400 500
Time, t(sec)

Figure 2. Comparison of Runge-Kutta methods of 1st, 2nd, and 4th order.

 6.4 Higher-Order Equations and Systems

Second-Order IVPs
An IVP
y  f (x, y, y) , y(x0 )  y0 , y(x0 )  u0 (1)
can be expressed by
y  u
u  f (x, y, u) (2)
Since y’(x0) = u0, then y(x0) = y0, u(x0) = u0.
Apply the Euler’s method to (2)
yn1  yn  hun
un1  un  hf (xn , yn , un ) (3)
whereas the RK4 method is applied
h
yn1  yn  (m1  2m2  2m3  m4 )
6 (4)
un1  u n  (k1  2k 2  2k3  k 4 )
h
6
where
m1  un k1  f (xn , yn , un )
m2  un  12 hk1 k2  f (xn  12 h, y n  12 hm1 , u n  12 hk1 )
m3  un  12 hk 2 k3  f (xn  12 h, y n  12 hm2 , u n  12 hk 2 )
m4  un  hk3 k 4  f (xn  h, yn  hm3 , un  hk3 )
In general,
y (n)  f (x, y, y', … , y (n1) )
 Example 1 Euler’s Method

Use the Euler’s method to obtain y(0.2), where

y  xy  y  0 , y(0)  1 , y(0)  2 (5)
Solution:
Let y’ = u, then (5) becomes
y  u
u  xu  y
From (3)
yn1  yn  hun
un1  un  h[xnun  yn ]
 Example 1 (2)

Using h = 0.1, y0 = 1, u0 = 2, we find

y1  y0  (0.1)u0  1 (0.1)2  1.2
u1  u0  (0.1)[x0u0  y0 ]  2  (0.1)[(0)(2) 1]  1.9
y2  y1  (0.1)u1  1.2  (0.1)(1.9)  1.3
u2  u1  (0.1)[x1u1  y1 ]  1.9  (0.1)[(0.1)(1.9) 1.2]  1.761
In the words, y(0.2) ≈ 1.39 and y’(0.2) ≈ 1.761.
Fig 6.4.1 shows the comparison of results between by
Euler’s method and by the RK4 method.
 1
C1  C 0  (k1  2k 2  2k 3  k 4 )h  8.1059 10 6 parts/m 3
6
C2  6.5705106 parts/m 3

 3t 

C(t)  110 7 e  50 
C(7)  6.5705 10 6 parts/m 3
 Numerical Solution of a System

The solution of a system of the form

dx1
 f1 (t, x1 , x2 , …, xn )
dt
dx2
 f 2 (t, x1 , x2 , …, xn )
dt
⁝ ⁝
dxn
 f n (t, x1 , x2 ,…, xn )
dt

can be approximated by numerical methods.

For example, by the RK4 method:
x  f (t, x, y)
y  g(t, x, y)
x(t0 )  x0 , y(t0 )  y 0
(6)

looks like this:

h
xn1  xn  xn  (m1  2m2  2m3  m4 )
6
h (7)
yn1  yn  (k1  2k 2  2k3  k 4 )
6
 where
m1  f (t n , xn , yn ) k1  g(t n , xn , yn )
m2  f (tn  12 h, xn  12 m1 , yn  12 k1 ) k2  g(tn  12 h, xn  12 m1 , yn  12 k1 )
m3  f (tn  12 h, xn  12 m2 , y n  12 k 2 ) k3  g(tn  12 h, xn  12 m2 , yn  12 k 2 )
m4  f (tn  h, xn  m3 , yn  k3 ) k4  g(tn  h, xn  m3 , yn  k3 )
(8)
 Example 3 RK4 Method

Consider x  2x  4 y
y  x  6 y
x(0)  1, y(0)  6
Use the RK4 method to approximate x(0.6) and y(0.6)
with h = 0.2 and h = 0.1.
Solution:
With h = 0.2 and the given data, from (8)
 Example 3 (2)

m1  f (t0 , x0 , y0 )  f (0, 1, 6)  2(1)  4(6)  22

k1  g(t0 , x0 , y0 )  g(0, 1, 6)  1(1)  6(6)  37
m2  f (t0  12 h, x 0  12 hm1 , y0  12 hk1 )  f (0.1,1.2, 9.7)  41.2
k2  g(t0  12 h, x0  12 hm1 , y0  12 hk1 )  g(0.1,1.2, 9.7)  57
m3  f (t0  12 h, x 0  12 hm2 , y0  12 hk 2 )  f (0.1, 3.12, 11.7)  53.04
k3  g(t0  12 h, x 0  12 hm2 , y 0  12 hk 2 )  g(0.1, 3.12, 11.7)  67.08
m4  f (t0  h, x0  hm3 , y0  hk3 )  f (0.2, 9.608,19.416)  96.88
k4  g(t0  h, x0  hm3 , y0  hk3 )  g(0.2, 9.608,19.416)  106.888
 Example 3 (3)

Therefore, from (7) we get

0.2
x1  x0  (m1  2m2  2m3  m4 )
6
0.2
 1 (2.2  2(41.2)  2(53.04)  96.88)  9.2453
6
0.2
y1  y 0  (k1  2k 2  2k3  k 4 )
6
0.2
 6 (37  2(57)  2(67.08) 106.888)  19.0683
6
See Fig 6.4.2 and Table 6.4.1, 6.4.2.