4: Differential Equations

We now discuss the numerical solution of ordinary differential equations. We will include the initial value problem and the boundary value problem.

Initial value problem

We begin with the simple Euler method, then discuss the more sophisticated RungeKutta methods, and conclude with the Runge-Kutta-Fehlberg method, as implemented in the MATLAB function ode $45 . \mathrm{m}$. Our differential equations are for $x=x(t)$, where the time $t$ is the independent variable.

Euler method

The Euler method is the most straightforward method to integrate a differential equation. Consider the first-order differential equation

$$\dot{x}=f(t, x), \nonumber$$

with the initial condition $x(0)=x_{0}$. Define $t_{n}=n \Delta t$ and $x_{n}=x\left(t_{n}\right)$. A Taylor series expansion of $x_{n+1}$ results in

$$\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) \\ &=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber$$

The Euler Method is therefore written as

$$x_{n+1}=x\left(t_{n}\right)+\Delta t f\left(t_{n}, x_{n}\right) . \nonumber$$

We say that the Euler method steps forward in time using a time-step $\Delta t$, starting from the initial value $x_{0}=x(0)$. The local error of the Euler Method is $\mathrm{O}\left(\Delta t^{2}\right)$. The global error, however, incurred when integrating to a time $T$, is a factor of $1 / \Delta t$ larger and is given by $\mathrm{O}(\Delta t)$. It is therefore customary to call the Euler Method a first-order method.

Modified Euler method

This method is a so-called predictor-corrector method. It is also the first of what we will see are Runge-Kutta methods. As before, we want to solve (4.1). The idea is to average the value of $\dot{x}$ at the beginning and end of the time step. That is, we would like to modify the Euler method and write

$$x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}\right)\right) . \nonumber$$

The obvious problem with this formula is that the unknown value $x_{n+1}$ appears on the right-hand-side. We can, however, estimate this value, in what is called the predictor step. For the predictor step, we use the Euler method to find

$$x_{n+1}^{p}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) \text {. } \nonumber$$

The corrector step then becomes

$$x_{n+1}=x_{n}+\frac{1}{2} \Delta t\left(f\left(t_{n}, x_{n}\right)+f\left(t_{n}+\Delta t, x_{n+1}^{p}\right)\right) . \nonumber$$

The Modified Euler Method can be rewritten in the following form that we will later identify as a Runge-Kutta method

$$\begin{align}\nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{1}\right), \\ &x_{n+1}=x_{n}+\frac{1}{2}\left(k_{1}+k_{2}\right) . \nonumber \end{align} \nonumber$$

Second-order Runge-Kutta methods

We now derive the complete family of second-order Runge-Kutta methods. Higherorder methods can be similarly derived, but require substantially more algebra.

We again consider the differential equation given by (4.1). A general second-order Runge-Kutta method may be written in the form

$$\begin{align} \nonumber &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}, \nonumber \end{align} \nonumber$$

with $\alpha, \beta, a$ and $b$ constants that define the particular second-order Runge-Kutta method. These constants are to be constrained by setting the local error of the second-order Runge-Kutta method to be $\mathrm{O}\left(\Delta t^{3}\right)$. Intuitively, we might guess that two of the constraints will be $a+b=1$ and $\alpha=\beta$.

We compute the Taylor series of $x_{n+1}$ directly, and from the Runge-Kutta method, and require them to be the same to order $\Delta t^{2}$. First, we compute the Taylor series of $x_{n+1}$. We have

$$\begin{aligned} x_{n+1} &=x\left(t_{n}+\Delta t\right) \\ &=x\left(t_{n}\right)+\Delta t \dot{x}\left(t_{n}\right)+\frac{1}{2}(\Delta t)^{2} \ddot{x}\left(t_{n}\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber$$

Now,

$$\dot{x}\left(t_{n}\right)=f\left(t_{n}, x_{n}\right) . \nonumber$$

The second derivative is more complicated and requires partial derivatives. We have

$$\begin{aligned} \ddot{x}\left(t_{n}\right) &\left.=\frac{d}{d t} f(t, x(t))\right]_{t=t_{n}} \\ &=f_{t}\left(t_{n}, x_{n}\right)+\dot{x}\left(t_{n}\right) f_{x}\left(t_{n}, x_{n}\right) \\ &=f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right) . \end{aligned} \nonumber$$

Therefore,

