# Transforming Stiffness and Chaos

Jan Scheffel  Electromagnetic Engineering and Fusion Science
KTH Royal Institute of Technology, Stockholm, Sweden

###### Abstract

Stiff and chaotic differential equations are challenging for time-stepping
numerical methods. For explicit methods, the required time step resolution
significantly exceeds the resolution associated with the smoothness of the
exact solution for specified accuracy. In order to improve efficiency, the
question arises whether transformation to asymptotically stable solutions can
be performed, for which neighbouring solutions converge towards each other at
a controlled rate. Employing the concept of local Lyapunov exponents, it is
demonstrated that chaotic differential equations can be successfully
transformed to obtain high accuracy, whereas stiff equations cannot. For
instance, the accuracy of explicit fourth order Runge-Kutta solution of the
Lorenz chaotic equations can be increased by two orders of magnitude.
Alternatively, the time step can be significantly extended with retained
accuracy.

###### keywords:

Time-spectral; stiff equation; time-averaged; smoothing algorithm; GWRM.

††journal: journal for publication

##  1 Introduction

When assessing the usefulness of a computational method for differential
equations, there are at least five properties of the problems addressed that
come into play: 1) stability (characterized by Lyapunov exponents
$\lambda_{i}<0$, with small magnitudes), 2) oscillations (number of extrema in
the computational interval), 3) steepness (normalized maximum rate of change),
4) stiffness (at least one $\lambda_{i}<0$, with large magnitude), and 5)
chaoticity (continuous, non-periodic solution in at least three variables [1,
2], with at least one $\lambda_{i}>0$). We will, in this work, focus on the
concepts of stiffness and chaoticity for explicit time-stepping methods,
making frequent use of the concept of local Lyapunov exponents.

The phenomenon of stiffness may severly limit the time steps of explicit
methods. The concept has been found difficult to characterize precisely but
substantial progress was made in the work by Cartwright [3]. It is generally
acknowledged that for stability a stiff differential equation requires step
lengths that are ”excessively small in relation to the smoothness of the exact
solution” in the interval [4]. The reason for this behaviour is that the
scales for decay of different parts of the solution are widely separated. To
handle stiffness within time-stepping methods one must usually resort to time-
consuming, less accurate implicit or semi-implicit methods. An alternative is
time-spectral methods [5].

The notion of Lyapunov exponents is helpful when distinguishing between
stiffness and chaoticity for nonlinear problems. A novel method for
approximate computation of local Lyapunov exponents (LLEs), without the
necessity of solving the differential equations, is presented in Section 2. By
”local” is here meant that the exponents are associated with the current time
step or time interval. We also present an improved approach for categorization
and formal separation of the properties of stiffness and chaoticity, in line
with that of [3]. LLEs are applied to examples of stiff linear and nonlinear
equations in Sections 3 and 4, illustrating their usefulness.

In Sections 5 and 6, the concept of LLEs is instrumental for proving that
chaotic differential equations can be transformed into differential equations
with asymptotically stable solutions, and solved efficiently at strongly
improved accuracy. It is also shown that stiff differential equations cannot
be transformed similarly since the difficulty to numerically resolve
neighbouring, rapidly diminishing solutions, is transferred to the numerical
resolution of the transformed functions. Notably, time-spectral methods [5]
remain unaffected by stiffness or chaos, as they are methodologically acausal,
or in a sense, implicit. Consequently, they often represent optimal choices
for addressing stiff or chaotic problems, a topic that will be elaborated upon
in a forthcoming publication.

A discussion and summary of the transformation method proposed in this work
can be found in Sections 7 and 8.

##  2 Determining local stiffness and chaoticity

It is of interest to find a formal and practical characterization for
determining stiffness and chaoticity of a given linear or nonlinear system of
differential equations, both locally and globally. The concept of Lyapunov
exponents is useful for separating the properties of stability, stiffness and
chaoticity. We here employ a standard definition for global Lyapunov
exponents, whereafter we move on to present some novel aspects for LLEs.

Assume that a coupled set of $N$ linear or nonlinear differential equations
are solved for two different initial conditions with initial, infinitesimal
separation $\bf{\delta{X}_{0}}$ (components $\delta X_{0i}$) in phase space.
The time-dependent separation of solution curves can be approximated within a
linear analysis;

$|\delta X_{i}(t)|\approx e^{\lambda_{i}t}|\delta X_{0i}|.$ (1)

This relation defines the spectrum of $N$ Lyapunov exponents $\lambda_{i}$,
determined by separations in $N$-dimensional space. The largest of these
Lyapunov exponents is termed the maximal Lyapunov exponent $\lambda_{M}$. For
a chaotic problem, $\lambda_{M}$ is positive, thus initially closely separated
orbits will diverge exponentially in time. For many problems, however,
$\lambda_{M}$ is instead moderately negative so that the solution is stable in
the sense that perturbed solutions remain close to the original solution, or
even asymptotically stable in which case original and perturbed solutions
converge towards each other over time [6]. The latter property saves
convergence and accuracy even for basic, low order explicit difference methods
with large time steps.

For stiff problems, the minimum Lyapunov exponents are strongly negative. The
associated rapid convergence may be difficult to resolve for finite difference
methods. Whether a problem is stiff or not requires an analysis of the local
variation of the solution, since stiffness implies that the local step size
must be much smaller than what is necessitated by numerical resolution of the
solution. For scalar functions $u=u(t)$, the absolute local curvature
$\kappa(t)$ is a suitable measure of the local variation;

$\kappa(t)=\frac{|u^{\prime\prime}(t)|}{(1+{u^{\prime}(t)}^{2})^{3/2}}\quad.$
(2)

A prime denotes differentiation with respect to $t$. Local Lyapunov exponents
(LLEs) [3] can now be defined for each phase space dimension $i$;

$\gamma_{i}(\tau;t)=\lim_{\sigma_{i}(\tau)\rightarrow
0}\frac{1}{\tau}ln\frac{\sigma_{i}(t+\tau)}{\sigma_{i}(t)}\quad.$ (3)

Here $\sigma_{i}(\tau)$ denotes the $N$ principal axes of a small ellipsoid in
phase space, representing the temporally varying separations. By comparing
with the standard definition of global Lyapunov exponents (compare Eq. (1))

$\lambda_{i}=\lim_{t\rightarrow\infty}\,\lim_{|\delta X_{0i}|\rightarrow
0}\,\frac{1}{t}ln\frac{|\delta X_{i}(t)|}{|\delta X_{0i}|}\quad,$ (4)

it is clear that local Lyapunov exponents, for any time $t$, estimate the
deviation of the original and perturbed solutions in an interval $[t,t+\tau]$,
or equivalently $t\in[t_{n},t_{n+1}]$ with $t_{n+1}-t_{n}=\tau$. Lyapunov
exponents are obtained from local Lyapunov exponents from the relation

$\lambda_{i}=\lim_{\tau\rightarrow\infty}\gamma_{i}(\tau;t)\quad.$ (5)

Following Cartwright [3], stiffness can be defined for nonlinear systems,
still allowing for chaoticity: A system is stiff in a given interval if in
that interval the most negative local Lyapunov exponent is large, while the
curvature of the trajectory is small. This criterion can be employed for stiff
equations using the measure

$R=\frac{|\gamma_{i_{min}}(\tau;t)|}{\kappa(t)}$ (6)

where $\gamma_{i_{min}}$ denotes the most negative of the $n$ local Lyapunov
exponents. Since at least one positive Lyapunov exponent is required for chaos
and at least one should be strongly negative for stiffness, stiff chaos
implies a large spread of Lyapunov exponents. The measure is preferable to
more pragmatic earlier definitions like ”Stiff equations are problems for
which explicit methods don’t work” [7]. It also extends the linear concept of
stiffness ratio [4] to nonlinear problems. Furthermore, the definition avoids
the assumption that stiff problems should be non-chaotic, as in Lambert’s
definition [4].

The definition enables distinction between, for example, systems that are
stiff but non-chaotic ($R$ large and $\gamma_{i_{max}}\leq 0$) and systems
that are both stiff and chaotic ($R$ large and $\gamma_{i_{max}}>0$). If the
curvature associated with stiffness $\kappa_{stiff}(t)$ is comparable to, or
exceeds, the curvature $\kappa(t)$ of $u(t)$, accuracy is determined by
stiffness rather than by the functional variation of $u(t)$.

The maximum step length ${\Delta t}_{max}$ for approximating a function $u(t)$
with local curvature $\kappa(t)$ at accuracy $\epsilon$ in
$t\in[t_{n},t_{n+1}]$ can be estimated as follows. A small segment of a local
curvature radius $1/\kappa(t)$ can be approximated to accuracy $\epsilon$ by a
secant with length $\Delta t$ as a function of $t$. A simple geometric
analysis straightforwardly yields the approximate relation

$(\Delta t)_{max}=2\sqrt{2}\sqrt{\frac{\epsilon}{\kappa(t^{*})}}$ (7)

Note that $(\Delta t)_{max}$ is a function of time $t^{*}=t-t_{n}$.

Next we determine the maximum step size $(\Delta t)_{stiff}$ that is allowed
by stiffness in $t\in[t_{n},t_{n+1}]$, for specified accuracy $\epsilon$. The
focus is on contributions to the solution, or neighbouring solutions, in this
interval, with local dependence of the form $\delta u(t)=\epsilon e^{\gamma
t^{*}}$ where $\gamma$ is negative for stiff problems. Since $\epsilon$ is the
desired accuracy, which should be maintained throughout the solution
algorithm, the choice of $\delta u(0)$ is reasonable. In the general case,
$\gamma=\gamma(t)$, but we will here assume that $\gamma$ is a constant, that
may however differ for each interval $[t_{n},t_{n+1}]$. If the differential
equation is stiff, terms like $\epsilon e^{\gamma t^{*}}$ with associated
curvature

$\kappa_{stiff}(t^{*})=\frac{\epsilon\gamma^{2}e^{\gamma
t^{*}}}{(1+\epsilon^{2}\gamma^{2}e^{2\gamma t^{*}})^{3/2}}$ (8)

must be resolved; it does not suffice to resolve $u(t)$ itself to maintain
$\epsilon$ accuracy overall. From Eq. (8), it is seen that the maximum value
of $\kappa_{stiff}$ is independent of $\epsilon$;
$\kappa_{stiff,max}=(2\sqrt{3}/9)|\gamma|$ and occurs for
$t^{*}_{max}=-ln(2\gamma^{2}\epsilon^{2})/(2\gamma)$.

At this time (for $t^{*}_{max}>0$) the maximum allowable time step for
accuracy $\epsilon$ can be obtained using Eq. (7);

$(\Delta t)_{stiff}=6\sqrt{\frac{\epsilon}{\sqrt{3}|\gamma|}}$ (9)

In most cases $2\gamma^{2}\epsilon^{2}<1$ however, and $t^{*}_{max}<0$, so
that the maximum value of $\kappa$ is obtained at the start $t^{*}=0$ of the
interval $[t_{n},t_{n+1}]$;

$\kappa_{stiff,max}=\frac{\epsilon\gamma^{2}}{(1+\epsilon^{2}\gamma^{2})^{3/2}}.$
(10)

For this case the maximum allowed step length is found at $t^{*}=0$ and is

$(\Delta
t)_{stiff}=\frac{2\sqrt{2}}{|\gamma|}(1+\epsilon^{2}\gamma^{2})^{3/4}.$ (11)

Thus we have estimated the maximum allowed step length $(\Delta t)_{stiff}$,
determined by the stiffness of the problem, to be compared with $(\Delta
t)_{max}$, required by the specific numerical method for resolving the
solution curve to desired accuracy. For $(\Delta t)_{max}\leq(\Delta
t)_{stiff}$, the differential equation should not be regarded as stiff. In
adaptive methods, $(\Delta t)_{max}$ and $(\Delta t)_{stiff}$ can be
continually evaluated for determining local stiffness or chaoticity. When, for
example, applying the transformation techniques described in Section 5, a more
intuitive measure of stiffness than the ratio $R$ of Eq. (6) is useful:

