08.03 - Runge-kutta 4th Order Method For Ode

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

-1–

Chapter 8-03 Runge-Kutta 4th Order Method for Ordinary Differential Equations Runge-Kutta 4th order method is a numerical technique to solve ordinary differential equation of the form dy = f ( x, y ), y (0) = y 0 dx So only first order ordinary differential equations can be solved by using Runge-Kutta 4th 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(3x ) − x 2 y 2 = , y (0) = 5 dx ey In this case 2 Sin(3x ) − x 2 y 2 f ( x, y ) = ey th Runge-Kutta 4 order method is based on the following

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

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

-2–

yi +1 = yi + (a1k1 + a2 k2 + a3k3 + a4 k4 )h (1) 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 3 dy 1 d2y 1 d3y 2 ( ) ( ) (xi +1 − xi ) yi +1 = yi + x − x + x − x + xi , y i i +1 i i +1 i 2 xi , yi 3 xi , y i 2! dx 3! dx dx (2) 4 1d y 4 + x , y ( xi +1 − xi ) 4! dx 4 i i dy = f ( x, y ) and xi +1 − xi = h Knowing that dx 1 1 1 (3) yi +1 = yi + f ( xi , yi )h + f ' ( xi , yi )h 2 + f '' ( xi , yi )h3 + f ''' (xi , yi )h 4 2! 3! 4! Based on equating Equation (2) and Equation (3), one of the popular solutions used is 1 yi +1 = yi + (k1 + 2k2 + 2k3 + k4 )h (4) 6 k1 = f ( xi , yi ) (5a) 1 1 ⎞ ⎛ k2 = f ⎜ xi + h, yi + k1h ⎟ 2 2 ⎠ ⎝ 1 1 ⎛ ⎞ k3 = f ⎜ xi + h, yi + k2 h ⎟ 2 2 ⎝ ⎠ k4 = f (xi + h, yi + k3h )

(5b)

(5c) (5d)

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 × 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. Solution dθ = −2.2067 × 10−12 θ 4 − 81 × 108 dt f (t ,θ ) = −2.2067 × 10−12 θ 4 − 81 × 108 1 θ i +1 = θ i + (k1 + 2k 2 + 2k 3 + k 4 )h 6 For i = 0 , t 0 = 0 , θ 0 = 1200 K

(

)

(

)

(

)

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

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

k1 = f (t0 ,θ0 ) = f (0,1200 ) = −2.2067 × 10−12 (12004 − 81 × 108 ) = −4.5579 1 1 ⎞ ⎛ k2 = f ⎜ t0 + h, θ0 + k1h ⎟ 2 2 ⎠ ⎝ 1 1 ⎛ ⎞ = f ⎜ 0 + (240),1200 + (− 4.5579)240 ⎟ 2 2 ⎝ ⎠ = f (120,653.05) = −2.2067 × 10−12 653.054 − 81 × 108 = −0.38347 1 1 ⎛ ⎞ k3 = f ⎜ t0 + h, θ0 + k2 h ⎟ 2 2 ⎝ ⎠ 1 1 ⎛ ⎞ = f ⎜ 0 + (240),1200 + (− 0.38347 )240 ⎟ 2 2 ⎝ ⎠ = f (120,1154.0 ) = 2.2067 × 10−12 1154.04 − 81 × 108 = −3.8954 1 1 ⎞ ⎛ k 4 = f ⎜ t 0 + h,θ 0 + k 3 h ⎟ 2 2 ⎝ ⎠ 1 1 ⎛ ⎞ = f ⎜ 0 + 240,1200 + (− 3.894)240 ⎟ 2 2 ⎝ ⎠ = f (120,265.10 ) = 2.2067 × 10−12 265.104 − 81 × 108 = 0.0069750 1 θ1 = θ 0 + (k1 + 2k 2 + 2k 3 + k 4 )h 6 1 = 1200 + (− 4.5579 + 2(− 0.38347 ) + 2(− 3.8954) + (0.069750))240 6 ( = 1200 + − 2.1848) 240 = 675.65 K θ1 is the approximate temperature at t = t1 = t0 + h = 0 + 240 = 240 θ1 = θ (240 ) ≈ 675.65 K For i = 1, t1 = 240,θ1 = 675.65 K k1 = f (t1 ,θ1 )

