08.04 - Runge-kutta 2nd Order Method For Ode

RUNGE-KUTTA 2ND ORDER METHOD

08.04.1

Chapter 08.04 Runge-Kutta 2nd Order Method for Ordinary Differential Equations Runge-Kutta 2nd order method is a numerical technique to solve ordinary differential equation of the form dy = f ( x, y ), y (0 ) = y 0 dx Only first order ordinary differential equations can be solved by using Runge-Kutta 2nd order method. In other sections, we have discussed how Euler and Runge-Kutta methods are used to solve higher order ordinary differential equations or coupled (simultaneous) differential equations. How does one write a first order differential equation in the above form? Example 1 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 Another example is given as follows Example 2 dy ey + x 2 y 2 = 2 Sin (3x ), y (0) = 5 dx is rewritten as dy 2 Sin(3 x) − x 2 y 2 = , y (0) = 5 dx ey In this case 2Sin(3x) − x 2 y 2 f ( x, y ) = ey Euler’s method is given by y i +1 = y i + f ( xi , y i )h where W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

(1)

RUNGE-KUTTA 2ND ORDER METHOD

08.04.2

x0 = 0 y0 = y0 h = xi +1 − xi To understand Runge-Kutta 2nd order method, we need to derive Euler’s method from Taylor series. 3 2 dy (xi +1 − xi ) + 1 d 2y (xi +1 − xi )2 + 1 d y3 (xi +1 − xi )3 + ... y i +1 = y i + dx xi , yi 2! dx x , y 3! dx x , y i

i

i

i

1 1 2 3 yi +1 = yi + f ( xi , yi ) + f ' ( xi , yi )( xi +1 − xi ) + f ' ' ( xi , yi )( xi +1 − xi ) + ... (2) 2! 3! As you can see the first two terms of the Taylor series y i +1 = yi + f ( xi , y i )h

are the Euler’s method and hence can be considered to be Runge-Kutta 1st order method. The true error in the approximation is given by f ′( xi , yi ) 2 f ′′( xi , yi ) 3 Et = h + h + ... (3) 2! 3! So how would a 2nd order method formula look like. It would include one more term of the Taylor series as follows. 1 (4) yi +1 = yi + f ( xi , yi )h + f ′( xi , yi )h 2 2! Let us take a generic example of a first ordinary differential equation dy = e − 2 x − 3 y , y (0) = 5 dx f (x, y ) = e−2 x − 3 y Now since y is a function of x, ∂f ( x, y ) ∂f ( x, y ) dy f ′( x, y ) = + (5) ∂x ∂y dx ∂ ∂ −2 x (e − 3 y ) (e −2 x − 3 y ) = (e − 2 x − 3 y ) + ∂x ∂y = −2e −2 x + ( −3)(e −2 x − 3 y )

[

]

= −5e −2 x + 9 y The 2nd order formula for the above example would be 1 y i +1 = y i + f (xi , y i )h + f ′(xi , y i )h 2 2! 1 yi +1 = yi + (e − 2 xi − 3 yi )h + (− 5e − 2 xi + 9 yi )h 2 2!

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

08.04.3

However, we already see the difficulty of having to find f ′( x, y ) in the above method. What Runge and Kutta did was write the 2nd order method as yi +1 = yi + (a1k1 + a 2 k 2 )h (6) where

k1 = f ( xi , yi ) k2 = f ( xi + p1h, yi + q11k1h )

(7)

This form allows one to take advantage of the 2nd order method without having to calculate f ′( x, y ) . So how do we find the unknowns a1 , a2 , p1 and q11 . Without proof, equating Equation (4) and (6) , gives three equations. a1 + a2 = 1 1 a2 p1 = 2 1 a2 q11 = 2 Since we have 3 equations and 4 unknowns, we can assume the value of one of the unknowns. The other three will then be determined from the three equations. Generally the value of a2 is chosen to evaluate the other three constants. The three 1 2 values generally used for a2 are , 1 and , and are known as Heun’s Method, 2 3 Midpoint method and Ralston’s method, respectively. Heun’s method 1 Here a2 = is chosen, giving 2 1 a1 = 2 p1 = 1 q11 = 1 resulting in 1 ⎞ ⎛1 yi +1 = yi + ⎜ k1 + k2 ⎟h 2 ⎠ ⎝2 where k1 = f ( x i , y i ) k 2 = f ( xi + h, y i + k1 h ) This method is graphically explained in Figure 1.

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

(8)

(9a) (9b)

RUNGE-KUTTA 2ND ORDER METHOD

08.04.4

Slope = f ( xi + h, yi + k1h )

y

Slope = f ( xi , yi )

yi+1, predicted

Average Slope =

1 [ f (xi + h, yi + k1h ) + f (xi , yi )] 2

yi

xi

xi+1

x

Figure 1. Runge-Kutta 2nd order method (Heun’s method) 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 ⎜ xi + h, yi + k1h ⎟ 2 2 ⎝ ⎠ Ralston’s method 2 Here a 2 = is chosen, giving 3 1 a1 = 3 3 p1 = 4 3 q11 = 4 resulting in

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

(10) (11a) (11b)

