Numerical methods for ODE

As we have seen time ad time again a useful tool from calculus for designing numerical methods is Taylor’s formula

f(z) = f(a) + \frac{x-a}{1!} f'(a) + \frac{(x-1)^2}{2!} f^{(2)}(a) + \ldots + \frac{(x-1)^{n-1}}{(n-1)!} + R_n,

where the reminder is

\frac{(x-1)^{n}}{n!} f^{(n)}(\xi), \ \ \xi \in [a,x].

Alternatively, by instead expanding $f(z)$ around $z+h$, we find that

f(z+h) = f(z) + h f'(z) + \frac{h^2}{2!} f''(z) + \mathcal{O} (h^3).

Now if we are looking to solve an ordinary differential equation

\dot{u}(t) = f(u(t),t),\ \ u(0) = u_0,

we can use the above formulas to approximate the time derivative on the left

\frac{u(t+\Delta t)-u(t)}{\Delta t} \approx \dot{u}(t) = f(u(t),t).

Rearranging the equation allows us to compute \(u_n \equiv u({t_n}) = u(n \Delta t)\) according to

\begin{eqnarray*} u_0 &=& u_0, \\ u_{n+1} &=& u_n + \Delta t \, f(u_n,t_n), \ \ n= 0,1,\ldots \end{eqnarray*}

The above method is called forward Euler.

Accuracy and stability

It turns out that if we compute the solution up to some fixed time the error decreases with order of accuracy of the above method is one, or in other words the error is proportional to \(\Delta t\).

This might be surprising as the error in the approximation of one time step is actually proportional to \((\Delta t)^2\). Roughly speaking the explanation goes like this, if we compute to some time \(T\) with a time step \(\Delta\) we have to take \(N_t = \frac{T}{\Delta t}\) steps, if we commit an error of order \((\Delta t)^2\) in each of those steps the final error can be as large as \(T \Delta t\) and hence of order of accuracy one.

In any case, any numerical method for solving ODEs has an error expansion which, assuming there is a maximum allowable error, sets an upper limit on \(\Delta t\). This is the accuracy limit. To illustrate the stability limit we consider the following exercise.

Exercise

There is also an upper limit on \(\Delta t\) that is dictated by stability considerations. Consider the special case when \(f(u(t),t) = \lambda u(t), \ \ \lambda < 0\) and \(u_0 = 1\)

  1. Find the solution to

    \dot{u}(t) = \lambda u(t), \ \ u(0) = 1.
  2. What is the steady state? Is the solution of the continuous problem stable?

  3. Use Euler’s method to compute the approximate solution at the first few time steps. What is the stability limit on \(\Delta t\)?

  4. Repeat the analysis (applied to \(f(u(t),t) = \lambda u(t)\)) for backward Euler

    \begin{eqnarray*} u_0 &=& u_0, \\ u_{n+1} &=& u_n + \Delta t \, f(u_{n+1},t_{n+1}), \ \ n= 0,1,\ldots \end{eqnarray*}

    When can this method be useful? What is the additional complication compared to forward Euler when applied to an ODE with a more complicated right hand side?