Taylor’s method#

Taylor’s method for solving systems of ordinary differential equations (ODEs) is conceptually simple. Given an ODE system in the general explicit from

(1)#\[\frac{d \boldsymbol{x}\left( t \right)}{dt}=\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right),\]

with initial conditions at \(t=t_0\)

(2)#\[\boldsymbol{x}_0=\boldsymbol{x}\left(t_0\right),\]

Taylor’s method approximates the value of the solution at \(t=t_1\) as the truncated Taylor series expansion of \(\boldsymbol{x}\left( t \right)\) around \(t=t_0\):

(3)#\[\boldsymbol{x}\left( t_1 \right) = \boldsymbol{x}_0 + \boldsymbol{x}'\left(t_0\right)h +\frac{1}{2}\boldsymbol{x}''\left(t_0\right)h^2+\ldots+\frac{\boldsymbol{x}^{\left( p \right)}\left(t_0\right)}{p!}h^p.\]

Here \(h=t_1-t_0\) is the integration timestep and \(p\) is the order of the Taylor method. Eq. (3) can be rewritten in more compact form as

(4)#\[\boldsymbol{x}\left( t_1 \right) = \sum_{n=0}^p \boldsymbol{x}^{\left[ n \right]} \left(t_0\right) h^n,\]

where we have defined

(5)#\[\boldsymbol{x}^{\left[ n \right]}\left( t \right) = \frac{1}{n!}\boldsymbol{x}^{\left( n \right)}\left( t \right)\]

as the normalised derivative of order \(n\) of \(\boldsymbol{x}\left( t \right)\).

The derivatives appearing in the Taylor polynomial (4) can be computed from the initial conditions (2) and the right-hand side of the ODE system (1). Thus, (4) effectively can be used as a time stepper whose precision and performance characteristics depend on the choice of the step size \(h\) and the Taylor order \(p\).

Taylor integrators need to be… tailored to the specific expression of \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\). That is, unlike other popular numerical integration methods (e.g., Runge-Kutta methods), Taylor’s method requires the user not only to provide \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\), but also to implement all the derivatives necessary to construct the Taylor polynomial (4). This task, if done by hand, can be extremely cumbersome, inefficient and error-prone, especially for large ODE systems and/or high-accuracy applications.

The main functionality provided by heyoka is the ability to automatically synthesise a complete Taylor integrator starting only from a symbolic representation of \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\) and a set of initial conditions. Specifically, given a differentiable symbolic expression for \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\), heyoka takes care of:

the computation of the high-order derivatives necessary to implement the time stepper (4) via a process of automatic differentiation,
the deduction of optimal values for the Taylor order \(p\) and the (adaptive) step size \(h\),
the propagation of the state of the system via the evaluation of the Taylor polynomial (4).

In order to represent symbolically \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\), heyoka relies on a small, self-contained symbolic expression system (similar to an extremely trimmed-down, bare-bones computer algebra system). The expression system is used to decompose \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\) into a sequence of elementary subexpression on which automatic differentiation rules are applied. The sequence of operations necessary to compute the high-order derivatives of \(\boldsymbol{F}\left(t, \boldsymbol{x}\left( t \right) \right)\) is then assembled and compiled just-in-time (via LLVM) to produce a time stepper function usable from regular C++ code.