$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$

# 3.3: The Runge-Kutta Method

• • William F. Trench
• Andrew G. Cowles Distinguished Professor Emeritus (Mathematics) at Trinity University
$$\newcommand{\vecs}{\overset { \rightharpoonup} {\mathbf{#1}} }$$ $$\newcommand{\vecd}{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$$$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\id}{\mathrm{id}}$$ $$\newcommand{\Span}{\mathrm{span}}$$ $$\newcommand{\kernel}{\mathrm{null}\,}$$ $$\newcommand{\range}{\mathrm{range}\,}$$ $$\newcommand{\RealPart}{\mathrm{Re}}$$ $$\newcommand{\ImaginaryPart}{\mathrm{Im}}$$ $$\newcommand{\Argument}{\mathrm{Arg}}$$ $$\newcommand{\norm}{\| #1 \|}$$ $$\newcommand{\inner}{\langle #1, #2 \rangle}$$ $$\newcommand{\Span}{\mathrm{span}}$$

In general, if $$k$$ is any positive integer and $$f$$ satisfies appropriate assumptions, there are numerical methods with local truncation error $$O(h^{k+1})$$ for solving an initial value problem

$\label{eq:3.3.1} y'=f(x,y),\quad y(x_0)=y_0.$

Moreover, it can be shown that a method with local truncation error $$O(h^{k+1})$$ has global truncation error $$O(h^k)$$. In Sections 3.1 and 3.2 we studied numerical methods where $$k=1$$ and $$k=2$$. We’ll skip methods for which $$k=3$$ and proceed to the Runge-Kutta method, the most widely used method, for which $$k=4$$. The magnitude of the local truncation error is determined by the fifth derivative $$y^{(5)}$$ of the solution of the initial value problem. Therefore the local truncation error will be larger where $$|y^{(5)}|$$ is large, or smaller where $$|y^{(5)}|$$ is small. The Runge-Kutta method computes approximate values $$y_1$$, $$y_2$$, …, $$y_n$$ of the solution of Equation \ref{eq:3.3.1} at $$x_0$$, $$x_0+h$$, …, $$x_0+nh$$ as follows: Given $$y_i$$, compute

\begin{align*} k_{1i}&=f(x_i,y_i),\\ k_{2i}&=f \left(x_i+{h\over2},y_i+{h\over2}k_{1i}\right),\\ k_{3i}&=f\left(x_i+{h\over2},y_i+{h\over2}k_{2i}\right),\\ k_{4i}&=f(x_i+h,y_i+hk_{3i}),\end{align*}

and

$y_{i+1}=y_i+{h\over6}(k_{1i}+2k_{2i}+2k_{3i}+k_{4i}).\nonumber$

The next example, which deals with the initial value problem considered in Examples and Example 3.3.1 , illustrates the computational procedure indicated in the Runge-Kutta method.

##### Example 3.3.1

Use the Runge-Kutta method with $$h=0.1$$ to find approximate values for the solution of the initial value problem

$\label{eq:3.3.2} y'+2y=x^3e^{-2x},\quad y(0)=1,$

at $$x=0.1,0.2$$.

Solution

Again we rewrite Equation \ref{eq:3.3.2} as

$y'=-2y+x^3e^{-2x},\quad y(0)=1, \nonumber$

which is of the form Equation \ref{eq:3.3.1}, with

$f(x,y)=-2y+x^3e^{-2x},\ x_0=0,\mbox{ and}\ y_0=1. \nonumber$

The Runge-Kutta method yields

\begin{aligned} k_{10} & = f(x_0,y_0) = f(0,1)=-2,\\ k_{20} & = f(x_0+h/2,y_0+hk_{10}/2)=f(.05,1+(.05)(-2))\\ &= f(.05,.9)=-2(.9)+(.05)^3e^{-.1}=-1.799886895,\\ k_{30} & = f(x_0+h/2,y_0+hk_{20}/2)=f(.05,1+(.05)(-1.799886895))\\ &= f(.05,.910005655)=-2(.910005655)+(.05)^3e^{-.1}=-1.819898206,\\ k_{40} & = f(x_0+h,y_0+hk_{30})=f(.1,1+(.1)(-1.819898206))\\ &=f(.1,.818010179)=-2(.818010179)+(.1)^3e^{-.2}=-1.635201628,\\ y_1&=y_0+{h\over6}(k_{10}+2k_{20}+2k_{30}+k_{40}),\\ &=1+{.1\over6}(-2+2(-1.799886895)+2(-1.819898206) -1.635201628)=.818753803,\\[4pt] k_{11} & = f(x_1,y_1) = f(.1,.818753803)=-2(.818753803))+(.1)^3e^{-.2}=-1.636688875,\\ k_{21} & = f(x_1+h/2,y_1+hk_{11}/2)=f(.15,.818753803+(.05)(-1.636688875))\\ &= f(.15,.736919359)=-2(.736919359)+(.15)^3e^{-.3}=-1.471338457,\\ k_{31} & = f(x_1+h/2,y_1+hk_{21}/2)=f(.15,.818753803+(.05)(-1.471338457))\\ &= f(.15,.745186880)=-2(.745186880)+(.15)^3e^{-.3}=-1.487873498,\\ k_{41} & = f(x_1+h,y_1+hk_{31})=f(.2,.818753803+(.1)(-1.487873498))\\ &=f(.2,.669966453)=-2(.669966453)+(.2)^3e^{-.4}=-1.334570346,\\ y_2&=y_1+{h\over6}(k_{11}+2k_{21}+2k_{31}+k_{41}),\\ &=.818753803+{.1\over6}(-1.636688875+2(-1.471338457)+2(-1.487873498)-1.334570346) \\&=.670592417.\end{aligned}

The Runge-Kutta method is sufficiently accurate for most applications.

##### Example 3.3.2

Table 3.3.1 shows results of using the Runge-Kutta method with step sizes $$h=0.1$$ and $$h=0.05$$ to find approximate values of the solution of the initial value problem

$y'+2y=x^3e^{-2x},\quad y(0)=1 \nonumber$

at $$x=0$$, $$0.1$$, $$0.2$$, $$0.3$$, …, $$1.0$$. For comparison, it also shows the corresponding approximate values obtained with the improved Euler method in Example 3.2.2, and the values of the exact solution

$y={e^{-2x}\over4}(x^4+4).\nonumber$

The results obtained by the Runge-Kutta method are clearly better than those obtained by the improved Euler method in fact; the results obtained by the Runge-Kutta method with $$h=0.1$$ are better than those obtained by the improved Euler method with $$h=0.05$$.

Improved Euler Runge-Kutta
x h=0.1 h=0.05 h=0.1 h-0.05 "exact"
0.0 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
0.1 0.820040937 0.819050572 0.818753803 0.818751370 0.818751221
0.2 0.672734445 0.671086455 0.670592417 0.670588418 0.670588174
0.3 0.552597643 0.550543878 0.549928221 0.549923281 0.549922980
0.4 0.455160637 0.452890616 0.452210430 0.452205001 0.452204669
0.5 0.376681251 0.374335747 0.373633492 0.373627899 0.373627557
0.6 0.313970920 0.311652239 0.310958768 0.310953242 0.310952904
0.7 0.264287611 0.262067624 0.261404568 0.261399270 0.261398947
0.8 0.225267702 0.223194281 0.222575989 0.222571024 0.222570721
0.9 0.194879501 0.192981757 0.192416882 0.192412317 0.192412038
1.0 0.171388070 0.169680673 0.169173489 0.169169356 0.169169104
Table 3.3.1 : Numerical solution of $$y'+2y=x^3e^{-2x},\ y(0)=1$$, by the Runge-Kuttta method and the improved Euler method.
##### Example 3.3.3

Table 3.3.2 shows analogous results for the nonlinear initial value problem

$y'=-2y^2+xy+x^2,\ y(0)=1. \nonumber$

We applied the improved Euler method to this problem in Example 3.2.3.

Improved Euler Runge-Kutta
x h=0.1 h=0.05 h=0.1 h-0.05 "exact"
0.0 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
0.1 0.840500000 0.838288371 0.837587192 0.837584759 0.837584494
0.2 0.733430846 0.730556677 0.729644487 0.729642155 0.729641890
0.3 0.661600806 0.658552190 0.657582449 0.657580598 0.657580377
0.4 0.615961841 0.612884493 0.611903380 0.611901969 0.611901791
0.5 0.591634742 0.588558952 0.587576716 0.587575635 0.587575491
0.6 0.586006935 0.582927224 0.581943210 0.581942342 0.581942225
0.7 0.597712120 0.594618012 0.593630403 0.593629627 0.593629526
0.8 0.626008824 0.622898279 0.621908378 0.621907553 0.621907458
0.9 0.670351225 0.667237617 0.666251988 0.666250942 0.666250842
1.0 0.730069610 0.726985837 0.726017378 0.726015908 0.726015790
Table 3.3.2 : Numerical solution of $$y'=-2y^2+xy+x^2,\ y(0)=1$$, by the Runge-Kuttta method and the improved Euler method.
##### Example 3.3.4

Tables 3.3.3 and 3.3.4 show results obtained by applying the Runge-Kutta and Runge-Kutta semilinear methods to to the initial value problem

$y'-2xy=1,\ y(0)=3, \nonumber$

which we considered in Examples 3.3.3 and 3.3.4

$$x$$ $$h=0.2$$ $$h=0.1$$ $$h=0.05$$ "Exact"
0.0 3.000000000 3.000000000 3.000000000 3.000000000
0.2 3.327846400 3.327851633 3.327851952 3.327851973
0.4 3.966044973 3.966058535 3.966059300 3.966059348
0.6 5.066996754 5.067037123 5.067039396 5.067039535
0.8 6.936534178 6.936690679 6.936700320 6.936700945
1.0 10.184232252 10.184877733 10.184920997 10.184923955
1.2 16.064344805 16.066915583 16.067098699 16.067111677
1.4 27.278771833 27.288605217 27.289338955 27.289392347
1.6 49.960553660 49.997313966 50.000165744 50.000377775
1.8 98.834337815 98.971146146 98.982136702 98.982969504
2.0 211.393800152 211.908445283 211.951167637 211.954462214
$$x$$ $$h=0.2$$ $$h=0.1$$ $$h=0.05$$ "Exact"
0.0 3.000000000 3.000000000 3.000000000 3.000000000
0.2 3.327853286 3.327852055 3.327851978 3.327851973
0.4 3.966061755 3.966059497 3.966059357 3.966059348
0.6 5.067042602 5.067039725 5.067039547 5.067039535
0.8 6.936704019 6.936701137 6.936700957 6.936700945
1.0 10.184926171 10.184924093 10.184923963 10.184923955
1.2 16.067111961 16.067111696 16.067111678 16.067111677
1.4 27.289389418 27.289392167 27.289392335 27.289392347
1.6 50.000370152 50.000377302 50.000377745 50.000377775
1.8 98.982955511 98.982968633 98.982969450 98.982969504
2.0 211.954439983 211.954460825 211.954462127 211.954462214
Table 3.3.4 : Numerical solution of $$y'-2xy=1,\ y(0)=3$$, by the Runge-Kutta semilinear method.

## The Case Where xₒ Is not The Left Endpoint

So far in this chapter we’ve considered numerical methods for solving an initial value problem

$\label{eq:3.3.3} y'=f(x,y),\quad y(x_0)=y_0$

on an interval $$[x_0,b]$$, for which $$x_0$$ is the left endpoint. We haven’t discussed numerical methods for solving Equation \ref{eq:3.3.3} on an interval $$[a,x_0]$$, for which $$x_0$$ is the right endpoint. To be specific, how can we obtain approximate values $$y_{-1}$$, $$y_{-2}$$, …, $$y_{-n}$$ of the solution of Equation \ref{eq:3.3.3} at $$x_0-h, \dots,x_0-nh$$, where $$h=(x_0-a)/n$$? Here’s the answer to this question:

##### Question

Consider the initial value problem

$\label{eq:3.3.4} z'=-f(-x,z),\quad z(-x_0)=y_0,$ on the interval $$[-x_0,-a]$$, for which $$-x_0$$ is the left endpoint. Use a numerical method to obtain approximate values $$z_1$$, $$z_2$$, …, $$z_n$$ of the solution of $$\eqref{eq:3.3.4}$$ at $$-x_0+h$$, $$-x_0+2h$$, …, $$-x_0+nh=-a$$. Then $$y_{-1}=z_1$$, $$y_{-2}=z_2$$, $$\dots$$, $$y_{-n}=z_n$$ are approximate values of the solution of $$\eqref{eq:3.3.3}$$ at $$x_0-h$$, $$x_0-2h$$, …, $$x_0-nh=a$$.

The justification for this answer is sketched in Exercise 3.3.23. Note how easy it is to make the change the given problem Equation \ref{eq:3.3.3} to the modified problem Equation \ref{eq:3.3.4} : first replace $$f$$ by $$-f$$ and then replace $$x$$, $$x_0$$, and $$y$$ by $$-x$$, $$-x_0$$, and $$z$$, respectively.

##### Example 3.3.5

Use the Runge-Kutta method with step size $$h=0.1$$ to find approximate values of the solution of

$\label{eq:3.3.5} (y-1)^2y'=2x+3,\quad y(1)=4$

at $$x=0$$, $$0.1$$, $$0.2$$, …, $$1$$.

Solution

We first rewrite Equation \ref{eq:3.3.5} in the form Equation \ref{eq:3.3.3} as

$\label{eq:3.3.6} y'={2x+3\over(y-1)^2},\quad y(1)=4.$

Since the initial condition $$y(1)=4$$ is imposed at the right endpoint of the interval $$[0,1]$$, we apply the Runge-Kutta method to the initial value problem

$\label{eq:3.3.7} z'={2x-3\over(z-1)^2},\quad z(-1)=4$

on the interval $$[-1,0]$$. (You should verify that Equation \ref{eq:3.3.7} is related to Equation \ref{eq:3.3.6} as Equation \ref{eq:3.3.4} is related to Equation \ref{eq:3.3.3}.) Table [table:3.3.5} shows the results. Reversing the order of the rows in Table [table:3.3.5} and changing the signs of the values of $$x$$ yields the first two columns of Table [table:3.3.6}. The last column of Table [table:3.3.6} shows the exact values of $$y$$, which are given by

$y=1+(3x^2+9x+15)^{1/3}.\nonumber$

(Since the differential equation in Equation \ref{eq:3.3.6} is separable, this formula can be obtained by the method of Section 2.2.)

$$x$$ $$z$$
-1.0 4.000000000
-0.9 3.944536474
-0.8 3.889298649
-0.7 3.834355648
-0.6 3.779786399
-0.5 3.725680888
-0.4 3.672141529
-0.3 3.619284615
-0.2 3.567241862
-0.1 3.516161955
0.0 3.466212070
Table 3.3.5 : Numerical solution of $$z'= {2x-3\over(z-1)^2},\ z(-1)=4$$, on $$[-1,0]$$.
$$x$$ $$y$$ Exact
0.00 3.466212070 3.466212074
0.10 3.516161955 3.516161958
0.20 3.567241862 3.567241864
0.30 3.619284615 3.619284617
0.40 3.672141529 3.672141530
0.50 3.725680888 3.725680889
0.60 3.779786399 3.779786399
0.70 3.834355648 3.834355648
0.80 3.889298649 3.889298649
0.90 3.944536474 3.944536474
1.00 4.000000000 4.000000000
Table 3.3.6 : Numerical solution of $$(y-1)^2y'=2x+3,\ y(1)=4$$, on $$[0,1]$$.

We leave it to you to develop a procedure for handling the numerical solution of Equation \ref{eq:3.3.3} on an interval $$[a,b]$$ such that $$a<x_0<b$$ (Exercises 3.3.26 and 3.3.27).