3.4: Runge-Kutta Methods
- Page ID
- 91060
As WE HAD SEEN IN THE LAST SECTION, we can use higher order Taylor methods to derive numerical schemes for solving
\[\dfrac{d y}{d t}=f(t, y), \quad y(a)=y_{0}, \quad t \in[a, b] \label{3.18} \]
using a scheme of the form
\[\begin{equation} \begin{aligned}\tilde{y}_{i+1} &=\tilde{y}_{i}+h T^{(n)}\left(t_{i}, \tilde{y}_{i}\right), \quad i=0,1, \ldots, N-1 \\ \tilde{y}_{0} &=y_{0} \end{aligned} \end{equation}\label{3.19} \]
where we have defined
\(T^{(n)}(t, y)=y^{\prime}(t)+\dfrac{h}{2} y^{\prime \prime}(t)+\cdots+\dfrac{h^{(n-1)}}{n !} y^{(n)}(t)\)
In this section we will find approximations of \(T^{(n)}(t, y)\) which avoid the need for computing the derivatives.
For example, we could approximate
\(T^{(2)}(t, y)=f(t, y)+\dfrac{h}{2} \dfrac{d f}{d t}(t, y)\)
by
\(T^{(2)}(t, y) \approx a f(t+\alpha, y+\beta)\)
for selected values of \(a, \alpha\), and \(\beta .\) This requires use of a generalization of Taylor’s series to functions of two variables. In particular, for small \(\alpha\) and \(\beta\) we have
\[\begin{equation} \begin{aligned}a f(t+\alpha, y+\beta)=& a\left[f(t, y)+\dfrac{\partial f}{\partial t}(t, y) \alpha+\dfrac{\partial f}{\partial y}(t, y) \beta\right. \left. +\dfrac{1}{2}\left(\dfrac{\partial^{2} f}{\partial t^{2}}(t, y) \alpha^{2}+2 \dfrac{\partial^{2} f}{\partial t \partial y}(t, y) \alpha \beta+\dfrac{\partial^{2} f}{\partial y^{2}}(t, y) \beta^{2}\right)\right] \\ &+\text { higher order terms. } \end{aligned} \end{equation}\label{3.20} \]
Furthermore, we need \(\dfrac{d f}{d t}(t, y) .\) Since \(y=y(t)\), this can be found using a generalization of the Chain Rule from Calculus III:
\(\dfrac{d f}{d t}(t, y)=\dfrac{\partial f}{\partial t}+\dfrac{\partial f}{\partial y} \dfrac{d y}{d t}\)
Thus,
\(T^{(2)}(t, y)=f(t, y)+\dfrac{h}{2}\left[\dfrac{\partial f}{\partial t}+\dfrac{\partial f}{\partial y} \dfrac{d y}{d t}\right]\)
Comparing this expression to the linear (Taylor series) approximation of \(a f(t+\alpha, y+\beta)\), we have
\[\begin{equation} \begin{aligned}
T^{(2)} & \approx a f(t+\alpha, y+\beta) \\
f+\dfrac{h}{2} \dfrac{\partial f}{\partial t}+\dfrac{h}{2} f \dfrac{\partial f}{\partial y} & \approx a f+a \alpha \dfrac{\partial f}{\partial t}+\beta \dfrac{\partial f}{\partial y}
\end{aligned} \end{equation} \label{3.21} \]
We see that we can choose
\(a=1, \quad \alpha=\dfrac{h}{2}, \quad \beta=\dfrac{h}{2} f.\)
This leads to the numerical scheme
\[\begin{equation} \begin{aligned}\tilde{y}_{i+1} &=\tilde{y}_{i}+h f\left(t_{i}+\dfrac{h}{2}, \tilde{y}_{i}+\dfrac{h}{2} f\left(t_{i}, \tilde{y}_{i}\right)\right), \quad i=0,1, \ldots, N-1, \\ \tilde{y}_{0} &=y_{0}, \end{aligned} \end{equation} \label{3.22} \]
The Midpoint or Second Order Runge-Kutta Method
This Runge-Kutta scheme is called the Midpoint Method, or Second Order, and it has order 2 if all second order derivatives of \(f(t, y)\) are bounded. The Midpoint or Second Order RungeKutta Method. Often, in implementing Runge-Kutta schemes, one computes the arguments separately as shown in the following MATLAB code snippet. (This code snippet could replace the Euler’s Method section in the code in the last section.)
% Midpoint Method
y(1)=1;
for i=2: N+1
k1=h/2*f(t(i-1),y(i-1));
k2=h*f(t(i-1)+h/2,
y(i-1)+k1);
y(i)=y(i-1)+k2;
t(i)=t(i-1)+h;
end
Compare the Midpoint Method with the 2nd Order Taylor’s Method for the problem
\[y^{\prime}=t^{2}+y, \quad y(0)=1, \quad t \in[0,1] \nonumber \]
.
The solution to this problem is \(y(t) = 3e^t − 2 − 2t − t^2\). In order to implement the 2nd Order Taylor’s Method, we need
\[\begin{equation} \begin{aligned}
&T^{(2)}=f(t, y)+\dfrac{h}{2} f^{\prime}(t, y) \\
&=t^{2}+y+\dfrac{h}{2}\left(2 t+t^{2}+y\right)
\end{aligned} \end{equation}\label{3.24} \]
The results of the implementation are shown in Table \(\PageIndex{1}\).
| Exact | Taylor | Error | Midpoint | Error |
|---|---|---|---|---|
| \(1.0000\) | \(1.0000\) | \(0.0000\) | \(1.0000\) | \(0.0000\) |
| \(1.1055\) | \(1.1050\) | \(0.0005\) | \(1.1053\) | \(0.0003\) |
| \(1.2242\) | \(1.2231\) | \(0.0011\) | \(1.2236\) | \(0.0006\) |
| \(1.3596\) | \(1.3577\) | \(0.0019\) | \(1.3585\) | \(0.0010\) |
| \(1.5155\) | \(1.5127\) | \(0.0028\) | \(1.5139\) | \(0.0016\) |
| \(1.6962\) | \(1.6923\) | \(0.0038\) | \(1.6939\) | \(0.0023\) |
| \(1.9064\) | \(1.9013\) | \(0.0051\) | \(1.9032\) | \(0.0031\) |
| \(2.1513\) | \(2.1447\) | \(0.0065\) | \(2.1471\) | \(0.0041\) |
| \(2.4366\) | \(2.4284\) | \(0.0083\) | \(2.4313\) | \(0.0053\) |
| \(2.7688\) | \(2.7585\) | \(0.0103\) | \(2.7620\) | \(0.0068\) |
| \(3.1548\) | \(3.1422\) | \(0.0126\) | \(3.1463\) | \(0.0085\) |
There are other way to approximate higher order Taylor polynomials. For example, we can approximate \(T^{(3)}(t, y)\) using four parameters by
\(T^{(3)}(t, y) \approx a f(t, y)+b f(t+\alpha, y+\beta f(t, y).\)
Expanding this approximation and using
\(T^{(3)}(t, y) \approx f(t, y)+\dfrac{h}{2} \dfrac{d f}{d t}(t, y)+\dfrac{h^{2}}{6} \dfrac{d f}{d t}(t, y)\)
we find that we cannot get rid of \(O\left(h^{2}\right)\) terms. Thus, the best we can do is derive second order schemes. In fact, following a procedure similar to the derivation of the Midpoint Method, we find that
\(a+b=1, \quad, \alpha b=\dfrac{h}{2}, \beta=\alpha\)
There are three equations and four unknowns. Therefore there are many second order methods. Two classic methods are given by the modified Euler method \(\left(a=b=\dfrac{1}{2}, \alpha=\beta=h\right)\) and Huen’s method \(\left(a=\dfrac{1}{4}, b=\dfrac{3}{4},\right.\), \(\left.\alpha=\beta=\dfrac{2}{3} h\right) .\)
The Fourth Order Runge-Kutta
The Fourth Order Runge-Kutta Method, which is most often used, is given the scheme
\[\begin{aligned}\tilde{y}_{0}=y_{0} \\k_{1}=h f\left(t_{i}, \tilde{y}_{i}\right) \\k_{2}=h f\left(t_{i}+\dfrac{h}{2}, \tilde{y}_{i}+\dfrac{1}{2} k_{1}\right) \\k_{3}=h f\left(t_{i}+\dfrac{h}{2}, \tilde{y}_{i}+\dfrac{1}{2} k_{2}\right) \\k_{4}=h f\left(t_{i}+h, \tilde{y}_{i}+k_{3}\right) \\k_{4}=h f\left(t_{i}+h, \tilde{y}_{i}+k_{3}\right) \end{aligned} \nonumber \]
\[\tilde{y}_{i+1}=\tilde{y}_{i}+\dfrac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \quad i=0,1, \ldots, N-1 . \nonumber \]
Again, we can test this on Example \(\PageIndex{1}\) with \(N=10\). The MATLAB implementation is given by
% Runge-Kutta 4th Order to solve dy/dt = f(t,y), y(a)=y0, on [a,b]
clear
a=0;
b=1;
N=10;
h=(b-a)/N;
% Slope function
f = inline(’t^2+y’,’t’,’y’);
sol = inline(’-2-2*t-t^2+3*exp(t)’,’t’);
% Initial Condition
t(1)=0;
y(1)=1;
% RK4 Method
y1(1)=1;
for i=2:N+1
k1=h*f(t(i-1),y1(i-1));
k2=h*f(t(i-1)+h/2,y1(i-1)+k1/2);
k3=h*f(t(i-1)+h/2,y1(i-1)+k2/2);
k4=h*f(t(i-1)+h,y1(i-1)+k3);
y1(i)=y1(i-1)+(k1+2*k2+2*k3+k4)/6;
t(i)=t(i-1)+h;
end
MATLAB has built-in ODE solvers, as do other software packages, like Maple and Mathematica. You should also note that there are currently open source packages, such as Python based NumPy and Matplotlib, or Octave, of which some packages are contained within the Sage Project.
MATLAB has built-in ODE solvers, such as ode45 for a fourth order Runge-Kutta method. Its implementation is given by
[t,y]=ode45(f,[0 1],1);
In this case \(f\) is given by an inline function like in the above RK 4 code. The time interval is enetered as \([0,1]\) and the 1 is the initial condition, \(y(0)=\)
However, ode 45 is not a straight forward RK \(_{4}\) implementation. It is a hybrid method in which a combination of 4 th and 5 th order methods are combined allowing for adaptive methods to handled subintervals of the integration region which need more care. In this case, it implements a fourth order Runge-Kutta-Fehlberg method. Running this code for the above example actually results in values for \(N=41\) and not \(N=10\). If we wanted to have the routine output numerical solutions at specific times, then one could use the following form
tspan=0:h:1;
[t,y]=ode45(f,tspan,1);
In Table \(\PageIndex{2}\) we show the solutions which results for Example \(\PageIndex{1}\) paring the \(R K_{4}\) snippet above with ode45. As you can see \(R K_{4}\) is much better than the previous implementation of the second order RK (Midpoint) Method. However, the MATLAB routine is two orders of magnitude better that \(\mathrm{RK}_{4}\).
| Exact | Taylor | Error | Midpoint | Error |
|---|---|---|---|---|
| \(1.0000\) | \(1.0000\) | o.0000 | \(1.0000\) | o.0000 |
| \(1.1055\) | \(1.1055\) | \(4.5894 \mathrm{e}-08\) | \(1.1055\) | \(-2.5083 \mathrm{e}-10\) |
| \(1.2242\) | \(1.2242\) | \(1.2335 \mathrm{e}-07\) | \(1.2242\) | \(-6.0935 \mathrm{e}-10\) |
| \(1.3596\) | \(1.3596\) | \(2.3850 \mathrm{e}-07\) | \(1.3596\) | \(-1.0954 \mathrm{e}-09\) |
| \(1.5155\) | \(1.5155\) | \(3.9843 \mathrm{e}-07\) | \(1.5155\) | \(-1.7319 \mathrm{e}-09\) |
| \(1.6962\) | \(1.6962\) | \(6.1126 \mathrm{e}-07\) | \(1.6962\) | \(-2.5451 \mathrm{e}-09\) |
| \(1.9064\) | \(1.9064\) | \(8.8636 \mathrm{e}-07\) | \(1.9064\) | \(-3.5651 \mathrm{e}-09\) |
| \(2.1513\) | \(2.1513\) | \(1.2345 \mathrm{e}-06\) | \(2.1513\) | \(-4.8265 \mathrm{e}-09\) |
| \(2.4366\) | \(2.4366\) | \(1.6679 \mathrm{e}-06\) | \(2.4366\) | \(-6.3686 \mathrm{e}-09\) |
| \(2.7688\) | \(2.7688\) | \(2.2008 \mathrm{e}-06\) | \(2.7688\) | \(-8.2366 \mathrm{e}-09\) |
| \(3.1548\) | \(3.1548\) | \(2.8492 \mathrm{e}-06\) | \(3.1548\) | \(-1.0482 \mathrm{e}-08\) |
There are many ODE solvers in MATLAB. These are typically useful if \(\mathrm{RK}_{4}\) is having difficulty solving particular problems. For the most part, one which is similar to odens but combining a second and third order scheme. Which is similar to ode 45 but combining a second and third order scheme. Applying the results to Example \(3.3\) we obtain the results in table \(\PageIndex{2}\) We are shown below.
% Second Order RK Method
y1(1)=1;
for i=2:N+1
k1=h*f(t(i-1),y1(i-1));
k2=h*f(t(i-1)+h/2,y1(i-1)+k1/2);
y1(i)=y1(i-1)+k2;
t(i)=t(i-1)+h;
end
tspan=0:h:1;
[t,y]=ode23(f,tspan,1);
| \(1.0000\) | \(1.0000\) | \(0.0000\) | \(1.0000\) | \(0.0000\) |
| \(1.1055\) | \(1.1053\) | \(0.0003\) | \(1.1055\) | \(2.7409 \mathrm{e}-06\) |
| \(1.2242\) | \(1.2236\) | \(0.0006\) | \(1.2242\) | \(8.7114 \mathrm{e}-06\) |
| \(1.3596\) | \(1.3585\) | \(0.0010\) | \(1.3596\) | \(1.6792 \mathrm{e}-05\) |
| \(1.5155\) | \(1.5139\) | \(0.0016\) | \(1.5154\) | \(2.7361 \mathrm{e}-05\) |
| \(1.6962\) | \(1.6939\) | \(0.0023\) | \(1.6961\) | \(4.0853 \mathrm{e}-05\) |
| \(1.9064\) | \(1.9032\) | \(0.0031\) | \(1.9063\) | \(5.7764 \mathrm{e}-05\) |
| \(2.1513\) | \(2.1471\) | \(0.0041\) | \(2.1512\) | \(7.8665 \mathrm{e}-05\) |
| \(2.4366\) | \(2.4313\) | \(0.0053\) | \(2.4365\) | \(0.0001\) |
| \(2.7688\) | \(2.7620\) | \(0.0068\) | \(2.7687\) | \(0.0001\) |
| \(3.1548\) | \(3.1463\) | \(0.0085\) | \(3.1547\) | \(0.0002\) |
We have seen several numerical schemes for solving initial value problems. There are other methods, or combinations of methods, which aim to refine the numerical approximations efficiently as if the step size in the current methods were taken to be much smaller. Some methods extrapolate solutions to obtain information outside of the solution interval. Others use one scheme to get a guess to the solution while refining, or correcting, this to obtain better solutions as the iteration through time proceeds. Such methods are described in courses in numerical analysis and in the literature. At this point we will apply these methods to several physics problems before continuing with analytical solutions.