# Taylor-Fourier integration

M.P. Calvo111IMUVA and Departamento de Matemática Aplicada, Facultad de
Ciencias, Universidad de Valladolid, Spain. E-mail<EMAIL_ADDRESS>, J.
Makazaga222Universidad del Pais Vasco (UPV/EHU), Donostia-San Sebastián,
Spain. E-mail<EMAIL_ADDRESS>and A. Murua333Universidad del Pais
Vasco (UPV/EHU), Donostia-San Sebastián, Spain. E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper we introduce an algorithm which provides approximate solutions
to semi-linear ordinary differential equations with highly oscillatory
solutions that after an appropriate change of variables can be written as a
non-autonomous system with $(2\pi/\omega)$-periodic dependence on $t$. The
proposed approximate solutions are written in closed form as functions
$X(t,\omega\,t)$ where $X(t,\theta)$ is, (i) a truncated Fourier series in
$\theta$ for fixed $t$, and (ii) a truncated Taylor series in $t$ for fixed
$\theta$ (that is the reason for the name of the proposed integrators). Such
approximations are intended to be uniformly accurate in $\omega$ (in the sense
that their accuracy is not deteriorated as $\omega\to\infty$). This feature
implies that Taylor-Fourier approximations become more efficient than the
application of standard numerical integrators for sufficiently high basic
frequency $\omega$. The main goal of the paper is to propose a procedure to
efficiently compute such approximations by combining power series arithmetic
techniques and the FFT algorithm. We present numerical experiments that
demonstrate the effectiveness of our approximation method through its
application to well-known problems of interest.

Keywords: Taylor series, Fourier methods, oscillatory problems, numerical
integrators.

##  1 Introduction

In this paper we consider differential systems of the form

$\frac{d}{dt}{\bf x}=\omega A{\bf x}+{\bf g}({\bf x}),$ (1)

where ${\bf g}:\mathbb{R}^{D}\to\mathbb{R}^{D}$ is a smooth map, $\omega$ is a
positive real parameter, and $A$ is a real $D\times D$ matrix whose
eigenvalues are integer multiples of the imaginary unit $i$. Hence, the matrix
exponential $e^{\theta A}$ is $(2\pi)$-periodic in $\theta$ and, therefore,
the solution ${\bf x}(t)=e^{t\omega A}{\bf x}(0)$ of (1) with ${\bf
g}\equiv{\bf 0}$ is $(2\pi/\omega)$-periodic in $t$. If $\omega$ is large, the
solution of the homogeneous problem associated to (1) is highly oscillatory.
Examples of differential systems (1) with these properties can be found, for
instance, when dealing with semi-discretized dispersive partial differential
equations [3, 9, 17], or in some formulations of the equations of motion in
satellite dynamics [16].

We focus on approximating the solution ${\bf x}(t)$ of (1) supplemented with
the initial condition

${\bf x}(0)={\bf x}_{0}\in\mathbb{R}^{D}.$ (2)

With the change of variables

${\bf x}(t)=e^{t\omega A}{\bf y}(t),$ (3)

the initial value problem (1)–(2) reads

$\frac{d}{dt}{\bf y}={\bf f}(\omega t,{\bf y}),\quad{\bf y}(0)={\bf x}_{0},$
(4)

with

${\bf f}(\theta,{\bf y})=e^{-\theta A}{\bf g}(e^{\theta A}{\bf y}),$ (5)

which is $(2\pi)$-periodic in $\theta$.

Provided that the basic frequency $\omega$ is high enough, solving numerically
(4) by a Runge-Kutta method demonstrates superior efficiency compared to the
direct application of the same Runge-Kutta scheme to the original equation
(1).

In the context of problems of the form (1) obtained from the space semi-
discretization of dispersive PDEs, the numerical integration with a Runge-
Kutta method of the transformed system (4), as a means to solve (1)–(2), leads
to the so-called Lawson’s generalized Runge-Kutta methods [11, 9].

A completely different context where problems of the considered form arise is
the numerical orbit propagation of artificial satellites. The equations of
motion of the satellite can be brought into the form (1), for instance, by
using the Kustaanheimo-Stiefel (KS) formulation [10]. In this context, (4) is
referred to as the VOP (variation of parameters) formulation of the system
(1). In particular, a VOP formulation of the KS equations was presented and
analyzed by Stiefel and Scheifele in [16]. Other VOP formulations (obtained by
using different sets of coordinates) of the equations of motion in satellite
dynamics are also of the form (4).

From now on we will focus on computing approximations to the solution ${\bf
y}(t)$ of the initial value problem (4) which are valid on an appropriate
interval around $t=0$. In the highly oscillatory case, for high enough basic
frequency $\omega$, the larger is the basic frequency $\omega$, the more
computational effort will be required by the Runge-Kutta method to achieve a
prescribed accuracy. In the present work, we aim at obtaining a closed form
approximation to the solution ${\bf y}(t)$ of (4) that exhibits uniform
accuracy across varying frequencies $\omega$. This implies that comparable
accuracy is attained with the same computational effort even for arbitrarily
high frequencies. The Taylor-Fourier approximations that we propose here draw
inspiration from previous works on the numerical integration of oscillatory
problems (see e.g. [4], [7], and references therein), and are related to the
methods proposed in [5] and in [2].

For a prescribed positive integer $M$, we propose to compute a sequence

$\left\\{{\bf y}^{[1]}(t),{\bf y}^{[2]}(t),\ldots,{\bf
y}^{[d]}(t),\ldots\right\\}$

of approximations to the solution ${\bf y}(t)$ of (4) of the form

${\bf y}^{[d]}(t)=\sum_{k=-M}^{M}e^{ik\omega t}\left(\sum_{j=0}^{d}t^{j}\,{\bf
y}_{k,j}^{[d]}\right),\quad d=1,2,\ldots,$ (6)

with appropriate coefficients ${\bf y}_{k,j}^{[d]}\in\mathbb{C}^{D}$, and
starting from ${\bf y}^{[0]}(t)\equiv{\bf x}_{0}$ (this means that we can use
(6) with $d=0$, setting ${\bf y}_{0,0}^{[0]}={\bf x}_{0}$). We will refer to
(6) as the (M,d)-Taylor-Fourier approximation to the solution of (4).

In order to compute the successive approximations (6) to the solution of (4),
for $d=0,1,\ldots$, we do the following:

  * •

Given ${\bf y}^{[d]}(t)$, we first insert it into the right hand side of (4),
and then we determine an approximation ${\bf z}^{[d]}(t)$ to ${\bf f}(\omega
t,{\bf y}^{[d]}(t))$ of the form

${\bf z}^{[d]}(t)=\sum_{k=-M}^{M}e^{ik\omega t}\left(\sum_{j=0}^{d}t^{j}\,{\bf
z}_{k,j}^{[d]}\right).$ (7)

This approximation, to be defined precisely in Section 2, is determined by
combining truncated Taylor expansions and trigonometric interpolation of
periodic functions. That is the reason for the name of the proposed
integrators.

  * •

Once ${\bf z}^{[d]}(t)$ has been obtained, we use it as the right-hand side of
(4) and integrate to obtain

$\begin{split}{\bf y}^{[d+1]}(t)&={\bf y}_{0}+\int_{0}^{t}{\bf
z}^{[d]}(s)ds\\\ &=\sum_{k=-M}^{M}e^{ik\omega
t}\left(\sum_{j=0}^{d+1}t^{j}\,{\bf y}_{k,j}^{[d+1]}\right).\end{split}$ (8)

It is straightforward to check that, by construction, ${\bf
y}_{k,d+1}^{[d+1]}=0$ for $k\neq 0$.

