
3.3: The Runge-Kutta Method

$$\newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} }$$

$$\newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}}$$

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}).$

The next example, which deals with the initial value problem considered in Examples and Example $$\PageIndex{1}$$, illustrates the computational procedure indicated in the Runge-Kutta method.

Example $$\PageIndex{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 $$\PageIndex{2}$$

Table $$\PageIndex{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 $$\PageIndex{2}$$:

Add text here. For the automatic number to work, you need to add the “AutoNum” template (preferably at3.2.2}, and the values of the exact solution

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

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$$.

Table $$\PageIndex{1}$$: Numerical solution of $$y'+2y=x^3e^{-2x},\ y(0)=1$$, by the Runge-Kuttta method and the improved Euler method.
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

Example $$\PageIndex{3}$$

Table $$\PageIndex{1}$$ 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 [exer:3.2.3}.

Table $$\PageIndex{2}$$: Numerical solution of $$y'=-2y^2+xy+x^2,\ y(0)=1$$, by the Runge-Kuttta method and the improved Euler method.
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

Example $$\PageIndex{4}$$

Tables $$\PageIndex{3}$$ and $$\PageIndex{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 $$\PageIndex{3}$$ and $$\PageIndex{4}$$

 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
 0 3 3 3 3 0.2 3.32785 3.32785 3.32785 3.32785 0.4 3.96606 3.96606 3.96606 3.96606 0.6 5.06704 5.06704 5.06704 5.06704 0.8 6.9367 6.9367 6.9367 6.9367 1 10.1849 10.1849 10.1849 10.1849 1.2 16.0671 16.0671 16.0671 16.0671 1.4 27.2894 27.2894 27.2894 27.2894 1.6 50.0004 50.0004 50.0004 50.0004 1.8 98.983 98.983 98.983 98.983 2 211.954 211.954 211.954 211.954

The Case Where xₒ Isn’t 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 [exer: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 $$\PageIndex{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}.$

(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.)

[table:3.3.5] Numerical solution of $$z'= {2x-3\over(z-1)^2},\ z(-1)=4$$, on $$[-1,0]$$.

 -1 4 -0.9 3.94454 -0.8 3.8893 -0.7 3.83436 -0.6 3.77979 -0.5 3.72568 -0.4 3.67214 -0.3 3.61928 -0.2 3.56724 -0.1 3.51616 0 3.46621

[table:3.3.6] Numerical solution of $$(y-1)^2y'=2x+3,\ y(1)=4$$, on $$[0,1]$$.

 0 3.46621 3.46621 0.1 3.51616 3.51616 0.2 3.56724 3.56724 0.3 3.61928 3.61928 0.4 3.67214 3.67214 0.5 3.72568 3.72568 0.6 3.77979 3.77979 0.7 3.83436 3.83436 0.8 3.8893 3.8893 0.9 3.94454 3.94454 1 4 4

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 [exer:3.3.26} and [exer:3.3.27}).