Out Of The Ordinary

| 0 Comments

Several years ago we saw how to use the trapezium rule to approximate integrals[1]. This works by dividing the interval of integration into a set of equally spaced values, evaluating the function being integrated, or integrand, at each of them and calculating the area under the curve formed by connecting adjacent points with straight lines to form trapeziums.
This was an improvement over an even more rudimentary scheme which instead placed rectangles spanning adjacent values with heights equal to the values of the function at their midpoints to approximate the area. Whilst there really wasn't much point in implementing this since it offers no advantage over the trapezium rule, it is a reasonable first approach to approximating the solutions to another type of problem involving calculus; ordinary differential equations, or ODEs.

ODE To Joy

The simplest ODEs, known as first order ODEs of the first degree, define the derivative of one variable with respect to another as a function of them both, rather than just the latter, so that we can't solve them by direct integration. For example
\[ \frac{\mathrm{d} y}{\mathrm{d} x} = -x \times y \]
which is solved by
\[ y = c \times e^{-\frac12 \times x^2} \]
for any constant \(c\), which we can confirm by calculating its derivative using the chain rule
\[ \frac{\mathrm{d} y}{\mathrm{d} x} = c \times e^{-\frac12 \times x^2} \times -\tfrac12 \times 2 \times x = y \times -x \]
If we know that the solution passes through a specific point \(\left(x^\ast, y^\ast\right)\) then we can tie down the value of \(c\) with
\[ \begin{align*} y^\ast &= c \times e^{-\frac12 \times x^{\ast^2}}\\ c &= y^\ast \times e^{\frac12 \times x^{\ast^2}}\\ y &= y^\ast \times e^{-\frac12 \times \left(x^2 - x^{\ast^2}\right)} \end{align*} \]
There are several tricks for solving ODEs with specific forms. To solve our example we note that
\[ \frac{1}{y} \times \mathrm{d} y = -x \times \mathrm{d} x \]
and so integrating both sides yields
\[ \begin{align*} \ln y &= -\tfrac12 \times x^2 + c\\ y &= e^c \times e^{-\frac12 \times x^2} \end{align*} \]
This is known as separation of variables and can be used to solve ODEs of the form
\[ f(y) \times \frac{\mathrm{d} y}{\mathrm{d} x} = g(x) \]
Another is for ODEs of the form
\[ \frac{\mathrm{d} y}{\mathrm{d} x} = f\left(\frac{y}{x}\right) \]
which are known as homogeneous ODEs and can be solved with the substitution
\[ \begin{align*} y &= x \times z\\ \frac{\mathrm{d} y}{\mathrm{d} x} &= z + x \times \frac{\mathrm{d} z}{\mathrm{d} x}\\ f\left(\frac{y}{x}\right) &= f(z) \end{align*} \]
so that
\[ \begin{align*} x \times \frac{\mathrm{d} z}{\mathrm{d} x} &= f(z) - z\\ \frac{1}{x} \times \mathrm{d} x &= \frac{1}{f(z) - z} \times \mathrm{d} z \end{align*} \]
Note that the degree of a first order ODE is the power to which the derivative is raised so that the general form of a first order ODE of the \(n^\mathrm{th}\) degree is
\[ \left(\frac{\mathrm{d} y}{\mathrm{d} x}\right)^n = f(x, y) \]
Whilst there are other tricks for exactly solving first order ODEs of the first degree they only apply to a few different forms and so we must typically resort to numerical approximation.

The Euler Method