The aim of the present paper is: (i) to introduce and give a precise
definition of Taylor-Fourier approximations, (ii) to propose a procedure to
compute them efficiently by making extensive use of the FFT algorithm, and
(iii) to illustrate the application of the proposed approximation method to
certain well known problems of interest. Our findings suggest that, provided
that the basic frequency $\omega$ is high enough, Taylor-Fourier
approximations can achieve comparable precision to standard numerical
integration methods applied to the transformed system (4), while demanding
less CPU time.

We want to stress that employing Taylor-Fourier approximations may be of
interest even in cases where the application of standard numerical integrators
to the transformed system (4) proves to be very efficient. Indeed, the
algorithm introduced in the present paper provides, for given $d$ and $M$,
closed expressions

${\bf X}(t)=e^{t\omega A}\left(\sum_{k=-M}^{M}e^{ik\omega
t}\sum_{j=0}^{d}t^{j}\,{\bf y}_{k,j}^{[d]}\right)$ (9)

of approximate solutions of (1)-(2), valid for a relatively large time
interval and that, once computed, can be efficiently evaluated at arbitrary
times $t$. This is in contrast with most time integrators commonly used, which
generate approximations to the exact solution of the problem either at a
prescribed final time $T$ or at a finite set $0<t_{1}<t_{2}<\cdots\leq T$ of
intermediate times. (In the later case, approximations to the solutions at
arbitrary times $t\in[0,T]$ are typically obtained with piecewise
interpolating polynomials.)

In our numerical experiments, we will compare the $(M,d)$-Taylor-Fourier
approximations (6) to the solution ${\bf y}(t)$ of (4) (equivalently, the
approximation (9) to the solution ${\bf x}(t)$ of (1)–(2)) with the
approximations obtained by numerically integrating (4) with highly efficient
Runge-Kutta methods (equivalently, numerically integrating (1)–(2) with
Lawson’s generalized Runge-Kutta methods). A state-of-the-art highly
performant library for the numerical integration of differential equations,
DifferentialEquation.jl [15], implemented in julia language [1], is freely
available and can be used to integrate the transformed system (4) with
optimized adaptive explicit Runge-Kutta integrators. In order to compare such
numerical solvers with the integrators proposed in this paper, we have carried
out a careful implementation of the Taylor-Fourier methods in julia language,
paying special attention to computational aspects such as memory management
and efficient power series arithmetic.

In Section 2 we give a precise definition of $(M,d)$-Taylor-Fourier
approximations, and demonstrate their applicability to problem classes
extending beyond real semi-linear problems of the form (1). Section 3 is
devoted to describe an algorithm designed to efficiently compute Taylor-
Fourier approximations. In Section 4 we present a couple of relevant numerical
examples where the performance of Taylor-Fourier integrators is analyzed.
Finally, in Section 5, we summarize our conclusions and outline avenues for
future research on Taylor-Fourier methods.

##  2 Taylor-Fourier methods

###  2.1 Definition of Taylor-Fourier approximations

So far, we have not specified how to determine the approximation (7) to ${\bf
f}(\omega t,{\bf y}^{[d]}(t))$ in the procedure introduced in Section 1. We
first note that (6) and (7) can be rewritten as

${\bf y}^{[d]}(t)=\sum_{j=0}^{d}t^{j}\,{\bf Y}_{j}^{[d]}(\omega t),\qquad{\bf
z}^{[d]}(t)=\sum_{j=0}^{d}t^{j}\,{\bf Z}_{j}^{[d]}(\omega t),$

respectively, where each component of the $D$-vector functions ${\bf
Y}_{j}^{[d]}(\theta)$ and ${\bf Z}_{j}^{[d]}(\theta)$ is a trigonometric
polynomial of degree $M$ in the angular variable $\theta$, that is,

${\bf Y}_{j}^{[d]}(\theta)=\sum_{k=-M}^{M}e^{ik\theta}{\bf
y}_{k,j}^{[d]},\qquad{\bf Z}_{j}^{[d]}(\theta)=\sum_{k=-M}^{M}e^{ik\theta}{\bf
z}_{k,j}^{[d]},$

with coefficients ${\bf y}_{k,j}^{[d]},{\bf z}_{k,j}^{[d]}\in\mathbb{C}^{D}$,
$0\leq j\leq d$, $-M\leq k\leq M$.

Given a positive integer $M$ and ${\bf Y}_{j}^{[d]}(\theta)$,
$j=0,1,\ldots,d$, we will determine ${\bf Z}_{j}^{[d]}(\theta)$,
$j=0,1,\ldots,d$ by imposing the following conditions:

  * •

For each $\theta_{n}=n\pi/M$, $n=0,1,2,\ldots,2M-1$,

${\bf f}\left(\theta_{n},\sum_{j=0}^{d}t^{j}\,{\bf
Y}_{j}^{[d]}(\theta_{n})\right)=\sum_{j=0}^{d}t^{j}\,{\bf
Z}_{j}^{[d]}(\theta_{n})+\mathcal{O}(t^{d+1}){\color[rgb]{1,0,0},}$ (10)

as $t\to 0$.

  * •

For $j=0,1,\ldots,d$,

${\bf z}_{-M,j}^{[d]}={\bf z}_{M,j}^{[d]}.$ (11)

The first condition is aimed at ensuring that

$\sum_{j=0}^{d}t^{j}\,{\bf Z}_{j}^{[d]}(\theta)\approx{\bf
f}\left(\theta,\sum_{j=0}^{d}t^{j}\,{\bf Y}_{j}^{[d]}(\theta)\right)$

for all $\theta\in[0,2\pi]$ and for sufficiently small $t\in\mathbb{R}$,
whereas the second one is included to ensure uniqueness.

To sum up, given a positive integer $M$, the successive approximations (6) to
the solution of (4) are uniquely determined by (10)–(11) and (8). We will
refer to the approximation (6) as the $(M,d)$-Taylor-Fourier approximation.

###  2.2 Application to more general equations

The definition of Taylor-Fourier approximations can be extended without
modifications to initial value problems of the form (4) in the complex domain.
In particular, Taylor-Fourier approximations can be successfully applied to
equations of the form (1) defined in the complex domain, provided that $A$ is
a complex $D\times D$ matrix whose eigenvalues are integer multiples of the
imaginary unit $i$, and the map ${\bf g}:\mathbb{C}^{D}\to\mathbb{C}^{D}$ is
such that ${\bf g}({\bf x})={\bf G}({\bf x},{\bf\bar{x}})$, where ${\bf
G}:\mathbb{C}^{2D}\to\mathbb{C}^{D}$ is a polynomial map.

Going back to equations in the real domain, Taylor-Fourier methods are not
limited to semi-linear equations of the form (1). They can be applied more
generally to problems of the form

$\frac{d}{dt}{\bf x}=\omega\,{\bf r}({\bf x})+{\bf g}({\bf x}),\quad{\bf
x}(0)={\bf x}_{0},$ (12)

where ${\bf r}:\mathbb{R}^{D}\to\mathbb{R}^{D}$ and ${\bf
g}:\mathbb{R}^{D}\to\mathbb{R}^{D}$ are smooth maps and $\omega$ is a positive
real parameter, provided that the following two assumptions hold:

  1. 1.

For each $t\in\mathbb{R}$ there exists a $t$-flow map
$\varphi_{t}:\mathbb{R}^{D}\to\mathbb{R}^{D}$ of the differential system
$\frac{d}{dt}{\bf x}={\bf r}({\bf x})$. That is, for each $t\in\mathbb{R}$ and
${\bf x}\in\mathbb{R}^{D}$,