$$\begin{align} \nonumber x_{n+1}=x_{n}+\Delta t f\left(t_{n}, x_{n}\right) & \\ &+\frac{1}{2}(\Delta t)^{2}\left(f_{t}\left(t_{n}, x_{n}\right)+f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{align} \nonumber$$

Second, we compute $x_{n+1}$ from the Runge-Kutta method given by (4.7). Combining (4.7) into a single expression, we have

$$x_{n+1}=x_{n}+a \Delta t f\left(t_{n}, x_{n}\right) \nonumber$$

$$+b \Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \nonumber$$

We Taylor series expand using

$$\begin{aligned} f\left(t_{n}+\alpha \Delta t, x_{n}+\beta \Delta t f\left(t_{n}, x_{n}\right)\right) & \\ &=f\left(t_{n}, x_{n}\right)+\alpha \Delta t f_{t}\left(t_{n}, x_{n}\right)+\beta \Delta t f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)+\mathrm{O}\left(\Delta t^{2}\right) . \end{aligned} \nonumber$$

The Runge-Kutta formula is therefore

$$\begin{aligned} x_{n+1}=x_{n}+(a+b) \Delta t f &\left(t_{n}, x_{n}\right) \\ &+(\Delta t)^{2}\left(\alpha b f_{t}\left(t_{n}, x_{n}\right)+\beta b f\left(t_{n}, x_{n}\right) f_{x}\left(t_{n}, x_{n}\right)\right)+\mathrm{O}\left(\Delta t^{3}\right) . \end{aligned} \nonumber$$

Comparing (4.9) and (4.10), we find

$$a+b=1, \quad \alpha b=1 / 2, \quad \beta b=1 / 2 . \nonumber$$

Since there are only three equations for four parameters, there exists a family of secondorder Runge-Kutta methods.

The modified Euler method given by (4.6) corresponds to $\alpha=\beta=1$ and $a=$ $b=1 / 2$. Another second-order Runge-Kutta method, called the midpoint method, corresponds to $\alpha=\beta=1 / 2, a=0$ and $b=1$. This method is written as

$$\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right), \\ &x_{n+1}=x_{n}+k_{2} . \end{aligned} \nonumber$$

Higher-order Runge-Kutta methods

The general second-order Runge-Kutta method was given by (4.7). The general form of the third-order method is given by

$$\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right), \\ &k_{2}=\Delta t f\left(t_{n}+\alpha \Delta t, x_{n}+\beta k_{1}\right), \\ &k_{3}=\Delta t f\left(t_{n}+\gamma \Delta t, x_{n}+\delta k_{1}+\epsilon k_{2}\right), \\ &x_{n+1}=x_{n}+a k_{1}+b k_{2}+c k_{3} . \end{aligned} \nonumber$$

The following constraints on the constants can be guessed: $\alpha=\beta, \gamma=\delta+\epsilon$, and $a+b+c=1$. Remaining constraints need to be derived.

The fourth-order method has a $k_{1}, k_{2}, k_{3}$ and $k_{4}$. The fifth-order method requires at least to $k_{6}$. The table below gives the minimum order of the method and the number of stages required.

order	2	3	4	5	6	7	8
minimum # stages	2	3	4	6	7	9	11

Because of the jump in the number of stages required between the fourth- and fifthorder methods, the fourth-order Runge-Kutta method has some appeal. The general fourth-order method with four stages has 13 constants and 11 constraints. A particularly simple fourth-order method that has been widely used in the past by physicists is given by

$$\begin{aligned} &k_{1}=\Delta t f\left(t_{n}, x_{n}\right) \\ &k_{2}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}\right) \\ &k_{3}=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}\right) \\ &k_{4}=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}\right) \\ &x_{n+1}=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right) . \end{aligned} \nonumber$$

Adaptive Runge-Kutta Methods

As in adaptive integration, it is useful to devise an ode integrator that automatically finds the appropriate $\Delta t$. The Dormand-Prince Method, which is implemented in MATLAB’s ode $45 . \mathrm{m}$, finds the appropriate step size by comparing the results of a fifth-order and fourth-order method. It requires six function evaluations per time step, and constructs both a fifth-order and a fourth-order method from the same function evaluations.

Suppose the fifth-order method finds $x_{n+1}$ with local error $\mathrm{O}\left(\Delta t^{6}\right)$, and the fourthorder method finds $x_{n+1}^{\prime}$ with local error $\mathrm{O}\left(\Delta t^{5}\right)$. Let $\varepsilon$ be the desired error tolerance of the method, and let $e$ be the actual error. We can estimate $e$ from the difference between the fifth- and fourth-order methods; that is,

$$e=\left|x_{n+1}-x_{n+1}^{\prime}\right| \text {. } \nonumber$$

Now $e$ is of $\mathrm{O}\left(\Delta t^{5}\right)$, where $\Delta t$ is the step size taken. Let $\Delta \tau$ be the estimated step size required to get the desired error $\varepsilon$. Then we have