By Taylor's theorem we have
\[ y(x+\Delta x) = y(x) + \Delta x \times y^\prime(x) + O\left(\Delta x^2\right) = y(x) + \Delta x \times f(x, y(x)) + O\left(\Delta x^2\right) \]
Derivation 1: The Euler Method's Approximation Error
Firstly, denote the error in \(y_i\) as \(\epsilon_i\) so that, by Taylor's theorem, we have
\[ \begin{align*} f\left(x_{i-1}, y_{i-1}\right) &= f\left(x_{i-1}, y\left(x_{i-1}\right) + \epsilon_{i-1}\right)\\ &= f\left(x_{i-1}, y\left(x_{i-1}\right)\right) + O\left(\epsilon_{i-1}\right) \end{align*} \]
and therefore
\[ \begin{align*} y_i &= y_{i-1} + \Delta x \times f\left(x_{i-1}, y_{i-1}\right)\\ &= y\left(x_{i-1}\right) + O\left(\epsilon_{i-1}\right) + \Delta x \times f\left(x_{i-1}, y_{i-1}\right)\\ &= y\left(x_i\right) + O\left(\epsilon_{i-1}\right) + O\left(\Delta x \times \epsilon_{i-1}\right) + O\left(\Delta x^2\right) \end{align*} \]
This implies that
\[ O\left(\epsilon_i\right) = O\left(\epsilon_{i-1}\right) + O\left(\Delta x^2\right) \]
and, noting that \(\epsilon_0\) equals zero, we have by induction
\[ O\left(\epsilon_i\right) = i \times O\left(\Delta x^2\right) \]
Finally, given that \(n\) is proportional to the reciprocal of \(\Delta x\), the order of the final error is
\[ O\left(\epsilon_n\right) = O\left(\frac{1}{\Delta x}\right) \times O\left(\Delta x^2\right) = O\left(\Delta x\right) \]

Derivation 2: The Euler Method's Rounding Error
Adding rounding to the order of the error in \(y_i\) yields
\[ \begin{align*} y_i &= \left(y_{i-1} + \Delta x \times f\left(x_{i-1}, y_{i-1}\right)\right) \times (1+O(\epsilon))\\ &= y\left(x_i\right) + O\left(\epsilon_{i-1}\right) + O\left(\Delta x^2\right) + O(\epsilon) \end{align*} \]
implying that
\[ O\left(\epsilon_i\right) = O\left(\epsilon_{i-1}\right) + O\left(\Delta x^2\right) + O(\epsilon) \]
By induction we now have
\[ \begin{align*} O\left(\epsilon_n\right) &= n \times \left(O\left(\Delta x^2\right) + O(\epsilon)\right)\\ &= O\left(\Delta x\right) + O\left(\frac{\epsilon}{\Delta x}\right) \end{align*} \]
for sufficiently small \(\Delta x\), where \(y^\prime\) is an alternative notation for the derivative of \(y\) and \(O\left(\Delta x^2\right)\) is an error term of order \(\Delta x^2\).
We can consequently approximate \(y\left(x_n\right)\) by starting from \(\left(x_0, y_0\right)\) and iteratively calculating
\[ y_i = y_{i-1} + \Delta x \times f\left(x_{i-1}, y_{i-1}\right) \]
with
\[ \Delta x = \frac{x_n-x_0}{n} \]
so that
\[ x_i = x_{i-1} + \Delta x \]
This is known as Euler's method and results in
\[ y_n = y\left(x_n\right) + O\left(\Delta x\right) \]
as proven in derivation 1.

Unfortunately this error analysis assumes infinitely precise real arithmetic and, since we'll be using floating point numbers, we must take into account rounding errors, much as we did when considering reasonable values for the step \(\delta\) in finite difference approximations of the derivatives of functions[2].
Adding rounding to it yields
\[ y_n = y\left(x_n\right) + O\left(\Delta x\right) + O\left(\frac{\epsilon}{\Delta x}\right) \]
where \(\epsilon\) is the difference between one and the smallest floating point number greater than one, as proven in derivation 2. The value of \(n\) that minimises the error will therefore satisfy
\[ O(\Delta x) = O(\sqrt{\epsilon}) \]
and, for no reason other than pragmatism, might as well be
\[ \begin{align*} |\Delta x| &= \frac{\left|x_n-x_0\right|}{n} \approx \left(1 + |y_0|\right) \times \sqrt{\epsilon}\\ n & = \Bigg\lceil\frac{\left|x_n-x_0\right|}{\left(1 + |y_0|\right) \times \sqrt{\epsilon}}\Bigg\rceil \end{align*} \]
where \(\approx\) means approximately equal to and the odd shaped brackets represent the ceiling of the term between them, being the smallest integer that is greater than or equal to it.
Whilst this puts a not wholly unreasonable upper bound upon the number of steps, in practice we would generally be willing to accept greater approximation errors for the sake of computational efficiency by choosing to take a much smaller number of them, at least having some idea of the errors that such choices might add to the result.

