A Kutta Above The Rest

| 0 Comments

We have recently been looking at ordinary differential equations, or ODEs, which relate the derivatives of one variable with respect to another to them with a function so that we cannot solve them with plain integration. Whilst there are a number of tricks for solving such equations if they have very specific forms, we typically have to resort to approximation algorithms such as the Euler method[1], with first order errors, and the midpoint method[2], with second order errors.
In this post we shall see that these are both examples of a general class of algorithms that can be accurate to still greater orders of magnitude.

A Fourth Order Approximation

We have specifically sought to approximate the solutions of first order ODEs of the form
\[ y^\prime(x) = f(x, y) \]
where \(y^\prime\) is the derivative of \(y\) with respect to \(x\), noting that \(n^\mathrm{th}\) order ODEs can be expressed as first order ODEs of \(n\) dimensional vectors with
\[ \mathbf{z}^\prime(x) = \mathbf{f}(x, \mathbf{z}) \]
where the \(i^\mathrm{th}\) element of \(\mathbf{z}\) is the \(i^\mathrm{th}\) derivative of \(y\) and
\[ f_i\left(x, \mathbf{z}(x)\right) = \begin{cases} f\left(x, \mathbf{z}(x)\right) & i = n-1\\ z_{i+1}(x) & \text{otherwise} \end{cases} \]
Exploiting Taylor's theorem Euler's method approximates the change in \(y\) resulting from a change in \(x\) with
\[ y(x + \Delta x) = y(x) + \Delta x \times f(x, y(x)) + O\left(\Delta x^2\right) \]
where \(O\left(\Delta x^2\right)\) is an error of order \(\Delta x^2\), and the midpoint method does so with
\[ y(x + \Delta x) = 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) \]
Defining
\[ \begin{align*} k_1 &= f(x, y(x))\\ k_2 &= f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12 \times \Delta x \times k_1\right)\\ k_3 &= f\left(x + \tfrac12 \times \Delta x, y(x) + \tfrac12 \times \Delta x \times k_2\right)\\ k_4 &= f\left(x + \Delta x, y(x) + \Delta x \times k_3\right) \end{align*} \]
we can do still better with
\[ y(x + \Delta x) = y(x) + \tfrac16 \times \Delta x \times (k_1 + 2 \times k_2 + 2 \times k_3 + k_4) + O\left(\Delta x^5\right) \]
and so the iteration
\[ \begin{align*} \Delta x &= \frac{x_n-x_0}{n}\\ k_{i,1} &= f(x_{i-1}, y_{i-1})\\ k_{i,2} &= f\left(x_{i-1} + \tfrac12 \times \Delta x, y_{i-1} + \tfrac12 \times \Delta x \times k_{i,1}\right)\\ k_{i,3} &= f\left(x_{i-1} + \tfrac12 \times \Delta x, y_{i-1} + \tfrac12 \times \Delta x \times k_{i,2}\right)\\ k_{i,4} &= f\left(x_{i-1} + \Delta x, y_{i-1} + \Delta x \times k_{i,3}\right)\\ y_i &= y_{i-1} + \tfrac16 \times \Delta x \times (k_{i,1} + 2 \times k_{i,2} + 2 \times k_{i,3} + k_{i,4}) \end{align*} \]
yields
\[ y_n = y\left(x_n\right) + O\left(\Delta x^4\right) \]

Worked Examples

Program 1 uses this to approximate the solution to
\[ y^\prime(x) = -x \times y(x) \]
given by
\[ y(x) = e^{-\frac12 \times x^2} \]
if \(y(0)\) equals one.

Program 1: Our Example First Order ODE

Here the errors are too small to allow us to visually distinguish between the exact solution, plotted in black, and the approximation, plotted in red, and so program 2 prints the differences between them for repeatedly doubled numbers of steps.

Program 2: First Order Approximation Errors

Whilst the errors reduce by the expected approximate factor of sixteen once the number of steps is sufficiently large, an important point to note is that the error increases on the final step because the floating point rounding errors start to dominate the approximation errors.

Program 3 applies the approximation to the second order ODE
\[ x^{\prime\prime}(t) = -2 \times x^\prime(t) - 101 \times x(t) \]
with initial conditions of
\[ \begin{align*} x(0) &= 1\\ x^\prime(0) &= 0 \end{align*} \]
which has the exact solution
\[ x(t) = e^{-t} \times \left(\cos(10 \times t) + \tfrac{1}{10} \times \sin(10 \times t)\right) \]
plotting the former in red and the latter in black.

Program 3: Our Example Second Order ODE

