Taylor's Method

Suppose the solution to the initial value problem:

$$ \frac{dy}{dt} = f(t,y), \quad a \leq t \leq b, \quad y(t_0) = y_0 \\ \label{eq:tm1} $$

has \((n+1)\) continuous derivatives. Expanding the solution \(y(t)\) in terms of its \(n^{th}\) order Taylor series at the point \(t_i\):

$$ f(t_i+h) = f(t_i) + \frac{h}{1!}\left.\frac{df}{dt}\right|_{t=t_i} + \frac{h^2}{2!}\left.\frac{d^2f}{dt^2}\right|_{t=t_i} + \dots + \frac{h^n}{n!}\left.\frac{d^nf}{dt^n}\right|_{t=t_i} + \frac{h^{(n+1)}}{(n+1)!}\left.\frac{d^{(n+1)}f}{dt^{(n+1)}}\right|_{t=t_i} \\ \label{eq:tm2} $$

Successive differentiation of the solution \(y(t)\) gives:

$$ \frac{dy}{dt} = f(t,y(t)) \\ \frac{d^2y}{dt^2} = \frac{dy}{dt} \\ \label{eq:tm3} $$

and in general

$$ \frac{d^ky}{dt^k} = \frac{d^{(k+1)}y}{dt^{(k+1)}} \\ \label{eq:tm4} $$

Substituting \eqref{eq:tm4} into \eqref{eq:tm2}:

$$ f(t_i+h) = f(t_i) + hf(t_i,y(t_i)) + \frac{h^2}{2!}\left.\frac{df}{dt}\right|_{t=t_i} + \dots + \frac{h^n}{n!}\left.\frac{d^{(n-1)}f}{dt^{(n-1)}}\right|_{t=t_i} + \frac{h^{(n+1)}}{(n+1)!}\left.\frac{d^nf}{dt^n}\right|_{t=t_i} \\ \label{eq:tm5} $$

The difference-equation method corresponding to equation \eqref{eq:tm5} may be obtained by omitting the remainder term \(\frac{h^{(n+1)}}{(n+1)!}\left.\frac{d^nf}{dt^n}\right|_{t=t_i}\) and noting that \(t_i + h = t_{i+1}\):

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ \omega_{i+1} & = \omega_i + hT^n(t_i,\omega_i), \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:tm6} \end{equation}

where:

$$ T^n(t_i,\omega_i) = \frac{h}{2!}\left.\frac{df}{dt}\right|_{t=t_i} + \frac{h^{(n-1)}}{n!}\left.\frac{d^{(n-1)}f}{dt^{(n-1)}}\right|_{t=t_i} \\ \label{eq:tm7} $$

The remainder term of \eqref{eq:tm5} is termed the local truncation error of the method. The local truncation error at a specified time step measures the amount by which the exact solution to the differential equations fails to satisfy the difference equation used for the approximation. For example, suppose that \(y_i = y(t_i)\) denotes the exact value to the solution of the initial value problem at \(t=t_i\) and that \(y_{i+1} = y(t_{i+1})\) denotes the approximation solution to the initial value problem at \(t=t_{i+1}\), generated by the difference equation \eqref{eq:tm6}. The local truncation error (\(y_{i+1} - y_i\) may therefore be written as:

\begin{equation} \begin{aligned} y_{i+1} - y_i &= y_{i+1} - \left(y_i + hf(t_i,y(t_i)) + \frac{h^2}{2!}\left.\frac{df}{dt}\right|_{t=t_i} + \dots + \frac{h^n}{n!}\left.\frac{d^{(n-1)}f}{dt^{(n-1)}}\right|_{t=t_i} + \frac{h^{(n+1)}}{(n+1)!}\left.\frac{d^nf}{dt^n}\right|_{t=\zeta_i} \right) \\ & = y_{i+1} - y_i - hT^n(t_i,y_i) - \frac{h^{(n+1)}}{(n+1)!}\left.\frac{d^nf}{dt^n}\right|_{t=\zeta_i} \\ & = y_{i+1} - y_i - hT^n(t_i,y_i) - \tau_{i+1}(h) \\ \end{aligned} \label{eq:tm8} \end{equation}

for some \(\zeta_i\) in \((t_i,t_{i+1})\). The local truncation error may therefore be written as:

$$ \tau_{i+1}(h) = \frac{y_{i+1} - y_i}{h} - T^n(t_i,y_i) = \frac{h^n}{(n+1)!}\left.\frac{d^nf}{dt^n}\right|_{t=\zeta_i} = O(h^n) \\ \label{eq:tm9} $$

for each \(i=1 \dots N-1\).

Runge-Kutta Methods

Runge-Kutta methods preserve the high order local truncation error of the Taylor methods but eliminate the need to determine and evaluate the derivatives of \(f(t,y)\). Suppose that \(f(t,y)\) and all its partial derivatives less than or equal to order \(n+1\) are continuous on \(D = \{(t,y) | a \leq t \leq b, c \lt y \lt d\}\) and that \((t,y) \in D\) and \((t+h,y+k) \in D\). Expressing \(f(t,y)\) as a Taylor series expansion about the point \(t,y)\):

