Kansas State University

Discussion

Euler's method and the improved Euler's method both try to approximate $y(x_0+h)$ by approximating the slope $m$ of the secant line from $(x_0,y(x_0))$ to $(x_0+h,y(x_0+h))$ and using the formula $y(x_0+h)=y(x_0)+mh$ (see figure 1).

Third Order Runge-Kutta Method

To approximate the solution of $$ \frac{dy}{dx}=f(x,y),\qquad y(x_0)=y_0 $$ compute $$ \begin{align} x_1&=x_0+h \\ k_1&=f(x_0,y_0) \\ k_2&=f(x_0+h,y_0+hk_1) \\ k_3&=f(x_0+h/2,y_0+(h/2)(k_1+k_2)/2) \\ y_1&=y_0+\frac{k_1+k_2+4k_3}{6}h \end{align} $$ Then $y(x_1)\approx y_1$.

This method uses the improved Euler method to find an approximate midpoint of the secant line and then takes a weighted average of the slopes at the left and right endpoints and the midpoint. Note that if $f(x,y)$ is just $f(x)$, a function of $x$ alone, then solving the differential equation $\displaystyle\frac{dy}{dx}=f(x)$ is just evaluating the integral $\displaystyle\int_{x_0}^{x_1} f(x)\quad dx$. In this case, the third order Runge-Kutta method is the same as Simpson's rule for numerical approximation of integrals from Calculus 2. Just as for Euler's method and the improved Euler's method, the RK3 method can be easily put in a spreadsheet. The first part of a spreadsheet for the RK3 method applied to $\displaystyle\frac{dy}{dx}=2xy$, $y(0)=1$ is shown below. Note that this is the problem for exercise 1, so you can use this spreadsheet to check your work.

Fourth Order Runge-Kutta method

Just as with Euler's method and the improved Euler's method, we can repeat the process to find an $x_2=x_1+h=x_0+2h$ and $y_2\approx y(x_2)=y(x_0+2h)$ and so on. The fourth order Runge-Kutta method listed above is powerful enough that it was a popular technique for practical computation (and not just for educational purposes). The steps are well suited to programming into a computer, just loop through the basic algorithm as many times as needed to get to the values you are interested in. The first part of a spreadsheet for the RK4 method applied to $\displaystyle\frac{dy}{dx}=2xy$, $y(0)=1$ is shown below. Note that this is the problem for exercise 2, so you can use this spreadsheet to check your work.

Exercises

1. Use the third order Runge-Kutta method with step sizes $h=0.1$ and $h=0.05$ to approximate $y(2)$ for the initial value problem $\displaystyle\frac{dy}{dx}=2xy$, $y(0)=1$. Find the error in both approximations.

2. Use the fourth order Runge-Kutta method with step sizes $h=0.1$ and $h=0.05$ to approximate $y(2)$ for the initial value problem $\displaystyle\frac{dy}{dx}=2xy$, $y(0)=1$. Find the error in both approximations.