(

)

(

)

(

)

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

-3–

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

= f (240,675.65)

(

)

(

)

(

)

-4–

= −2.2067 × 10−12 675.654 − 81 × 108 = −0.44199 1 1 ⎞ ⎛ k2 = f ⎜ t1 + h, θ1 + k1h ⎟ 2 2 ⎠ ⎝ 1 1 ⎛ ⎞ = f ⎜ 240 + (240),675.65 + (− 0.44199)240 ⎟ 2 2 ⎝ ⎠ ( ) = f 360,622.61 = −2.2067 × 10−12 622.614 − 81 × 108 = −0.31372 1 1 ⎛ ⎞ k3 = f ⎜ t1 + h, θ1 + k2 h ⎟ 2 2 ⎝ ⎠ 1 1 ⎛ ⎞ = f ⎜ 240 + (240),675.65 + (− 0.31372)240 ⎟ 2 2 ⎝ ⎠ ( ) = f 360,638.00 = −2.2067 × 10−12 638.004 − 81 × 108 = −0.34775 1 ⎞ ⎛ k 4 = f ⎜ t1 + h,θ 1 + k 3 h ⎟ 2 ⎝ ⎠ 1 ⎛ ⎞ = f ⎜ 240 + 240,675.65 + (− 0.34775)240 ⎟ 2 ⎝ ⎠ ( ) = f 360,592.19 = 2.2067 × 10−12 (592.194 − 81 × 108 ) = −0.25351 1 θ 2 = θ1 + (k1 + 2k 2 + 2k 3 + k 4 )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 t = t2 t2 = t1 + h = 240 + 240 = 480 θ 2 = θ (10 ) ≈ 594.91K Figure 1 compares the exact solution with the numerical solution using Runge-Kutta 4th order method step size of h=240.

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

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

-5–

Temperature, θ(K)

1600 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 for different step sizes

Table 1 and Figure 2 shows the effect of step size on the value of the calculated temperature at t=480 seconds. Table 1. Value of temperature at time, t=480s for different step sizes

Step size, h 480 240 120 60 30

Et |εt | % θ (480 ) -90.278 737.85 113.94 594.91 52.660 8.1319 646.16 1.4122 0.21807 647.54 0.033626 0.0051926 647.57 0.00086900 0.00013419

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

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

-6–

Temperature, θ(480)

800 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

In Figure 3, we are comparing the exact results with Euler’s method (Runge-Kutta 1st order method), Heun’s method (Runge-Kutta 2nd order method) and Runge-Kutta 4th order method.

Temperature, θ(K)

1400 1200

4th order

1000 800

Exact

600

Heun

400

Euler

200 0

0

100

200

300

400

500

Time, t(sec)

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

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

RUNGE-KUTTA 4TH ORDER METHOD - ORDINARY DIFFERENTIAL EQUATIONS

-7–

The formula described in this chapter was developed by Runge. This formula is same 1 as Simpson’s rd rule, if f ( x, y ) is only a function of x . There are other versions of 3 the fourth order method just like there are several versions of the second order methods. Formula developed by Kutta is 1 (6) yi +1 = yi + (k1 + 3k 2 + 3k 3 + k 4 )h 8 where k1 = f ( xi , yi ) (7a) 1 1 ⎛ ⎞ k 2 = f ⎜ xi + h, yi + hk1 ⎟ 3 3 ⎝ ⎠ 2 1 ⎛ ⎞ k 3 = f ⎜ xi + h, yi − hk1 + hk 2 ⎟ 3 3 ⎝ ⎠ k 4 = f ( xi + h, yi + hk1 − hk 2 + hk 3 ) This formula is the same as the Simpson’s

(7b) (7c) (7d) 3 th rule, if f ( x, y ) is only a function of 8

x.

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