Ordinary Differential Equations

Euler’s Method
y

dy
 f  x, y , y 0  y 0 True value
dx
y1, Predicted
Rise value
Slope  x0,y0 Φ
Run
y1  y 0
 Step size, h
x1  x0
 f  x0 , y 0  x

y1  y 0  f  x0 , y 0  x1  x0 
Figure 1 Graphical interpretation of the first step of Euler’s method
 y0  f  x0 , y0 h

2
 Euler’s Method
y

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

3
 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
dy
 2 y  1.3e  x , y 0   5
dx
is rewritten as
dy
 1.3e  x  2 y, y 0   5
dx
In this case

f  x, y   1.3e  x  2 y
4
 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
 2.2067  10 12  4  81  10 8 ,  0   1200 K
dt

Find the temperature at t  480 seconds using Euler’s method. Assume a step size of
h  240 seconds.

5
 Solution
Step 1:
d
 2.2067  10 12  4  81  10 8 
dt
f t ,   2.2067  10 12  4  81 108 
 i 1   i  f ti ,  i h
1   0  f t 0 ,  0 h
 1200  f 0,1200 240
 1200   2.2067  10 12 1200 4  81  10 8 240
 1200   4.5579 240
 106.09 K
1 is the approximate temperature at t  t1  t0  h  0  240  240
 240  1  106.09K
6
 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 .32 K

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

 480   2  110.32K

7
 Solution Cont

The exact solution of the ordinary differential equation is given by the

solution of a non-linear equation as
  300
0.92593 ln  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.57 K

8
 Comparison of Exact and
Numerical Solutions
1400

1200
Temperature, θ(K)

1000
Exact Solution
800

600

400
h=240
200

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

Figure 3. Comparing exact and Euler’s method

9
 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.57 K (exact)

10
 Comparison with exact results
1500

1000 Exact solution

Temperature, θ(K)

500
h=120
h=240
0
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

11
 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.

12
 Errors in Euler’s Method
It can be seen that Euler’s method has large errors. This can be illustrated using
Taylor series.
dy 1 d2y 1 d3y
y i 1  y i  xi 1  xi   2
 x i 1  x i 2
 3
 x i 1  x i 3
 ...
dx xi , yi 2! dx x , y 3! dx x , y
i i i i

1 1
y i 1  y i  f ( xi , y i ) xi 1  xi   f ' ( xi , y i ) xi 1  xi   f ' ' ( xi , y i ) xi 1  xi   ...
2 3

2! 3!
As you can see the first two terms of the Taylor series

yi 1  yi  f xi , yi h are the Euler’s method- Runge Kutta 1st Order Method

The true error in the approximation is given by

f  xi , yi  2 f  xi , yi  3 Et  h 2
Et  h  h  ...
2! 3!
13
Ordinary differential Equations:
Runge-Kutta 2nd Order Method
 Runge-Kutta 2nd Order Method
dy
For  f ( x, y ), y (0)  y0
dx

Runge Kutta 2nd order method is given by

yi 1  yi  a1k1  a2 k 2 h a1+a2=1
a2p1=1/2
where
a2q11=1/2
k1  f xi , yi  3 equations 4 unknowns- assume 1 and
other 3 then can be determined
k2  f  xi  p1h, yi  q11k1h 

Generally a2 is chosen1/2, 1 and 2/3

15
 Heun’s Method
Heun’s method
y
Slope  f  xi  h, yi  k1h 
Here a2=1/2 is chosen
1 yi+1, predicted
a1 
2 Slope  f  xi , yi 
p1  1
q11  1 Average Slope 
1
 f xi  h, yi  k1h   f xi , yi 
2
yi

resulting in
1 1 
yi 1  yi   k1  k2 h
2 2  xi xi+1
x

where
Figure 1 Runge-Kutta 2nd order method (Heun’s method)
k1  f  x i , y i 
k 2  f  x i  h, y i  k 1 h 

16
 Midpoint Method
Here a2  1 is chosen, giving

a1  0
1
p1 
2
1
q11 
2
resulting in
yi 1  yi  k2h

where

k1  f xi , yi 
 1 1 
k 2  f  x i  h , y i  k1 h 
 2 2 

17
 Ralston’s Method
Here a2  2 is chosen, giving
3
1
a1 
3
3
p1 
4
3
q11 
4
resulting in
1 2 
yi 1  yi   k1  k 2 h
3 3 
where
k1  f xi , yi 
 3 3 
k 2  f  x i  h , y i  k1 h 
 4 4 
18
 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
dy
 2 y  1.3e  x , y 0   5
dx
is rewritten as
dy
 1.3e  x  2 y, y 0   5
dx
In this case

f  x, y   1.3e  x  2 y
19
 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
 2.2067  10 12  4  81  10 8 , 0  1200 K
dt

Find the temperature at t  480 seconds using Heun’s method. Assume a step size of
h  240 seconds.
d
 2.2067  10 12  4  81 10 8 
dt
f t ,   2.2067  10 12  4  81 10 8 
1 1 
 i 1   i   k1  k 2 h
2 2 
20
 Solution
Step 1: i  0, t0  0,  0   (0)  1200 K

k1  f t0 ,  o  k 2  f t0  h,  0  k1h 

 f 0,1200  f 0  240,1200   4.5579 240 
 2.2067 1012 12004  81108   f 240,106.09 
 4.5579  2.2067  10 12 106.09 4  81 108 
 0.017595
1 1 
1   0   k1  k 2 h
2 2 
1 1 
 1200    4.5579   0.017595 240
2 2 
 1200   2.2702 240
 655.16 K
21
 Solution Cont
