3.3: The Runge-Kutta Method

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),\\[4pt] k_{2i}&=f \left(x_i+{h\over2},y_i+{h\over2}k_{1i}\right),\\[4pt] k_{3i}&=f\left(x_i+{h\over2},y_i+{h\over2}k_{2i}\right),\\[4pt] 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,\\[4pt] k_{20} & = f(x_0+h/2,y_0+hk_{10}/2)=f(.05,1+(.05)(-2))\\[4pt] &= f(.05,.9)=-2(.9)+(.05)^3e^{-.1}=-1.799886895,\\[4pt] k_{30} & = f(x_0+h/2,y_0+hk_{20}/2)=f(.05,1+(.05)(-1.799886895))\\[4pt] &= f(.05,.910005655)=-2(.910005655)+(.05)^3e^{-.1}=-1.819898206,\\[4pt] k_{40} & = f(x_0+h,y_0+hk_{30})=f(.1,1+(.1)(-1.819898206))\\[4pt] &=f(.1,.818010179)=-2(.818010179)+(.1)^3e^{-.2}=-1.635201628,\\[4pt] y_1&=y_0+{h\over6}(k_{10}+2k_{20}+2k_{30}+k_{40}),\\[4pt] &=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,\\[4pt] k_{21} & = f(x_1+h/2,y_1+hk_{11}/2)=f(.15,.818753803+(.05)(-1.636688875))\\[4pt] &= f(.15,.736919359)=-2(.736919359)+(.15)^3e^{-.3}=-1.471338457,\\[4pt] k_{31} & = f(x_1+h/2,y_1+hk_{21}/2)=f(.15,.818753803+(.05)(-1.471338457))\\[4pt] &= f(.15,.745186880)=-2(.745186880)+(.15)^3e^{-.3}=-1.487873498,\\[4pt] k_{41} & = f(x_1+h,y_1+hk_{31})=f(.2,.818753803+(.1)(-1.487873498))\\[4pt] &=f(.2,.669966453)=-2(.669966453)+(.2)^3e^{-.4}=-1.334570346,\\[4pt] y_2&=y_1+{h\over6}(k_{11}+2k_{21}+2k_{31}+k_{41}),\\[4pt] &=.818753803+{.1\over6}(-1.636688875+2(-1.471338457)+2(-1.487873498)-1.334570346) \\[4pt]&=.670592417.\end{aligned} \nonumber \]

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

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.

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

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.

	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

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

Table 3.3.4 : Numerical solution of \(y'-2xy=1,\ y(0)=3\), by the Runge-Kutta semilinear method.

\(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

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

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

\(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.6 : Numerical solution of \((y-1)^2y'=2x+3,\ y(1)=4\), on \([0,1]\).

\(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

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