$$e / \varepsilon=(\Delta t)^{5} /(\Delta \tau)^{5}, \nonumber$$

or solving for $\Delta \tau$,

$$\Delta \tau=\Delta t\left(\frac{\varepsilon}{e}\right)^{1 / 5} \nonumber$$

On the one hand, if the actual error is less that the desired error, $e<\varepsilon$, then we accept $x_{n+1}$ and do the next time step using the larger value of $\Delta \tau$. On the other hand, if the actual error is greater than the desired error, $e>\varepsilon$, then we reject the integration step and redo the time step using the smaller value of $\Delta \tau$. In practice, one usually increases the time step slightly less and decreases the time step slightly more to prevent the wastefulness of too many failed time steps.

System of differential equations

Our numerical methods can be easily adapted to solve higher-order differential equations, or equivalently, a system of differential equations. First, we show how a secondorder differential equation can be reduced to two first-order equations. Consider

$$\ddot{x}=f(t, x, \dot{x}) . \nonumber$$

This second-order equation can be rewritten as two first-order equations by defining $u=\dot{x}$. We then have the system

$$\dot{x}=u, \quad \dot{u}=f(t, x, u) . \nonumber$$

This trick also works for higher-order equations. For example, the third-order equation

$$\dddot x=f(t, x, \dot{x}, \ddot{x}), \nonumber$$

can be written as

$$\dot{x}=u, \quad \dot{u}=v, \quad \dot{v}=f(t, x, u, v) . \nonumber$$

We can generalize the Runge-Kutta method to solve a system of differential equations. As an example, consider the following system of two odes,

$$\dot{x}=f(t, x, y), \quad \dot{y}=g(t, x, y), \nonumber$$

with the initial conditions $x(0)=x_{0}$ and $y(0)=y_{0}$. The generalization of the commonly used fourth-order Runge-Kutta method would be

$$\begin{aligned} k_{1} &=\Delta t f\left(t_{n}, x_{n}, y_{n}\right), \\ l_{1} &=\Delta t g\left(t_{n}, x_{n}, y_{n}\right), \\ k_{2} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ l_{2} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{1}, y_{n}+\frac{1}{2} l_{1}\right), \\ k_{3} &=\Delta t f\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ l_{3} &=\Delta t g\left(t_{n}+\frac{1}{2} \Delta t, x_{n}+\frac{1}{2} k_{2}, y_{n}+\frac{1}{2} l_{2}\right), \\ k_{4} &=\Delta t f\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ l_{4} &=\Delta t g\left(t_{n}+\Delta t, x_{n}+k_{3}, y_{n}+l_{3}\right), \\ x_{n+1} &=x_{n}+\frac{1}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \\ y_{n+1} &=y_{n}+\frac{1}{6}\left(l_{1}+2 l_{2}+2 l_{3}+l_{4}\right) . \end{aligned} \nonumber$$

Boundary value problems

Shooting method

We consider the general ode of the form

$$\frac{d^{2} y}{d x^{2}}=f(x, y, d y / d x), \nonumber$$

with two-point boundary conditions $y(0)=A$ and $y(1)=B$. We will first formulate the ode as an initial value problem. We have

$$\frac{d y}{d x}=z, \quad \frac{d z}{d x}=f(x, y, z) . \nonumber$$

The initial condition $y(0)=A$ is known, but the second initial condition $z(0)=b$ is unknown. Our goal is to determine $b$ such that $y(1)=B$.

In fact, this is a root-finding problem for an appropriately defined function. We define the function $F=F(b)$ such that

$$F(b)=y(1)-B, \nonumber$$

where $y(1)$ is the numerical value obtained from integrating the coupled first-order differential equations with $y(0)=A$ and $z(0)=b$. Our root-finding routine will want to solve $F(b)=0$. (The method is called shooting because the slope of the solution curve for $y=y(x)$ at $x=0$ is given by $b$, and the solution hits the value $y(1)$ at $x=1$. This looks like pointing a gun and trying to shoot the target, which is $B$.)

To determine the value of $b$ that solves $F(b)=0$, we iterate using the secant method, given by

$$b_{n+1}=b_{n}-F\left(b_{n}\right) \frac{b_{n}-b_{n-1}}{F\left(b_{n}\right)-F\left(b_{n-1}\right)} . \nonumber$$

We need to start with two initial guesses for $b$, solving the ode for the two corresponding values of $y(1)$. Then the secant method will give us the next value of $b$ to try, and we iterate until $|y(1)-B|<$ tol, where tol is some specified tolerance for the error.