$Q(t)\equiv\frac{(\Delta t)_{max}}{(\Delta t)_{stiff}}>1,$ (12)

A solution algorithm for which $Q(t)<1$ is thus not limited by stiffness. Note
that ${(\Delta t)}_{max}$ is given by Eq. (7) and ${(\Delta t)}_{stiff}$ is
found from either of Eqs. (9), for $2\gamma^{2}\epsilon^{2}\geq 1$, or Eq.
(11), for $2\gamma^{2}\epsilon^{2}<1$. This formalism will now be applied to
stiff and chaotic problems.

##  3 A stiff linear problem

The success of numerical integration techniques relies to a large extent on
that neighbouring solutions are asymptotically stable, converging towards each
other over time and damping out numerical errors. Stiffness appears, however,
when this convergence becomes extremely rapid, producing strong local
gradients that upset explicit schemes unless very small time steps are used.
We now solve a typical, stiff linear problem [6] and use the formalism of
Section 2 to diagnose its stiffness properties;

$\frac{du}{dt}=-au+at+a+1,\quad u(0)=1.$ (13)

The general solution is $u(t)=1+t+ce^{-at}$, with $c$ being an arbitrary
constant. The solution satisfying the given initial condition is $u(t)=1+t$.
For $a\gg 1$ and large $t$ the general solution becomes a straight line,
independent of the initial condition. For $u(0)=1$ the solution is smooth
(since $c=0$), but for neighbouring initial conditions (other values of $c$)
the solution changes very rapidly for small $t$, being problematic for
explicit methods. The cause of this behaviour is the strongly negative
Lyapunov exponent, due to the large value of $a$, in combination with the
extremely small trajectory curvature ([3, 6]). For this problem, we have
access to the exact solution, which shows that even small numerical errors
will make the equation stiff.

###  3.1 Approximate method for estimating local stiffness

We now introduce an approximative formal method to determine the character,
being it stable, stiff or chaotic, of linear and non-linear initial-value
problems, using LLEs. Other methods for determining LLEs are proposed in the
literature [8], but often feature higher computational cost due to the need
for computing long system trajectories.

Denoting the exact solution with $u(t;0)$ and an $\epsilon$-perturbed
($\epsilon\ll 1$) solution $u(t;\epsilon)$, obtained by solving Eq. (13) for
$u(0)=1+\epsilon$ instead, the deviation of the solutions becomes

$\delta u(t)=u(t;\epsilon)-u(t;0).$ (14)

Subtracting the two versions of Eq. (13) for $u(t;\epsilon)$ and $u(t;0)$,
there results

$\frac{d\delta u}{dt}=-a\delta u$ (15)

The solution to Eq. (15), for each time interval, is

$\delta u(t)=\delta u(t_{n})e^{-at^{*}}$ (16)

where $t\in[t_{n},t_{n+1}]$, with $t_{n}$ arbitrary and $t^{*}=t-t_{n}$. Hence
the solution to Eq. (13) is asymptotically stable for moderately large $a>0$
and stiff for $a\gg 1$. To study stiffness, we perturb the initial condition
$u(0;\epsilon)=1+\epsilon=1.05$ and set $a=300$. The perturbed solution
$u(t;\epsilon)=1+t+0.05e^{-300t}$ deviates marginally, but rapidly from
$u(t;0)=1+t$ at low $t$.

Figure 1: Comparison, for solution of Eq. (13) at $\epsilon=0.001$ accuracy,
between maximum time step $(\Delta t)_{max}$ for numerical resolution and
maximum time step $(\Delta t)_{stiff}$ for resolving stiffness of perturbed
solution $u(t;\epsilon)=1+t+0.05e^{-300t}$. It is seen that the ODE is locally
stiff for $t<0.004$, where $Q>1$.

The ODE is not stiff everywhere; see Figure 1. For $\delta
u(t_{0})=\epsilon=0.001$ there results $2\epsilon^{2}a^{2}=0.18<1$, indicating
that the maximal stiffness obtains at the beginning of the interval, and
$\kappa_{stiff,max}\approx 90$. This renders $(\Delta t)_{stiff}=0.0094$.
Using Eqs. (7) and (11), $(\Delta t)_{max}$ is computed as function of $t$.
The problem is not stiff when the curvature of $u(t;\epsilon)$ is sufficiently
large, so that $Q(t)=(\Delta t)_{max}/(\Delta t)_{stiff}<1$, which occurs for
$t>0.004$, as can be seen in Figure 1. Apparently the condition $Q(t)>1$ for
stiffness is an intuitively transparent and straightforward diagnostic.

Somewhat surprisingly, fourth order Runge-Kutta (RK4) solves this stiff
problem well. Stability can be reached already using as few as 12 time steps,
although at a relative error of $0.03$. An accuracy of $0.001$ is obtained
using 25 time steps. Solution of Eq. (13) using an implicit time-stepping
method (second order Trapezoidal method) [9], requires some 40 time steps for
this accuracy.

In summary we have, using a simple model problem, seen how LLEs can be used to
define and compute the time-dependent diagnostics $\Delta t_{max}$, $\Delta
t_{stiff}$ and $Q(t)$ for assessing local stiffness of an ODE.

##  4 Stiff nonlinear problems

We next consider nonlinear stiff problems. Two example problems will be
studied; a differential equation for modelling explosive combustion and a
problem in kinetic analysis of chemical reactions.

###  4.1 Explosive combustion

The dynamic evolution of a flame can, in spherically symmetric geometry, be
modelled [10] with the equation

$\frac{du}{dt}=u^{2}-u^{3},\quad u(0)=d.$ (17)