RUNGE-KUTTA 2ND ORDER METHOD

08.04.5

1 2 yi +1 = yi + ( k1 + k 2 )h 3 3

where

(12)

k1 = f ( xi , yi )

(13a)

3 3 ⎛ ⎞ (13b) k 2 = f ⎜ x i + h , y i + k1 h ⎟ 4 4 ⎝ ⎠ Example 3 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 × 108 ) dt Find the temperature at t = 8 minutes using Runge-Kutta 2nd order method. Assume a step size of h = 4 minutes. Solution dθ = −2.2067 × 10 −12 θ 4 − 81 × 10 8 dt f (t , θ ) = −2.2067 × 10 −12 θ 4 − 81 × 10 8 Per Heun’s method given by Equations (8) and (9) 1 ⎞ ⎛1 θ i +1 = θ i + ⎜ k1 + k 2 ⎟h 2 ⎠ ⎝2 k1 = f (ti , θ i )

(

)

(

)

k 2 = f (ti + h, θ i + k1h ) i = 0, t0 = 0, θ 0 = θ (0) = 1200 k1 = f (t 0 ,θ o )

= f (0,1200) = −2.2067 × 10 −12 (1200 4 − 81 × 10 8 ) = −4.5579 k 2 = f (t 0 + h,θ 0 + k1 h ) = f (0 + 240,1200 + (− 4.5579)240) = f (240,106.09) = −2.2067 × 10 −12 (106.09 4 − 81 × 10 8 ) = 0.017595 1 ⎞ ⎛1 θ 1 = θ 0 + ⎜ k1 + k 2 ⎟ h 2 ⎠ ⎝2

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

1 ⎛1 ⎞ = 1200 + ⎜ (− 4.5579 ) + (0.017595 )⎟240 2 ⎝2 ⎠ = 1200 + (− 2.2702)240 = 655.16K i = 1, t1 = t 0 + h = 0 + 240 = 240, θ1 = 655.16 K

08.04.6

k1 = f (t1 ,θ1 ) = f (240,655.16) = −2.2067 × 10 −12 (655.16 4 − 81 × 10 8 ) = −0.38869 k 2 = f (t1 + h, θ1 + k1 h ) = f (240 + 240,655.16 + (− 0.38869)240) = f (480,561.87 ) = −2.2067 × 10 −12 (561.87 4 − 81 × 10 8 ) = −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 θ 2 = θ (480) = 584.27 K The results from Heun’s method are compared with exact results in Figure 2. 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 is θ ( 480) = 647.57 K

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

08.04.7

Temperature, θ(K)

1200 h=120

Exact

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

Using smaller step size would increases the accuracy of the result as given in Table 1 and Figure 3 below. Table 1. Effect of step size for Heun’s method Step size, h θ (480) Et ∈t %

480 240 120 60 30

-393.87 584.27 651.35 649.91 648.21

1041.4 63.304 -3.7762 -2.3406 -0.63219

160.82 9.7756 0.58313 0.36145 0.097625

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

08.04.8

Temperature, θ(480)

800 600 400 200 0 0 -200

100

200 300 Step size, h

400

500

-400

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

In Table 2, the Euler’s method and Runge-Kutta 2nd order method results are shown as a function of step size, Table 2. Comparison of Euler and the Runge-Kutta methods θ ( 480) Step size, 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

while in Figure 4, the comparison is shown over the range of time.

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

08.04.9

1200

Temperature,

θ(K)

1100

Midpoint

1000

Ralston

900

Heun

800

Analytical

700

600

Euler

500 0

100

200

300

400

500

600

Time, t (sec)

Figure 4. Comparison of Euler and Runge Kutta methods with exact results over time.

How do these three methods compare with results obtained if we found f ′( x, y ) directly? Of course, we know that since we are including first three terms in the series, if the solution is a polynomial of order two or less (that is, quadratic, linear or constant), any of the three methods are exact. But for any other case the results will be different. Let us take the example of dy = e − 2 x − 3 y , y (0 ) = 5 . dx If we directly find the f ′( x, y ) , the first three terms of Taylor series gives 1 yi +1 = yi + f ( xi , yi )h + f ′( xi , yi )h 2 2! where f (x, y ) = e−2 x − 3 y

f ′(x, y ) = −5e−2 x + 9 y For a step size of h = 0.2 , using Heun’s method, we find y (0.6) = 1.0930

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc

RUNGE-KUTTA 2ND ORDER METHOD

08.04.10

The exact solution y ( x ) = e−2 x + 4e−3 x gives y (0.6) = e−2(0.6 ) + 4e−3(0.6 ) = 0.96239 Then the absolute relative true error is 0.96239 − 1.0930 ∈t = × 100 0.96239 = 13.571% For the same problem, the results from the Euler and the three Runge-Kutta method are given below Table 3. Comparison of Euler’s and Runge-Kutta 2nd order methods y(0.6) Exact Euler Direct 2nd Heun Midpoint Ralston Value 0.96239 0.4955 1.0930 1.1012 1.0974 1.0994 48.514 13.571 14.423 14.029 14.236 ∈t %

W:\mws\gen\08ode\mws_gen_ode_txt_runge2nd.doc