$$ f(t+h,y+k) = f(t,y) + \left[h\left.\frac{\partial f}{\partial t}\right|_{(t,y)} + k\left.\frac{\partial f}{\partial y}\right|_{(t,y)} \right] + \left[\frac{h^2}{2}\left.\frac{\partial^2 f}{\partial t^2}\right|_{(t,y)} + hk\left.\frac{\partial^2 f}{\partial t \partial y}\right|_{(t,y)} + \frac{k^2}{2}\left.\frac{\partial^2 f}{\partial y^2}\right|_{(t,y)}\right] + \dots \\ \label{eq:rk1} $$

The first step in deriving a Runge-Kutta method is to determine values for \(a_1\), \(\alpha_1\) and \(\beta_1\) with the property that \(a_1f(t_i+\alpha_1,y+\beta_1)\) approximates:

$$ T^2(t,y) = f(t,y) + \frac{h}{2}\left.\frac{df}{dt}\right|_{(t,y)} \\ \label{eq:rk2} $$

with error no greater than \(O(h^2)\), the local truncation error of the Taylor method of order two. Since (from the multivariable Chain rule):

$$ \frac{df}{dt} = \frac{\partial f}{\partial t} + \frac{\partial f}{\partial y} \frac{dy}{dt} \\ \label{eq:rk3} $$

and

$$ \frac{dy}{dt} = f(t,y) \\ \label{eq:rk4} $$

then

$$ T^2(t,y) = f(t,y) + \frac{h}{2}\left.\frac{\partial f}{\partial t}\right|_{(t,y)} + \frac{h}{2}\left.\frac{\partial f}{\partial y}\right|_{(t,y)} f(t,y) \\ \label{eq:rk5} $$

Expanding \(a_1f(t_i+\alpha_1,y+\beta_1)\) as a Taylor series of order one about the point \((t,y)\):

$$ a_1f(t_i+\alpha_1,y+\beta_1) = a_1f(t,y) + a_1\left[\alpha_1 \left.\frac{\partial f}{\partial t}\right|_{(t,y)} + \beta_1 \left.\frac{\partial f}{\partial y}\right|_{(t,y)} \right] + a_1R_1(t + \alpha_1, y + \beta_1) \\ \label{eq:rk6} $$

where:

$$ R_1(t + \alpha_1, y + \beta_1) = \left[\frac{{\alpha_1}^2}{2} \left.\frac{\partial^2 f}{\partial t^2}\right|_{(\zeta,\mu)} + \alpha_1\beta_1 \left.\frac{\partial^2 f}{\partial t \partial y}\right|_{(\zeta,\mu)} + \frac{{\beta_1}^2}{2} \left.\frac{\partial^2 f}{\partial y^2}\right|_{(\zeta,\mu)}\right] \\ \label{eq:rk7} $$

for some \(\zeta\) between \(t\) and \(t+\alpha_1\) and \(\mu\) between \(y\) and \(y+\beta_1\). Equating \(T^2(t,y)\) and \(a_1f(t_i+\alpha_1,y+\beta_1)\) and matching coefficients of\(f(t,y)\) and its partial derivatives:

$$ f(t,y) : a_1 = 1 \quad \left.\frac{\partial f}{\partial t}\right|_{(t,y)} : a_1\alpha_1 = \frac{h}{2} \quad \left.\frac{\partial f}{\partial y}\right|_{(t,y)} : a_1\beta_1 = \frac{h}{2}f(t,y) \\ \label{eq:rk8} $$