Program 1 demonstrates Euler's method by using it to approximate the solution to our example ODE, plotting its results in red and the exact solution in black

Program 1: Our Example First Order ODE

and program 2 shows its errors for progressively doubled numbers of steps.

Program 2: First Order Approximation Errors

Here we can see that the error roughly halves as the number of steps doubles, at least once there are four or more steps so that \(\Delta x\) is less than one.

Bring On A Higher Order

Euler's method can easily be extended to first order vector ODEs such as
\[ \mathbf{y}^\prime(x) = \mathbf{f}(x, \mathbf{y}(x)) \]
where \(\mathbf{y}^\prime\) is the vector of the derivatives of the elements of \(\mathbf{y}\) and \(\mathbf{f}\) is a vector valued function, with
\[ \mathbf{y}_i = \mathbf{y}_{i-1} + \Delta x \times \mathbf{f}\left(x_{i-1}, \mathbf{y}_{i-1}\right) \]
This is particularly useful for approximating the solutions of \(n^\mathrm{th}\) order ODEs
\[ y^{(n)}(x) = f\left(x, y(x), y^\prime(x), y^{\prime\prime}(x), \dots, y^{(n-1)}(x)\right) \]
where \(y^{(n)}\) is the \(n^\mathrm{th}\) derivative of \(y\) with respect to \(x\). To do so we first define an \(n\) dimensional vector \(\mathbf{z}\) having the elements
\[ z_i(x) = y^{(i)}(x) \]
with the convention that \(y^{(0)}\) is \(y\), so that for \(i\) less than \(n-1\) we have
\[ z^\prime_i(x) = z_{i+1}(x) \]
We can consequently express an \(n^\mathrm{th}\) order ODE as
\[ \mathbf{z}^\prime(x) = \mathbf{f}(x, \mathbf{z}(x)) \]
where the elements of \(\mathbf{f}\) are
\[ 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} \]
For a worked example we shall use damped harmonic motion which occurs when a particle is acted upon by two forces, one toward the origin and proportional to the distance from it and the other in opposition to the velocity and proportional to it. The particle's acceleration therefore satisfies
\[ \begin{align*} F(t) &= -F_v \times x^\prime(t) - F_x \times x(t) = m \times x^{\prime\prime}(t)\\ x^{\prime\prime}(t) &= -\frac{F_v}{m} \times x^\prime(t) - \frac{F_x}{m} \times x(t) \end{align*} \]
There's a trick for solving ODEs of this form too. Specifically, given
\[ a \times x^{\prime\prime}(t) + b \times x^\prime(t) + c \times x(t) = 0 \]
make the substitution
\[ x(t) = e^{z \times t} \]
to yield
\[ a \times z^2 \times e^{z \times t} + b \times z \times e^{z \times t} + c \times e^{z \times t} = 0 \]
dividing through by \(e^{z \times t}\) yields a quadratic equation in \(z\) whose, possibly complex, roots \(z_0\) and \(z_1\) give us the general form of the solution
\[ x(t) = A \times e^{z_0 \times t} + B \times e^{z_1 \times t} \]
We can express the constants \(A\) and \(B\) in terms of the initial position and velocity by noting that
\[ \begin{align*} x(0) &= A \times e^{0} + B \times e^{0} = A + B\\ x^\prime(0) &= A \times z_0 \times e^0 + B \times z_1 \times e^0 = A \times z_0 + B \times z_1 \end{align*} \]
and solving for them.

