Taylor's method for ODE initial value problems in Julia

Jorge A. Pérez and Luis Benet

Instituto de Ciencias Físicas

Universidad Nacional Autónoma de México

Thanks: DGAPA-UNAM: IG-100616, and CONACYT.


Taylor's method provides a high accuracy method for the integration of initial value problems.

The error in each step of the integration can be made close to the round-off error.

99942 Apophis


Credit: Osservatorio Astronomico Sormano via Wikipedia

Semi-major axis: 0.92228 AU; eccentricity: 0.19108

Aphelion: 1.09851 AU; perihelion: 0.74605 AU

Important close approach April 13th, 2029.

Taylor's method

Consider the initial value problem

\begin{equation*} \dot{x} = f(x, t), \qquad x(0) = x_0. \end{equation*}

Approximate locally the solution by a (Taylor) polynomial in $t$

\begin{equation*} x(t) = \sum_{n=0}\, x_{[n]}\, t^n. \end{equation*}

Using the Taylor expansion of $f(x(t),t)$ around $t=0$

\begin{equation*} f(x(t), t) = \sum_{n=0}\, f_{[n]}\, t^n, \end{equation*}

we arrive at the recursion relation

\begin{equation*} x_{[n+1]} = \frac{f_{[n]}}{n+1}, \qquad x_{[0]} = x_0. \end{equation*}

We use TaylorSeries.jl to obtain the expansion of $f(x(t),t)$ around a given value of $t$.

To assure convergence of the series, the coefficients $x_{[n]}$ must be computed to a large enough order.

The time step is determined through the last computed coefficients, $ | x_{[n]} h^n | \leq \epsilon$, then

\begin{equation*} h \leq \min\Big\{ \Big(\frac{\epsilon}{x_{[n]}}\Big)^{1/n}, \,\Big(\frac{\epsilon}{x_{[n-1]}}\Big)^{1/(n-1)}\Big\}. \end{equation*}

Note that higher order of the Taylor expansion provides a larger time step.

A blackboard example

Let's consider the equation

\begin{equation*} \dot{x} = -x, \qquad x(0) = 1 \end{equation*}

whose solution is $x(t) = \exp(-t)$.

Clearly, $x_{[0]}=1$ and $f_{[n]} = -x_{[n]}$. Then,

\begin{eqnarray*} x_{[k+1]} & = & \frac{f_{[k]}}{k+1} = -\frac{x_{[k]}}{k+1} = (-1)^2\frac{x_{[k-1]}}{k (k+1)}\\ & = & \dots = (-1)^{k+1} \frac{x_{[0]}}{(k+1)!} = \frac{(-1)^{k+1}}{(k+1)!}, \end{eqnarray*}

which is the (k+1)-Taylor coefficient of $x(t) = \exp(-t)\,$.

Jet transport

Jet transport is a tool that allows the propagation under the flow of a small neighborhood in phase space around a given initial condition, instead of propagating single initial conditions only.

To compute the propagation of $\mathbf{x}_0+\delta\mathbf{x}$, where $\delta\mathbf{x}$ are independent small displacements in phase space, one has to solve high-order variational equations. The idea is to treat $\mathbf{x}_0+\delta\mathbf{x}$ as a (truncated) polynomial in the $\delta\mathbf{x}$ variables.

Jet transport works with any ordinary differential equations integrator, provided the computations can be done using polynomial algebra.

An example

From D. Pérez-Palau et al, Celest. Mech, Dyn. Astron. 123, 239-262 (2015).

Let us consider the differential equations for the harmonic oscillator:

\begin{eqnarray*} \dot{x} & = & y, \\ \dot{y} & = & -x, \end{eqnarray*}

with the initial condition $\mathbf{x}_0=[x_0, y_0]^T$.

Euler's method corresponds to:

\begin{equation*} \mathbf{x}_{n+1} = \mathbf{x}_n + h \,\mathbf{f}(\mathbf{x}_n). \end{equation*}

Instead of the initial condition, we consider the polynomial

\begin{equation*} P_{0,\mathbf{x}_0}\,(\delta\mathbf{x}) = [x_0+\delta x, y_0 + \delta y]^T. \end{equation*}


\begin{eqnarray*} \mathbf{x}_1 &=& P_{h, \mathbf{x}_0}(\delta\mathbf{x}) = P_{0,\mathbf{x}_0}(\delta\mathbf{x}) + h\, \mathbf{f}(P_{0,\mathbf{x}_0}(\delta\mathbf{x}))\\ & = & \left( \begin{array}{c} x_0 + h y_0 \\ y_0 - h x_0 \end{array} \right) + \left( \begin{array}{cc} 1 & h \\ -h & 1 \end{array} \right) \left( \begin{array}{c} \delta x\\ \delta y \end{array} \right). \end{eqnarray*}

At each step of Taylor's method, we can write

\begin{equation*} \phi(t_n+h; \mathbf{x}_n) = \sum_{i=0}^p \mathbf{x}^{(i)}(\mathbf{x}_n, t_n) \,h^i. \end{equation*}

The Lorenz system

\begin{eqnarray*} \dot{x_1} & = & σ (x_2-x_1), \\ \dot{x_2} & = & x_1 (ρ-x_3)-x_2, \\ \dot{x_3} & = & x_1 x_2-β x_3. \end{eqnarray*}
In [1]:

# Some experimental things used are here
# Pkg.checkout("TaylorIntegration", "jp/poincare")
using TaylorIntegration, BenchmarkTools

using Plots
INFO: Recompiling stale cache file /Users/Jorge/.julia/lib/v0.6/TaylorSeries.ji for module TaylorSeries.
In [2]:
#Lorenz system parameters and ODE:
const σ, β, ρ = 16.0, 4.0, 45.92

function lorenz!(t, x, dx)
    dx[1] = σ*(x[2]-x[1])
    dx[2] = x[1]*(ρ-x[3])-x[2]
    dx[3] = x[1]*x[2]-β*x[3]
lorenz! (generic function with 1 method)
In [3]:
const x0 = [19.0, 20.0, 50.0] #the initial condition

const t0 = 0.0   #the initial time

const tmax = 100.0 #final time of integration
In [4]:
@time tv, xv = taylorinteg(lorenz!, x0, t0, tmax, 28, 1e-20; 
  0.679648 seconds (3.37 M allocations: 341.860 MiB, 6.05% gc time)
In [5]:
plot(xv[:,1], xv[:,2], xv[:,3])
title!("Lorenz attractor")