Zoe Farmer

import numpy as np
import matplotlib.pyplot as plt

%pylab inline

Populating the interactive namespace from numpy and matplotlib

1. Derive the Jacobian $D_{\vec{x} }\vec{F}$ for the Lorenz system:

\[\vec{F}(\vec{x}, a, r, b) = \left[\begin{array}{c} \dot{x}\\ \dot{y}\\ \dot{z} \end{array}\right] = \left[\begin{array}{c} a(y - x)\\ rx - y - xz\\ xy - bz \end{array}\right]\]

The Jacobian of this system is defined as

\[\left[\begin{array}{ccc} \partial_x \dot{x} & \partial_y \dot{x} & \partial_z \dot{x}\\ \partial_x \dot{y} & \partial_y \dot{y} & \partial_z \dot{y}\\ \partial_x \dot{z} & \partial_y \dot{z} & \partial_z \dot{z}\\ \end{array}\right]\]

Which is equal to

\[\left[\begin{array}{ccc} -a & a & 0\\ r - z & -1 & -x\\ y & x & -b\\ \end{array}\right]\]

2. Write down the associated variational system $\dot{\delta} = D_{\vec{x} } \vec{F} \delta$. The product of the Jacobian matrix $D_{\vec{x} } \vec{F}$ and the $n \times n$ matrix of variations

\[\delta = \left[\begin{array}{ccc} \delta_{xx} & \delta_{yx} & \delta_{zx}\\ \delta_{xy} & \delta_{yy} & \delta_{zy}\\ \delta_{xz} & \delta_{yz} & \delta_{zz}\\ \end{array}\right]\]

yields an $n \times n$ matrix of the derivatives of the variation $\dot{\delta}$.

Put more simply this is

\[\left[\begin{array}{ccc} -a & a & 0\\ r - z & -1 & -x\\ y & x & -b\\ \end{array}\right] \left[\begin{array}{ccc} \delta_{xx} & \delta_{yx} & \delta_{zx}\\ \delta_{xy} & \delta_{yy} & \delta_{zy}\\ \delta_{xz} & \delta_{yz} & \delta_{zz}\\ \end{array}\right]\]

Multiplying it out we get

\[\left[ \begin{array}{ccc} a \delta _{\text{xy} }-a \delta _{\text{xx} } & a \delta _{\text{yy} }-a \delta _{\text{yx} } & a \delta _{\text{zy} }-a \delta _{\text{zx} } \\ (r-z) \delta _{\text{xx} }-\delta _{\text{xy} }-x \delta _{\text{xz} } & (r-z) \delta _{\text{yx} }-\delta _{\text{yy} }-x \delta _{\text{yz} } & (r-z) \delta _{\text{zx} }-\delta _{\text{zy} }-x \delta _{\text{zz} } \\ y \delta _{\text{xx} }+x \delta _{\text{xy} }-b \delta _{\text{xz} } & y \delta _{\text{yx} }+x \delta _{\text{yy} }-b \delta _{\text{yz} } & y \delta _{\text{zx} }+x \delta _{\text{zy} }-b \delta _{\text{zz} } \\ \end{array} \right]\]

3. A combination of this variational derivative and the original system derivative can be used to integrate the $(n^2 + n)$-dimensional variational equation

\[\left\{ \begin{array}{c} \dot{\vec{x} }\\ \dot{\delta}\\ \end{array} \right\} = \left\{ \begin{array}{c} \vec{F}\\ D_{\vec{x} } \vec{F} \delta \end{array} \right\}\]

from the initial condition

\[\left\{\begin{array}{c} \vec{x_0}\\ I \end{array}\right\}\]

with $t = t_0$. The time evolution of the first $n$ elements of this set of initial conditions follows the trajectory $\Phi_t (\vec{x_0}, t_0)$. The row and column sums of the matrix formed by the next $n^2$ elements are different ways to look at the evolved versions of the inital variations: the first column sum, for example, gives the $x$ component of the evolved variation, while the first row sum tells you what the $x$-piece of the original variation has grown into.