Here $d\in]0,1[$ is an arbitrary constant. This ODE is stiff for $t>1/d$ when
$d$ is small. Solutions for $d=0.1$, $d=0.05$ and $d=0.01$ are shown in Fig.
3.

Figure 2: GWRM solutions $u(t)$ of Eq. (17) for $d=0.1$, $0.05$ and $0.01$.

How can the claim that the equation is stiff be justified? We employ again the
method described in the previous section. Assuming the existence of an
$\epsilon$-perturbed solution to Eq. (17), we label the two solutions $u(t;0)$
and $u(t;\epsilon)$, letting $\delta u(t)=u(t;\epsilon)-u(t;0)$, and obtain

$\frac{d\delta u(t)}{dt}=\frac{du(t;\epsilon)}{dt}-\frac{du(t;0)}{dt}=-\delta
uf(u(t;0))$ (18)

to first order in $\delta u$, where $f$ is a function to be found. By
subtraction of the two versions of Eq. (17) we indeed have

$\frac{d\delta u}{dt}=-\delta u(3u^{2}(t;0)-2u(t;0)).$ (19)

As $u(t;0)\rightarrow 1$, the solution to Eq. (19) becomes $\delta u=\delta
u_{0}e^{-(t-t_{0})}=\delta u_{0}e^{-\tau}$, which explains stiffness when
$t\rightarrow 2/d>>1$ for small $d$ in this interval. If a constant step
length $\Delta t$ is employed throughout, it will for $d\leq 0.05$ be decided
by $(\Delta t)_{max}$ (obtained from Eqs. (7) and (17)) near $t=1/d$ and by
stiffness through $(\Delta t)_{stiff}$ (from Eqs. (11) and (19)) for $t>1/d$.

Figure 3: Illustration of stiffness of Eq. (17) for $d=0.01$. Eq. (19) is
solved in the two intervals $[105,110]$ and $[115,120]$. The time step for
$\epsilon=0.001$ accuracy becomes limited by stiffness sufficiently beyond
$t=1/d=100$ where the solution curve $u(t)$ approaches a straight line; thus
$(\Delta t)_{stiff}<(\Delta t)_{max}$, or $Q>1$ here.

The stiffness of Eq. (17), in the domain beyond $t=1/d$, is demonstrated for
$d=0.01$ in Figure 3. We employ $u(t;0)=1$ and $\delta u_{0}=\epsilon$. The
time step for $\epsilon=0.001$ accuracy becomes limited by stiffness ($(\Delta
t)_{stiff}<(\Delta t)_{max}$) beyond $t=1/d$ where the solution curve $u(t;0)$
approaches a straight line.

Explicit solvers are unsuitable when $d<<1$, that is for strong stiffness.
Time-spectral ODE solvers, like the GWRM [10], are however very efficient; for
$d=10^{-4}$, 69 time intervals provide an accuracy of $10^{-4}$ with an
efficiency comparable to optimized, implicit Matlab solvers.

###  4.2 Strong stiffness - autocatalytical chemical reactions

When performing a kinetic analysis of autocatalytic reactions, the following
system of equations emerge [11]:

$\begin{split}\frac{dx}{dt}=-ax+byz\\\ \frac{dy}{dt}=ax-byz-cy^{2}\\\
\frac{dz}{dt}=cy^{2}\end{split}$ (20)

with $a=0.04$, $b=10^{4}$, and $c=3\cdot 10^{7}$. The initial conditions are
$(x,y,z)=(1,0,0)$. Clearly, the time scales of the reactions are very
different. This generally causes no problems for short intervals, but for very
long intervals, like $[10^{-6},10^{6}]$ explicit methods are likely to fail
because of the stiffness caused by the large difference between the (reaction
rate) constants. As can be seen from the time-adaptive GWRM solution in Figure
4 (note the logarithmic time axis), $y(t)$ features a strong transient
initially , whereafter smooth solutions hold for all three variables
$x(t),y(t)$ and $z(t)$.

Figure 4: Solution of the Robertson equations (20), employing accuracy
$\epsilon=0.001$. Note that $y(t)$ has been magnified a factor $10^{4}$. The
solution was obtained by a GWRM solver.

In order to formally study stiffness, we employ the method earlier described.
We thus let $\delta u_{1}\equiv
u_{1}(t;\epsilon)-u_{1}(t;0)=x(t;\epsilon)-x(t;0)\equiv\delta x$ and proceed
similarly for the $y$\- and $z$\- components. Subtraction of the unperturbed
and perturbed equations yields, to lowest order in $\delta u_{i}$;

$\begin{split}\frac{d\delta x}{dt}=-a\delta x+b(y\delta z+z\delta y)\\\
\frac{d\delta y}{dt}=a\delta x-b(y\delta z+z\delta y)-2cy\delta y\\\
\frac{d\delta z}{dt}=2cy\delta y\end{split}$ (21)

The coefficients for $\delta u_{i}$, including $x,y,z$, are assumed constant
in each interval $[t_{n},t_{n+1}]$ that is solved for. They do, of course,
take on different values in each interval. We now have a system of ODE’s, so
it is not apparent how the Lyapunov exponents can be determined for deciding
stiffness. Since the system is linear, the following theorem may be applied.
Assume the system of ODE’s

$\frac{d\bf u}{dt}=A\bf u$ (22)

where $\textbf{u}\in R^{n}$ and $\bf{A}$ is a constant, quadratic,
diagonalizable matrix with distinct eigenvalues ${\gamma}_{k}\in C$ and
eigenvectors $\textbf{c}_{k}\in C^{n}$. The solution to Eq. (22) is

$\textbf{u}(t)=\sum_{k=1}^{n}a_{k}e^{\gamma_{k}}\textbf{c}_{k}$ (23)

where $a_{k}$ are arbitrary constants. Comparing with the system of Eqs. (21)
we have $\textbf{u}=\delta\textbf{u}$ and $A=J$, where $J$ is the Jacobian of
the system. Since we are mainly interested in local Lyapunov exponents for the
intervals $[t_{n},t_{n+1}]$, we designate the eigenvalues with $\gamma_{k}$,
rather than as $\lambda_{k}$. These are obtained from the characteristic
equation

$|J-\gamma I|=0$ (24)

where $I$ is the identity matrix. In this case

$J=\begin{pmatrix}-a&bz&by\\\ a&-bz-2cy&-by\\\ 0&2cy&0\end{pmatrix}$

We may solve Eq. (24) for the eigenvalues numerically. For the initial values,
the Robertson equations are not stiff; we obtain $\gamma_{1}=0$,
$\gamma_{2}=0$ and $\gamma_{3}=-0.04$ from Eq. (24). But already at
$y=10^{-6}$ (maintaining $x=1$ and $z=0$) the system becomes stiff
($\kappa_{i}$ are large); $\gamma_{1}=0$, $\gamma_{2}=-0.05$ and
$\gamma_{3}=-60$. Extreme stiffness is obtained for $t>10^{-4}$, as
$y>0.4\cdot 10^{-4}$ and $\gamma_{3}<-2400$. At this point, $(\Delta
t)_{stiff}\approx 2\sqrt{2}/|\gamma_{3}|<10^{-3}$, decreasing further for
longer times, and the time step is completely determined by the stiffness of
the solution.

An RK4 algorithm with adaptive step size has been employed to solve the
problem (20). After over 61 000 time steps the algorithm stagnates near
$t=96$. Similar behaviour was encountered for other parameters. It is clear
that explicit time-stepping methods are not suited for stiff problems of this
type.

Thus an implicit standard trapezoid algorithm with adaptive step size was
applied. For an accuracy $\epsilon=0.001$, a solution was obtained, employing
118 steps for the whole interval, with initial step length 0.1. Employing a
similar step size adaption algorithm for the same accuracy, the GWRM time-
spectral method found the solution about 80 times faster, using 49 time
intervals.

##  5 Transforming stiff and chaotic differential equations

Local Lyapunov stability analysis shows that properties like stiffness and
chaoticity depend on the form of the differential equations. Thus, we may ask
if it is possible to transform the differential equations to produce
asymptotically stable solutions. These could then be efficiently and
accurately integrated numerically, whereafter a back-transformation is made to
the sought solutions. We will find that chaotic equations can be successfully
transformed, whereas stiff equations cannot. In the following, this assertion
will be demonstrated, employing a general formulation.

For simplicity we wish to find asymptotically stable solutions to the
following single ODE in a time interval $[0,T]$:

$\frac{du}{dt}=f(t,u),\quad u(0)=u_{0}$ (25)

where $f$ is an arbitrary function of the independent variable $t$ and the
variable $u(t)$, time-normalized so that $T$ is of order unity. We assume Eq.
(25) to be either stiff (with neighbouring solutions $\propto
exp(\kappa_{0}t)$, with $\kappa_{0}<0$ and $|\kappa_{0}|>>1$) or chaotic
(neighbouring solutions $\delta u(t)\propto exp(\kappa_{0}t)$, with
$\kappa_{0}>0$). Strictly, chaotic ODEs feature at least three degrees of
freedom, and are thus represented by at least three ODEs [Sprott], but, as
will become apparent, the same line of reasoning holds.

In a local Lyapunov stability analysis we can assume that $\kappa$ is
approximately constant in each solution interval $[t_{n},t_{n+1}]$. For an
asymptotically stable solution, it should hold that $\kappa<0$ and that
$|\kappa t|>>1$ does not hold.

###  5.1 General transformation

A general transformation for solving Eq.(25) can be written

$u(t)\equiv g(t,z(t))$ (26)

where $g$ may be chosen arbitrarily. In the following, the main focus will be
on stiffness. Thus we may, using Eq. (12), assume that $Q_{u}(t)>>1$ at least
in some solution intervals, that is for a number of time steps. The goal is
now to find an ODE for $z(t)$ for which $Q_{z}(t)\leq 1$ everywhere,
guaranteeing asymptotic stability of the solution. From Eqs. (25) and (26) we
find

$\frac{dz}{dt}=-\frac{\partial g}{\partial t}/\frac{\partial g}{\partial
z}+{f(t,g(t,z))}/\frac{\partial g}{\partial z}\quad.$ (27)

In order to determine

$\frac{d\delta z}{dt}=\frac{d(z+\delta z)}{dt}-\frac{dz}{dt}$ (28)

we use the following linearizations

$\frac{\partial g}{\partial t}(t,z+\delta z)=\frac{\partial g}{\partial
t}(t,z)+\frac{\partial^{2}g}{\partial z\partial t}(t,z)\delta z$ (29)

$\frac{\partial g}{\partial z}(t,z+\delta z)=\frac{\partial g}{\partial
z}(t,z)+\frac{\partial^{2}g}{\partial z^{2}}(t,z)\delta z$ (30)

$f(t,g(t,z+\delta z))=f(t,g(t,z))+\frac{\partial f}{\partial g}\frac{\partial
g}{\partial z}\delta z\quad.$ (31)

Expanding $(\frac{\partial g}{\partial z}(t,z+\delta z))^{-1}$ there results,
to first order,

$\frac{d\delta z}{dt}=[(\frac{\partial g}{\partial
t}-f)\frac{\nicefrac{{\partial^{2}g}}{{\partial z^{2}}}}{(\nicefrac{{\partial
g}}{{\partial z}})^{2}}-\frac{\nicefrac{{\partial^{2}g}}{{\partial t\partial
z}}}{\nicefrac{{\partial g}}{{\partial z}}}+\frac{\partial f}{\partial
g}]\delta z\equiv\kappa_{g}\delta z\quad.$ (32)

with the approximate dependence, for $\kappa_{g}$ $\approx$ constant within a
time interval,

$\delta z(t)\sim e^{\kappa_{g}t}$ (33)

We may also write

$\frac{d\delta u}{dt}=\frac{\partial f}{\partial u}\delta
u\equiv\kappa_{f}\delta u$ (34)

Now, since

$\frac{\partial g}{\partial t}+\frac{\partial g}{\partial z}\frac{dz}{dt}=f,$
(35)

we find from Eq. (32), for $\partial^{2}g/\partial z^{2}\neq 0$, using
$\partial f/\partial g=\partial f/\partial u=\kappa_{f}$,

$\frac{dz}{dt}=(\kappa_{f}-\kappa_{g})\frac{\nicefrac{{\partial g}}{{\partial
z}}}{\nicefrac{{\partial^{2}g}}{{\partial
z^{2}}}}-\frac{\nicefrac{{\partial^{2}g}}{{\partial t\partial
z}}}{\nicefrac{{\partial^{2}g}}{{\partial z^{2}}}}\quad.$ (36)

For an asymptotically stable solution $z(t)$, it holds that $\kappa_{g}<0$ and
$|\kappa_{g}\lesssim 1|$. Assuming that Eq. (25) is stiff, $\kappa_{f}<0$ and
$|\kappa_{f}>>1|$. Thus $\kappa_{g}-\kappa_{f}>>1$ and the ODE in Eq. (36) may
be stiff, or strongly chaotic, depending on the choice of $g(t,z(t))$.

###  5.2 General transformation, first order

The case $\partial^{2}g/\partial z^{2}=0$ implies that $g$ is linear in $z$;

$g(t,z(t))=A(t)z(t)+B(t)\enskip.$ (37)

Then we, from Eqs. (32) and (34), obtain

$\frac{-dA/dt}{A}+\frac{\partial f}{\partial
g}=\frac{-dA/dt}{A}+\kappa_{f}=\kappa_{g},$ (38)

with formal solution

$A(t)=A(0)e^{(\kappa_{f}-\kappa g)t}.$ (39)

Since $-(\kappa_{f}-\kappa_{g})\approx-\kappa_{f}>>1$, $A(t)$ features, as
indicated by Eq. (34), a strong local curvature similar to that characterizing
the stiffness of Eq. (25). Consequently the ODE for $z(t)=(u(t)-B(t))/A(t)$
will be correspondingly difficult to resolve numerically. The result of the
transformation (37) is simply a transition from a stiff ODE for $u(t)$, with
typical time step $(\Delta t)^{u}_{stiff}$, to an asymptotically stable ODE
for $z(t)$, with typical time step $(\Delta t)^{z}_{max}$, for which it
however holds that

$\mathcal{O}((\Delta t)^{z}_{max})=\mathcal{O}((\Delta t)^{u}_{stiff})$ (40)

This implies that, using the linear transformation (37), neither the linear
nor the nonlinear stiff ODEs of Sections 3 and 4 can be transformed to ODEs
that can be more efficiently solved for explicit finite difference schemes.
Asymptotic stability is attained at the cost of degraded numerical resolution
of the sought solution. For the stiff linear ODE in Section 3 we have
$\partial f/\partial u=\kappa_{f}=-a$ so that $A(t)\approx A(0)e^{-at}$. For
the stiff nonlinear ODE in Section 4, modelling explosive combustion, the
time-normalized ($t=(2/d)\tilde{t}$) relation $\partial f/\partial
u=(2/d)(2u-3u^{2})\approx-2/d$ holds in the stiff interval $\tilde{t}>1/2$
where $u\approx 1$ so that $A(\tilde{t})\approx A(0)e^{-(2/d)\tilde{t}}$. In
both cases $A(t)$, and thus $z(t)$, require high temporal resolution for given
accuracy.

###  5.3 General transformation, second order

Solution of $z(t)$ from Eq. (36) requires that $g(t,z(t))$ is expanded to
second order. We obtain

$\frac{\partial g}{\partial z}\approx A+B(t-t_{0})+C(z-z_{0})$ (41)

where

$A\equiv\frac{\partial g}{\partial
z}(t_{0},z_{0}),B\equiv\frac{\partial^{2}g}{\partial t\partial
z}(t_{0},z_{0}),C\equiv\frac{\partial^{2}g}{\partial z^{2}}(t_{0},z_{0}).$
(42)

Considering individual time intervals $[t_{n},t_{n+1}]$, we can let
$t_{0}=t_{n}=0$ without losing generality. The notation $z_{0}=z(t_{0})=z(0)$
is also used. Inserting Eq. (41) into Eq. (36), the following solution is
obtained, exact to second order;

$z(t)=\frac{A}{C}(e^{(\kappa_{f}-\kappa_{g})t}-1)-\frac{B}{C}t+z_{0}.$ (43)

Recognizing that $-(\kappa_{f}-\kappa_{g})>>1$, the temporal resolution
required is again approximately the same as that required due to the stiffness
of the problem.

In order to also study the effect of explicit $t$-dependence in $g(t,z(t))$,
the following transformation is useful:

$g(t,z(t))=g_{0}e^{\kappa_{z}t}z^{\alpha}$ (44)

where $g_{0}$, $\kappa_{z}$ are constants and $\alpha>1$. Linearizing in $t$,
using Eqs. (41)-(43), there obtains

$z(t)\approx z(0)e^{((\kappa_{f}-\kappa_{g}-\kappa_{z})/(\alpha-1))t},$ (45)

Choosing $\kappa_{z}$ so that $\kappa_{f}-\kappa_{g}-\kappa_{z}$ becomes
small, requires that $\kappa_{z}\approx\kappa_{f}$, which, according to Eq.
(44), violates the assumption that $u(t)=g(t,z(t))$ should be a smooth
function of $t$.

That explicit time dependence in $g(t,z(t))$ cannot compensate for the effect
of the term $\kappa_{f}-\kappa_{g}$ in $z(t)$ can also be seen by setting
$dz/dt\approx 0$ in Eq. (36). It follows that

$(\kappa_{f}-\kappa_{g})\frac{\partial g}{\partial
z}\approx\frac{\partial^{2}g}{\partial t\partial z}$ (46)

The solution to this equation is obtained by integration in $t$ and $z$;

$g(t,z(t))=e^{(\kappa_{f}-\kappa_{g})t}(\int_{0}^{t}e^{-(\kappa_{f}-\kappa_{g})t^{\prime}}T(t^{\prime})dt^{\prime}+Z(z))$
(47)

where $T(t)$ and $Z(z)$ are arbitrary functions. Since we may assume $Z(z)\neq
0$, the last term will be non-smooth and again the problem requires a
numerical time resolution amounting to that needed for resolving stiffness.

If the problem (25) on the other hand were chaotic, so that $\kappa_{f}>0$, we
could maintain the asymptotically stable conditions $\kappa_{g}<0$ and
$|\kappa_{g}\lesssim 1|$ for the numerical solution $z(t)$ by a suitable
choice of $\kappa_{z}$. It should be noted that the function $g(t,z(t))$ must
be assumed given in order to solve for $z(t)$.

##  6 Removing numerical chaoticity

The question thus arises how the numerical errors caused by chaotic
differential equations can be mitigated locally, using variable
transformations. The idea is to, via a suitable transformation, temporarily
remove the exponential divergence between neighbouring solutions in parameter
space for the new variable. The sensitivity to initial conditions at each new
time step is thus alleviated. Due to the nonlinearity of the equations, the
net result is substantially higher accuracy than what can be attained by
standard time-stepping solvers like explicit RK4. We will now illustrate how
accuracy can be substantially improved for a well known set of chaotic ODEs.

###  6.1 The Lorenz 1984 chaotic equations

The Lorenz 1984 model [12] is a set of three non-linear ordinary differential
equations that features chaotic, rather than stiff, behaviour related to
meteorological systems. It is one of the simplest models of Hadley
circulation. Thus, whereas the Lorenz 1984 model is not an accurate numerical
weather prediction (NWP) model, it is of interest for numerical analysis of
simple, yet non-linear chaotic behaviour that may occur in NWP. The equations
are

$\displaystyle\frac{dx}{dt}$ $\displaystyle=-y^{2}-z^{2}-ax+aF,$ (48)
$\displaystyle\frac{dy}{dt}$ $\displaystyle=xy-bxz-y+G,$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=bxy+xz-z$

Further details of the model are found in [12]. We note here that the
variables $x(t)$, $y(t)$, and $z(t)$ represent certain meteorological systems
such as wind currents and large-scale eddies. The coefficients $a$, $b$, $F$,
and $G$ are chosen within certain limits to act as damping, coupling, and
amplification of the physical processes of the system. In Figure 5, a solution
for $a=0.25$, $b=4.0$, $F=8.0$ and $G=1.0$ is shown. The initial conditions
are $(x,y,z)=(0.96,-1.1,0.5)$.

Figure 5: Solution of the Lorenz equations (48) with parameters $a=0.25$,
$b=4.0$, $F=8.0$ and $G=1.0$. The initial conditions are
$(x,y,z)=(0.96,-1.1,0.5)$. From top to bottom; $x(t)+8$, $y(t)+4$, $z(t)$. A
high accuracy, time-spectral GWRM code was used.

This set of equations is difficult to solve for any explicit time-stepping
method, since chaoticity causes any error in the solutions at time $t$ to be
strongly amplified in the interval $[t_{n},t_{n+1}]$. To see this, we analyze
the eigenvalues as in the previous section, thus obtaining the system

$\begin{split}\frac{d\delta x}{dt}=-2y\delta y-2z\delta z-a\delta x\\\
\frac{d\delta y}{dt}=x\delta y+y\delta x-bx\delta z-bz\delta x-\delta y\\\
\frac{d\delta z}{dt}=bx\delta y+by\delta x+x\delta z+z\delta x-\delta
z\end{split}$ (49)

The corresponding Jacobian becomes

$J=\begin{pmatrix}-a&-2y&-2z\\\ y-bz&x-1&-bx\\\ by+z&bx&x-1\end{pmatrix}$ (50)

For the initial values, the Lorenz equations are chaotic. By solving the
characteristic equation $|J-\gamma I|=0$ for $t=0$ there obtains
$\gamma_{1}=1.9$, $\gamma_{2}=-1.1+4.5i$ and $\gamma_{3}=-1.1-4.5i$. It is
seen from Eq. (49) that the parameters $F$ and $G$ do not affect the character
of the ODE’s. The system is found to be robust in the sense that chaoticity is
always manifested by at least one eigenvalue being positive for any value of
$b$. For $a\geq 3.1$, however, the system will initially be non-chaotic and
may also locally be non-chaotic, for example for $t\approx 13-14$. The
equations are not stiff for any values of $a$ and $b$. Stiffness and
chaoticity are indeed different properties.

It is of interest to determine the local eigenvalues of the characteristic
equation for the full solution interval $t\in[0,30]$. In Figure 6, these have
been computed at 400 equidistant times. It is seen that at least one local
Lyapunov exponent is positive, indicating chaoticity, in the entire interval
except for in small intervals near $t=14$.

Figure 6: Eigenvalues $\gamma$ of the characteristic equation $|J-\gamma I|=0$
for the Lorenz local Lyapunov exponent system (49) in the time interval
$[0,30]$. The top curve represents the real root, whereas the bottom curve
images the real part of the remaining two, complex conjugated roots.

###  6.2 Transforming the Lorenz equations to eliminate numerical chaoticity

Recognizing the exponential dependence of chaoticity we choose the following,
simple transformations for the Lorenz equations, where $A_{i}(t)$ are explicit
functions of $t$, and $\epsilon_{i}$ as well as $\mu_{i}$ are constants:

$\begin{split}x(t)=A_{1}(t)z_{1}(t)=\epsilon_{1}e^{\mu_{1}t}z_{1}(t)\\\
y(t)=A_{2}(t)z_{2}(t)=\epsilon_{2}e^{\mu_{2}t}z_{2}(t)\\\
z(t)=A_{3}(t)z_{3}(t)=\epsilon_{3}e^{\mu_{3}t}z_{3}(t)\end{split}$ (51)

Solving for the $z$ variables, Eq. (48) becomes

$\begin{split}\frac{dz_{1}}{dt}=-\mu_{1}z_{1}-\frac{\epsilon_{2}^{2}}{\epsilon_{1}}e^{(2\mu_{2}-\mu_{1})t}z_{2}^{2}-\frac{\epsilon_{3}^{2}}{\epsilon_{1}}e^{(2\mu_{3}-\mu_{1})t}z_{3}^{2}-az_{1}+\frac{a}{\epsilon_{1}}e^{-\mu_{1}t}F\\\
\frac{dz_{2}}{dt}=-\mu_{2}z_{2}+\epsilon_{1}e^{\mu_{1}t}z_{1}z_{2}-b\frac{\epsilon_{1}\epsilon_{3}}{\epsilon_{2}}e^{(\mu_{1}-\mu_{2}+\mu_{3})t}z_{1}z_{3}-z_{2}+\frac{G}{\epsilon_{2}}e^{-\mu_{2}t}\\\
\frac{dz_{3}}{dt}=-\mu_{3}z_{3}+b\frac{\epsilon_{1}\epsilon_{2}}{\epsilon_{3}}e^{(\mu_{1}+\mu_{2}-\mu_{3})t}z_{1}z_{2}+\epsilon_{1}e^{\mu_{1}}z_{1}z_{3}-z_{3}\\\
\end{split}$ (52)

This system will be solved with a standard RK4 algorithm, using $N$
equidistant time steps for $t\in[0,T]$. In order to avoid excessively large
numerical exponential factors, the time domain is separated into $K$ intervals
with $N/K$ time steps in each. In all cases to follow $N=600$. At the end of
each interval, new initial conditions for $z_{i}(t)$ are generated using Eq.
(51). For simplicity, we here set $\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=1$.

####  6.2.1 Method 1

The problem is now to optimize the choice of $\mu_{i}$ in order to obtain
asymptotically stable, rather than chaotic, solutions $z_{i}(t)$. Given the
knowledge of the local Lyapunov exponents displayed in Figure 6, one may be
tempted to simply set all $\mu_{i}=\mu_{0}$, where $\mu_{0}$ is a positive
constant, for example in the interval $[2,3]$. This approach may not be
effective, since optimal choices of $\mu_{i}$ are expected to be
interdependent. A number of RK4 code runs were run manually, individually
varying $\mu_{i}$. To compute the local Lyapunov exponents for the transformed
equations (52), we determine, as before, the system

$\begin{split}\frac{d\delta z_{1}}{dt}=-\mu_{1}\delta z_{1}-2z_{2}\delta
z_{2}-2z_{3}\delta z_{3}-a\delta z_{1}\\\ \frac{d\delta
z_{2}}{dt}=-\mu_{2}\delta z_{2}+z_{1}\delta z_{2}+z_{2}\delta
z_{1}-bz_{1}\delta z_{3}-bz_{3}\delta z_{1}-\delta z_{2}\\\ \frac{d\delta
z_{3}}{dt}=-\mu_{3}\delta z_{3}+bz_{1}\delta z_{2}+bz_{2}\delta
z_{1}+z_{1}\delta z_{3}+z_{3}\delta z_{1}-\delta z_{3}\end{split}$ (53)

We make the approximation $t=0$ in Eq. (52), corresponding to short interval
lengths $t_{n+1}-t_{n}$. The eigenvalues (local Lyapunov exponents) $\gamma$
are thus obtained from the determinant $|J^{*}-\gamma I|=0$ where

$J^{*}=\begin{pmatrix}-a-\mu_{1}&-2z_{2}&-2z_{3}\\\
z_{2}-bz_{3}&z_{1}-1-\mu_{2}&-bz_{1}\\\
bz_{2}+z_{3}&bz_{1}&z_{1}-1-\mu_{3}\end{pmatrix}$ (54)

In Figure 7 the eigenvalues $\gamma$ for the choice $\mu_{1}=2.592$,
$\mu_{2}=1.944$ and $\mu_{3}=1.539$ (constant throughout the computations) are
shown. It is seen that the local Lyapunov exponents have been reduced
substantially, so that the system (52), solved with retained explicit time
dependence for accuracy, becomes non-chaotic over a large fraction of the time
domain $[0,30]$. The ratio $N/K=10$ (10 time steps in each time interval) was
chosen. We term this approach Method 1. Three additional methods are discussed
in the following, whereafter all methods will be compared in Figure 8.

Figure 7: Eigenvalues $\gamma$ of the characteristic equation $|J^{*}-\gamma
I|=0$ for the Lorenz local Lyapunov exponent system (53) with
$\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=1$, using $\mu_{1}=2.592$,
$\mu_{2}=1.944$ and $\mu_{3}=1.539$ in the time interval $[0,30]$. The top
curve represents the real root, whereas the bottom curve images the real part
of the remaining two, complex conjugated roots. A large fraction of the
temporal domain is non-chaotic.

####  6.2.2 Method 2

It may be noted that, in the determinant $|J^{*}-\gamma I|=0$, the eigenvalues
and the transformation parameters will occur as terms $\mu_{i}+\gamma$ only.
This shows that i) the chosen transformation (51) features a high degree of
simplicity and interpretability and ii) that $\mu_{i}$ typically should be
chosen sufficiently positive to reduce the real parts of the eigenvalues
$\gamma$ towards negative values, corresponding to asymptotically stable
solutions. A complication is that the interrelation between the constants
$\mu_{i}$ and the eigenvalues is unclear, since for each choice of
$(\mu_{1},\mu_{2},\mu_{3})$ there generally exists three eigenvalues $\gamma$.