$\frac{\partial}{\partial t}\varphi_{t}({\bf x})={\bf r}(\varphi_{t}({\bf
x})){\color[rgb]{1,0,0}.}$ (13)

In that case, problem (12) is transformed with the change of variables

${\bf x}(t)=\varphi_{t\omega}({\bf y}(t))$

into (4), where

${\bf f}(\theta,{\bf y})=\varphi^{\prime}_{\theta}({\bf y})^{-1}{\bf
g}(\varphi_{\theta}({\bf y})).$ (14)

Here, $\varphi^{\prime}_{\theta}({\bf y})$ denotes the Jacobian matrix of the
map $\varphi_{\theta}$ evaluated at the vector ${\bf y}$.

  2. 2.

The function ${\bf f}(\theta,{\bf y})$ defined in (14) is $(2\pi)$-periodic in
$\theta$.

Under these two assumptions, the solution ${\bf x}(t)$ of (12) can be
approximated by

${\bf x}^{[d]}(t)=\varphi_{\omega\,t}\left({\bf y}^{[d]}(t)\right),$ (15)

where

${\bf y}^{[d]}(t)=\sum_{k=-M}^{M}e^{ik\omega t}\sum_{j=0}^{d}t^{j}\,{\bf
y}_{k,j}^{[d]}$

is the $(M,d)$-Taylor-Fourier approximation to the solution of (4)–(14).

In particular, the above-mentioned two assumptions hold for systems of the
form

$\begin{split}\frac{d}{dt}\phi_{j}&=k_{j}\omega+g_{j}({\bf x}),\quad
j=1,\ldots,d,\\\ \frac{d}{dt}a_{j}&=g_{d+j}({\bf x}),\quad j=1,\ldots,D-d,\\\
\end{split}$ (16)

where ${\bf x}=(\phi_{1},\ldots,\phi_{d},a_{1},\ldots,a_{D-d})$,
$k_{j}\in\mathbb{Z}$, and for each $j$, $g_{j}$ is $2\pi$-periodic in each of
the angle variables $\phi_{j}$. In that case,

$\varphi_{t}(\phi_{1},\ldots,\phi_{d},a_{1},\ldots,a_{D-d})=\left(\begin{matrix}\phi_{1}+k_{1}\,t\\\
\vdots\\\ \phi_{d}+k_{d}\,t\\\ a_{1}\\\ \vdots\\\ a_{D-d}\\\
\end{matrix}\right),$

and the change of variables ${\bf x}(t)=\varphi_{t\omega}({\bf y}(t))$, with
${\bf
y}=(\tilde{\phi}_{1},\ldots,\tilde{\phi}_{d},\tilde{a}_{1},\ldots,\tilde{a}_{D-d})$,
transforms the system (16) into

$\displaystyle\frac{d}{dt}\tilde{\phi}_{j}$
$\displaystyle=g_{j}(\tilde{\phi}_{1}+k_{1}\,\omega\,t,\ldots,\tilde{\phi}_{d}+k_{d}\,\omega\,t,\tilde{a}_{1},\ldots,\tilde{a}_{D-d}),\quad
j=1,\ldots,d,$ $\displaystyle\frac{d}{dt}\tilde{a}_{j}$
$\displaystyle=g_{d+j}(\tilde{\phi}_{1}+k_{1}\,\omega\,t,\ldots,\tilde{\phi}_{d}+k_{d}\,\omega\,t,\tilde{a}_{1},\ldots,\tilde{a}_{D-d}),\quad
j=1,\ldots,D-d.$

##  3 Efficient implementation of Taylor-Fourier integrators

We will next show, for given positive integers $M$ and $d$, and given vector
coefficients ${\bf y}_{k,j}^{[d]}\in\mathbb{C}^{D}$, $0\leq j\leq d$, $-M\leq
k\leq M$, how to compute efficiently the new coefficients ${\bf
y}_{k,j}^{[d+1]}\in\mathbb{C}^{D}$, $0\leq j\leq d+1$, $-M\leq k\leq M$. This
computation is carried out in four sequential steps (for simplicity of
notation, we omit superscripts $[d]$ in ${\bf Y}_{j}^{[d]}$ and ${\bf
Z}_{j}^{[d]}$, when they are not necessary). We recall that
$\theta_{n}=n\pi/M,\quad n=0,1,\ldots,2M-1$.

###  3.1 Step 1: Evaluation of trigonometric polynomials

For $n=0,1,\ldots,2M-1$, we calculate the vector coefficients

${\bf Y}_{n,j}={\bf Y}_{j}(\theta_{n})=\sum_{k=-M}^{M}e^{ikn\pi/M}{\bf
y}_{k,j}^{[d]},\quad j=0,1,\ldots,d,$ (17)

by applying the inverse DFT

${\bf Y}_{n,j}=\frac{1}{2M}\sum_{k=0}^{2M-1}e^{ikn\pi/M}\hat{\bf
Y}_{k,j},\quad n=0,1,\ldots,2M-1,$

where for $j=0,1,\ldots,d$,

$\begin{array}[]{l}\hat{\bf Y}_{0,j}=2M\,{\bf y}_{0,j}^{[d]},\quad\hat{\bf
Y}_{M,J}=2M({\bf y}_{M,j}^{[d]}+{\bf y}_{-M,j}^{[d]}),\\\\[8.0pt] \hat{\bf
Y}_{k,j}=2M\,{\bf y}_{k,j}^{[d]},\quad\hat{\bf Y}_{2M-k,j}=2M{\bf
y}_{-k,j}^{[d]},\quad k=1,\ldots,M-1.\end{array}$

###  3.2 Step 2: Truncated Taylor expansions

For each $n=0,1,\ldots,2M-1$, we compute ${\bf Z}_{n,j}$, $j=0,1,\ldots,d$,
such that

${\bf f}\left({\theta_{n},}\sum_{j=0}^{d}t^{j}\,{\bf
Y}_{n,j}\right)=\sum_{j=0}^{d}t^{j}\,{\bf Z}_{n,j}+{\cal O}(t^{d+1}).$

This means that we obtain the Taylor polynomial of degree $d$ of the left hand
side of the previous equality. The computation of the required Taylor
expansion can be done in a very efficient way when appropriate tools to
perform power series arithmetic are available.

###  3.3 Step 3: Trigonometric interpolation

For each $j=0,1,\ldots,d$, we determine the coefficients ${\bf z}_{k,j}$,
$-M\leq k\leq M$, of the trigonometric polynomial

${\bf Z}_{j}(\theta)=\sum_{k=-M}^{M}e^{ik\theta}{\bf z}_{k,j},$

that satisfy the interpolatory conditions

${\bf Z}_{j}(\theta_{n})={\bf Z}_{n,j},\quad n=0,1,2,\ldots,2M-1,$

together with ${\bf z}_{-M,j}={\bf z}_{M,j}$.

This task can be efficiently done by computing the discrete Fourier Transform
(DFT)

$\hat{\bf Z}_{k,j}=\sum_{n=0}^{2M-1}e^{-ikn\pi/M}\,{\bf Z}_{n,j},\quad
k=0,1,\ldots,2M-1,$ (18)

and setting

$\begin{array}[]{l}{\bf z}_{0,j}=\dfrac{1}{2M}\hat{\bf Z}_{0,j},\quad{\bf
z}_{M,j}={\bf z}_{-M,j}=\dfrac{1}{4M}\hat{\bf Z}_{M,J},\\\\[8.0pt] {\bf
z}_{k,j}=\dfrac{1}{2M}\hat{\bf Z}_{k,j},\quad{\bf
z}_{-k,j}=\dfrac{1}{2M}\hat{\bf Z}_{2M-k,j},\quad k=1,\ldots,M-1.\end{array}$
(19)