Once again the differences are too small to see in this chart and so program 4 prints them for doubling numbers of steps, showing that the errors ultimate reduce by the expected order of magnitude.

Program 4: Second Order Approximation Errors

Further Generalisation

This scheme is often referred to as the Runge-Kutta method[3][4], although it's more accurate to call it a Runge-Kutta method since it is one of many, including the Euler and midpoint methods, that take the form
\[ \begin{align*} k_i &= \begin{cases} f(x, y(x)) & i=1\\ f\left(x + c_i \times \Delta x, y(x) + \Delta x \times \left(\sum_{j=1}^{i-1} a_{i,j} \times k_j\right)\right) & \text{otherwise} \end{cases}\\ y(x + \Delta x) &\approx y(x) + \Delta x \times \sum_{i=1}^n b_i \times k_i \end{align*} \]
where \(\sum\) is the summation sign. Specific examples are typically represented with a Butcher tableau[5]
\[ \begin{array}{c|ccccc} 0\\ c_2 & a_{2,1}\\ c_3 & a_{3,1} & a_{3,2}\\ \vdots & \vdots & \vdots & \ddots\\ c_n & a_{n,1} & a_{n,2} & \dots &a_{n,n-1}\\ \hline & b_1 & b_2 & \dots & b_{n-1} & b_n \end{array} \]
which has a strictly lower triangular matrix of the coefficients \(a_{i,j}\), meaning that only the elements in columns with lower indices than the rows are non-zero, a vector of the weights \(b_i\) and another of the nodes \(c_i\).
The method described above has the tableau
\[ \begin{array}{c|cccc} 0\\ \tfrac12 & \tfrac12\\ \tfrac12 & 0 & \tfrac12\\ 1 & 0 & 0 & 1\\ \hline & \tfrac16 & \tfrac13 & \tfrac13 & \tfrac16 \end{array} \]
whilst the Euler method is given by
\[ \begin{array}{c|c} 0\\ \hline & 1 \end{array} \]
and the midpoint method by
\[ \begin{array}{c|cc} 0\\ \tfrac12 & \tfrac12\\ \hline & 0 & 1 \end{array} \]
Note that the elements of a Butcher tableau are not arbitrary but must be carefully chosen if it is to represent a valid method for approximating the solutions of ODEs. Furthermore, whilst it is frequently the case that the nodes of a Butcher tableau satisfy
\[ c_i = \sum_{j=1}^{i-1} a_{i,j} \]
it is not a requirement, although it is necessary that the sum of the weights equals one.

The Implementation

The ak.rungeKutta function defined in listing 1 creates a Runge-Kutta method from the derivative f, a step length dx and the coefficients, weights and nodes a, b and c, first validating the inputs before returning a function that implements the specified method.

Listing 1: ak.rungeKuttaODE
ak.rungeKuttaODE = function(f, dx, a, b, c) {
 var tc = ak.nativeType(c);
 var s, alb, k;

 dx = Math.abs(dx);
 checkArgTypes(f, dx, a, b, c, tc);
 checkTableau(a, b, c, tc);

 if(tc===ak.UNDEFINED_T) c = makeNodes(a);
 s = b.reduce(function(s,x){return s+x;}, 0);
 if(!isFinite(s) || s===0) {
  throw new Error('invalid total weight in ak.rungeKuttaODE');
 }

 alb = makeLB(a);
 a = a.map(function(r){return r.slice(0);});
 b = b.slice(0);
 if(tc!==ak.UNDEFINED_T) c = c.slice(0);
 k = new Array(a.length+1);

 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.rungeKuttaODE');
  }
  n = ak.ceil(Math.abs(x1-x0)/dx);
  if(n>ak.INT_MAX) throw new Error('too many steps in ak.rungeKuttaODE');
  return ak.nativeType(y0)===ak.NUMBER_T
       ? numberRungeKuttaODE(f, n, x0, x1, y0, a, b, c, s, alb, k)
       : generalRungeKuttaODE(f, n, x0, x1, y0, a, b, c, s, alb, k);
 };
};

Note that rather than use the step length directly the returned function uses the number of steps n that yields the closest step length that is not less than it. It also assumes that the weights are relative and so will use their sum s to normalise them by dividing the final weighted sum by it.

The first validity test of the arguments, given in listing 2, simply checks that f is a function, dx is a non-zero number and that a, b and c are array of the correct sizes, albeit allowing the latter to be undefined which we shall take to imply that its elements are the sums of those in the rows of a.