Rather than use this trick to solve a specific ODE of this form, we shall take the rather easier route of starting from a valid solution and deducing the parameters of the ODE that it solves
\[ \begin{align*} x(t) &= e^{-t} \times \left(\cos(10 \times t) + \tfrac{1}{10} \times \sin(10 \times t)\right)\\ x^\prime(t) &= -e^{-t} \times \left(\cos(10 \times t) + \tfrac{1}{10} \times \sin(10 \times t)\right) + e^{-t} \times \left(-10 \times \sin(10 \times t) + \tfrac{1}{10} \times 10 \times \cos(10 \times t)\right)\\ &= -x(t) + e^{-t} \times \left(-10 \times \sin(10 \times t) + \cos(10 \times t)\right)\\ x^{\prime\prime}(t) &= -x^\prime(t) -e^{-t} \times \left(-10 \times \sin(10 \times t) + \cos(10 \times t)\right) + e^{-t} \times \left(-100 \times \cos(10 \times t) - 10 \times \sin(10 \times t)\right)\\ &= -x^\prime(t) - \left(x^\prime(t) + x(t)\right) - 100 \times x(t)\\ &= -2 \times x^\prime(t) - 101 \times x(t) \end{align*} \]
The initial conditions are therefore
\[ \begin{align*} x(0) &= e^0 \times \left(\cos 0 + \tfrac{1}{10} \times \sin 0\right) = 1\\ x^\prime(0) &= -x(0) + e^0 \times \left(-10 \times \sin 0 + \cos 0\right) = -1 + 1 = 0 \end{align*} \]
The Euler method approximates the solution of this ODE with
\[ \begin{align*} x_i &= x_{i-1} + \Delta t \times x^\prime_{i-1}\\ x^\prime_i &= x^\prime_{i-1} + \Delta t \times \left(-2 \times x^\prime_{i-1} - 101 \times x_{i-1}\right) \end{align*} \]
starting from \(x_0\) equal to one and \(x^\prime_0\) equal to zero. Program 3 plots its results in red, showing that it has the same decaying oscillatory behaviour as the exact solution, plotted in black, but consistently overshoots its minima and maxima.

Program 3: Our Example Second Order ODE

Although the error is significantly larger than it was for our first order example, after a sufficient number of steps have been taken it once again roughly halves as they double, as demonstrated by program 4.

Program 4: Second Order Approximation Errors

If you have the patience then you can increase the value of n in program 3 by a factor of ten to see for yourself how much more accurate the Euler method becomes as the number of steps increases.

ak.eulerODE

When it comes to our implementation ak.eulerODE, which is defined in EulerODE.js, we initialise it with the step length rather than the number of steps, as shown by listing 1, since the former is directly related to the approximation error.

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

 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.eulerODE');
  }
  n = ak.ceil(Math.abs(x1-x0)/dx);
  if(n>ak.INT_MAX) throw new Error('too many steps in ak.eulerODE');
  return ak.nativeType(y0)===ak.NUMBER_T ? numberEulerODE(f, n, x0, x1, y0)
                                         : generalEulerODE(f, n, x0, x1, y0);
 };
};

Here we're checking that the right hand side of the ODE is a function and that the step length is a non-zero number before returning an approximation that first checks that the bounds and step length are valid.

In the univariate case when y0 is a number we forward to the numberEulerODE function which is implemented in listing 2

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

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

and demonstrated by program 5 which applies it to the first order example ODE.

Program 5: A First Order ak.eulerODE

Otherwise we forward to the generalEulerODE helper function, defined in listing 3, which supports any type that has arithmetic operators defined the ak library.

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

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

In particular, we can use our ak.vector type to approximate higher order ODEs, as used for the second order example ODE by program 6.

Program 6: A Second Order ak.eulerODE

References

[1] Trapezium Artistry, www.thusspakeak.com, 2015.

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

Leave a comment