Remark: For systems of the form (4) in the real domain, the symmetries present
in the complex vector resulting from the application of DFT to a real vector
can be used to compute it more efficiently (with the real FFT algorithm) than
in the general case of complex vectors. In particular, this implies that each
${\bf z}_{-k,j}$ coincides with the complex conjugate of ${\bf z}_{k,j}$,
which allows us to halve the number of coefficients that must be computed and
saved.

###  3.4 Step 4: Computing quadratures

Given

${\bf z}(t)=\sum_{k=-M}^{M}e^{ik\omega t}\left(\sum_{j=0}^{d}t^{j}\,{\bf
z}_{k,j}\right),$

we want to write

${\bf y}^{[d+1]}(t)={\bf y}_{0}+\int_{0}^{t}{\bf z}(s)ds{.}$ (20)

To simplify the notation, we will drop the superscript $[d+1]$ from the vector
coefficients ${\bf y}_{k,j}^{[d+1]}$ in the remaining of the present
subsection, so that we will write

${\bf y}^{[d+1]}(t)=t^{d+1}\,{\bf y}_{0,d+1}+\sum_{k=-M}^{M}e^{ik\omega
t}\left(\sum_{j=0}^{d}t^{j}\,{\bf y}_{k,j}\right).$

First, we notice that ${\bf z}(t)$ can be written as

${\bf z}(t)={\bf p}(t)+\sum_{j=0}^{d}t^{j}\,\tilde{{\bf Z}}_{j}(\omega t),$
(21)

where

${\bf p}(t)=\sum_{j=0}^{d}t^{j}\,{\bf z}_{0,j},\quad\tilde{{\bf
Z}}_{j}(\theta)=\sum_{k=1}^{M}(e^{ik\theta}{\bf z}_{k,j}+e^{-ik\theta}{\bf
z}_{-k,j}),\quad j=0,\ldots,d.$ (22)

Therefore,

${\bf y}(t)={\bf y}_{0}+{\bf q}(t)+\sum_{j=0}^{d}t^{j}\,\tilde{{\bf
Y}}_{j}(\omega t),$ (23)

with

$\displaystyle\tilde{{\bf Y}}_{j}(\theta)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{M}(e^{ik\theta}{\bf y}_{k,j}+e^{-ik\theta}{\bf
y}_{-k,j}),\quad j=0,\ldots,d,$ (24) $\displaystyle{\bf q}(t)$
$\displaystyle=$ $\displaystyle-\tilde{{\bf Y}}_{0}(0)+\int_{0}^{t}{\bf
p}(s)ds.$ (25)

The coefficients ${\bf y}_{k,j},{\bf y}_{-k,j}$, $j=0,1,\ldots,{d+1}$,
$k=1,\ldots,M$, can be efficiently obtained as follows:

  * •

Set $\tilde{{\bf Y}}_{d+1}(\theta)\equiv{\bf 0}$, and notice that taking
derivatives in (23) with respect to $\theta=\omega t$, and using (20) and (21)
one gets

$\omega\frac{d}{d\theta}\tilde{{\bf Y}}_{j}(\theta)=\tilde{{\bf
Z}}_{j}(\theta)-(j+1)\tilde{{\bf Y}}_{j+1}(\theta),\quad j=0,1,\ldots,d.$

This means that setting ${\bf y}_{k,d+1}={\bf 0}$ for $k\neq 0$, we can
compute recursively for $j=d,d-1,\ldots,1,0$,

${\bf y}_{k,j}=\frac{1}{ik\omega}\left({\bf z}_{k,j}-(j+1){\bf
y}_{k,j+1}\right).$

  * •

For $k=0$, one gets from (25) and (22)

${\bf y}_{0,j+1}=\frac{{\bf z}_{0,j}}{j},\quad j=0,\ldots,d,$

and,

${\bf y}_{0,0}={\bf y}_{0}-\sum_{k=1}^{M}({\bf y}_{k,0}+{\bf y}_{-k,0}).$

##  4 Numerical illustrations

In this section we include two numerical examples to illustrate the
performance of the Taylor-Fourier approximations proposed in this paper, when
compared with state-of-the-art methods. All codes we have used to carry out
the numerical experiments are available in
https://github.com/jmakazaga/Taylor-Fourier.

###  4.1 The cubic nonlinear Schrödinger equation

The first example is the cubic nonlinear Schrödinger (NLS) equation

$iu_{t}+u_{xx}+|u|^{2}u=0,\quad x\in[0,2\pi],$ (26)

with periodic boundary conditions and prescribed initial data

$u(x,0)=u^{0}(x).$ (27)

For the numerical solution of (26)-(27), we use the method of lines and then,
we first discretize in space with the spectral collocation method as follows.

We fix an integer $J$, set $h=\pi/J$, introduce the uniform grid of points
$x_{j}=(j-1)h,1\leq j\leq 2J$, and denote ${\bf
U}(t)=[U_{1}(t),\ldots,U_{2J}(t)]^{T}$, where $U_{j}(t)$ is an approximation
to $u(x_{j},t)$, for $1\leq j\leq 2J$. We also consider the vector $\hat{\bf
U}(t)=[\hat{U}_{1}(t),\ldots,\hat{U}_{2J}(t)]^{T}$ obtained by applying the
discrete Fourier transform to ${\bf U}(t)$, that is, for $k=1,\ldots,2J$,

$\hat{U}_{k}(t)=\frac{1}{2J}\sum_{j=1}^{2J}U_{j}(t)e^{-i(k-1)x_{j}}.$

In matrix notation, we write it as $\hat{\bf U}={\cal F}{\bf U}$. The vector
${\bf U}$ can be recovered from $\hat{\bf U}$ by applying the inverse Fourier
transform, i.e., ${\bf U}={\cal F}^{-1}\hat{\bf U}$. More explicitly, for
$j=1,\ldots,2J$,

$U_{j}(t)=\frac{1}{2J}\sum_{k=1}^{2J}\hat{U}_{k}(t)e^{i(k-1)x_{j}},$

or equivalently,

$U_{j}(t)=\sum_{k=-J}^{J}\hat{u}_{k}(t)e^{ikx_{j}},$

where

$\begin{array}[]{l}\hat{u}_{0}(t)=\dfrac{1}{2J}\hat{U}_{1}(t),\quad\hat{u}_{J}(t)=\hat{u}_{-J}(t)=\dfrac{1}{4J}\hat{U}_{J+1}(t),\\\\[8.0pt]
\hat{u}_{k}(t)=\dfrac{1}{2J}\hat{U}_{k+1}(t),\quad\hat{u}_{-k}(t)=\dfrac{1}{2J}\hat{U}_{2J-k+1}(t),\quad
k=1,\ldots,J-1.\end{array}$ (28)

Notice that

$\tilde{u}(x,t)=\sum_{k=-J}^{J}\hat{u}_{k}(t)e^{ikx}$

is a trigonometric polynomial satisfying the interpolation conditions

$\tilde{u}(x_{j},t)=U_{j}(t),\quad j=1,\ldots,2J.$

Both, the discrete Fourier transform and its inverse, can be efficiently
computed by applying the fast Fourier transform (FFT) algorithm. This is
particularly efficient when $J=2^{l}$.

The spectral collocation semi-discretization of the nonlinear Schrödinger
equation is obtained by requiring that

$i\tilde{u}_{t}(x_{j},t)+\tilde{u}_{xx}(x_{j},t)+|\tilde{u}(x_{j},t)|^{2}\tilde{u}(x_{j},t)=0,\quad\quad
j=1,\ldots,2J.$

By taking into account that

