Runge-Kutta method
The formula for the fourth order Runge-Kutta method (RK4) is given below. Consider the problem y = f (t, y) y(t0 ) = α Define h to be the time step size and ti = t0 + ih. Then the following formula w0 = α k1 = hf (ti , wi ) k2 = hf k3 = hf k1 h ti + , wi + 2 2 h k2 ti + , wi + 2 2

k4 = hf (ti + h, wi + k3 ) 1 wi+1 = wi + (k1 + 2k2 + 2k3 + k4 ) 6 computes an approximate solution, that is wi ≈ y(ti ). Let us look at an example: y = y − t2 + 1 y(0) = 0.5
1 The exact solution for this problem is y = t2 + 2t + 1 − 2 et , and we are interested in the value of y for 0 ≤ t ≤ 2.

1. 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: t0 = 0, t1 = 0.5, t2 = 1, t3 = 1.5, t4 = 2. Step 0 t0 = 0, w0 = 0.5. Step 1 t1 = 0.5 k1 k2 K3 K4 w1 Step 2 t2 = 1 k1 k2 K3 K4 w2 = hf (t1 , w1 ) = 0.5f (0.5, 1.425130208333333) = 1.087565104166667 = hf (t1 + h/2, w1 + k1 /2) = 0.5f (0.75, 1.968912760416667) = 1.203206380208333 = hf (t1 + h/2, w1 + k2 /2) = 0.5f (0.75, 2.0267333984375) = 1.23211669921875 = hf (t1 + h, w1 + K3 ) = 0.5f (1, 2.657246907552083) = 1.328623453776042 = w1 + (k1 + 2k2 + 2k3 + k4 )/6 = 2.639602661132812 1 = hf (t0 , w0 ) = 0.5f (0, 0.5) = 0.75 = hf (t0 + h/2, w0 + k1 /2) = 0.5f (0.25, 0.875) = 0.90625 = hf (t0 + h/2, w0 + k2 /2) = 0.5f (0.25, 0.953125) = 0.9453125 = hf (t0 + h, w0 + K3 ) = 0.5f (0.5, 1.4453125) = 1.09765625 = w0 + (k1 + 2k2 + 2k3 + k4 )/6 = 1.425130208333333

Step 3 t3 = 1.5 k1 k2 K3 K4 w3 = hf (t2 , w2 ) = 0.5f (1, 2.639602661132812) = 1.319801330566406 = hf (t2 + h/2, w2 + k1 /2) = 0.5f (1.25, 3.299503326416016) = 1.368501663208008 = hf (t2 + h/2, w2 + k2 /2) = 0.5f (1.25, 3.323853492736816) = 1.380676746368408 = hf (t2 + h, w2 + K3 ) = 0.5f (1.5, 4.020279407501221) = 1.385139703750610 = w2 + (k1 + 2k2 + 2k3 + k4 )/6 = 4.006818970044454

Step 4 t4 = 2 k1 k2 K3 K4 w4 = hf (t3 , w3 ) = 0.5f (1.5, 4.006818970044454) = 1.378409485022227 = hf (t3 + h/2, w3 + k1 /2) = 0.5f (1.75, 4.696023712555567) = 1.316761856277783 = hf (t3 + h/2, w3 + k2 /2) = 0.5f (1.75, 4.665199898183346) = 1.301349949091673 = hf (t3 + h, w3 + K3 ) = 0.5f (2, 5.308168919136127) = 1.154084459568063 = w3 + (k1 + 2k2 + 2k3 + k4 )/6 = 5.301605229265987

