## 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{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}

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{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}

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{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}

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{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}

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{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}

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{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}

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