$\displaystyle i\,\tilde{u}_{xx}(x_{j},t)$ $\displaystyle=$
$\displaystyle\sum_{k=-J}^{J}-i\,k^{2}\,\hat{u}_{k}(t)e^{ikx_{j}}$
$\displaystyle=$
$\displaystyle\frac{1}{2J}\sum_{k=1}^{J}\lambda_{k}\,\hat{U}_{k}(t)e^{i(k-1)x_{j}},$

where

$\lambda_{k}=\left\\{\begin{array}[]{rcl}-i(k-1)^{2}&\mbox{ for
}&k=1,\ldots,J+1,\\\ -i(2J-k+1)^{2}&\mbox{ for
}&k=J+2,\ldots,2J,\end{array}\right.$

we arrive at the semi-discretized NLS equation

$\frac{d}{dt}{\bf U}=A{\bf U}+{\bf g}({\bf U}),$ (29)

where

${\bf g}({\bf U})=i\,\left(\begin{matrix}|U_{1}|^{2}U_{1}\\\ \vdots\\\
|U_{2J}|^{2}U_{2J}\end{matrix}\right),$ (30)

$A={{\cal F}}^{-1}D{{\cal F}}$, and $D$ is the diagonal matrix with diagonal
entries

$(\lambda_{1},\ldots,\lambda_{2J})=(0,-i,-2^{2}i,\ldots,-(J-1)^{2}i,-J^{2}i,-(J-1)^{2}i,\ldots,-2^{2}i,-i).$

The time-dependent change of variables ${\bf U}=e^{tA}{\bf W}$ transforms (29)
into

$\frac{d}{dt}{\bf W}={\bf f}(t,{\bf W}),$ (31)

where ${\bf f}(t,{\bf W})$ is $(2\pi)$-periodic in $t$ and is defined by

${\bf f}(\theta,{\bf W})=e^{-\theta A}{\bf g}\left(e^{\theta A}{\bf
W}\right).$

Products of a vector in $\mathbb{C}^{2J}$ by matrices $e^{\theta A}$ and
$e^{-\theta A}$ can be computed efficiently by performing two FFT in
$\mathbb{C}^{2J}$, since $e^{\theta A}={\cal F}^{-1}e^{\theta D}{\cal F}$, and
$e^{\theta D}$ is the diagonal matrix with diagonal entries
$(e^{\lambda_{1}\theta},\ldots,e^{\lambda_{2J}\theta})$. On the other hand,
due to the particular form of the nonlinearity ${\bf g}$ in (30), the Taylor
expansions required to apply the Taylor-Fourier integrator are cheaply
available since sums and products of power series in $t$ can be computed in a
very efficient way.

In the numerical experiments we have considered the NLS equation with initial
data $u(x,0)=\eta(x)$ equal to the step function

$\eta(x)=\left\\{\begin{array}[]{l}-1\quad\mbox{ if }x\in[0,\pi),\\\
\hphantom{-}1\quad\mbox{ if }x\in[\pi,2\pi],\end{array}\right.$ (32)

a rough initial data that has widely appeared in the recent literature [13,
14, 6]. For instance, in [6] the initial condition (32) is used to illustrate
the Talbot effect.

Figure 1: Numerical solution of (26)-(32) at $t=0.3$ (top left), $t=0.31$ (top
right), $t=0.314$ (bottom left) and $t=\pi/10$ (bottom right) computed with
$J=2^{9}$, and the $(2^{11},3)$-Taylor-Fourier approximation. Real part is
plotted with blue solid line and imaginary part with brown dashed line.

In the first experiment, we take $J=2^{7}$, $2^{8}$, $2^{9}$, $2^{10}$,
$T=\pi/10$, $d\leq 5$, and for each $J$, $M=J$, $2J$, $4J$. Once computed the
coefficients of the $(M,d)$-Taylor-Fourier approximation, we have evaluated
${\bf W}^{[d]}(t)$ at $t=0.3$, $0.31$, $0.314$ and $\pi/10$, as in [6]. The
results for $J=2^{9}$, $M=4J$ and $d=3$ are shown in Figure 1, to be compared
with those in Figures 5 and 6 in [6]. Different from [6], we have included in
the same plot the real and the imaginary part of the numerical solutions.

Secondly, to illustrate how the choice of $d$ and $M$ affects the numerical
solution computed with the Taylor-Fourier method, we include in Figure 2 the
numerical approximations to the solution of (26)-(32) at $t=\pi/10$ when
$J=2^{9}$, and the combinations $M=4J$, $d=1$ (left), $M=J$, $d=3$ (center)
and $M=J$, $d=1$ (right) are used. Plots in Figure 2 must be compared with the
bottom right plot in Figure 1 where the $(4J,3)$-Taylor-Fourier approximation
is displayed. Comparing the plots corresponding to the same value of $M$ and
reducing $d$ from 3 to 1, we observe that, although the overall shape of the
numerical solutions is rather similar, the maximum and minimum values of both
solutions are clearly different. Looking now at the solutions obtained with
the same value of $d$ we see that when decreasing $M$ from $4J$ to $J$ small
oscillations appear.

Figure 2:  $(M,d)$-Taylor-Fourier approximations to the solution of (26)-(32)
at $t=\pi/10$ for $J=2^{9}$, computed with $M=4J$, $d=1$ (left), $M=J$, $d=3$
(center) and $M=J$, $d=1$ (right). Real part is plotted with blue solid line
and imaginary part with brown dashed line.

In order to test the effect of changing the basic frequency $\omega$, we
consider (26) with initial data

$u(0,x)=\epsilon\,\eta(x)$ (33)

and $t\in[0,\epsilon^{-2}\pi/10]$, with different values of $\epsilon>0$. This
is equivalent to considering $u(t,x)=\epsilon\,v(\epsilon^{2}\,t,x)$, where
$v(\tau,x)$ is the solution of

$i\,v_{\tau}=-\omega\,v_{xx}-|v|^{2}v,\quad
x\in[0,2\pi],\quad\tau\in[0,\pi/10],$

with $\omega=\epsilon^{-2}$.

In the abstract and in Section 1 we claim that the accuracy of Taylor-Fourier
approximations does not deteriorate as the basic frequency $\omega$ increases.
We next check whether this holds true in the particular case of the ODE system
obtained after the space discretization with $J=2^{6}$ of (26) for the
$(2^{11},7)$-Taylor-Fourier approximation, with initial data (33) for
different values of $\epsilon>0$. In Figure 3, we represent (in logarithmic
scale) the error $|\tilde{U}_{j}(t)-U_{j}(t)|$ ($j=1,\ldots,2J$) at time
$t=\epsilon^{-2}\pi/10$ versus $x_{j}$, $j=1,\ldots,2J$ for $\epsilon=2^{-m}$,
$m=1,2,3,4$. Here, $\tilde{U}_{j}(t)$ denotes the reference solution of the
semi-discretized system computed using the Lawson integrator based on the
implementation DP5 in the julia library DifferentialEquation.jl of the 5th
order embedded explicit Runge-Kutta pair proposed in [8], with tolerance for
absolute and relative errors equal to $10^{-12}$. We observe that the errors
do not increase as $\omega=\epsilon^{-2}$ increases, but rather decrease
linearly with $\epsilon=\omega^{-1/2}$ in this case.

Figure 3: Errors of the $(2^{11},7)$-Taylor-Fourier approximation at
$t=\epsilon^{-2}\pi/10$ for $\epsilon=2^{-m}$, $m=1,2,3,4$ with $J=2^{6}$ and
initial data (33).

Finally, we have tested the efficiency (time-discretization error
$|\tilde{U}_{j}(t)-U_{j}(t)|$ versus CPU time) of Taylor-Fourier
approximations to the solution $u(t,x)$ of (26) with initial data (33). We
have considered different space discretization parameters depending on the
value of $\epsilon$: $J=2^{9}$ for $\epsilon=1$, and $J=2^{6}$ for
$\epsilon=2^{-4}$. We have plotted in Figure 4 the results corresponding to
$(M,d)$-Taylor-Fourier integrations with different choices of the parameters
$M$ and $d$.

The reference solution of the space-discretized problem used to measure time-
discretization errors has been computed as indicated above. For comparison
purposes, we have also included errors obtained with the same scheme for
tolerances $10^{-3}$, $\ldots$, $10^{-6}$ when $\epsilon=1$, and for
tolerances $10^{-3}$, $\ldots$, $10^{-7}$ for $\epsilon=2^{-4}$. In both
plots, a black dashed line representing the spatial error for the chosen value
of $J$ (estimated by a precise time-integration of the semi-discretized
problem with $2J$ spatial nodes) has been included. Notice that although $J$
is larger for $\epsilon=1$ than for $\epsilon=2^{-4}$, the spatial error is
also larger for $\epsilon=1$ than for $\epsilon=2^{-4}$. This is due to the
factor $\epsilon$ that multiplies the rough initial data (32). In both cases
we have used $2\leq d\leq 5$.

Figure 4: Errors against CPU time when $J=2^{9}$, the initial data is (32) and
$T=\pi/10$ (left), and when $J=2^{6}$, the initial data is (33) with
$\epsilon=2^{-4}$ and $T=\epsilon^{-2}\pi/10$ (right). Diamonds correspond to
the Lawson method based on DP5 and circles to the $(M,d)$-Taylor-Fourier
approximations

.

We observe in the left plot of Figure 4 that the Taylor-Fourier integrations
are considerably less efficient than the Lawson-DP5 integrations for the case
$\epsilon=1$. In the right plot, for $\epsilon=2^{-4}$, we see that the
performance of both integrators is rather similar. The low efficiency of the
Taylor-Fourier methods for $\epsilon=1$ is due to the high dimension of the
semi-discretized problem ($J=2^{9}$), and to the fact that in order to get
temporal errors not larger than the spatial errors, the parameter $M$ varies
from $2^{8}$ to $2^{11}$, increasing the computational cost of the method. For
$\epsilon=2^{-4}$, the dimension of the semi-discretized problem has been
divided by $8$, and the problem is more difficult to be solved using standard
integrators because its solution is really highly oscillatory. Then, although
$M$ ranges between $2^{6}$ and $2^{9}$, the advantage of Lawson-DP5 with
respect to Taylor-Fourier is much smaller.

###  4.2 A perturbed Kepler problem

The second example is a perturbed Kepler problem (the main problem of
satellite theory), which accounts for the oblateness of the earth. More
precisely, it is the Kepler problem perturbed by the $J_{2}$-term,

$\displaystyle\ddot{\bf q}$ $\displaystyle=-\frac{\mu}{r^{3}}{\bf q}-\nabla
V({\bf q}),$ (34) $\displaystyle{\bf q}(t_{0})$ $\displaystyle={\bf
q}_{0},\quad\dot{\bf q}(t_{0})=\dot{\bf q}_{0},$ (35)

where

${\bf q}=\left(\begin{matrix}x\\\ y\\\ z\end{matrix}\right),\quad r=\|{\bf
q}\|=\sqrt{x^{2}+y^{2}+z^{2}},\quad\dot{\bf q}=\frac{d}{dt}{\bf
q},\quad\ddot{\bf q}=\frac{d^{2}}{dt^{2}}{\bf q},$

and

$V({\bf
q})=J_{2}\,\frac{\mu\,R_{e}^{2}}{2\,r^{3}}\,\left(3\,\sin^{2}\theta-1\right),\quad\sin\theta=\frac{z}{r}.$
(36)

Here, $J_{2}$ is the zonal harmonic coefficient ($J_{2}=1.08262668\times
10^{-3}$), $\mu$ is the standard gravitational parameter of the Earth
($\mu=398600.44189\,km^{3}/seg^{2}$) and $R_{e}$ is the equatorial radius of
the Earth ($R_{e}=6378.137\,km$).

From now on, we will use the notation $\epsilon:=J_{2}\,\mu\,R_{e}^{2}$. From
(36), we get that

$\nabla
V(x,y,z)=\frac{3\,\epsilon}{2\,r^{5}}\left(\begin{matrix}(1-5\,\sin^{2}\theta)\,x\\\
(1-5\,\sin^{2}\theta)\,y\\\ (3-5\,\sin^{2}\theta)\,z\end{matrix}\right).$

####  4.2.1 KS-formulation of the equations of motion

We will rewrite the equations by making use of the Kustaanheimo-Stiefel (KS)
transformation.

The position variables ${\bf q}=(x,y,z)\in\mathbb{R}^{3}$ are replaced by new
variables ${\bf u}=(u_{1},u_{2},u_{3},u_{4})\in\mathbb{R}^{4}$ satisfying the
relation ${\bf q}=L({\bf u}){\bf u}$, where

$L({\bf u})=\left(\begin{matrix}u_{1}&-u_{2}&-u_{3}&u_{4}\\\
u_{2}&u_{1}&-u_{4}&-u_{3}\\\ u_{3}&u_{4}&u_{1}&u_{2}\\\ \end{matrix}\right).$

Clearly, it holds that $\|{\bf u}\|^{2}=\|{\bf q}\|=r$.

Instead of the physical time $t$, a new independent variable (the fictitious
time)

$\tau=\int_{0}^{t}\frac{1}{\|{\bf q}(t)\|}\,dt$

will be used (we will denote with a prime the derivative with respect to
$\tau$), and the physical time $t$ is considered as an additional dependent
variable of the equations of motion.

Following [16], the solution curve ${\bf q}(t)$ of the initial value problem
(34)–(35) can be parametrized as ${\bf q}=L({\bf u}(\tau)){\bf u}(\tau)$,
$t=t(\tau)$, where $({\bf u}(\tau),t(\tau))$ is the solution of the initial
value problem

$\displaystyle{\bf u}^{\prime\prime}$ $\displaystyle=-\frac{h}{2}{\bf
u}-\nabla R({\bf u}),$ (37) $\displaystyle{\bf u}(0)$ $\displaystyle={\bf
u}_{0},\quad{\bf u}^{\prime}(0)={\bf u}^{\prime}_{0},$ (38) $\displaystyle
t^{\prime}$ $\displaystyle=\|{\bf u}\|^{2},$ (39) $\displaystyle t(0)$
$\displaystyle=t_{0},$ (40)

where

$h=\frac{\mu}{\|{\bf q}_{0}\|}-\frac{1}{2}\,\langle\dot{\bf q}_{0},\dot{\bf
q}_{0}\rangle-V({\bf q}_{0}),$

and

$R({\bf u})=\frac{1}{2}\,\|{\bf u}\|^{2}\,V(L({\bf u}){\bf
u})=\frac{\epsilon}{4\,\|{\bf u}\|^{4}}\,\left(3\,\sin^{2}\theta-1\right),$

with $\sin\theta=2\,(u_{1}u_{3}+u_{2}u_{4})/\|{\bf u}\|^{2}$.

The initial condition ${\bf u}_{0}\in\mathbb{R}^{4}$ is chosen in such a way
that

${\bf q}_{0}=L({\bf u}_{0}){\bf u}_{0},$ (41)

and ${\bf u}_{0}^{\prime}\in\mathbb{R}^{4}$ is determined as

${\bf u}_{0}^{\prime}=\frac{1}{2}\,L({\bf u}_{0})^{T}\dot{\bf q}_{0}.$ (42)

There are infinitely many choices of ${\bf u}_{0}\in\mathbb{R}^{4}$ satisfying
(41) for a given ${\bf q}_{0}=(x_{0},y_{0},z_{0})$. Among them, we will
uniquely determine ${\bf u}_{0}=(u_{1,0},u_{2,0},u_{3,0},u_{4,0})$ by imposing
the following conditions:

$\begin{array}[]{ll}\mbox{if }x_{0}\geq 0,\quad u_{1,0}=u_{4,0}\geq 0,\\\
\mbox{if }x_{0}<0,\quad u_{2,0}=u_{3,0}\geq 0.\end{array}$ (43)

This leads to

$u_{1,0}=u_{4,0}=\frac{1}{2}\,\sqrt{r_{0}+x_{0}},\quad
u_{2,0}=\frac{y_{0}u_{1,0}+z_{0}u_{4,0}}{r_{0}+x_{0}},\quad
u_{3,0}=\frac{z_{0}u_{1,0}-y_{0}u_{4,0}}{r_{0}+x_{0}},$

if $x_{0}\geq 0$, and

$u_{2,0}=u_{3,0}=\frac{1}{2}\,\sqrt{r_{0}-x_{0}},\quad
u_{1,0}=\frac{y_{0}u_{2,0}+z_{0}u_{3,0}}{r_{0}-x_{0}},\quad
u_{4,0}=\frac{z_{0}u_{2,0}-y_{0}u_{3,0}}{r_{0}-x_{0}},$

if $x_{0}<0$, where $r_{0}=\|{\bf q}_{0}\|$.

The gradient of $R({\bf u})$ can be computed as

$\nabla R({\bf u})=\ \frac{\epsilon}{\|{\bf
u}\|^{6}}\,(1-6\sin^{2}\theta)\left(\begin{matrix}u_{1}\\\ u_{2}\\\ u_{3}\\\
u_{4}\end{matrix}\right)+\frac{3\,\epsilon}{\|{\bf
u}\|^{6}}\,\sin\theta\left(\begin{matrix}u_{3}\\\ u_{4}\\\ u_{1}\\\
u_{2}\end{matrix}\right).$

Clearly, the system of second order ODEs (37) is of the form (1) with

$\omega=\sqrt{h/2},\quad{\bf x}=\left(\begin{matrix}{\bf u}\\\ {\bf
u}^{\prime}\end{matrix}\right),\quad
A=\left(\begin{matrix}0&\omega^{-1}I_{4}\\\
-\omega\,I_{4}&0\end{matrix}\right),\quad g({\bf x})=\left(\begin{matrix}{\bf
0}\\\ -\nabla R({\bf u})\end{matrix}\right),$

where $I_{4}$ denotes the $4\times 4$ identity matrix.

####  4.2.2 Stiefel and Scheifele’s formulation

From now on, we assume that $h>0$, and denote $\omega=\sqrt{h/2}$. The
solution ${\bf u}(\tau)$ of the initial value problem (37)–(38) can be written
as [16]

${\bf
u}(\tau)=\cos(\omega\tau)\bm{\alpha}(\tau)+\omega^{-1}\sin(\omega\tau)\bm{\beta}(\tau),$

where $\bm{\alpha}(\tau)$ and $\bm{\beta}(\tau)$ are the solutions of the
8-dimensional ODE system

$\begin{split}\bm{\alpha}^{\prime}&=\omega^{-1}\sin(\omega\tau)\nabla
R(\cos(\omega\tau)\bm{\alpha}+\omega^{-1}\sin(\omega\tau)\bm{\beta}),\\\
\bm{\beta}^{\prime}&=-\cos(\omega\tau)\nabla
R(\cos(\omega\tau)\bm{\alpha}+\omega^{-1}\sin(\omega\tau)\bm{\beta}),\end{split}$
(44)

supplemented with the initial conditions

$\bm{\alpha}(0)={\bf u}_{0},\quad\bm{\beta}(0)={\bf u}^{\prime}_{0}.$ (45)

Actually, (44) is the VOP formulation (4) corresponding to (37).

The solution $t(\tau)$ of (39)–(40) can be written as

$t(\tau)=t_{0}+\int_{0}^{\tau}\|\cos(\omega\sigma)\bm{\alpha}(\sigma)+\omega^{-1}\sin(\omega\sigma)\bm{\beta}(\sigma)\|^{2}\,d\sigma,$

which can be explicitly computed once appropriate $(M,d)$-Taylor-Fourier
approximations

$\bm{\alpha}(\tau)\approx\sum_{k=-M}^{M}\sum_{j=0}^{d}e^{ik\omega\tau}\,\tau^{j}\,\bm{\alpha}_{k,j},\quad\bm{\beta}(\tau)\approx\sum_{k=-M}^{M}\sum_{j=0}^{d}e^{ik\omega\tau}\,\tau^{j}\,\bm{\beta}_{k,j},$

of $\bm{\alpha}(\tau)$ and $\bm{\beta}(\tau)$ are obtained.

####  4.2.3 High precision numerical propagation of a geostationary orbit

We now consider the initial state of a geostationary satellite [12, pg. 116].
The initial position coordinates (in kilometers) and velocity vectors (in
kilometers per second) are

$\displaystyle{\bf q}_{0}$ $\displaystyle=(42149.1336,0,0),$
$\displaystyle\dot{\bf q}_{0}$
$\displaystyle=(0,3.075823259987749,0.0010736649055318406).$

Recall from the discussion above that the initial conditions (45) for the
equations (44) can be computed from (41), (42) and (43).

We have computed several $(M,d)$-Taylor-Fourier approximations to the solution
of the VOP formulation (44), and found that $M=8$ and $d=8$ is a good choice
in that particular initial value problem. For comparison purposes, we have
solved (44) numerically (with absolute and relative tolerance set to
$10^{-13}$) for 400 revolutions with an explicit RK method of order 9 due to
Verner [18], efficiently implemented with adaptive time-stepping as algorithm
Vern9 in DifferentialEquations.jl. We have estimated the errors of both
numerical approximations by comparing them with a reference solution computed
(in high precision floating point arithmetic) with Vern9 with absolute and
relative tolerance set to $10^{-20}$.

Figure 5:  Relative errors in position (left) and absolute errors in physical
time (right) versus fictitious time $\tau$ (scaled by $P=2\pi/\omega$) of the
numerical integration with Vern9 with tol=$10^{-13}$ (top) and the
$(8,8)$-Taylor-Fourier approximation (bottom).

In Figure 5 we display the relative errors in position ${\bf q}$ (left
subplot) and the absolute errors in physical time $t$ (right subplot) versus
the fictitious time $\tau$ scaled by the period $P=2\pi/\omega$ of the
satellite orbits. The upper (resp. lower) plots correspond to errors of Vern9
(resp. Taylor-Fourier). Observe that the absolute errors in physical time are
smaller for the Taylor-Fourier approximation. Moreover, while similar relative
errors in position are obtained for both approximations at the end of the
integration interval, computing the Taylor-Fourier approximation required
considerably less CPU time than Vern9 (by a factor of 20 on our laptop).

We emphasize that while Vern9 exhibits the expected pattern of increasing
position errors due to the propagation of local truncation and round-off
errors, the Taylor-Fourier approximation maintains positional errors below a
threshold of approximately $3\times 10^{-15}$ (corresponding to the errors due
to truncation of Fourier modes) for approximately 380 revolutions. Using
$M=16$ (not shown here) instead of $M=8$, that threshold lowers to round-off
levels. Observe that after $380$ revolutions, these errors start to grow. We
have checked that error growth is approximately proportional to $t^{d+1}$. The
validity of the Taylor-Fourier approximation can be extended to longer time
intervals by increasing the degree $d$. For instance, with $d=9$, the errors
in position remain below the above-mentioned threshold for up to 500
revolutions.

####  4.2.4 High precision numerical propagation of a highly eccentric orbit

We now make similar computations with the initial state of a highly eccentric
orbit (with eccentricity $e=0.7679436$) obtained from [12, pg. 51]. In
particular,

$\displaystyle{\bf q}_{0}$
$\displaystyle=(11959.886901183693,-16289.448826603336,-5963.757695165331),$
$\displaystyle\dot{\bf q}_{0}$
$\displaystyle=(4.724300951633136,-1.1099935305609756,-0.3847854410416176).$

Figure 6:  Relative errors in position (left) and absolute errors in physical
time (right) versus fictitious time $\tau$ (scaled by $P=2\pi/\omega$) of the
numerical integration with Vern9 with tol=$10^{-13}$ (top) and the
$(128,14)$-Taylor-Fourier approximation (bottom).

We have computed several $(M,d)$-Taylor-Fourier approximations to the solution
of the corresponding initial value problem (44), and found that $M=128$ and
$d=14$ gives a precision similar to that obtained with $M=8$ and $d=8$ for the
geostationary orbit. This is due to the eccentricity of the current orbit
which is much higher than the eccentricity $e\approx 10^{-3}$ of the
geostationary orbit, and makes the initial value problem (44) more difficult
to be solved numerically (in the sense of requiring more computational effort
to achieve similar accuracy). An additional implication of the higher
eccentricity is that the interval of validity of the Taylor-Fourier
approximations is reduced.

Figure 6 is similar to Figure 5, the only differences being that we display
the errors for up to 40 revolutions (instead of 400) and that $(M,d)=(128,14)$
(instead of $(M,d)=(8,8)$). Vern9 is still applied with absolute and relative
tolerances set to $10^{-13}$.

As in the previous example, the absolute errors in physical time are smaller
for the Taylor-Fourier approximation. In this second example, the relative
errors in position are also smaller for the Taylor-Fourier approximation,
while the computation of the Taylor-Fourier approximation required
considerably less CPU time than Vern9 (by a factor of $2.6$ on our laptop).

In this example, the $(128,14)$-Taylor-Fourier approximation maintains
position errors below a threshold of approximately $8\times 10^{-16}$ up to
$35$ revolutions, and beyond that, exhibits an error growth that is closely
proportional to $t^{d+1}$. We have checked that this threshold is not lowered
by increasing $M$, unless a higher precision floating point arithmetic is
employed.

##  5 Conclusions

We have proposed approximate solutions to semi-linear ordinary differential
equations with highly oscillatory solutions that after an appropriate change
of variables can be written as a non-autonomous system with
$(2\pi/\omega)$-periodic dependence on $t$. These approximate solutions, that
can be written in closed form combining truncated Fourier and Taylor series,
are valid for a relatively large time interval and, once computed, can be
efficiently evaluated at arbitrary times in that interval. The resulting
approximate solutions are uniformly accurate in the basic frequency $\omega$,
which makes them more efficient than standard numerical integrators when
$\omega$ is large enough. In addition, we have provided a detailed procedure
to compute efficiently such approximate solutions making use of power series
arithmetic techniques and the FFT algorithm.

We have also reported numerical results to show the performance of the Taylor-
Fourier approximations when compared with state-of-the-art methods. When
applied to a perturbed Kepler problem, the proposed approximations become more
efficient than the numerical integration of the transformed system using the
adaptive algorithm Vern9, even for highly eccentric orbits. The application of
the Taylor-Fourier integration procedure to the semi-discretized cubic
nonlinear Schrödinger equation is not so advantageous from the computational
point of view when compared with the Lawson integrator based on the variable-
step implementation of DP5. This is mainly due to the high dimension of the
semi-discretized problem. However, the uniform accuracy of Taylor-Fourier
approximations with respect to the basic frequency $\omega$ implies that they
are more efficient than standard numerical methods (which typically experience
accuracy deterioration as the frequency $\omega$ increases) for high enough
values of $\omega$.

A detailed error analysis of Taylor-Fourier approximations, and their
connection with other related topics such as stroboscopic averaging and the
solution of the transport equation will be the subject of future work.

Acknowledgements. All the authors have received funding by the Spanish State
Research Agency through project PID2022-136585NB-C22,
MCIN/AEI/10.13039/501100011033, European Union. M.P.C. was also supported by
project VA169P20 (Junta de Castilla y León, ES) cofinanced by FEDER funds
(EU). J.M. and A.M. were partially supported by the Department of Education of
the Basque Government through the Consolidated Research Group MATHMODE
(ITI456-22).

## References

  * [1] J. Bezanson, A. Edelman, S. Karpimski & V.B. Shah, Julia: A Fresh Approach to Numerical computing, SIAM Review 59, 65–98 (2017).
  * [2] L. Brugnano, M. Calvo, J.I. Montijano & L. Rández, Fourier methods for oscillatory differential problems with a constant high frequency, AIP Conference Proc. 1863, 020003 (2017).
  * [3] B. Cano & A. González-Pachón, Exponential time integration of solitary waves of cubic Schrödinger equation, Appl. Numer. Math. 91, 26–45 (2015)
  * [4] P. Chartier, A. Murua & J.M. Sanz-Serna, A formal series approach to averaging: exponentially small error estimates, Discrete and Continuous Dynamical Systems 32, 3009–3027 (2012).
  * [5] P. Chartier, N. Crouseilles, M. Lemou & F. Méhats, Uniformly accurate numerical schemes for highly oscillatory Klein–-Gordon and nonlinear Schrödinger equations, Numer. Math. 129, 211–250 (2015).
  * [6] G. Chen & P.J. Olver, Numerical simulation of nonlinear dispersive quantization, Discrete and Continuous Dynamical Systems 34, 991–1008 (2014).
  * [7] D. Cohen, E. Hairer & Ch. Lubich, Modulated Fourier Expansions of Highly Oscillatory Differential Equations, Found. Comput. Math. 3, 327–345 (2003).
  * [8] J.R. Dormand & P.J. Prince, A family of embedded Runge-–Kutta formulae, J. Comput. Appl. Math. 6, 19–26 (1980).
  * [9] M. Hochbruck & A. Ostermann, Exponential Integrators, Acta Numerica 19, 209–286 (2010).
  * [10] P. Kustaanheimo & E. Stiefel, Perturbation Theory of Kepler Motion Based on Spinor Regularization, J. Reine Angew. Math. 218, 204-–219 (1965).
  * [11] J.D. Lawson, Generalized Runge–Kutta processes for stable systems with large Lipschitz constants, SIAM J. Numer. Anal. 4, 372–380 (1967).
  * [12] O. Montenbruck & E. Gill, Satellite Orbits. Models, Methods and Applications, Springer, Berlin (2000).
  * [13] A. Ostermann & K. Schratz, Low regularity exponential-type integrators for semilinear Schrödinger equations, Foundations of Computational Mathematics 18, 731–-755 (2018).
  * [14] A. Ostermann & F. Yao, A fully discrete low-regularity integrator for the nonlinear Schrödinger equation, J. Sci. Comput. 91, 9 (2022).
  * [15] C. Rackauckas & Q. Nie, Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia, Journal of Open Research Software 5(1):15 (2017).
  * [16] E. Stiefel & G. Scheifele, Linear and Regular Celestial Mechanics, Springer Verlag, Berlin (1971).
  * [17] Ll.N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia (2000).
  * [18] J. Verner, https://www.sfu.ca/ jverner/.