Listing 2: Checking The Arguments
function checkArgTypes(f, dx, a, b, c, tc) {
 if(ak.nativeType(f)!==ak.FUNCTION_T) {
  throw new Error('invalid function in ak.rungeKuttaODE');
 }
 if(isNaN(dx) || dx===0) {
  throw new Error('invalid step length in ak.rungeKuttaODE');
 }
 if(ak.nativeType(a)!==ak.ARRAY_T) {
  throw new Error('invalid matrix in ak.rungeKuttaODE');
 }
 if(ak.nativeType(b)!==ak.ARRAY_T) {
  throw new Error('invalid weights in ak.rungeKuttaODE');
 }
 if(b.length!==a.length+1) {
  throw new Error('matrix/weights size mismatch in ak.rungeKuttaODE');
 }

 if(tc!==ak.UNDEFINED_T) {
  if(tc!==ak.ARRAY_T) {
   throw new Error('invalid nodes in ak.rungeKuttaODE');
  }
  if(c.length!==a.length) {
   throw new Error('matrix/nodes size mismatch in ak.rungeKuttaODE');
  }
 }
}

Here we're assuming that a doesn't include the first row of empty coefficients and c doesn't include the first unused element of the nodes which we shall have to take into account during the calculation.

The checkTableau function defined in listing 3 ensures that the elements of a, b and c are finite numbers, additionally checking that the former is strictly lower triangular.

Listing 3: Checking The Tableau
function checkTableau(a, b, c, tc) {
 var ai, i, j;

 for(i=0;i<a.length;++i) {
  ai = a[i];

  if(ak.nativeType(b[i])!==ak.NUMBER_T || !isFinite(b[i])) {
   throw new Error('invalid weight in ak.rungeKuttaODE');
  }
  if(ak.nativeType(ai)!==ak.ARRAY_T || ai.length!==i+1) {
   throw new Error('invalid matrix row in ak.rungeKuttaODE');
  }
  for(j=0;j<=i;++j) {
   if(ak.nativeType(ai[j])!==ak.NUMBER_T || !isFinite(ai[j])) {
    throw new Error('invalid matrix element in ak.rungeKuttaODE');
   }
  }
  if(tc!==ak.UNDEFINED_T
  && (ak.nativeType(c[i])!==ak.NUMBER_T || !isFinite(c[i]) || c[i]===0)) {
   throw new Error('invalid node in ak.rungeKuttaODE');
  }
 }
 if(ak.nativeType(b[i])!==ak.NUMBER_T || !isFinite(b[i])) {
  throw new Error('invalid weight in ak.rungeKuttaODE');
 }
}

If c is defined then we additionally require that its elements are non-zero. Note that this doesn't check that the tableau represents a valid Runge-Kutta method because, frankly, that's really really difficult.

If the nodes are undefined then the makeNodes function defined in listing 4 creates nodes that are the sums of the elements in the rows of a, additionally checking that they are finite and non-zero.

Listing 4: Making The Nodes
function makeNodes(a) {
 var c = new Array(a.length);
 var i;

 for(i=0;i<a.length;++i) {
  c[i] = a[i].reduce(function(s,x){return s+x;}, 0);
  if(!isFinite(c[i]) || c[i]===0) {
   throw new Error('invalid node in ak.rungeKuttaODE');
  }
 }
 return c;
}

Noting that Runge-Kutta methods frequently have leading zeros in their matrices of coefficients the makeLB function defined in listing 5 records the indices of the first non-zero elements in each of their rows so that we can skip over zero terms.

Listing 5: The First Non-Zero Matrix Row Elements
function makeLB(a) {
 var alb = new Array(a.length);
 var ai, i, j;

 for(i=0;i<a.length;++i) {
  ai = a[i];

  j = 0;
  while(j<=i && ai[j]===0) ++j;
  alb[i] = j;
 }
 return alb;
}

Finally, we come to the implementation of Runge-Kutta methods. If y0 is a number then we use the inbuilt JavaScript arithmetic operators to implement the steps of the iterations, as shown by listing 6.

Listing 6: The Scalar Implementation
function numberRungeKuttaODE(f, n, x0, x1, y0, a, b, c, s, alb, k) {
 var kn = k.length;
 var dx = (x1-x0)/n;
 var i, j, l, dy, aj, kj;

 for(i=0;i<n;++i) {
  k[0] = f(x0, y0);
  dy = b[0]*k[0];
  for(j=1;j<kn;++j) {
   aj = a[j-1];
   l = alb[j-1];

   kj = y0;
   if(l!==j) {
    kj += k[l]*aj[l]*dx;
    while(++l<j) kj += k[l]*aj[l]*dx;
   }
   kj = f(x0+c[j-1]*dx, kj);
   dy += b[j]*kj;
   k[j] = kj;
  }
  x0 += dx;
  y0 += dy*dx/s;
 }
 return y0;
}

