Finding The Middle Ground


Last time[1] we saw how we can use Euler's method to approximate the solutions of ordinary differential equations, or ODEs, which define the derivative of one variable with respect to another as a function of them both, so that they cannot be solved by direct integration. Specifically, it uses Taylor's theorem to estimate the change in the first variable that results from a small step in the second, iteratively accumulating the results for steps of a constant length to approximate the value of the former at some particular value of the latter.
Unfortunately it isn't very accurate, yielding an accumulated error proportional to the step length, and so this time we shall take a look at a way to improve it.

A Cheerful Symmetry

Given the first order ODE
\[ y^\prime(x) = f\left(x, y(x)\right) \]
where \(y^\prime\) is the derivative of \(y\) with respect to \(x\), we have from Taylor's theorem
\[ \begin{align*} y(x + \Delta x) &= y(x) + \Delta x \times y^\prime(x) + O\left(\Delta x^2\right)\\ &= y(x) + \Delta x \times f\left(x, y(x)\right) + O\left(\Delta x^2\right) \end{align*} \]
where \(O\left(\Delta x^2\right)\) is an error term of order \(\Delta x^2\). If we instead take a lead from the symmetric finite difference[2] approximation of the derivative with
\[ \begin{align*} y\left(x - \tfrac12 \times \Delta x\right) &= y(x) - \tfrac12 \times \Delta x \times y^\prime(x) + \tfrac14 \times \Delta x^2 \times y^{\prime\prime}(x) + O\left(\Delta x^3\right)\\ y\left(x + \tfrac12 \times \Delta x\right) &= y(x) + \tfrac12 \times \Delta x \times y^\prime(x) + \tfrac14 \times \Delta x^2 \times y^{\prime\prime}(x) + O\left(\Delta x^3\right) \end{align*} \]
\[ \begin{align*} y\left(x + \tfrac12 \times \Delta x\right) - y\left(x - \tfrac12 \times \Delta x\right) &= \Delta x \times y^\prime(x) + O\left(\Delta x^3\right)\\ y\left(x + \tfrac12 \times \Delta x\right) = y\left(x - \tfrac12 \times \Delta x\right) &+ \Delta x \times y^\prime(x) + O\left(\Delta x^3\right) \end{align*} \]
then a shift in \(x\) of \(\tfrac12 \times \Delta x\) yields
\[ \begin{align*} y\left(x + \Delta x\right) &= y\left(x\right) + \Delta x \times y^\prime\left(x + \tfrac12 \times \Delta x\right) + O\left(\Delta x^3\right)\\ &= y\left(x\right) + \Delta x \times f\left(x + \tfrac12 \times \Delta x, y\left(x + \tfrac12 \times \Delta x\right)\right) + O\left(\Delta x^3\right) \end{align*} \]
giving us an extra order of magnitude of accuracy.

Now I'm sure you've spotted the flaw in this scheme; we don't know the value of \(y\) at \(x + \tfrac12 \times \Delta x\). We do, however, know how to approximate it, giving us
\[ \begin{align*} f\left(x + \tfrac12 \times \Delta x, y\left(x + \tfrac12 \times \Delta x\right)\right) &= f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12 \times \Delta x \times f\left(x, y(x)\right) + O\left(\Delta x^2\right)\right)\\ &= f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12 \times \Delta x \times f\left(x, y(x)\right)\right) + O\left(\Delta x^2\right) \end{align*} \]
with the last step coming from applying Taylor's theorem to \(f\). This yields the midpoint method for approximating the change in \(y\)
\[ \begin{align*} y\left(x + \Delta x\right) &= y\left(x\right) + \Delta x \times \left(f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12\Delta x \times f\left(x, y(x)\right)\right) + O\left(\Delta x^2\right)\right) + O\left(\Delta x^3\right)\\ &= y\left(x\right) + \Delta x \times f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12 \times \Delta x \times f\left(x, y(x)\right)\right) + O\left(\Delta x^3\right) \end{align*} \]
We can turn this into an iterative algorithm with
\[ \begin{align*} \Delta x &= \frac{x_n-x_0}{n}\\ y_i &= y_{i-1} + \Delta x \times f\left(x_{i-1} + \tfrac12 \times \Delta x, y_{i-1} + \tfrac12 \times \Delta x \times f\left(x_{i-1}, y_{i-1}\right)\right) \end{align*} \]
\[ y_n = y\left(x_n\right) + O\left(\Delta x^2\right) \]
by the same argument that we made for the Euler method's cumulative error.

Similarly following our previous argument, taking rounding errors into account yields
\[ y_n = y\left(x_n\right) + O\left(\Delta x^2\right) + O\left(\frac{\epsilon}{\Delta x}\right) \]
where \(\epsilon\) is the difference between one and the smallest floating number greater than one, representing the precision of floating point arithmetic. The optimal choice of \(\Delta x\) will consequently satisfy
\[ O(\Delta x) = O\left(\epsilon^\frac13\right) \]
and, following the same justification of pragmatism, we can guess at an upper bound for \(n\) of
\[ \begin{align*} |\Delta x| &= \frac{\left|x_n-x_0\right|}{n} \approx \left(1 + |y_0|\right) \times \epsilon^\frac13\\ n & = \Bigg\lceil\frac{\left|x_n-x_0\right|}{\left(1 + |y_0|\right) \times \epsilon^\frac13}\Bigg\rceil \end{align*} \]
where the vertical bars represent the magnitude of the term between them and the strangely shaped brackets represent the ceiling of the term that they enclose, being the smallest integer that is greater than it, although we should again almost certainly prefer a much smaller value for the sake of efficiency.

