Runge-Kutta method

1 Runge-Kutta methodThe formula for the fourth order Runge-Kutta method (RK4) is given below. Consider theproblem{y =f(t,y)y(t0) = Definehto be the time step size andti=t0+ih. Then the following formulaw0= k1=hf(ti, wi)k2=hf(ti+h2, wi+k12)k3=hf(ti+h2, wi+k22)k4=hf(ti+h, wi+k3)wi+1=wi+16(k1+ 2k2+ 2k3+k4)computes an approximate solution, that iswi y(ti).Let us look at an example:{y =y t2+ 1y(0) = exact solution for this problem isy=t2+ 2t+ 1 12et, and we are interested in the value ofyfor0 t We first solve this problem using RK4 withh= 0tot= 2with step sizeh= , it takes 4 steps:t0= 0,t1= ,t2= 1,t3= ,t4= 0t0= 0,w0= 1t1= (t0,w0) = (0, ) = (t0+h/2,w0+k1/2) = ( , ) = (t0+h/2,w0+k2/2) = ( , ) = (t0+h,w0+K3) = ( , ) = + (k1+ 2k2+ 2k3+k4)/6 = 2t2= 1k1=hf(t1,w1) = ( , ) = (t1+h/2,w1+k1/2) = ( , ) = (t1+h/2,w1+k2/2) = ( , ) = (t1+h,w1+K3) = (1, ) = + (k1+ 2k2+ 2k3+k4)/6 = 3t3= (t2,w2) = (1, ) = (t2+h/2,w2+k1/2) = ( , ) = (t2+h/2,w2+k2/2) = ( , ) = (t2+h,w2+K3) = ( , ) = + (k1+ 2k2+ 2k3+k4)/6 = 4t4= 2k1=hf(t3,w3) = ( , ) = (t3+h/2,w3+k1/2) = ( , ) = (t3+h/2,w3+k2/2) = ( , ) = (t3+h,w3+K3) = (2, ) = + (k1+ 2k2+ 2k3+k4)/6 = let s compare what we got with the exact solutiontiExact solutiony(ti)Numerical solutionwiError|wi y(ti)| this can be done by using Matlab:function rungekuttah = ;t = 0;w = ;fprintf( Step 0: t = % , w = % \n , t, w).}}

2 For i=1:4k1 = 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 = % , w = % \n , i, t, w);end%%%%%%%%%%%%%%%%%%function v = f(t,y)v = y-t 2+1;22. Solve the problem using RK4 withh= you need to do is to replaceh = ;andfor i=1:4in the above Matlab program intoh = i=1:10. Then we havetiExact solutiony(ti)Numerical solutionwiError|wi y(ti)| Solve the problem using RK4 withh= , all need to do is to seth = i=1:40. Then we havetiExact solutiony(ti)Numerical solutionwiError|wi y(ti)| Comparing the resultsBy comparing the above results, we can see thatforh= ,the error att= ,using4stepsforh= ,the error att= ,using10stepsforh= ,the error att= ,using40stepsThis brings out a question.

3 How to find the most appropriate step size, if we want to make sure theerror is less than a given , for example, = Adaptive step size control and the Runge-Kutta -Fehlberg methodThe answer is, we willuse adaptive step size control during the computation. The idea is to start with a moderate stepsize. When we detect the expected error is larger than , reduce the step size and recalculate thecurrent step. When we detect the expected error is less than , keep the current step and slightlyenlarge the step size in the next step. This requires us to have a good estimation of the expectederror .3 Here s the formula for the Runge-Kutta -Fehlberg method (RK45).w0= k1=hf(ti, wi)k2=hf(ti+h4, wi+k14)k3=hf(ti+3h8, wi+332k1+932k2)k4=hf(ti+12h13, wi+19322197k1 72002197k2+72962197k3)k5=hf(ti+h, wi+439216k1 8k2+3680513k3 8454104k4)k6=hf(ti+h2, wi 827k1+ 2k2 35442565k3+18594104k4 1140k5)wi+1=wi+25216k1+14082565k3+219741 04k4 15k5 wi+1=wi+16135k1+665612825k3+2856156430k4 950k5+255k6R=1h| wi+1 wi+1| = ( R)1/4ifR keepwas the current step solutionand move to the next step with step size hifR > recalculate the current step with step size hThe Matlab program is given below:function rk45epsilon = ;h = ;t = 0;w = ;i = 0;fprintf( Step %d: t = % , w = % \n , i, t, w);while t<2h = min(h, 2-t);k1 = h*f(t,w);4k2 = h*f(t+h/4, w+k1/4);k3 = h*f(t+3*h/8, w+3*k1/32+9*k2/32);k4 = h*f(t+12*h/13, w+1932*k1/2197-7200*k2/2197+7296*k3/2197 );k5 = h*f(t+h, w+439*k1/216-8*k2+3680*k3/513-845*k4/410 4).

4 K6 = h*f(t+h/2, w-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104 -11*k5/40);w1 = w + 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5 ;w2 = w + 16*k1/135+6656*k3/12825+28561*k4/56430-9 *k5/50+2*k6/55;R = abs(w1-w2)/h;delta = *(epsilon/R) (1/4);if R<=epsilont = t+h;w = w1;i = i+1;fprintf( Step %d: t = % , w = % \n , i, t, w);h = delta*h;elseh = delta*h;endend%%%%%%%%%%%%%%%%%%function v = f(t,y)v = y-t 2+1;Run this program, we see that starting withh= , the program outputsStep 0: t = , w = 1: t = , w = 2: t = , w = 3: t = , w = 4: t = , w = 5: t = , w = 6: t = , w = 7: t = , w = 8: t = , w = the error is|y(2) w8|=