Otherwise, if y0 is one of our arithmetic types then we use our overloaded arithmetic operations to do the same, as shown in listing 7.

Listing 7: The General Implementation
function generalRungeKuttaODE(f, n, x0, x1, y0, a, b, c, s, alb, k) {
 var kn = k.length;
 var dx = (x1-x0)/n;
 var i, j, l, dy, aj, kj;

 for(i=0;i<n;++i) {
  k[0] = f(x0, y0);
  dy = ak.mul(b[0], k[0]);
  for(j=1;j<kn;++j) {
   aj = a[j-1];
   l = alb[j-1];

   kj = y0;
   if(l!==j) {
    kj = ak.add(kj, ak.mul(k[l], aj[l]*dx));
    while(++l<j) if(aj[l]!==0) kj = ak.add(kj, ak.mul(k[l], aj[l]*dx));
   }
   kj = f(x0+c[j-1]*dx, kj);
   if(b[j]!==0) dy = ak.add(dy, ak.mul(b[j], kj));
   k[j] = kj;
  }
  x0 += dx;
  y0 = ak.add(y0, ak.mul(dy, dx/s));
 }
 return y0;
}

Since the parameters of a Runge-Kutta method must be fine tuned, rather than expect the user to provide them, the ak library provides implementations of specific examples of various orders, as shown by listing 8.

Listing 8: Classic Runge-Kutta Algorithms
ak.eulerRK1ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [], [1]);
};

ak.midpointRK2ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[0.5]], [0,1]);
};

ak.heunRK2ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[1]], [0.5,0.5]);
};

ak.ralstonRK2ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[2/3]], [0.25,0.75]);
};

ak.kuttaRK3ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[0.5],[-1,2]], [1,4,1]);
};

ak.heunRK3ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[1/3],[0,2/3]], [1,0,3]);
};

ak.ralstonRK3ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[0.25],[0,0.75]], [2,3,4]);
};

ak.classicRK4ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[0.5],[0,0.5],[0,0,1]], [1,2,2,1]);
};

ak.kuttaRK4ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx, [[1/3],[-1/3,1],[1,-1,1]], [1,3,3,1]);
};

Program 5 compares first, second and fourth order approximations of our example second order ODE, plotting the exact solution in black and the approximations in green, blue and red respectively.

Program 5: A Trio Of Runge-Kutta Methods

Note that there are Runge-Kutta algorithms of even higher orders of accuracy, as illustrated by listing 9.

Listing 9: Some Fifth Order Algorithms
ak.fehlbergRK5ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx,
  [[1/4]
   [3/32,       9/32],
   [1932/2197, -7200/2197,  7296/2197],
   [439/216,   -8,          3680/513,  -845/4104],
   [-8/27,      2,         -3544/2565,  1859/4104, -11/40]],
  [33440,0,146432,142805,-50787,10260]);
};

ak.cashKarpRK5ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx,
  [[1/5],
   [3/40,        9/40],
   [3/10,       -9/10,     6/5],
   [-11/54,      5/2,     -70/27,     35/27],
   [1631/55296,  175/512,  575/13824, 44275/110592, 253/4096]],
  [9361,0,38500,20125,0,27648]);
};

ak.dormandPrinceRK5ODE = function(f, dx) {
 return ak.rungeKuttaODE(f, dx,  
  [[1/5],
   [3/40,        9/40],
   [44/45,      -56/15,      32/9],
   [19372/6561, -25360/2187, 64448/6561, -212/729],
   [9017/3168,  -355/33,     46732/5247,  49/176, -5103/18656]],
  [12985,0,64000,92750,-45927,18656]);
};

These aren't generally used as they are but play an important part in a somewhat more robust scheme that we shall take a look at next time. In the meantime you can find the implementations of all of these functions in ak.rungeKuttaODE.js.

References

[1] Out Of The Ordinary, www.thusspakeak.com, 2021.

[2] Finding The Middle Ground, www.thusspakeak.com, 2021.

[3] Runge, C., Über die numerische Auflösung von Differentialgleichungen, Mathematische Annalen, Vol. 46, No. 2, Springer, 1895.

[4] Kutta, M., Beitrag zur näherungsweisen Integration totaler Differentialgleichungen, Zeitschrift für Mathematik und Physik, Vol. 46, 1901.

[5] Butcher, J., Coefficients for the study of Runge-Kutta integration processes, Journal of the Australian Mathematical Society, Vol. 3, No. 2, 1963.

Leave a comment