An option would be to simply set all eigenvalues equal (multiple roots), to a
suitable negative value for asymptotic stability and solve for
$(\mu_{1},\mu_{2},\mu_{3})$. The corresponding, multiple root solutions would
then be of the form $\delta z_{i}=(c_{2i}t^{2}+c_{1i}t+c_{0i})e^{\gamma t}$,
where the $c$’s are constants. Assuming short time intervals, $\delta
z_{i}\approx c_{0i}e^{\gamma t}$ where the time dependence is governed by the
exponential (Lyapunov-like) term alone. This approach is, however, only
moderately successful with respect to minimizing the error for longer times
($t>10$). This is due to the rapid time variation of the eigenvalues, as can
be seen in Figs. 6 and 7; knowledge of the LLEs is not sufficiently helpful
for providing asymptotic stability in the subsequent time interval.

Method 2 is an algorithm, that instead compensates for local chaoticity. The
maximum real part of the three eigenvalues $\gamma_{max}(t_{n})$ of
$|J^{*}-\gamma I|=0$ is computed in each time interval $[t_{n},t_{n+1}]$, and
the values $\mu_{i}=q\gamma_{max}(t_{n})$, where $q$ is a constant, are used
for the subsequent interval. For the case shown here, the initial values were
$\mu_{1}=2.0,\mu_{2}=2.0,\mu_{3}=2.0$ and $q=1.0$ was used for $N/K=10$.