Step 2: i  1, t1  t 0  h  0  240  240,1  655.16 K
k1  f t1 ,1  k 2 f t1  h,1  k1h 
 f 240,655.16  f 240  240,655.16   0.38869240
 f 480,561.87
 2.2067 1012 655.164  81108 
 2.2067 1012 561.87 4  81108 
 0.38869
 0.20206

1 1 
 2  1   k1  k 2 h
2 2 
1 1 
 655.16    0.38869    0.20206 240
2 2 
 655.16   0.29538 240
 584.27 K
22
 Solution Cont

The exact solution of the ordinary differential equation is given by the

solution of a non-linear equation as
  300
0.92593 ln  1.8519 tan 1 0.0033333   0.22067  10 3 t  2.9282
  300

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

 (480)  647.57 K

23
 Comparison with exact results
1200

Exact h=120
Temperature, θ(K)

800

h=240
400
h=480

0
0 100 200 300 400 500

-400
Time, t(sec)

Figure 2. Heun’s method results for different step sizes

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

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

480 −393.87 1041.4 160.82

240 584.27 63.304 9.7756
120 651.35 −3.7762 0.58313
60 649.91 −2.3406 0.36145
30 648.21 −0.63219 0.097625

 (480)  647.57 K (exact)

25
 Effects of step size on Heun’s
Method
800
Temperature, θ(480)

600

400

200

0
0 100 200 300 400 500
-200 Step size, h

-400

Figure 3. Effect of step size in Heun’s method

26
 Comparison of Euler and Runge-
Kutta 2nd Order Methods
Table 2. Comparison of Euler and the Runge-Kutta methods

Step size, (480)

h
Euler Heun Midpoint Ralston
480 −987.84 −393.87 1208.4 449.78
240 110.32 584.27 976.87 690.01
120 546.77 651.35 690.20 667.71
60 614.97 649.91 654.85 652.25
30 632.77 648.21 649.02 648.61

 (480)  647.57 K (exact)

27
 Comparison of Euler and Runge-
Kutta 2nd Order Methods

Table 2. Comparison of Euler and the Runge-Kutta methods

Step size, t %
h
Euler Heun Midpoint Ralston
480 252.54 160.82 86.612 30.544
240 82.964 9.7756 50.851 6.5537
120 15.566 0.58313 6.5823 3.1092
60 5.0352 0.36145 1.1239 0.72299
30 2.2864 0.097625 0.22353 0.15940

 (480)  647.57 K (exact)

28
 Comparison of Euler and Runge-
Kutta 2nd Order Methods
1200

1100
Temperature, θ(K)

1000 Midpoint

900 Ralston

800 Heun

700
Analytical
600
Euler
500
0 100 200 300 400 500 600

Time, t (sec)

Figure 4. Comparison of Euler and Runge Kutta 2nd order methods with
exact results.
29
Ordinary differential Equations:
Runge-Kutta 4th Order Method
 Runge-Kutta 4th Order Method
For dy
 f ( x, y ), y (0)  y0
dx
Runge Kutta 4th order method is given by
1
yi 1  yi  k1  2k2  2k3  k 4 h
6
where
k1  f  xi , yi 
 1 1 
k 2  f  xi  h, yi  k1h 
 2 2 
 1 1 
k3  f  xi  h, yi  k 2 h 
 2 2 

k4  f xi  h, yi  k3h 
31
 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
dy
 2 y  1.3e  x , y 0   5
dx
is rewritten as
dy
 1.3e  x  2 y, y 0   5
dx
In this case

f  x, y   1.3e  x  2 y
32
 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
 2.2067  10 12  4  81  10 8 , 0  1200 K
dt

Find the temperature at t  480 seconds using Runge-Kutta 4th order method.

Assume a step size of h  240 seconds.

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

1
 i 1   i  k1  2k 2  2k 3  k 4 h
6
33
 Solution
Step 1: i  0, t0  0,  0   (0)  1200

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

 1 1   1 1 
k 2  f  t0  h,  0  k1h   f  0  240,1200   4.5579 240 
 2 2   2 2 
 f 120,653.05  2.2067 10 12 653.054  81108   0.38347
 1 1   1 1 
k3  f  t0  h,  0  k 2 h   f  0  240,1200   0.38347 240 
 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  k3 h   f 0  240,1200   3.984240

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

34
 Solution Cont
1
1   0  k1  2k2  2k3  k4 h
6
1
 1200   4.5579  2 0.38347   2 3.8954  0.069750240
6
1
 1200   2.1848240
6
 675.65K

1 is the approximate temperature at

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

 240  1  675.65K

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

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

 1 1   1 1 
k 2  f  t1  h, 1  k1h   f  240  240 , 675.65   0.44199 240 
 2 2   2 2 
 f 360,622.61  2.2067  10 12 622.614  81108   0.31372
 1 1   1 1 
k3  f  t1  h, 1  k 2 h   f  240  240 ,675.65   0.31372240 
 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  k3h   f 240  240,675.65   0.34775240 

 f 480,592.19  2.2067 10 12 592.19 4  81108   0.25351

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

t2  t1  h  240  240  480

 480   2  594.91K

37
 Solution Cont

The exact solution of the ordinary differential equation is given by the

solution of a non-linear equation as
  300
0.92593 ln  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.57 K

38
 Comparison with exact results
1600
Temperature, θ(K)

1200

h=120
800 Exact
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

39
 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.57 K (exact)
40
 Effects of step size on Runge-
Kutta 4th Order Method
800
Temperature, θ(480)

600

400

200

0
0 100 200 300 400 500
-200
Step size, h

Figure 2. Effect of step size in Runge-Kutta 4th order method

41
 Comparison of Euler and Runge-
Kutta Methods
1400
Temperature, θ(K)

1200
4th order
1000

800
Exact
600
Heun
400

200 Euler

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

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

42