Transcript Runge-Kutta 2nd Order Method for Solving Ordinary

Runge 2nd Order Method
Mechanical Engineering Majors
Authors: Autar Kaw, Charlie Barker
http://numericalmethods.eng.usf.edu
Transforming Numerical Methods Education for STEM
Undergraduates
7/20/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 solid steel shaft at room temperature of 27°C is needed to be
contracted so it can be shrunk-fit into a hollow hub. It is placed in
a refrigerated chamber that is maintained at −33°C. The rate of
change of temperature of the solid shaft q is given by
6 4
5 3
3 2


3
.
69

10
θ

2
.
33

10
θ

1
.
35

10
θ 
dθ
6

θ  33
 5.3310 
2

dt
  5.4210 θ  5.588

q 0   27C
Find the temperature of the steel shaft after 24 hours. Take
a step size of h = 43200 seconds.


dθ
 5.33 10 6  3.69  10 6 θ 4  2.33 10 5 θ 3  1.35 10 3 θ 2  5.42 10  2 θ  5.588 θ  33
dt


f t,θ   5.33106  3.69106 θ 4  2.33105 θ 3  1.35103 θ 2  5.42102 θ  5.588 θ  33
1
2
1
2


qi 1  qi   k1  k2 h
8
http://numericalmethods.eng.usf.edu
Step 1: For i  0, t0  0,
Solution
q0  27
4
3
6
5









3
.
69

10
27

2
.
33

10
27
6 



k1  f t0 , q o   f 0, 27   5.3310
27  33
2
  1.35103 27  5.42102 27  5.588






 0.0020893
k 2  f t0  h, q 0  k1h   f 0  43200, 27   0.002089343200  f 43200,63.278


  3.69106  63.2784  2.33105  63.2783


6






  5.3310

63
.
278

33
  1.35103  63.2782  5.4210 2  63.278  5.588






 0.0092607
1 
1
1
 0.0020893  1  0.009260743200
2 
2
2
2

 27   0.005675043200 218.16C
q1  q 0   k1  k2 h  27  
9
http://numericalmethods.eng.usf.edu
Solution Cont
Step 2: i  1, t1  43200, q1  218.16C
k1  f t1 ,q1   f 43200,218.16
4
3
6
5









3
.
69

10

218
.
16

2
.
33

10

218
.
16
6 





  5.3310

218
.
16

33
  1.35103  218.162  5.4210 2  218.16  5.588






 8.4304
k 2 f t1  h, q1  k1h   f 43200 43200,218.16   8.430443200  f 86400,364410


  3.69106 (364410) 4  2.33105 (364410)3


6




  5.3310
 364410 33
2

3

2
  1.3510  364410  5.4210  364410  5.588






 1.26381017


1 
1
1
 8.4304  1  1.26381017 43200
2 
2
2
2

 218.16   6.31901016 43200 2.72981021C
q 2  q1   k1  k2 h  218.16  

10

http://numericalmethods.eng.usf.edu
Solution Cont
The solution to this nonlinear equation at t=86400s is
q 86400  26.099C
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
86400
43200
21600
10800
5400
q 86400
|t | %
Et
−58466
58440
−2.7298×1021 2.7298×1021
−24.537
−1.5619
−25.785
−0.31368
−26.027
−0.072214
223920
1.0460×1011
5.9845
1.2019
0.27670
q 86400  26.099C (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
q (86400)
Step size,
h
Euler
Heun
Midpoint
Ralston
86400
43200
21600
10800
5400
−153.52
−464.32
−29.541
−27.795
−26.958
−58466
−2.7298×1021
−24.537
−25.785
−26.027
−774.64
−0.33691
−24.069
−25.808
−26.039
−12163
−19.776
−24.268
−25.777
−26.032
q 86400  26.099C
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
16
t %
Step
size,
h
Euler
Heun
Midpoint
Ralston
86400
43200
21600
10800
5400
448.34
737.97
14.426
7.0957
3.5755
12216011
5.7292
1.1993
0.27435
1027.6
76.360
7.0508
1.0707
0.22604
23844
42.571
6.6305
1.2135
0.25776
q (86400)  25.217C
(exact)
7.906410
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