####  6.2.3 Methods 3 and 4

In Method 3, $\mu_{i}$ values for subsequent time intervals are computed from
accumulated averages $<\gamma_{max}>(t_{n})$ of the maximum eigenvalues
$\gamma_{max}(t_{n})$. High accuracy is obtained for the choice
$\mu_{1}(t_{n})=-1.5<\gamma_{max}>(t_{n}),\mu_{2}(t_{n})=-0.66<\gamma_{max}>(t_{n}),\mu_{3}(t_{n})=-0.5<\gamma_{max}>(t_{n})$,
employing $\mu_{1}(t_{1})=2.16,\mu_{2}(t_{1})=1.62,\mu_{3}(t_{1})=1.28$.

High accuracy also results from Method 4, where averages
$<\gamma_{max}>(t_{n})$ are taken from the last two previous time intervals
only. Both methods employ 40 time steps in each time interval, that is
$N/K=40$.

####  6.2.4 Method comparisons

In Figure 8 the absolute errors for $x(t)$ of Eq. (48), computed by the
transformation Methods 1-4, are compared with that of standard RK4 solution.
It is seen that radical improvements in accuracy is obtained using the
transformations. Methods 3 and 4 are about 100 times more accurate in the
entire interval than standard RK4, and display only a weak deterioration with
time. In Figure 9 the exact $x(t)$ solution, with same parameters as in Figure
5, is compared to that obtained from standard RK4 and Method 3, for the
numerical parameters given earlier.

It should be noted from Figure 8 that although RK4 accuracy is substantially
improved by the four transformation methods discussed, asymptotic stability
(error tending to zero) is not achieved since the Lorenz 1984 equations are
inherently chaotic; any error in the solution will be amplified for higher
times even for infinite numerical accuracy. The effect of the transformations
is to render the $z(t)$ solutions asymptotically stable, thereby nonlinearly
mitigating errors in the numerical, explicit time-stepping algorithm.

Figure 8: Absolute errors for RK4 solution of $x(t)$ from the Lorenz equations
(48), with same parameters as in Figure 5, using standard RK4 (top, black
points) and for 4 transformation methods, as described in the text: Method 1
(orange), Method 2 (brown), Method 3 (blue), Method 4 (green). Figure 9:
$x(t)$ obtained from the Lorenz equations (48), with Figure 5 parameters:
exact solution (points), standard RK4 (squares), and Method 3 (blue).

####  6.2.5 Time step extension

We now seek to investigate whether our objective has been achieved, namely,
that the influence of chaotic behavior has been reduced to the point where the
achieved level of accuracy is primarily constrained by numerical approximation
thresholds rather than by time step limitations caused by chaotic dynamics. To
investigate this, we compare the local $(\Delta t)_{max}(t)$, given by Eq. (7)
and the curvature of the exact solution, with the step size
$\Delta=30/600=0.05$, employed in the present computations.

The curvature is computed from the most non-smooth solution $Z(t)$, rather
than from $X(t)$ or $Y(t)$. We use the maximum error $(\epsilon=0.012)$ of
Method 3 in the time interval $[0,30]$ for computing $(\Delta t)_{max}(t)$. It
is clear, from Figure 10, that $\Delta$, being fixed in the interval, could
not have been longer, since then the requirement $(\Delta
t)_{max}(t)\geq\Delta$ would be locally violated. Inserted is also $(\Delta
t)_{max}(t)$ as obtained from standard RK4, for which $\epsilon=1.0$. Clearly,
the limited accuracy for standard RK4 is determined by chaoticity, rather than
by $(\Delta t)_{max}(t)$. Thus we may conclude that the transformations of
Methods 3 and 4 remove the time step limitations that are associated with the
chaoticity of the Lorenz equation (48).

This conclusion is further supported by Figure 11, where Method 3 is applied
using a larger number of time steps $N=1620$ ($\Delta=0.0185$) with $K=60$,
overall obtaining an accuracy $\epsilon=0.001$. It is again seen that longer
$\Delta$ would violate the maximum allowable time step for numerical
resolution $(\Delta t)_{max}(t)$ at this accuracy.

Figure 10: Comparison of $(\Delta t)_{max}(t)$, given by Eq. (7) and obtained
from the curvature of the exact solution for the overall accuracy
$\epsilon=0.012$ obtained by Method 3, with step size $\Delta=30/600=0.05$.
The curvature is computed from $Z(t)$, being the most non-smooth solution.
Inserted is also $(\Delta t)_{max}(t)$ (in blue) obtained by standard RK4,
where $\epsilon=1.0$.  Figure 11: $(\Delta t)_{max}(t)$ as in Figure 10, but
for Method 3 only, employing overall accuracy $\epsilon=0.001$. Obtained for
$N=1620$ ($\Delta=0.0185$, also displayed in the figure) with $K=60$.

##  7 Discussion

Chaoticity, when appearing in systems of differential equations, is inherent
and cannot be transformed away. But the numerical errors introduced by the
discreteness of time-stepping methods can be mitigated by transformations to
systems that produce asymptotically stable solutions. In these systems, the
tendency for neighbouring solutions to converge to each other compensates for
the sensitivity to initial conditions that characterize chaotic systems. In
Section 6.2 we suggest an analytically simple transformation method along
these lines. Based on local Lyapunov exponents, it is applicable to time
intervals consisting of a single or a number of time steps. At the end of each
interval, a back transformation to the original variables are used, producing
new initial conditions.

The back transformation involves an exponential amplification of the numerical
errors but, as we have seen by solving the chaotic Lorenz 1984 ODEs, the net
result is substantially higher accuracy than what can be attained by standard
time-stepping solvers like explicit fourth order Runge-Kutta. Employing four
different methods for determining the transformation parameters $\mu_{i}$, it
was found that the time step length was no longer determined by chaoticity but
by discrete, numerical resolution of the exact solution at the given accuracy.

Interesting topics for further research is to find robust methods for
determining parameters $\mu_{i}$ that result in accurate solutions for a
variety of chaotic problems. The transformation of Section 6.2 should be
generalizable also to PDEs. Thus, problems in areas like numerical weather
prediction and turbulence, which for long have been plagued by time step
limitations, could be addressed.

Turning to stiff systems, it is no major surprise that these problems cannot
be remedied by transformations. In the formal analysis presented it also
becomes clear, by analysis of the ODEs governing the local Lyapunov exponents,
that the original difficulty to numerically resolve neighbouring, fast
declining solutions, is transferred to the numerical resolution of the
transformed functions. The required analytical tools $(\Delta t)_{stiff}$ and
$(\Delta t)_{max}$ were introduced in Section 2.

Fixed time step lengths are used in this work in order to isolate the effect
of the transformations alone. It is expected that adaptive step length
algorithms, in combination with suitable methods for choosing transformation
parameters (like $\mu_{i}$) could help to optimize the effect of the
transformation, thus enabling even longer time steps for solution of chaotic
ODEs and PDEs with retained accuracy.

##  8 Conclusion

Transform methods for improving accuracy of explicit time-stepping ODE solvers
are introduced. It is shown that the numerical accuracy of solutions to
chaotic ODEs can be improved by at least two orders of magnitude, at retained
step length, by transforming to ODEs with asymptotically stable solutions. The
chaotic Lorenz 1984 equations serves as an example.

For stiff problems, however, accuracy remains a problem when transforming to
problems with asymptotically stable solutions. Both linear and nonlinear
examples of stiff systems are studied. The results rest on the concept of
local Lyapunov exponents (LLEs), for which we have defined a measure $Q(t)$
for diagnosing asymptotic stability, stiffness and chaos. Applications to
stiff linear and nonlinear ODEs are given.

## References

  * [1] M. W. Hirsch, S. Smale, Differential Equations, Dynamical Systems and Linear Algebra, Academic, New York, 1974.
  * [2] J. C. Sprott, Some simple chaotic jerk functions, Am. J. Phys. 65 (6) (1997) 537–543.
  * [3] J. H. Cartwright, Nonlinear stiffness, Lyapunov exponents, and attractor dimension, Physics Letters A 264 (1999) 298–302.
  * [4] J. Lambert, Numerical Methods for Ordinary Differential Systems, New York: Wiley, 1992.
  * [5] J. Scheffel, A Spectral Method in Time for Initial-Value Problems, American Journal of Computational Mathematics 02 (2012) 173–193.  doi:10.4236/ajcm.2012.23023.
  * [6] M. T. Heath, Scientific Computing - an Introductory Survey, McGraw-Hill, 2005.
  * [7] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential - Algebraic Problems, Springer, 2nd ed. 1996.
  * [8] D. Ayers, J. Lau, J. Amezcua, A. Carrassi, V. Ojha, Supervised machine learning to estimate instabilities in chaotic systems: estimation of local Lyapunov exponents, arXiv, (2022).  doi:10.48550/ARXIV.2202.04944.
URL https://arxiv.org/abs/2202.04944

  * [9] J. Scheffel, K. Lindvall, H. F. Yik, A time-spectral approach to numerical weather prediction, Computer Physics Communications 226 (2017) 127–135.  doi:10.1016/j.cpc.2018.01.010.
  * [10] J. Scheffel, Time-spectral solution of initial-value problems, in: Partial Differential Equations: Theory, Analysis and Applications, Nova Science Publishers, Inc., 2011, pp. 1–49.
  * [11] H. H. Robertson, The solution of a set of reaction rate equations, in: J. Walsh (Ed.), Numerical Analysis: An Introduction, Academic Press, London, 1966, pp. 178–182.
  * [12] E. N. Lorenz, Irregularity: a fundamental property of the atmosphere, Tellus A 36 (1984) 98–110.  doi:10.3402/tellusa.v36i2.11473.


# Global Vegetation Modeling with
Pre-trained Weather Transformers

