Download Range kutta method for solving differential equation and more Lecture notes Inorganic Chemistry in PDF only on Docsity!
Runge-Kutta method
The formula for the fourth order Runge-Kutta method (RK4) is given below. Consider the problem (^) {
y′^ = f (t, y) y(t 0 ) = α
Define h to be the time step size and ti = t 0 + ih. Then the following formula
w 0 = α
k 1 = hf (ti, wi)
k 2 = hf
ti +
h 2
, wi +
k 1 2
k 3 = hf
ti +
h 2
, wi +
k 2 2
k 4 = hf (ti + h, wi + k 3 )
wi+1 = wi +
(k 1 + 2k 2 + 2k 3 + k 4 )
computes an approximate solution, that is wi ≈ y(ti). Let us look at an example: (^) { y′^ = y − t^2 + 1 y(0) = 0. 5
The exact solution for this problem is y = t^2 + 2t + 1 − 12 et, and we are interested in the value of y for 0 ≤ t ≤ 2.
- We first solve this problem using RK4 with h = 0. 5. From t = 0 to t = 2 with step size h = 0. 5 , it takes 4 steps: t 0 = 0, t 1 = 0. 5 , t 2 = 1, t 3 = 1. 5 , t 4 = 2. Step 0 t 0 = 0, w 0 = 0. 5. Step 1 t 1 = 0. 5
k 1 = hf (t 0 , w 0 ) = 0. 5 f (0, 0 .5) = 0. 75 k 2 = hf (t 0 + h/ 2 , w 0 + k 1 /2) = 0. 5 f (0. 25 , 0 .875) = 0. 90625 K 3 = hf (t 0 + h/ 2 , w 0 + k 2 /2) = 0. 5 f (0. 25 , 0 .953125) = 0. 9453125 K 4 = hf (t 0 + h, w 0 + K 3 ) = 0. 5 f (0. 5 , 1 .4453125) = 1. 09765625 w 1 = w 0 + (k 1 + 2k 2 + 2k 3 + k 4 )/6 = 1. 425130208333333 Step 2 t 2 = 1
k 1 = hf (t 1 , w 1 ) = 0. 5 f (0. 5 , 1 .425130208333333) = 1. 087565104166667 k 2 = hf (t 1 + h/ 2 , w 1 + k 1 /2) = 0. 5 f (0. 75 , 1 .968912760416667) = 1. 203206380208333 K 3 = hf (t 1 + h/ 2 , w 1 + k 2 /2) = 0. 5 f (0. 75 , 2 .0267333984375) = 1. 23211669921875 K 4 = hf (t 1 + h, w 1 + K 3 ) = 0. 5 f (1, 2 .657246907552083) = 1. 328623453776042 w 2 = w 1 + (k 1 + 2k 2 + 2k 3 + k 4 )/6 = 2. 639602661132812
Step 3 t 3 = 1. 5
k 1 = hf (t 2 , w 2 ) = 0. 5 f (1, 2 .639602661132812) = 1. 319801330566406 k 2 = hf (t 2 + h/ 2 , w 2 + k 1 /2) = 0. 5 f (1. 25 , 3 .299503326416016) = 1. 368501663208008 K 3 = hf (t 2 + h/ 2 , w 2 + k 2 /2) = 0. 5 f (1. 25 , 3 .323853492736816) = 1. 380676746368408 K 4 = hf (t 2 + h, w 2 + K 3 ) = 0. 5 f (1. 5 , 4 .020279407501221) = 1. 385139703750610 w 3 = w 2 + (k 1 + 2k 2 + 2k 3 + k 4 )/6 = 4. 006818970044454
Step 4 t 4 = 2
k 1 = hf (t 3 , w 3 ) = 0. 5 f (1. 5 , 4 .006818970044454) = 1. 378409485022227 k 2 = hf (t 3 + h/ 2 , w 3 + k 1 /2) = 0. 5 f (1. 75 , 4 .696023712555567) = 1. 316761856277783 K 3 = hf (t 3 + h/ 2 , w 3 + k 2 /2) = 0. 5 f (1. 75 , 4 .665199898183346) = 1. 301349949091673 K 4 = hf (t 3 + h, w 3 + K 3 ) = 0. 5 f (2, 5 .308168919136127) = 1. 154084459568063 w 4 = w 3 + (k 1 + 2k 2 + 2k 3 + k 4 )/6 = 5. 301605229265987
Now let’s compare what we got with the exact solution
ti Exact solution y(ti) Numerical solution wi Error |wi − y(ti)| 0.0 0.5 0.5 0 0.5 1.425639364649936 1.425130208333333 0. 1.0 2.640859085770477 2.639602661132812 0. 1.5 4.009155464830968 4.006818970044454 0. 2.0 5.305471950534675 5.301605229265987 0.
All this can be done by using Matlab:
function rungekutta h = 0.5; t = 0; w = 0.5; fprintf(’Step 0: t = %12.8f, w = %12.8f\n’, t, w); for i=1: k1 = hf(t,w); k2 = hf(t+h/2, w+k1/2); k3 = hf(t+h/2, w+k2/2); k4 = hf(t+h, w+k3); w = w + (k1+2k2+2k3+k4)/6; t = t + h; fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w); end
function v = f(t,y) v = y-tˆ2+1;
Here’s the formula for the Runge-Kutta-Fehlberg method (RK45).
w 0 = α
k 1 = hf (ti, wi)
k 2 = hf
ti +
h 4
, wi +
k 1 4
k 3 = hf
ti +
3 h 8
, wi +
k 1 +
k 2
k 4 = hf
ti +
12 h 13
, wi +
k 1 −
k 2 +
k 3
k 5 = hf
ti + h, wi +
k 1 − 8 k 2 +
k 3 −
k 4
k 6 = hf
ti +
h 2
, wi −
k 1 + 2k 2 −
k 3 +
k 4 −
k 5
wi+1 = wi +
k 1 +
k 3 +
k 4 −
k 5
w ˜i+1 = wi +
k 1 +
k 3 +
k 4 −
k 5 +
k 6
R =
h
| w˜i+1 − wi+1|
δ = 0. 84
( (^) ε R
if R ≤ ε keep w as the current step solution and move to the next step with step size δh
if R > ε recalculate the current step with step size δh
The Matlab program is given below:
function rk epsilon = 0.00001; h = 0.2; t = 0; w = 0.5; i = 0; fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w); while t< h = min(h, 2-t);
k1 = h*f(t,w);
k2 = hf(t+h/4, w+k1/4); k3 = hf(t+3h/8, w+3k1/32+9k2/32); k4 = hf(t+12h/13, w+1932k1/2197-7200k2/2197+7296k3/2197); k5 = hf(t+h, w+439k1/216-8k2+3680k3/513-845k4/4104); k6 = hf(t+h/2, w-8k1/27+2k2-3544k3/2565+1859k4/4104-11k5/40); w1 = w + 25k1/216+1408k3/2565+2197k4/4104-k5/5; w2 = w + 16k1/135+6656k3/12825+28561k4/56430-9k5/50+2k6/55; R = abs(w1-w2)/h; delta = 0.84(epsilon/R)ˆ(1/4);
if R<=epsilon t = t+h; w = w1; i = i+1; fprintf(’Step %d: t = %6.4f, w = %18.15f\n’, i, t, w); h = deltah; else h = deltah; end end
function v = f(t,y) v = y-tˆ2+1;
Run this program, we see that starting with h = 0. 2 , the program outputs
Step 0: t = 0.0000, w = 0. Step 1: t = 0.2000, w = 0. Step 2: t = 0.4353, w = 1. Step 3: t = 0.6766, w = 1. Step 4: t = 0.9264, w = 2. Step 5: t = 1.1902, w = 3. Step 6: t = 1.4806, w = 3. Step 7: t = 1.8537, w = 4. Step 8: t = 2.0000, w = 5.
Notice the error is |y(2) − w 8 | = 0. 00001486603807077103