The parameters \(a_1\), \(\alpha_1\) and \(\beta_1\) are uniquely determined to be:

$$ a_1 = 1 \quad \alpha_1 = \frac{h}{2} \quad \beta_1 = \frac{h}{2}f(t,y) \\ \label{eq:rk9} $$

therefore:

$$ T^2(t,y) = f\left(t + \frac{h}{2}, y + \frac{h}{2}f(t,y)\right) - R_1\left(t + \frac{h}{2}, y + \frac{h}{2}f(t,y)\right) \\ \label{eq:rk10} $$

where:

$$ R_1\left(t + \frac{h}{2}, y + \frac{h}{2}f(t,y)\right) = \left[\frac{h^2}{8} \left.\frac{\partial^2 f}{\partial t^2}\right|_{(\zeta,\mu)} + \frac{h^2}{4} \left.\frac{\partial^2 f}{\partial t \partial y}\right|_{(\zeta,\mu)} + \frac{h^2}{8} [f(t,y)]^2 \left.\frac{\partial^2 f}{\partial y^2}\right|_{(\zeta,\mu)} \right] \\ \label{eq:rk11} $$

If all the second order partial derivatives are bounded, \(R_1\) is \(O(h)\).

The difference equation method resulting from replacing \(T^2(t,y)\) in Taylor's method of order two by \(a_1f(t_i+\alpha_1,y+\beta_1)\) is a specific Runge-Kutta method known as the Midpoint method. The difference equation may be written as:

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ \omega_{i+1} & = \omega_i + hf\left(t_i + \frac{h}{2}, \omega_i + \frac{h}{2}f(t_i,\omega_i)\right), \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:rk12} \end{equation}

If instead of replacing \(T^2(t,y)\) in Taylor's method of order two by \(a_1f(t_i+\alpha_1,y+\beta_1)\), an altternative expression of the form \(a_1f(t,y) + a_2f(t+\alpha_2,y+\delta_2f(t,y))\) is used, then there an infinite number of choices of the values of the four parameters. If \(a_1 = a_2 = 0.5\) and \(\alpha_2 = \delta_2 = h\) then the Modified Euler method is obtained with a difference equation which may be written as:

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ \omega_{i+1} & = \omega_i + \frac{h}{2}\left[f(t,\omega_i) + f(t_{i+1}, \omega_i + hf(t_i,\omega_i)) \right], \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:rk13} \end{equation}

If \(a_1 = \frac{1}{4}\), \(a_2 = \frac{3}{4}\) and \(\alpha_2 = \delta_2 = \frac{2}{3}h\) then Heun's method is obtained with a difference equation which may be written as:

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ \omega_{i+1} & = \omega_i + \frac{h}{4}\left[f(t,\omega_i) + 3f\left(t_{i+1}+\frac{2}{3}h, \omega_i + \frac{2}{3}hf(t_i,\omega_i)\right) \right], \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:rk14} \end{equation}

The Midpoint, Modified Euler and Heun's methods are classified as Runge-Kutta methods of order two due to the order of their local truncation errors.

The most common Runge-Kutta method in use is the Runge-Kutta method of order four. It is derived by replacing \(T^4(t,y)\) in Taylor's method of order four by a function of the form \(f\left(t+\alpha_1,y+\delta_1f\left(t+\alpha_2,y+\delta_2f(t,y)\right)\right)\). In difference equation form, it may be written as:

\begin{equation} \begin{aligned} \omega_0 & = \alpha \\ k_1 & = hf(t_i,\omega_i) \\ k_2 & = hf\left(t_i+\frac{h}{2}, \omega_i + \frac{1}{2}k_1\right) \\ k_3 & = hf\left(t_i + \frac{h}{2}, \omega_i + \frac{1}{2}k_2\right) \\ k_4 & = hf\left(t_{i+1}, \omega_i + k_3\right) \\ \omega_{i+1} & = \omega_i + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4), \quad i=0,\dots,N-1 \\ \end{aligned} \label{eq:rk15} \end{equation}

This method has local truncation error \(O(h^4)\), provided the solution \(y(t)\) has five continuous derivatives.