3. Show that approximating $y(b)$ using one step of the third order Runge-Kutta method for the initial value problem $\displaystyle\frac{dy}{dx}=f(x)$, $y(a)=y_0$ gives the same formula as Simpson's rule for $\displaystyle\int_a^b f(x)\quad dx.$ Simpson's rule is covered in your calculus textbook (or, if you don't have your book any more, either Wikipedia or Wolfram MathWorld have reasonable explanations).

Discussion

Fehlberg's improvement to the Runge-Kutta method is based on two assumptions.

1. The error in a single step of the improved Euler's method is about $C'h^3$ and the error in a single step of the third order Runge-Kutta method is about $C''h^4$ where $C'$ and $C''$ are constants that depend on the problem but not the step size. (The error in a single step is the error in $y(x_0+h)$ computed with a step-size of $h$. It is also called the local error.)
2. $C''h^4$ is much smaller than $C'h^3$.

To see how Fehlberg used these assumptions to compute a good step-size, consider the following equations where $y_1'$ is the approximation from the improved Euler's method, $y_1''$ is the approximation from the third order Runge-Kutta method, and $y(h)$ is the true value. $$\begin{align} y_1'&\approx y(h)+C'h^3 \\ y_1''&\approx y(h)+C''h^4 \\ |y_1'-y_1''|&\approx|C'h^3-C''h^4|\approx |C'|h^3 \end{align} $$ where the last line uses the assumption that $C''h^4$ is very small compared to $C'h^3$. Now we can solve this last line for the constant $C'$ to get $$ |C'|\approx \frac{|y_1'-y_1''|}{h^3} $$ Once we have this approximation for $C'$, we can pick a step-size $h_1$ to get the local error of the size we want. If we want the local error to be about size $T$, we just take a step-size $h_{\text{new}}$ where $$ h_{\text{new}}=h\left(\frac{T}{|y_1'-y_1''|}\right)^{1/3} $$ and then the error should be about $$ \begin{align} |\text{error}|&\approx |C'h_{\text{new}}^3| \\ &\approx\frac{|y_1'-y_1''|}{h^3}h^3\frac{T}{|y_1'-y_1''|} \\ &\approx T \end{align} $$ You might be a little worried about how all the errors in the different approximations mount up as we carry out all these computations to get our new step-size $h_{\text{new}}$. This is a serious consideration and is dealt with by introducing a chicken factor, usually taken to be 0.9. We actually use a step-size $$ h_1=0.9h\left(\frac{T}{|y_1'-y_1''|}\right)^{1/3}. $$ This is a bit smaller than the approximation we computed above, and thus a bit safer too.

The Runge-Kutta-Fehlberg 2(3) method uses exactly this technique to pick the right step-size. Suppose the initial value problem we want to solve is $$ \frac{dy}{dx}=f(x,y),\qquad y(x_0)=y_0 $$ We have an initial step-size $h$ (taken to be whatever value you fancy, we will update it automatically as needed). We compute the improved Euler's and RK3 estimates in the usual fashion. $$\begin{align} k1&=f(x_0,y_0) \\ k2&=f(x_0+h,y_0+h*k_1) \\ k3&=f(x_0+h/2,y_0+h*(k1+k2)/4) \\ y_1'&=y_0+h*(k1+k2)/2 \\ y_1''&=y_0+h\frac{k1+k2+4k3}{6} \end{align} $$ Next we compute the estimated error $$ |\text{error}|\approx|y_1'-y_1''| $$ If this error is small enough, say within a tolerance of $T=.001*\max(|y_0|,1)$ (which is the tolerance used by the matlab program ode23), then we accept this step-size for the current step and let $$\begin{align} x_1&=x_0+h\\ y_1&=y_1'' \end{align}$$ If the error is greater than $T$, we reject this step-size for the current step and leave $x_0$ and $y_0$ as they are. In either case, we choose a new step-size $$ h_{\text{new}}=.9h\left(\frac{T}{|y_1'-y_1''|}\right)^{1/3} $$ We then either compute the next step with the new step-size (if our error was less than $T$) or we repeat the current step with the new step-size $h_{\text{new}}$ (if the error was greater than $T$) and try again to find $x_1$ and $y_1$.

Example

We need to pick an initial step-size to get things started. We have to go from $x_0=0$ to $x=1$ so why not go for it all in one shot and guess $h=1$ initially. Actually, the reason not to use that large a step-size for the beginning is that it is usually too large and often takes a couple of iterations before you get the step-size down to a reasonable size and actually get the process off and running. Choosing a step-size of about the length of the interval divided by 16 or 32 is more typical. But I wanted to be sure I had to reject the estimate sometime in the course of the example so I decided to start off wrong with too large a step size to be sure that happened. We carry out the following computations. $$\begin{align} k1&=0\\ k2&=1\\ k3&=.75\\ y_1'&=.5\\ y_1''&=.666666\ldots\\ |y_1'-y_1''|&=.166666\ldots>.01 \end{align}$$ The estimated error is greater than the tolerance $.01$ so we reject the initial step-size of $h=1$. We compute a new step-size and try again. $$ \begin{align} h_{\text{new}}&=.9*1*(.01/.166666\ldots)^{1/3}=.3523380877\\ k1&=0\\ k2&=.35233880877\\ k3&=.2072045759\\ y_1'&=.062071064\\ y_1''&=.069361064\\ |y_1'-y_1''|&=.00729<.01 \end{align}$$ This time the estimated error is less than the tolerance so we accept the step-size and estimate and compute $$\begin{align} x_1&=x_0+h=.3523380877\\ y_1&=y_1''=.069361064 \end{align}$$ We now compute a new step-size $h_{\text{new}}$ and go on to the next step (Twice in this problem, we compute a new step-size and it turns out to exactly equal the old step-size. This is a freak accident. I couldn't have planned that to happen if I tried.) $$ \begin{align} h_{\text{new}}&=.9*.3523380877*(.01/.00729\ldots)^{1/3}=.3523380877\\ k1&=.4216991517\\ k2&=.9276179121\\ k3&=.7162817215\\ y_2'&=.3061881158\\ y_2''&=.3165523026\\ |y_2'-y_2''|&=.0103641868>.01 \end{align}$$ This time the estimated error is greater than the tolerance so we reject the step-size and try again with the same $x_1$ and $y_1$. $$ \eqalign{ h_{\text{new}}&=.9*.3523380877*(.01/.0103641868)^{1/3}=.3133456655\cr k1&=.4216991515\cr k2&=.8671824185\cr k3&=.6793383478\cr y_2'&=.2712937907\cr y_2''&=.2785837907\cr |y_2'-y_2''|&=.00729<.01\cr }$$ The estimated error is less than the tolerance so we accept the step-size and compute $$\eqalign{ x_2&=x_1+h=.6656837532\cr y_2&=y_2''=.2785837907\cr }$$ We now compute the new step-size for the next step and repeat the process. $$ \eqalign{ h_{\text{new}}&=.9*.3133456655*(.01/.00729)^{1/3}=.3133456655\cr k1&=.9442675439\cr k2&=1.553495351\cr k3&=1.296606171\cr y_3'&=.669915379\cr y_3''&=.679849358\cr |y_3'-y_3''|&=.0099695568<.01\cr }$$ The estimated error is less than the tolerance so we accept the estimate and compute $$\eqalign{ x_3&=x_2+h=.9790294187\cr y_3&=y_3''=.679849358\cr }$$ We now compute the step-size for the next step $$h_{\text{new}}=.9*.3133456655*(.01/.0099695568)^{1/3} =.2822978588$$ But this step-size is too large since $x_3+h=1.261327277>1$ and so it would put us past our final value for $x$. Therefore we shrink $h$ to hit $x=1$ exactly. $$ \eqalign{ h_{\text{new}}&=1-x_3=.0209705813\cr k1&=1.658914354\cr k2&=1.714673334\cr k3&=1.687086169\cr y_4'&=.7152579833\cr y_4''&=.7152620701\cr |y_4'-y_4''|&=.00000408681<.01\cr }$$ The estimated error is less than the tolerance so we accept this estimate and get the final computations $$\eqalign{ x_4&=x_3+h=1\cr y_4&=y_4''=.7152620701\cr }$$ So our approximation for $y(1)=y(x_4)\approx y_4=.7152620701$. The true value is $y(1)=.7182818285$ and the error is .0030197584 or a percent error of 0.42%. You may be wondering how we got a relative error of less than 1% when we had a tolerance of 1% error in each individual step. The trick is that we estimate the error using the less accurate improved Euler's method but use the estimates from the more accurate RK3 method. That way we have a safety margin for our estimate of the error. There are many different Runge-Kutta-Fehlberg methods, all of which involve comparing two different Runge-Kutta approximations to get an estimated error and step-size. The RK2(3) method is a low order method which gives moderately accurate results with only 3 function evaluations at each step. If high accuracy is important, you should use a higher order method. On the other hand, if you just need an answer to within 1% for a problem that doesn't cover too large a range of $x$ values, the RK2(3) method is quick, cheap and effective (with the recommended tolerance of $T=.001$).

Exercises

4. In the lab you discovered halving the step size reduced the error in the improved Euler method by a factor of about 4. Explain why this fits with the global error being proportional to $h^2$ for the improved Euler method. Note, this is why it is called a second order method and denoted with the 2 in naming the RKF23 method)

5. In the ingredients for the Runge-Kutta-Fehlberg method, we assume the local error of the improved Euler's method is about $C'h^3$. How does this fit with the error of the approximation for $y(1)$ being about $k'h^2$ (which is what it means for the error to be proportional to $h^2$)? Hint: How many steps of size $h$ does it take to get from $x=0$ to $x=1$?

6. Explain why you are guaranteed that if the estimated error is greater than the tolerance, so you reject the estimate, the revised step-size will be smaller than the previous step-size (i.e. $ h_{\text{new}}<h$). If the estimated error is less than the tolerance, so you accept the estimate, will the revised step-size necessarily increase for the next step?

For the next two problems you may want to program the RK2(3) method into a calculator or computer. Unfortunately, unlike the first extra credit assignment, this method isn't well suited to a spreadsheet because you keep having to test the step-size and possibly change it (perhaps more than once) at each step, though some students have used a spreadsheet successfully anyway.
7. Approximate $y(2)$ using the RK2(3) method for $y'=2xy$, $y(0)=1$ and a tolerance of $T=.001$. You can use whatever initial step-size you want. Then compute the exact answer for the initial value problem and find the percent error in the RK2(3) method.

8. Approximate $y(2)$ using the RK2(3) method for $y'=x+y$, $y(0)=0$ with a tolerance of $T=.001$.