The variational equations

In this tutorial, we will show how it is possible to exploit heyoka.py’s expression system to implement first-order variational equations for a system of ODEs.

First, let’s start with a brief recap on the variational equations. For simplicity, we will consider an autonomous ODE system in the variables \(x\) and \(y\) (the extension to more variables just complicates the notation):

\[\begin{split} \begin{cases} \frac{dx}{dt} = f\left(x, y \right)\\ \frac{dy}{dt} = g\left(x, y \right) \end{cases}, \end{split}\]

with initial conditions

\[\begin{split} \begin{cases} x\left( 0 \right) = x_0\\ y\left( 0 \right) = y_0 \end{cases}. \end{split}\]

The solutions of this ODE system will be the functions \(x\left( t; x_0, y_0 \right)\) and \(y\left( t; x_0, y_0 \right)\). Our objective is to determine the evolution in time of the first-order derivatives of these solutions with respect to the initial conditions \(\left( x_0, y_0 \right)\):

\[\begin{split} \begin{aligned} \left(\frac{\partial x}{\partial x_0}\right)\left(t; x_0, y_0\right),\\ \left(\frac{\partial x}{\partial y_0}\right)\left(t; x_0, y_0\right),\\ \left(\frac{\partial y}{\partial x_0}\right)\left(t; x_0, y_0\right),\\ \left(\frac{\partial y}{\partial y_0}\right)\left(t; x_0, y_0\right).\\ \end{aligned} \end{split}\]

These functions tell us how the the solutions of the ODE system react to changes in the initial values \(\left( x_0, y_0 \right)\). One way of determining numerically the evolution in time of the partial derivatives is to introduce them as additional state variables in the original ODE system, together with an additional set of differential equations (i.e., the variational equations). Let’s see how this is done for, e.g., \(\frac{\partial x}{\partial x_0}\). We take the time derivative of \(\frac{\partial x}{\partial x_0}\),

\[ \frac{d \left( \frac{\partial x}{\partial x_0} \right)}{dt}, \]

and expand it, via elementary calculus rules, to

\[ \frac{d \left( \frac{\partial x}{\partial x_0} \right)}{dt} = \frac{\partial \left( \frac{dx}{dt} \right) }{\partial x_0}=\frac{\partial f\left(x, y\right)}{\partial x}\frac{\partial x}{\partial x_0}+\frac{\partial f \left(x, y\right)}{\partial y}\frac{\partial y}{\partial x_0}. \]

We can write similar ODEs for the other partial derivatives, append them below the original ODE system and finally obtain an augmented ODE system containing the variational equations:

\[\begin{split} \begin{cases} \frac{dx}{dt} = f\left(x, y \right)\\ \frac{dy}{dt} = g\left(x, y \right)\\ \frac{d \left( \frac{\partial x}{\partial x_0} \right)}{dt} = \frac{\partial f\left(x, y\right)}{\partial x}\frac{\partial x}{\partial x_0}+\frac{\partial f\left(x, y\right)}{\partial y}\frac{\partial y}{\partial x_0}\\ \frac{d \left( \frac{\partial x}{\partial y_0} \right)}{dt} = \frac{\partial f\left(x, y\right)}{\partial x}\frac{\partial x}{\partial y_0}+\frac{\partial f\left(x, y\right)}{\partial y}\frac{\partial y}{\partial y_0}\\ \frac{d \left( \frac{\partial y}{\partial x_0} \right)}{dt} = \frac{\partial g\left(x, y\right)}{\partial x}\frac{\partial x}{\partial x_0}+\frac{\partial g\left(x, y\right)}{\partial y}\frac{\partial y}{\partial x_0}\\ \frac{d \left( \frac{\partial y}{\partial y_0} \right)}{dt} = \frac{\partial g\left(x, y\right)}{\partial x}\frac{\partial x}{\partial y_0}+\frac{\partial g\left(x, y\right)}{\partial y}\frac{\partial y}{\partial y_0} \end{cases}. \end{split}\]

The introduction of the variational equations in an ODE system of \(n\) equations results in additional \(n^2\) equations and variables. The initial conditions for the new variables are:

\[\begin{split} \begin{cases} \frac{\partial x}{\partial x_0} \left( 0 \right) = 1\\ \frac{\partial x}{\partial y_0} \left( 0 \right) = 0\\ \frac{\partial y}{\partial x_0} \left( 0 \right) = 0\\ \frac{\partial y}{\partial y_0} \left( 0 \right) = 1\\ \end{cases}. \end{split}\]

Let us now try to implement the variational equations for the simple pendulum:

import heyoka as hy

# Create the state variables x and v.
x, v = hy.make_vars("x", "v")

# Create the symbolic variables for the variational equations.
x_x0, x_y0, y_x0, y_y0 = hy.make_vars("x_x0", "x_y0", "y_x0", "y_y0")

# Create the right-hand side of the ODE system.
f = v
g = -9.8 * hy.sin(x)
x_x0_p = hy.diff(f, x) * x_x0 + hy.diff(f, v) * y_x0
x_y0_p = hy.diff(f, x) * x_y0 + hy.diff(f, v) * y_y0
y_x0_p = hy.diff(g, x) * x_x0 + hy.diff(g, v) * y_x0
y_y0_p = hy.diff(g, x) * x_y0 + hy.diff(g, v) * y_y0

# Create the integrator object.
ta = hy.taylor_adaptive(
                        # Definition of the ODE system.
                        sys = [
                            (x, f),
                            (v, g),
                            (x_x0, x_x0_p),
                            (x_y0, x_y0_p),
                            (y_x0, y_x0_p),
                            (y_y0, y_y0_p),
                        ],
                        # Initial conditions.
                        state = [0.05, 0.025, 1., 0., 0., 1.])

ta

C++ datatype            : double
Tolerance               : 2.220446049250313e-16
High accuracy           : false
Compact mode            : false
Taylor order            : 20
Dimension               : 6
Time                    : 0
State                   : [0.05, 0.025, 1, 0, 0, 1]

Here we used the diff() function to perform the symbolic differentiation of f and g with respect to the state variables.

Let us now propagate the state of the system (including the variational equations) up to \(t=50\):

ta.propagate_until(50.)
ta.state

array([ 0.03744787,  0.10667026,  0.80315068, -0.17679098,  1.82916222,
        0.84245788])

Next, we are going to create another integrator for the simple pendulum, this time without the variational equations and perturbing the initial value of \(x\) by \(10^{-8}\) with respect to the original initial conditions:

ta_dx = hy.taylor_adaptive(
                        # Definition of the ODE system.
                        sys = [
                            (x, f),
                            (v, g)
                        ],
                        # Initial conditions.
                        state = [0.05 + 1e-8, 0.025])

ta_dx.propagate_until(50.)

(<taylor_outcome.time_limit: -4294967299>,
 0.20208180443700685,
 0.2182277296610072,
 240,
 None,
 None)

Because of the perturbation on the initial value of \(x\), the value of \(x\) at \(t=50\) will differ in the two integrators:

ta.state[0] - ta_dx.state[0]

-8.031506786021492e-09

However, because in the variational integrator we now have the value of \(\frac{\partial x}{\partial x_0}\) at \(t=50\), we can compute an approximation of the \(\Delta x\) induced by the perturbation on the initial state as

\[ \Delta x \approx \frac{\partial x}{\partial x_0}\left( 50 \right) \cdot 10^{-8}. \]

Indeed:

(ta.state[0] + ta.state[2]*1e-8) - ta_dx.state[0]

2.0816681711721685e-17