Integrate the Lorenz variational equation using RK4 (NOT ARK4) from the following initial conditions for 100 steps. Use $a=16$, $r = 45$, $b = 4$, and a timestep of $0.001$. In each case give the components of the evolved matrix $\delta$ and the evolved variations (the columns sums of $\delta$ at the endpoint of the trajectory. Use $t_0 = 0$. You need only turn in these twelve numbers for each question.

Let’s first get our RK4 solver in here.

def _step_rk4(f, t0, y0, dt):
    """
    Perform one Fourth-Order Runge-Kutta Step
    
    :param f: function
                ODE function with input f(t, y)
    :param t0: float
                Initial starting time
    :param y0: numpy.ndarray()
                Initial ODE state (vector)
    :param dt: float
                Timestep
    """
    k1 = f(t0, y0)
    k2 = f(t0 + (dt / 2), y0 + (k1 * dt / 2))
    k3 = f(t0 + (dt / 2), y0 + (k2 * dt / 2))
    k4 = f(t0 + dt, y0 + (k3 * dt))
    y = y0 + dt * ((k1 / 6) + (k2 / 3) +
                   (k3 / 3) + (k4 / 6))
    return (t0 + dt,
            (y0 + dt * ((k1 / 6) + (k2 / 3) +
                (k3 / 3) + (k4 / 6))))



def mrk4(f, t0, y0, dt, n, writecsv=''):
    """
    Fixed-Step Fourth-Order Runge-Kutta ODE Solver
    
    :param f: function
        ODE function with input f(t, y)
    :param t0: float
        Initial starting time
    :param y0: numpy.ndarray()
        Initial ODE state (vector)
    :param dt: float
        Timestep
    :param n: int
        Number of iterations (steps) to perform
    :param writecsv: bool
        :default: False
        Write to csv file?
    """
    dim = y0.size

    # Establish blank solution trajectory
    # [[y00, ..., y0n, t0],
    #  [y10, ..., y1n, t1],
    # ...]
    traj = np.zeros((n + 1, dim + 1), dtype=np.float64)
    
    # Set initial position
    traj[0, 0:dim] = y0
    traj[0, -1]    = t0
    
    # Iterate
    for i in range(1, n + 1):
        (traj[i, -1],
         traj[i, 0:dim]) = _step_rk4(f,
                                    traj[i - 1, -1],
                                    traj[i - 1, 0:dim],
                                    dt)

    if writecsv != '':
        with open(writecsv, 'w') as f:
            csvwriter = csv.writer(f)
            [csvwriter.writerow(line) for line in m]
    return traj

(a)

\[\left[ \begin{array}{c} x\\ y\\ z\\ \delta_{xx}\\ \delta_{xy}\\ \delta_{xz}\\ \delta_{yx}\\ \delta_{yy}\\ \delta_{yz}\\ \delta_{zx}\\ \delta_{zy}\\ \delta_{zz}\\ \end{array} \right] = \left[ \begin{array}{c} 0\\ 1\\ 2\\ 1\\ 0\\ 0\\ 0\\ 1\\ 0\\ 0\\ 0\\ 1\\ \end{array} \right]\]

First we generate our lorenz system.

def gen_lorenz(a, r, b):
    lorenz = lambda t, y: np.array([a * (y[1] - y[0]),
                                    r * y[0] - y[1] - y[0] * y[2],
                                    y[0] * y[1] - b * y[2],
                                    -a * y[3] + a * y[6],
                                    -a * y[4] + a * y[7],
                                    -a * y[5] + a * y[8],
                                    (r - y[2]) * y[3] - y[6] - y[0] * y[9],
                                    (r - y[2]) * y[4] - y[7] - y[0] * y[10],
                                    (r - y[2]) * y[5] - y[8] - y[0] * y[11],
                                    y[1] * y[3] + y[0] * y[6] - b * y[9],
                                    y[1] * y[4] + y[0] * y[7] - b * y[10],
                                    y[1] * y[5] + y[0] * y[8] - b * y[11]])
    return lorenz

Now we can run it through our solver since we’ve established our variational equation.

initial_state = np.array([0, 1, 2, 1, 0, 0, 0, 1, 0, 0, 0, 1],
                         dtype=np.float64)

lorenz = gen_lorenz(16, 45, 4)

final = mrk4(lorenz, 0, initial_state, 0.001, 100)[-1]

np.array([[final[3:6], final[6:9], final[9:12]]])

array([[[ 2.40394187,  1.91887689, -0.02967512],
        [ 5.17512111,  4.20716588, -0.08930062],
        [ 0.48078252,  0.37413678,  0.66511534]]])

(b)

\[\left[ \begin{array}{c} x\\ y\\ z\\ \delta_{xx}\\ \delta_{xy}\\ \delta_{xz}\\ \delta_{yx}\\ \delta_{yy}\\ \delta_{yz}\\ \delta_{zx}\\ \delta_{zy}\\ \delta_{zz}\\ \end{array} \right] = \left[ \begin{array}{c} 10\\ -5\\ 2\\ 1\\ 0\\ 0\\ 0\\ 1\\ 0\\ 0\\ 0\\ 1\\ \end{array} \right]\]

initial_state = np.array([10, -5, 2, 1, 0, 0, 0, 1, 0, 0, 0, 1],
                         dtype=np.float64)

lorenz = gen_lorenz(16, 45, 4)

final = mrk4(lorenz, 0, initial_state, 0.001, 100)[-1]

np.array([[final[3:6], final[6:9], final[9:12]]])

array([[[ 2.13103438,  1.63916327, -0.48580854],
        [ 3.86002906,  3.00832887, -1.11058482],
        [ 3.12194096,  2.53248319, -0.01806069]]])

(c)

\[\left[ \begin{array}{c} x\\ y\\ z\\ \delta_{xx}\\ \delta_{xy}\\ \delta_{xz}\\ \delta_{yx}\\ \delta_{yy}\\ \delta_{yz}\\ \delta_{zx}\\ \delta_{zy}\\ \delta_{zz}\\ \end{array} \right] = \left[ \begin{array}{c} 0\\ -1\\ 2\\ 1\\ 0\\ 0\\ 0\\ 1\\ 0\\ 0\\ 0\\ 1\\ \end{array} \right]\]

initial_state = np.array([0, -1, 2, 1, 0, 0, 0, 1, 0, 0, 0, 1],
                         dtype=np.float64)

lorenz = gen_lorenz(16, 45, 4)

final = mrk4(lorenz, 0, initial_state, 0.001, 100)[-1]

np.array([[final[3:6], final[6:9], final[9:12]]])

array([[[ 2.40394187,  1.91887689,  0.02967512],
        [ 5.17512111,  4.20716588,  0.08930062],
        [-0.48078252, -0.37413678,  0.66511534]]])

(d) Look carefully at the evolved matrices of variations and describe some of their interesting features. From which point (a, b, or c) do the variations grow fastest? In which direction? Do you notice any symmetries of gross differences between the different points?

These three matrices are listed.

\[\begin{cases} (a) = \left[\begin{array}{ccc} 2.40394187&1.91887689&-0.02967512\\ 5.17512111&4.20716588&-0.08930062\\ 0.48078252&0.37413678&0.66511534 \end{array}\right]\\ (b) = \left[\begin{array}{ccc} 2.13103438 &1.63916327&-0.48580854\\ 3.86002906 &3.00832887&-1.11058482\\ 3.12194096 &2.53248319& -0.01806069 \end{array}\right]\\ (c) = \left[\begin{array}{ccc} 2.40394187& 1.91887689& 0.02967512\\ 5.17512111& 4.20716588& 0.08930062\\ -0.48078252& -0.37413678& 0.66511534 \end{array}\right] \end{cases}\]

We’ll note that matrices $(a)$ and $(c)$ are very similar, which makes a lot of sense as the initial states were almost identical.

We can also examine the column sums and obtain the following results:

\[\begin{cases} a \approx \{8, 6.5, -0.6\}\\ b \approx \{9, 7.1, -1.5\}\\ c \approx \{7, 5.7, 0.76\} \end{cases}\]

Using this result we see that the set of $(b)$ points have the fastest growing variations, and this growth is in the $x$ direction.