Transcript Runge-Kutta 2nd Order Method for Solving Ordinary

Runge 2nd Order Method
Civil Engineering Majors
Authors: Autar Kaw, Charlie Barker
http://numericalmethods.eng.usf.edu
Transforming Numerical Methods Education for STEM
Undergraduates
7/21/2015
http://numericalmethods.eng.usf.edu
1
Runge-Kutta 2nd Order Method
http://numericalmethods.eng.usf.edu
Runge-Kutta 2nd Order Method
For
dy
 f ( x, y ), y (0)  y0
dx
Runge Kutta 2nd order method is given by
yi 1  yi  a1k1  a2 k2 h
where
k1  f xi , yi 
k2  f xi  p1h, yi  q11k1h
3
http://numericalmethods.eng.usf.edu
Heun’s Method
Heun’s method
Slope f xi  h, yi  k1h 
y
Here a2=1/2 is chosen
1
a1 
2
yi+1, predicted
Slope f xi , yi 
p1  1
q11  1
resulting in
1 
1
yi 1  yi   k1  k2 h
2 
2
where
k1  f xi , yi 
Average Slope 
yi
xi
1
 f xi  h, yi  k1h   f xi , yi 
2
xi+1
x
Figure 1 Runge-Kutta 2nd order method (Heun’s method)
k 2  f xi  h, yi  k1h
4
http://numericalmethods.eng.usf.edu
Midpoint Method
Here a2  1 is chosen, giving
a1  0
p1 
1
2
q11 
1
2
resulting in
yi 1  yi  k2h
where
k1  f xi , yi 
1
1


k 2  f  xi  h, yi  k1h 
2
2


5
http://numericalmethods.eng.usf.edu
Ralston’s Method
Here a 2  2 is chosen, giving
3
1
3
3
p1 
4
3
q11 
4
resulting in
a1 
2 
1
yi 1  yi   k1  k2 h
3 
3
where
k1  f xi , yi 
3
3


k 2  f  xi  h, yi  k1h 
4
4


6
http://numericalmethods.eng.usf.edu
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
7
http://numericalmethods.eng.usf.edu
Example
A polluted lake with an initial concentration of a bacteria is 107
parts/m3, while the acceptable level is only 5×106 parts/m3. The
concentration of the bacteria will reduce as fresh water enters the
lake. The differential equation that governs the concentration C of
the pollutant as a function of time (in weeks) is given by
dC
 0.06C  0, C (0)  10 7
dt
Find the concentration of the pollutant after 7 weeks. Take a
step size of 3.5 weeks. dC
dt
 0.06C
f t , C   0.06C
8
1 
1
Ci 1  Ci   k1  k 2 h
2 
2
http://numericalmethods.eng.usf.edu
Solution
7
Step 1: i  0, t0  0, C0  10


 
k1  f t0 , C0   f 0,107  0.06 107  600000

k2  f t0  h, C0  k1h   f 0  3.5, 107   6000003.5





 f 3.5, 7.9 106  0.06 7.9 106  474000
1 
1
C1  C0   k1  k 2 h
2 
2
1
1

 107    600000   4740003.5
2
2

 107   5370003.5
C1 is the approximate
concentration of bacteria at
t  t1  t0  h  0  3.5  3.5 weeks
C3.5  C1  8.1205106 parts/m3
 8.1205106 parts/m3
9
http://numericalmethods.eng.usf.edu
Solution Cont
6
3
Step 2: i  1, t1  t0  h  0  3.5  3.5, C1  8.120510 parts/m




k1  f t1, C1   f 3.5, 8.1205106  0.06 8.1205106  487230


k 2 f t1  h, C1  k1h  f 3.5  3.5, 8.1205106   4872303.5
 f 7, 6415200  0.066415200  384910
1 
1
C2  C1   k1  k 2 h
2 
2
1
1

 8.1205 106    487230   3849103.5
2
2

 8.1205 106   4360703.5
 6.5943 106 part s/m3
C2 is the approximate
concentration of bacteria at
t  t2  t1  h  3.5  3.5  7 weeks
C7  C2  6.5943106 parts/m3
10
http://numericalmethods.eng.usf.edu
Solution Cont
The exact solution of the ordinary differential equation is
given by the solution of a non-linear equation as
 3t 


7  50 
C(t )  110 e
The solution to this nonlinear equation at t=7 weeks is
C(7)  6.5705106 parts/m3
11
http://numericalmethods.eng.usf.edu
Comparison with exact results
Figure 2. Heun’s method results for different step sizes
12
http://numericalmethods.eng.usf.edu
Effect of step size
Table 1. Effect of step size for Heun’s method
Step size, h
C 7 
Et
|t | %
7
3.5
1.75
0.875
0.4375
6.6520×106
6.5943×106
6.5760×106
6.5718×106
6.5708×106
−111530
−23784
−5489.1
−1318.8
−323.24
1.6975
0.36198
0.083542
0.020071
0.0049195
C(7)  6.5705106 (exact)
13
http://numericalmethods.eng.usf.edu
Effects of step size on Heun’s
Method
Figure 3. Effect of step size in Heun’s method
14
http://numericalmethods.eng.usf.edu
Comparison of Euler and RungeKutta 2nd Order Methods
Table 2. Comparison of Euler and the Runge-Kutta methods
Step size,
h
7
3.5
1.75
0.875
0.4375
C (7 )
Euler
Heun
Midpoint
Ralston
5.8000×106
6.2410×106
6.4160×106
6.4960×106
6.5340×106
6.6820×106
6.5943×106
6.5760×106
6.5718×106
6.5708×106
6.6820×106
6.5943×106
6.5760×106
6.7518×106
6.5708×106
6.6820×106
6.5943×106
6.5760×106
6.5718×106
6.5708×106
C(7)  6.5705106
15
(exact)
http://numericalmethods.eng.usf.edu
Comparison of Euler and RungeKutta 2nd Order Methods
Table 2. Comparison of Euler and the Runge-Kutta methods
Step size,
h
7
3.5
1.75
0.875
0.4375
t %
Euler
Heun
Midpoint
Ralston
11.726
5.0144
2.3447
1.1362
0.55952
1.6975
0.36198
0.083542
0.020071
0.0049195
1.6975
0.36198
0.083542
0.020071
0.0049195
1.6975
0.36198
0.083542
0.020071
0.0049195
C(7)  6.5705106 (exact)
16
http://numericalmethods.eng.usf.edu
Comparison of Euler and RungeKutta 2nd Order Methods
Figure 4. Comparison of Euler and Runge Kutta 2nd order
methods with exact results.
17
http://numericalmethods.eng.usf.edu
Additional Resources
For all resources on this topic such as digital audiovisual
lectures, primers, textbook chapters, multiple-choice
tests, worksheets in MATLAB, MATHEMATICA, MathCad
and MAPLE, blogs, related physical problems, please
visit
http://numericalmethods.eng.usf.edu/topics/runge_kutt
a_2nd_method.html
THE END
http://numericalmethods.eng.usf.edu