Now let’s compare what we got with the exact solution ti 0.0 0.5 1.0 1.5 2.0 Exact solution y(ti ) 0.5 1.425639364649936 2.640859085770477 4.009155464830968 5.305471950534675 Numerical solution wi 0.5 1.425130208333333 2.639602661132812 4.006818970044454 5.301605229265987 Error |wi − y(ti )| 0 0.000509156316603 0.001256424637665 0.002336494786515 0.003866721268688

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:4 k1 = h*f(t,w); k2 = h*f(t+h/2, w+k1/2); k3 = h*f(t+h/2, w+k2/2); k4 = h*f(t+h, w+k3); w = w + (k1+2*k2+2*k3+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; 2

2. Solve the problem using RK4 with h = 0.2. All you need to do is to replace h = 0.5; and for i=1:4 in the above Matlab program into h = 0.2 and for i=1:10. Then we have ti 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 Exact solution y(ti ) 0.5 0.829298620919915 1.214087651179365 1.648940599804746 2.127229535753766 2.640859085770477 3.179941538631726 3.732400016577663 4.283483787802442 4.815176267793527 5.305471950534674 Numerical solution wi 0.5 0.829293333333333 1.214076210666667 1.648922017041600 2.127202684947944 2.640822692728752 3.179894170232231 3.732340072854980 4.283409498318406 4.815085694579435 5.305363000692655 Error |wi − y(ti )| 0 0.000005287586582 0.000011440512698 0.000018582763146 0.000026850805823 0.000036393041726 0.000047368399496 0.000059943722683 0.000074289484036 0.000090573214092 0.000108949842019

3. Solve the problem using RK4 with h = 0.05. Again, all need to do is to set h = 0.05 and for i=1:40. Then we have ti 0.0 0.05 0.10 ··· 1.80 1.85 1.90 1.95 2.00 Exact solution y(ti ) Numerical solution wi 0.5 0.5 0.576864451811988 0.576864446614583 0.657414540962176 0.657414530368210 ··· ··· 4.815176267793529 4.815175898599096 4.942590238699086 4.942589852008494 5.067052778860367 5.067052374183828 5.188156209705356 5.188155786548850 5.305471950534677 5.305471508400809 Error |wi − y(ti )| 0 0.000000005197405 0.000000010593966 ··· 0.000000369194433 0.000000386690592 0.000000404676539 0.000000423156505 0.000000442133868

4. Comparing the results By comparing the above results, we can see that for h = 0.5, for h = 0.2, for h = 0.05, the error at t = 2 is 0.003866721268688, using 4 steps the error at t = 2 is 0.000108949842019, using 10 steps the error at t = 2 is 0.000000442133868, using 40 steps

This brings out a question. How to find the most appropriate step size, if we want to make sure the error is less than a given ε, for example, ε = 0.00001? 5. Adaptive step size control and the Runge-Kutta-Fehlberg method The answer is, we will use adaptive step size control during the computation. The idea is to start with a moderate step size. When we detect the expected error is larger than ε, reduce the step size and recalculate the current step. When we detect the expected error is less than ε, keep the current step and slightly enlarge the step size in the next step. This requires us to have a good estimation of the “expected error”.

3

Here’s the formula for the Runge-Kutta-Fehlberg method (RK45). w0 = α k1 = hf (ti , wi ) k2 = hf k3 = hf k4 = hf k5 = hf k6 = hf wi+1 = wi + k1 h ti + , wi + 4 4 ti + ti + 3h 3 9 , wi + k1 + k2 8 32 32 12h 1932 7200 7296 , wi + k1 − k2 + k3 13 2197 2197 2197 439 3680 845 k1 − 8k2 + k3 − k4 216 513 4104

ti + h, wi +

8 3544 1859 11 h k3 + k4 − k5 ti + , wi − k1 + 2k2 − 2 27 2565 4104 40

25 1408 2197 1 k1 + k3 + k4 − k5 216 2565 4104 5 16 6656 28561 9 2 wi+1 = wi + ˜ k1 + k3 + k4 − k5 + k6 135 12825 56430 50 55 1 ˜ R = |wi+1 − wi+1 | h δ = 0.84 if R ≤ ε if R > ε ε R
1/4

keep w as the current step solution and move to the next step with step size δh recalculate the current step with step size δh

The Matlab program is given below: function rk45 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