Pascal Janetzky, Florian Gallusser, Simon Hentschel, Andreas Hotho, Anna
Krause
Data Science Chair, Center for Artificial Intelligence and Data Science
(CAIDAS),
University of Würzburg
{janetzky,gallusser,hentschel,hotho<EMAIL_ADDRESS>

###### Abstract

Accurate vegetation models can produce further insights into the complex
interaction between vegetation activity and ecosystem processes. Previous
research has established that long-term trends and short-term variability of
temperature and precipitation affect vegetation activity. Motivated by the
recent success of Transformer-based Deep Learning models for medium-range
weather forecasting, we adapt the publicly available pre-trained FourCastNet
to model vegetation activity while accounting for the short-term dynamics of
climate variability. We investigate how the learned global representation of
the atmosphere’s state can be transferred to model the normalized difference
vegetation index (NDVI). Our model globally estimates vegetation activity at a
resolution of $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ while relying only
on meteorological data. We demonstrate that leveraging pre-trained weather
models improves the NDVI estimates compared to learning an NDVI model from
scratch. Additionally, we compare our results to other recent data-driven NDVI
modeling approaches from machine learning and ecology literature. We further
provide experimental evidence on how much data and training time is necessary
to turn FourCastNet into an effective vegetation model. Code and models will
be made available upon publication.

##  1 Introduction

Environmental changes affect the dynamics of terrestrial vegetation, which is
involved in controlling water, energy and CO2 fluxes (Richardson et al.,
2013), and is thus crucial for providing ecosystem services such as food,
fiber and fuel (Piao et al., 2020). Hence, a profound understanding of the
complex interplay of climate system variables and vegetation changes is
desirable to achieve sustainable ecological management.

Previous studies have shown that observed changes in vegetation can be
attributed to both long-term and short-term changes in temperature and
precipitation, i.e., climate change and climate variability (Burrell et al.,
2020; Chen et al., 2019; Higgins et al., 2023; Liu et al., 2022; Seddon et
al., 2016; Zhu et al., 2016). While the spatial arrangement of vegetation on a
large scale is primarily dictated by climatic factors, the interplay between
climate variability and the short-term dynamics of vegetation introduces a
higher level of complexity (Papagiannopoulou et al., 2017; Pelletier et al.,
2015). Different Machine Learning (ML) approaches have been suggested to
capture the complex nonlinear dynamics of those short-term dynamics. However,
the employed models are either limited to a specific region (Robin et al.,
2022; Smith et al., 2023) or use a coarse global resolution with one pixel
covering at least $0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$($\approx$
$55\text{\,}\mathrm{km}$) (Chen et al., 2021; Kraft et al., 2019). While there
are statistical approaches that globally quantify the effect of climate
variability on vegetation change on a finer spatial resolution up to
$0.083\text{\,}\mathrm{\SIUnitSymbolDegree}$, they only consider
meteorological data on a coarse time scale, e.g., one data point per month
(Burrell et al., 2020; Seddon et al., 2016). Nonetheless, the availability of
long-term weather reanalysis datasets such as ERA5 (Hersbach et al., 2020),
which comprises hourly high-resolution measurements of
$0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ ($\approx$
$27\text{\,}\mathrm{km}$) per pixel, provide the opportunity to model
dependencies of short-term changes in meteorological variables on vegetation
activity on a fine spatial and temporal resolution.

Recently, Deep Learning (DL) models have demonstrated the capability to
efficiently parse and exploit those vast amounts of meteorological data in the
context of medium-range weather forecasting. Architectural improvements and
increased compute availability have led to DL-based weather models that now
perform on par with commonly used numerical weather systems (Bi et al., 2023;
Lam et al., 2023; Pathak et al., 2022). These approaches learn a spatial
representation of the atmosphere’s state by forecasting future atmospheric
states. Previous studies have already shown that their trained atmospheric
models can be finetuned to effectively solve other climate-related tasks such
as statistical downscaling and climate projections (Lessig et al., 2023;
Nguyen et al., 2023).

Based on these advances, this work investigates how the pre-trained weather
forecasting model FourCastNet (FCN) (Pathak et al., 2022) can be adapted for
globally modeling the normalized difference vegetation index (NDVI) (Tucker &
Sellers, 1986; Vermote, 2019), a commonly used index for approximating
vegetation activity (Ferchichi et al., 2022). We outline an approach building
upon a state-of-the-art DL architecture for processing spatio-temporal data,
which enables global modeling of the NDVI at a high spatial
($0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$) and temporal (daily) resolution
with a single model. We investigate how to utilize FCN’s atmospheric knowledge
by comparing a finetuned model versus a model trained from scratch.
Additionally, we analyze the training time and data needed to make FCN an
effective vegetation model in three ablation studies.

##  2 Pre-trained weather models for vegetation modeling

### Dataset

For our study, we use daily global weather data from ERA5 (Hersbach et al.,
2020) at a resolution of $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ ($720$ x
$1440$ pixel) from the years 1982 to 2013. We use the same $20$ predictor
variables as Pathak et al. (2022): zonal and meridional wind velocity
($10\text{\,}\mathrm{m}$ above ground, at $1000\text{\,}\mathrm{hPa}$,
$850\text{\,}\mathrm{hPa}$ and $500\text{\,}\mathrm{hPa}$), temperature
($2\text{\,}\mathrm{m}$ above ground, at $850\text{\,}\mathrm{hPa}$ and
$500\text{\,}\mathrm{hPa}$), geopotential (at $1000\text{\,}\mathrm{hPa}$,
$850\text{\,}\mathrm{hPa}$, $500\text{\,}\mathrm{hPa}$, and
$50\text{\,}\mathrm{hPa}$), relative humidity (at $850\text{\,}\mathrm{hPa}$
and $500\text{\,}\mathrm{hPa}$) surface pressure, mean sea level pressure, and
total column water vapor. One sample has the dimensionality ($20$ x $720$ x
$1440$). The NDVI data (Vermote, 2019) is our target variable, regridded
linearly from originally $0.083\text{\,}\mathrm{\SIUnitSymbolDegree}$ to
ERA5’s $0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$ resolution. This vegetation
index is computed from satellite observations as the normalized difference
between the spectral reflectances in the near-infrared and red wavebands
(Tucker & Sellers, 1986). It ranges from $-1$ to $1$. Negative values indicate
water, positive values around zero indicate barren land, and values close to
one indicate dense vegetation. NDVI data after 2013 was not considered for
this study, as our analysis (cf. Fig. 5) shows a noticeable data shift
beginning in 2014

### Method

To investigate the applicability of pre-trained weather models for globally
modelling vegetation activity, we use the FCN Deep Learning model, whose pre-
trained weights are publicly available (Pathak et al., 2023). FCN is a
comparatively lightweight weather model (cf. Bi et al. (2023); Chen et al.
(2023); Lessig et al. (2023); Nguyen et al. (2023)) with a total of
$73\text{\,}\mathrm{\text{million}}$ parameters distributed over $8$
Transformer-like (Vaswani et al., 2017) encoder blocks, see Fig. 3 for an
architectural overview. Each of these blocks has
$5\text{\,}\mathrm{\text{million}}$ parameters and uses an Adaptive Fourier
Neural Operator layer (Guibas et al., 2021) replacing the attention mechanism
(Bahdanau et al., 2015).

We adopt the FCN to the NDVI modeling task by replacing the original weather
prediction head with a randomly initialized fully-connected layer with the
$\tanh$ activation function. For modeling the effects of short-term climate
variability on vegetation activity, FCN is trained on modeling the NDVI for
the same timestep as the daily input weather variables. For finetuning, we
initialize the original FCN model with the pre-trained weights, while for
training FCN from scratch, we freshly initialize all model weights.

### Comparison models

As a simple baseline, we designed and hyperparameter-optimized a convolutional
neural network (CNN) with the details given in Appendix C. We further compare
us with two recent data-driven models from ecology literature:

The first approach is a global long short-term memory (LSTM) (Hochreiter &
Schmidhuber, 1997) model by Kraft et al. (2019), trained on single-location
time-series of meteorological variables at a $15\text{\,}\mathrm{d}$ temporal
resolution and a $0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ spatial
resolution, half of our resolution, with globally shared weights. This
approach reflects the so-called memory effect of vegetation, i.e., that
preceding vegetation states can have longer-lasting effects on vegetation
activity (De Keersmaecker et al., 2015). The second data-driven approach
trains separate local, weekly state space model (SSM) on $100$ locations
across the globe, guided by equations describing the interplay between the
used climate-forcing data (Higgins et al., 2023).

### Experimental setup

The dataset is split into a training (1982-2010), validation (2011-2012), and
test (2013) set. Training details for FCN can be found in Appendix B. We
perform three individual ablation studies:. In study I, we vary the number of
finetuning epochs between $1$ to $200$ epochs. The number of frozen parameters
during finetuning is varied in study II, where we freeze between one and eight
of FCN’s Transformer blocks in ascending order. Lastly, in study III, the
number of finetuning data is modified by selecting a random
$10\text{\,}\mathrm{\char 37\relax}$ to $90\text{\,}\mathrm{\char 37\relax}$
subset of the training years.

### Evaluation

We evaluate all models by computing the root mean squared error (RMSE) and
$\text{R}^{2}$ score on the test set. The $\text{R}^{2}$ score measures the
goodness of fit of the model proportional to the temporal variation of the
target NDVI values at a given pixel and ranges from $1$ (best) to $-\infty$.
For global evaluation in comparison with the LSTM results reported by Kraft et
al. (2019), we match their evaluation scheme, thereby removing the same noisy
pixels and accounting for varying pixel areas through latitude-weighting (cf.
Appendix D). Also, our models’ outputs and target values are aggregated to
$15$-day averages to match their temporal resolution. For local evaluation in
comparison with the SSMs results provided by Higgins et al. (2023), we
evaluate our model locally at the same $100$ locations. At these locations,
FCN model outputs and target values are aggregated to $7$-day averages as in
Higgins et al. (2023), while evaluating the SSMs exclusively on the 2013 test
year.

##  3 Results and discussion

Table 1: Test year results. Left: Latitude weighted global evaluation of the finetuned FCN, FCN trained from scratch, our baseline CNN and the LSTM (Kraft et al., 2019). Right: Unweighted averages for local evaluation on $100$ locations for comparison with Higgins et al. (2023). ${\dagger}$: local model with global weight-sharing with different variables at $0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ resolution. ${\ddagger}$: local models with different variables. See Section 2 and Appendix C for details. Model | FCN finetune | FCN scratch | CNN | LSTM${\dagger}$ | FCN finetune | SSM${\ddagger}$
---|---|---|---|---|---|---
Evaluation | global, 15-daily | local, 7-daily
RMSE | $0.0403$ | $0.0512$ | $0.0431$ | $0.017$ | $0.0547$ | $0.0548$
$\text{R}^{2}$ | $0.6331$ | $0.4977$ | $0.6061$ | $0.904$ | $0.5151$ | $0.4038$

Figure 1: Results for ablation studies I & II. Left: varying number of
finetuning epochs. Right: Varying amount of training data. Results reported in
Table 1 are highlighted in both plots.  Figure 2: Results for ablation study
III: varying number of frozen Transformer blocks during finetuning. Runtimes
are averaged over five epochs.

Finetuning the learned atmospheric representation of FCN for vegetation
modeling outperforms an NDVI model trained from scratch, as Table 1 shows.
Here, the scratch model reaches an $\text{R}^{2}$ of $0.4977$ (RMSE:
$0.0512$). Finetuning the same model strongly improves NDVI modeling
performance up to an $\text{R}^{2}$ of $0.6331$, which is higher than the
strong hyperparameter-optimized CNN baseline with an $\text{R}^{2}$ of
$0.6061$ (RMSE: $0.0431$). For the FCN, a parameter search was not considered
in this study but could enhance the global ecosystem model further.

To contextualise our performance, the pixel-wise LSTM model by Kraft et al.
(2019) reaches both the highest $\text{R}^{2}$ of $0.904$ and the lowest RMSE
of $0.017$ but was trained at both lower temporal and spatial resolution than
our approach. Also, as described in Section 2, the time-series approach uses
past meteorological data which allows it to model the so-called memory effect
of vegetation. Considering this memory effect seems to be important for
modeling the NDVI from weather data and incorporating the respective
information into our model might close the observed performance gap compared
to the LSTM model.

The average performance on the $100$ locations selected by Higgins et al.
(2023) shows that our finetuned FCN reaches a $28\text{\,}\mathrm{\char
37\relax}$ higher $\text{R}^{2}$ score than the SSMs. A closer analysis of
these locations (see Appendix Table 2) shows that a single global model can
generally learn biome-specific vegetation patterns, but performance is higher
in forested regions than in regions with mainly barren land (cf. Appendix Fig.
4). Here, low- to mid-latitude ranges in the Northern Hemisphere are generally
well-modelled, while the performance in the Southern Hemisphere and high-
latitude regions is worse. This diminished performance may stem from limited
data availability towards the poles.

In ablation study I, performance most strongly rises until $80$ finetuning
epochs, as shown in Fig. 1. Afterwards, performance stagnates at around an
$\text{R}^{2}$ of $0.62$, as more finetuning epochs do not lead to further
improvements.

The results for study II in Fig. 1 show that finetuning FCN on more data
improves modeling performance. The largest jump occurs when doubling the
amount of finetuning data from $10\text{\,}\mathrm{\char 37\relax}$ to
$20\text{\,}\mathrm{\char 37\relax}$. Beyond, further increases still lead to
performance improvements albeit at a slower rate. Extrapolating Fig. 1, this
trend suggests that additional data may still enhance performance.

Freezing up to three Transformer blocks in the FCN model results in only minor
performance loss, as the results for study III in Fig. 2 show. When more
blocks are frozen, the $\text{R}^{2}$ scores can drop below the CNN baseline.
However, more frozen blocks reduce the average per-epoch runtime. It drops
from $355\text{\,}\mathrm{s}$ for the full model to $195\text{\,}\mathrm{s}$
when finetuning only the newly added vegetation-modeling head. These results
suggest that selectively freezing a moderate number of Transformer blocks can
provide a speedup compared to training the full model, while the performance
decreases marginally.

##  4 Conclusion

In this work, we investigated how a pre-trained weather model can be adapted
for globally modeling vegetation activity as measured by the normalized
difference vegetation index. We finetuned FourCastNet to model the NDVI from
$20$ meteorological variables from the ERA5 dataset and reach a globally
averaged test set $\text{R}^{2}$ of $0.6331$. This indicates that a weather
model finetuned for modeling vegetation activity from high temporal and
spatial resolution meteorological data can capture substantial amounts of the
NDVI’s variability. Our results further show that training from scratch
performs worse than finetuning a pre-trained weather model. This suggests that
during its pre-training phase, the weather model acquires structural knowledge
about the atmosphere, which is beneficial for vegetation modeling and which
probably is not attained when training from scratch.

While meteorological data partially reflects the impact of climate variability
on vegetation, other factors like atmospheric carbon dioxide, soil-related
properties and the so-called memory effect are known to be part of the complex
interplay between environmental driving forces and vegetation activity De
Keersmaecker et al. (2015); Piao et al. (2020). Hence, incorporating further
relevant variables into the model while preserving the information in the pre-
trained weather models about atmospheric dynamics is an area of future work.
Lastly, we want to highlight that explainable artificial intelligence
techniques allow examining attributions of the model’s input to its output,
such that Deep Learning models can contribute to enhancing our understanding
of how globally changing environmental factors affect local ecosystems.

## Acknowledgement

We are thankful to the HPC team at FAU for providing the computational
resources necessary for our research. This research was conducted in the
BigData@Geo 2.0 project which is co-financed by the European Regional
Development Fund (ERDF). The authors gratefully acknowledge the scientific
support and HPC resources provided by the Erlangen National High Performance
Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-
Nürnberg (FAU) under the NHR project ID b214cb. NHR funding is provided by
federal and Bavarian state authorities. NHR@FAU hardware is partially funded
by the German Research Foundation (DFG) – 440719683.

## References

  * Bahdanau et al. (2015) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio.  Neural machine translation by jointly learning to align and translate.  In Yoshua Bengio and Yann LeCun (eds.), _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * Bi et al. (2023) Kaifeng Bi, Lingxi Xie, Hengheng Zhang, Xin Chen, Xiaotao Gu, and Qi Tian.  Accurate medium-range global weather forecasting with 3d neural networks.  _Nature_ , 619(7970):533–538, 2023.
  * Burrell et al. (2020) A. L. Burrell, J. P. Evans, and M. G. De Kauwe.  Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification.  _Nature Communications_ , 11(1):3853, July 2020\.  ISSN 2041-1723.  doi: 10.1038/s41467-020-17710-7.
  * Chen et al. (2019) Chen Chen, Bin He, Wenping Yuan, Lanlan Guo, and Yafeng Zhang.  Increasing interannual variability of global vegetation greenness.  _Environmental Research Letters_ , 14(12):124005, November 2019.  ISSN 1748-9326.  doi: 10.1088/1748-9326/ab4ffc.
  * Chen et al. (2023) Kang Chen, Tao Han, Junchao Gong, Lei Bai, Fenghua Ling, Jing-Jia Luo, Xi Chen, Leiming Ma, Tianning Zhang, Rui Su, et al.  Fengwu: Pushing the skillful global medium-range weather forecast beyond 10 days lead.  _arXiv preprint arXiv:2304.02948_ , 2023.
  * Chen et al. (2021) Zhiting Chen, Hongyan Liu, Chongyang Xu, Xiuchen Wu, Boyi Liang, Jing Cao, and Deliang Chen.  Modeling vegetation greenness and its climate sensitivity with deep-learning technology.  _Ecology and Evolution_ , 11(12):7335–7345, 2021\.  ISSN 2045-7758.  doi: 10.1002/ece3.7564.
  * De Keersmaecker et al. (2015) Wanda De Keersmaecker, Stef Lhermitte, Laurent Tits, Olivier Honnay, Ben Somers, and Pol Coppin.  A model quantifying global vegetation resistance and resilience to short-term climate anomalies and their relationship with vegetation cover.  _Global Ecology and Biogeography_ , 24(5):539–548, 2015.  ISSN 1466-8238.  doi: 10.1111/geb.12279.
  * Ferchichi et al. (2022) Aya Ferchichi, Ali Ben Abbes, Vincent Barra, and Imed Riadh Farah.  Forecasting vegetation indices from spatio-temporal remotely sensed data using deep learning-based approaches: A systematic literature review.  _Ecological Informatics_ , 68:101552, May 2022.  ISSN 1574-9541.  doi: 10.1016/j.ecoinf.2022.101552.
  * Guibas et al. (2021) John Guibas, Morteza Mardani, Zongyi Li, Andrew Tao, Anima Anandkumar, and Bryan Catanzaro.  Adaptive fourier neural operators: Efficient token mixers for transformers.  _arXiv preprint arXiv:2111.13587_ , 2021.
  * Henighan et al. (2020) Tom Henighan, Jared Kaplan, Mor Katz, Mark Chen, Christopher Hesse, Jacob Jackson, Heewoo Jun, Tom B Brown, Prafulla Dhariwal, Scott Gray, et al.  Scaling laws for autoregressive generative modeling.  _arXiv preprint arXiv:2010.14701_ , 2020.
  * Hersbach et al. (2020) Hans Hersbach, Bill Bell, Paul Berrisford, Shoji Hirahara, András Horányi, Joaquín Muñoz-Sabater, Julien Nicolas, Carole Peubey, Raluca Radu, Dinand Schepers, Adrian Simmons, Cornel Soci, Saleh Abdalla, Xavier Abellan, Gianpaolo Balsamo, Peter Bechtold, Gionata Biavati, Jean Bidlot, Massimo Bonavita, Giovanna De Chiara, Per Dahlgren, Dick Dee, Michail Diamantakis, Rossana Dragani, Johannes Flemming, Richard Forbes, Manuel Fuentes, Alan Geer, Leo Haimberger, Sean Healy, Robin J. Hogan, Elías Hólm, Marta Janisková, Sarah Keeley, Patrick Laloyaux, Philippe Lopez, Cristina Lupu, Gabor Radnoti, Patricia de Rosnay, Iryna Rozum, Freja Vamborg, Sebastien Villaume, and Jean-Noël Thépaut.  The era5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.  doi: https://doi.org/10.1002/qj.3803.
  * Higgins et al. (2023) Steven I Higgins, Timo Conradi, and Edward Muhoko.  Shifts in vegetation activity of terrestrial ecosystems attributable to climate trends.  _Nature Geoscience_ , 16(2):147–153, 2023.
  * Hochreiter & Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Hoffmann et al. (2022) Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, et al.  An empirical analysis of compute-optimal large language model training.  _Advances in Neural Information Processing Systems_ , 35:30016–30030, 2022.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kraft et al. (2019) Basil Kraft, Martin Jung, Marco Körner, Christian Requena Mesa, José Cortés, and Markus Reichstein.  Identifying dynamic memory effects on vegetation state using recurrent neural networks.  _Frontiers in big Data_ , 2:31, 2019.
  * Lam et al. (2023) Remi Lam, Alvaro Sanchez-Gonzalez, Matthew Willson, Peter Wirnsberger, Meire Fortunato, Ferran Alet, Suman Ravuri, Timo Ewalds, Zach Eaton-Rosen, Weihua Hu, et al.  Learning skillful medium-range global weather forecasting.  _Science_ , 382(6677):1416–1421, 2023.
  * Lessig et al. (2023) Christian Lessig, Ilaria Luise, Bing Gong, Michael Langguth, Scarlet Stadler, and Martin Schultz.  Atmorep: A stochastic model of atmosphere dynamics using large scale representation learning.  _arXiv preprint arXiv:2308.13280_ , 2023.
  * Liu et al. (2022) Cuiyan Liu, Jianyu Liu, Qiang Zhang, Hui Ci, Xihui Gu, and Aminjon Gulakhmadov.  Attribution of NDVI Dynamics over the Globe from 1982 to 2015\.  _Remote Sensing_ , 14(11):2706, January 2022.  ISSN 2072-4292.  doi: 10.3390/rs14112706.
  * Loshchilov & Hutter (2016) Ilya Loshchilov and Frank Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  _arXiv preprint arXiv:1608.03983_ , 2016.
  * Nguyen et al. (2023) Tung Nguyen, Johannes Brandstetter, Ashish Kapoor, Jayesh K Gupta, and Aditya Grover.  Climax: A foundation model for weather and climate.  _arXiv preprint arXiv:2301.10343_ , 2023.
  * Papagiannopoulou et al. (2017) Christina Papagiannopoulou, Diego G. Miralles, Stijn Decubber, Matthias Demuzere, Niko E. C. Verhoest, Wouter A. Dorigo, and Willem Waegeman.  A non-linear Granger-causality framework to investigate climate–vegetation dynamics.  _Geoscientific Model Development_ , 10(5):1945–1960, May 2017.  ISSN 1991-959X.  doi: 10.5194/gmd-10-1945-2017.
  * Pathak et al. (2022) Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.  _arXiv preprint arXiv:2202.11214_ , 2022.
  * Pathak et al. (2023) Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.  FourCastNet: pretrained weights.  https://github.com/NVlabs/FourCastNet?tab=readme-ov-file#version-notes, 2023\.  [Online; last accessed February 2024].
  * Pelletier et al. (2015) Jon D. Pelletier, A. Brad Murray, Jennifer L. Pierce, Paul R. Bierman, David D. Breshears, Benjamin T. Crosby, Michael Ellis, Efi Foufoula-Georgiou, Arjun M. Heimsath, Chris Houser, Nick Lancaster, Marco Marani, Dorothy J. Merritts, Laura J. Moore, Joel L. Pederson, Michael J. Poulos, Tammy M. Rittenour, Joel C. Rowland, Peter Ruggiero, Dylan J. Ward, Andrew D. Wickert, and Elowyn M. Yager.  Forecasting the response of Earth’s surface to future climatic and land use changes: A review of methods and research needs.  _Earth’s Future_ , 3(7):220–251, 2015.  ISSN 2328-4277.  doi: 10.1002/2014EF000290.
  * Piao et al. (2020) Shilong Piao, Xuhui Wang, Taejin Park, Chi Chen, Xu Lian, Yue He, Jarle W. Bjerke, Anping Chen, Philippe Ciais, Hans Tømmervik, Ramakrishna R. Nemani, and Ranga B. Myneni.  Characteristics, drivers and feedbacks of global greening.  _Nature Reviews Earth & Environment_, 1(1):14–27, January 2020.  ISSN 2662-138X.  doi: 10.1038/s43017-019-0001-x.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Richardson et al. (2013) Andrew D Richardson, Trevor F Keenan, Mirco Migliavacca, Youngryel Ryu, Oliver Sonnentag, and Michael Toomey.  Climate change, phenology, and phenological control of vegetation feedbacks to the climate system.  _Agricultural and Forest Meteorology_ , 169:156–173, 2013\.
  * Robin et al. (2022) Claire Robin, Christian Requena-Mesa, Vitus Benson, Lazaro Alonso, Jeran Poehls, Nuno Carvalhais, and Markus Reichstein.  Learning to forecast vegetation greenness at fine resolution over Africa with ConvLSTMs, November 2022.
  * Seddon et al. (2016) Alistair W. R. Seddon, Marc Macias-Fauria, Peter R. Long, David Benz, and Kathy J. Willis.  Sensitivity of global terrestrial ecosystems to climate variability.  _Nature_ , 531(7593):229–232, March 2016.  ISSN 1476-4687.  doi: 10.1038/nature16986.
  * Smith et al. (2023) Michael J Smith, Luke Fleming, and James E Geach.  Earthpt: a foundation model for earth observation.  _arXiv preprint arXiv:2309.07207_ , 2023.
  * Tucker & Sellers (1986) Compton J Tucker and PJ Sellers.  Satellite remote sensing of primary production.  _International journal of remote sensing_ , 7(11):1395–1416, 1986.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Vermote (2019) Eric Vermote.  Noaa climate data record (cdr) of avhrr normalized difference vegetation index (ndvi), version 5, 2019.  Downloaded from https://www.ncei.noaa.gov/data/land-normalized-difference-vegetation-index/access/, last accessed February 2024.
  * Zhu et al. (2016) Zaichun Zhu, Shilong Piao, Ranga B. Myneni, Mengtian Huang, Zhenzhong Zeng, Josep G. Canadell, Philippe Ciais, Stephen Sitch, Pierre Friedlingstein, Almut Arneth, Chunxiang Cao, Lei Cheng, Etsushi Kato, Charles Koven, Yue Li, Xu Lian, Yongwen Liu, Ronggao Liu, Jiafu Mao, Yaozhong Pan, Shushi Peng, Josep Peñuelas, Benjamin Poulter, Thomas A. M. Pugh, Benjamin D. Stocker, Nicolas Viovy, Xuhui Wang, Yingping Wang, Zhiqiang Xiao, Hui Yang, Sönke Zaehle, and Ning Zeng.  Greening of the Earth and its drivers.  _Nature Climate Change_ , 6(8):791–795, August 2016.  ISSN 1758-6798.  doi: 10.1038/nclimate3004.

##  Appendix A Extended results and supplementary figures

Embedding Input Spatial Mixing Fourier Space Channel Mixing Linear Layer NDVIL
timesTransformer Block

Figure 3: Overview of the used architecture based on the FourCastNet model Pathak et al. (2022). We initialize the model from pre-trained weights Pathak et al. (2023) and replace the weather-specific head with a linear head for modelling the normalized difference vegetation index (NDVI). Figure 4: Global visualization of the $\text{R}^{2}$ score on the entire test set. $\text{R}^{2}$ scores below $0$ are clipped to 0 for ease of visualization. Performance is the strongest for continental Europe and North America and decreases towards higher latitude regions.  Table 2: Extended results for local evaluation. Unweighted average of RMSE and $R^{2}$ score per biome of the $100$ location from Higgins et al. (2023). The finetuned model improves the averaged $R^{2}$ score across all biomes compared to the SSMs except for boreal forests and tundra. Our global modeling approach thus captures biome-specific NDVI dynamics.  | RMSE | R2 | Samples
---|---|---|---
Biome | FCN | SSM | FCN | SSM |
Boreal forest | $0.065\,645$ | $0.083\,442$ | $0.724\,402$ | $0.832\,095$ | 16
Grassland | $0.041\,639$ | $0.044\,908$ | $0.471\,556$ | $0.448\,353$ | 14
Mediterranean-type | $0.036\,379$ | $0.039\,334$ | $0.192\,081$ | $-0.600\,848$ | 5
Tropical forest | $0.078\,943$ | $0.050\,779$ | $0.140\,531$ | $-0.051\,524$ | 16
Savanna | $0.051\,218$ | $0.051\,586$ | $0.715\,146$ | $0.662\,785$ | 18
Shrubland | $0.047\,838$ | $0.044\,490$ | $0.302\,104$ | $0.196\,631$ | 16
Temperate forest | $0.059\,774$ | $0.069\,294$ | $0.772\,599$ | $0.524\,717$ | 12
Tundra | $0.042\,269$ | $0.060\,027$ | $0.831\,170$ | $0.914\,885$ | 9
Figure 5: Distribution of the normalized difference vegetation index (NDVI)
data from 1982 to 2023. NDVI data was only used until 2013 due to the
noticeable data shift afterwards.

##  Appendix B FourCastNet training details

We train/finetune the FCN with the Adam optimizer (Kingma & Ba, 2014) and a
learning rate of $0.0001$ with cosine annealing (Loshchilov & Hutter, 2016)
for $80$ epochs (except when varying the number of training epochs in ablation
study I, see Section 2) using an $l_{2}$ loss. Using a binary mask, the loss
is only computed for locations with valid NDVI observations. The weights of
the model with the lowest validation loss are kept. We trained our models
using a single node equipped with 8 NVIDIA A100 40GB GPUs.

##  Appendix C Baseline and comparison models details

Table 3: Hyperparameter search ranges for the baseline CNN network. Hyperparameter | Search space | Step size
---|---|---
n_layers | 3 to 8 | $1$
learning rate |  $1\mathrm{e}{-5}$ to $1\mathrm{e}{-5}$ | log uniform
epochs | 50 to 100 | $10$
n_filters | 16 to 512 | $16$
kernel_size | 3 to 7 | $2$

We use a convolutional neural network (CNN) as a baseline model. The CNN’s
hyperparameters were optimized over $60$ trials on the validation set with the
search ranges given in Table 3. The hyperparameter-optimized CNN consists of
six convolution layers with $64$, $16$, $512$, $128$, $512$, $128$ kernels,
respectively. The used kernel sizes are $3$, $5$, $5$, $5$, $3$, $5$, each
with a stride of $1$. The CNN’s last layer is a fully connected dense layer
with $64$ neurons, whose output is reshaped to the target resolution of $720$
x $1440$. The CNN was trained with a learning rate of $0.00006$ for $80$
epochs.

##  Appendix D Evaluation setting

### Global evaluation

For global evaluation in comparison to LSTM models by Kraft et al. (2019), we
compute 15-days averages of our model output and target values to match their
temporal resolution. To provide a fair comparison to the reported results of
Kraft et al. (2019), we replicated their evaluation setting to the best of our
knowledge, since their code is not publicly available. To remove noisy pixels
as defined by Kraft et al. (2019), we remove pixels with
$50\text{\,}\mathrm{\char 37\relax}$ missing data in the time dimension,
pixels with less than $20\text{\,}\mathrm{\char 37\relax}$ land mass, and
barren-land pixels, which together removes coastal, high-latitude and desert
regions. Further, to account for the varying size of pixels across different
latitudes as Kraft et al. (2019), we use latitude-weighted RMSE and
$\text{R}^{2}$ scores, with the latitude weighting factors $w_{1},...,w_{I}$
for $I$ evaluated pixels given by

$w_{i}=\frac{\cos(\text{lat}(i))}{\frac{1}{N_{\text{lat}}}\sum_{j=1}^{N_{\text{lat}}}\cos(\text{lat}(j))}\quad\forall
i\in\\{1,...,I\\}\quad.$

We assume that our latitude-weighting is identical to their employed area
weighting scheme, i.e. given the latitude weights $w_{1},...,w_{I}$, the
corresponding pixel areas $\mathcal{A}_{1},...,\mathcal{A}_{I}$ and their
total area $A$

$\frac{\mathcal{A}_{i}}{A}=\frac{w_{i}}{\sum_{j=1}^{I}w_{j}}\quad\forall
i\in\\{1,...,I\\}\quad.$

With this assumption, we can show that the reported biome-weighted
$\text{RMSE}_{\text{global}}$ for the biomes $A^{1},...,A^{B}$ in Kraft et al.
(2019) is identical to our latitude-weighted RMSE:

RMSE
$\displaystyle=\frac{1}{\sum_{j=1}^{I}w_{j}}\sum_{i=1}^{I}w_{i}\cdot\text{RMSE}_{i}$
$\displaystyle=\frac{1}{A}\sum_{i=1}^{I}\mathcal{A}_{i}\cdot\text{RMSE}_{i}$
$\displaystyle=\frac{1}{A}\sum_{b=1}^{B}\sum_{k\in
I^{b}}\mathcal{A}^{b}_{k}\cdot\text{RMSE}^{b}_{k}$
$\displaystyle=\frac{1}{A}\sum_{b=1}^{B}A^{b}\cdot(\sum_{k\in
I^{b}}\frac{\mathcal{A}^{b}_{k}}{A^{b}}\cdot\text{RMSE}^{b}_{k})$
$\displaystyle=\frac{1}{A}\sum_{b=1}^{B}A^{b}\cdot\text{RMSE}^{b}$
$\displaystyle=\text{RMSE}_{\text{global}}$

The above equations also hold for the biome-weighted $R^{2}_{\text{global}}$
in Kraft et al. (2019) and our latitude-weighted $R^{2}$ score. However, we do
not perform the same 10-fold spatio-temporal cross-validation due to the long
training time of our model. Additionally, note that the LSTM model is
evaluated at a coarser resolution of
$0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$than our FCN model
($0.25\text{\,}\mathrm{\SIUnitSymbolDegree}$) and was trained on a different
set of variables. Those are six dynamic meteorological variables and $21$
static variables including water capacity, water table depth and land cover
fractions.

### Local evaluation

For local evaluation in comparison to the state space models (SSM) results
provided by Higgins et al. (2023), we average our model output and prediction
at the same $100$ locations to weekly resolution. We then compute unweighted
average RMSE and $R^{2}$ scores across these locations. Note that an
individual SSM is trained per location, using location-specific air
temperature $2\text{\,}\mathrm{m}$ above the surface, soil temperature, soil
moisture, surface solar radiation, and atmospheric carbon dioxide at
$0.083\text{\,}\mathrm{\SIUnitSymbolDegree}$ spatial resolution as climate-
forcing data. The target data is the NDVI at a weekly temporal resolution.

##  Appendix E Data scaling in ecosystem modeling

The results for ablation study II, visualized in Fig. 1, show that over the
evaluated range of number of finetuning data, performance increases as more
data is used. A similar behaviour was observed for earth observation data by
Smith et al. (2023), who trained a pixel-wise autoregressive Transformer model
(Radford et al., 2019) on satellite-derived earth observation data. They noted
that scaling their training data also scales model performance. In lieu of a
scaling “law” that applies to earth observation data, they assume their model
performance to follow the scaling law observed for large language models:

$N\sim 20\;D,$

where $N$ is the number of model parameters and $D$ is the number of training
tokens (Hoffmann et al., 2022). In the following, we assume that this tendency
also applies to our vegetation-modeling approach.

Our training dataset contains $29$ training years ($1982$ to $2010$) with
$365$ days per years, leading to $29*365=10585$ training samples. One sample
has dimensionality 720 x 1440, and the patch size used to tokenize this image-
like input is $8$. For one training sample, this leads to $720/8=90$ tokens
covering the latitudinal direction, and $1440/8=180$ tokens for the
longitudinal direction, or $90*180=16200$ total tokens for one sample. Over
the entire training data, we thus have $16200*10585=171477000$ tokens in
total.

The FourCastNet we use has $73\text{\,}\mathrm{\text{million}}$ parameters.
Assuming that the mentioned scaling law applies, we would thus need to train
FCN on $20\;*$
$73\text{\,}\mathrm{\text{million}}$$\;=\;$$1.46\text{\,}\mathrm{\text{billion}}$
tokens for optimal performance111Or, alternatively, we should scale down the
model to $171477000$ total tokens $/\;20=8.573.850$ parameters, which would
roughly be a two-Transformer block model. The observed results in Fig. 2
currently indicate that this configuration does not show the best performance,
at least for a setting where these blocks are the only ones being trained, and
further blocks – and thus parameters – exist but are frozen.  However, this
calculation assumes the validity of the large language model scaling law. For
autoregressive, image-like data-generating models, different data-scaling
might be better suited, such as the ones proposed by Henighan et al. (2020).
Applying these laws to ecosystem models is thus an open research direction.