Two Worked Examples

To demonstrate the midpoint method we shall use the same example ODEs that we did last time. Firstly
\[ y^\prime(x) = -x \times y(x) \]
which, given an initial condition of \(y(0)\) equal to one, has the solution
\[ y(x) = e^{-\frac12 \times x^2} \]
Program 1 plots the exact solution in black, the midpoint method approximation in red and the Euler method approximation in green.

Program 1: Our Example First Order ODE

Clearly the midpoint method is significantly more accurate that Euler's method, which program 2 confirms is a consequence of its quadratic, rather than linear, approximation error by demonstrating that it eventually reduces by a factor of four, rather than two, as the number of steps doubles.

Program 2: First Order Approximation Errors

Our second example was the second order ODE
\[ x^{\prime\prime}(t) = -2 \times x^\prime(t) - 101 \times x(t) \]
which, for the initial conditions
\[ \begin{align*} x(0) &= 1\\ x^\prime(0) &= 0 \end{align*} \]
is solved by
\[ x(t) = e^{-t} \times \left(\cos(10 \times t) + \tfrac{1}{10} \times \sin(10 \times t)\right) \]
Recall that the trick that we used to apply the Euler method to this example was to recast it as a pair of first order ODEs for
\[ \begin{align*} z_0(t) &= x(t)\\ z_1(t) &= x^\prime(t) \end{align*} \]
\[ \begin{align*} z^\prime_0(t) &= x^\prime(t) = z_1(t)\\ z^\prime_1(t) &= x^{\prime\prime}(t) = -2 \times x^\prime(t) - 101 \times x(t) = -2 \times z_1(t) - 101 \times z_0(t) \end{align*} \]
Program 3 uses this same trick and again plots the exact solution in black, the midpoint approximation in red and the Euler approximation in green.

Program 3: Our Example Second Order ODE

Here we also find that the midpoint method is far more accurate than Euler's method and program 4 demonstrates that its error is similarly quadratic.

Program 4: Second Order Approximation Errors


The implementation of the midpoint method is unsurprisingly similar to that of ak.eulerODE, as shown by ak.midpointODE in listing 1.

Listing 1: ak.midpointODE
ak.midpointODE = function(f, dx) {
 if(ak.nativeType(f)!==ak.FUNCTION_T) {
  throw new Error('invalid function in ak.midpointODE');
 dx = Math.abs(dx);
 if(isNaN(dx) || dx===0) {
  throw new Error('invalid step length in ak.midpointODE');

 return function(x0, x1, y0) {
  var n;
  if(ak.nativeType(x0)!==ak.NUMBER_T || !isFinite(x0)
  || ak.nativeType(x1)!==ak.NUMBER_T || !isFinite(x1)) {
   throw new Error('invalid interval in ak.midpointODE');
  n = ak.ceil(Math.abs(x1-x0)/dx);
  if(n>ak.INT_MAX) throw new Error('too many steps in ak.midpointODE');
  return ak.nativeType(y0)===ak.NUMBER_T ? numberMidpointODE(f, n, x0, x1, y0)
                                         : generalMidpointODE(f, n, x0, x1, y0);

This follows the same approach of first checking that f is a function and the step length dx is a non-zero number before defining a function to calculate the approximation. After checking that its start and end points are finite numbers it calculates the number of steps and then forwards to either the numberMidpointODE or the generalMidpointODE helper function depending upon whether the starting value is a number or not.
The latter case is supported since there is nothing in the definition of an ODE that requires \(y\) to be a scalar, just that it supports the required arithmetic operations.
The former is given in listing 2

Listing 2: The Scalar Implementation
function numberMidpointODE(f, n, x0, x1, y0) {
 var dx = (x1-x0)/n;
 var dx2 = dx/2;
 var i;

 for(i=0;i<n;++i) {
  y0 += dx * f(x0+dx2, y0+dx2*f(x0, y0));
  x0 += dx;
 return y0;

and is demonstrated by program 5.

Program 5: A First Order ak.midpointODE

The latter is given in listing 3

Listing 3: The General Implementation
function generalMidpointODE(f, n, x0, x1, y0) {
 var dx = (x1-x0)/n;
 var dx2 = dx/2;
 var i;

 for(i=0;i<n;++i) {
  y0 = ak.add(y0, ak.mul(dx, f(x0+dx2, ak.add(y0, ak.mul(dx2, f(x0, y0))))));
  x0 += dx;
 return y0;

and demonstrated by program 6.

Program 6: A Second Order ak.midpointODE

Here we're using our ak.vector type to represent the pair of variables that we used to transform the second order ODE into two first order ODEs, exploiting the fact that vector operations yield identical results to those of the scalar operations upon them.

Next time we shall see how we can generalise the midpoint method to yield higher orders of accuracy, but until then you can find the implementation of ak.midpointODE in MidpointODE.js.


[1] Out Of The Ordinary,, 2021.

[2] You're Going To Have To Think! Why Finite Differences Won't Cure Your Calculus Blues,, 2014.

Leave a comment