```python
from resources.workspace import *
from IPython.display import display
from scipy.integrate import odeint
import copy

%matplotlib inline
```

# Instability, stability and chaos in dynamical systems

Chaos is commonly understood by the <b>butterfly effect</b>; "A buttefly that flaps its wings in Brazil can "cause" a hurricane in Texas". As opposed to the opinions of Descartes/Newton/Laplace, chaos effectively means that even in a deterministic (non-stochastic) universe, we can only predict "so far" into the future.  We will introduce two very typical "toy" models that exhibit these features.

## Examples of models chaotic systems

### The "Lorenz-95" model

The [Lorenz 95/ 96 system](http://eaps4.mit.edu/research/Lorenz/Predicability_a_Problem_2006.pdf) is a one dimensional model, designed to simulate atmospheric convection.  Each variable <span style='font-size:1.25em'>$x^j$ </span> can be considered some atmospheric quantity in one of $m$ sectors along a single lattitude.  The differential equation for <span style='font-size:1.25em'>$x^j$ </span> reads,
<h3>$$
\frac{{\rm d} x^j}{{\rm d} t} \triangleq -x^{j-1} x^{j-2} + x^{j-1}x^{j+1} - x^j + F,
$$</h3>
where all indices $j$ are taken modulo $m$. 

There are **no accurate physics** represented in this model.  Rather, the model only seeks to capture qualitative features of the atmosphere, in that:
<ul>
    <li> there is external forcing, determined by a parameter $F$;</li>
    <li> there is internal dissipation, simulated by linear terms;</li>
    <li> there is advection, simulated by quadratic terms.</li>
</ul>

The number of sectors $m$ is assumed to be at least $m=4$ but typically it is taken that $m=40$. See the following link for
[further description](resources/DA_intro.pdf#page=23).

**Exc 4.2**: The system above has an easy to find fixed point, i.e., a point <span style='font-size:1.25em'>$\mathbf{x}_0$</span> such that
<h3>$$ \frac{{\rm d}}{{\rm d} t} \mathbf{x}_0 \equiv 0$$ </h3>

Can you identify one?


```python
# Example solution

# show_answer('fixed_point')
```

A fixed point is **stable** if for any perturbations sufficiently small, a trajectory evolved from this perturbation must be attracted to the fixed point.  That means, all nearby solutions will settle to a solution that doesn't change in time.  A classification of various fixed point dynamics is illustrated in the figure below.

<div style='width:900px'>

</div>

**By Freesodas (Gimp) [<a href="http://www.gnu.org/copyleft/fdl.html">GFDL</a> or <a href="https://creativecommons.org/licenses/by-sa/4.0">CC BY-SA 4.0</a>], <a href="https://commons.wikimedia.org/wiki/File:Stability_Diagram.png">via Wikimedia Commons</a>**

In the Lorenz-95 model, different values for $F$ will produce different behaviors in the model.  For some values of $F$, the fixed point from **Exc 4.2** is stable.  For some values of $F$, perturbations won't be drawn to the fixed point, but will settle to another kind of **steady behavior**.  For some values of $F$, small perturbations will behave wildy, with growth and decay that is difficult to predict.

**Exc 4.4**: Run the code below to interactively plot the behavior of the Lorenz-95 sytem.  The figure on the left hand side below plots a time lapse of the values of <span style='font-size:1.25em'>$x^j$</span> on the $y$-axis, while the $x$-axis varies each sector $j$, modulo $m=10$.  The time variable is given by "T" below. 

The figure on the right hand side below plots the time series of the total engergy of the system, defined by
<h3>$$\begin{align}
\mathbf{E}(\mathbf{x}) &= \frac{1}{2} \sum_{j=1}^m \left(\mathbf{x}^j\right)^2.
\end{align}$$
</h3>
Note that for $T>30$, we only plot the time series in the interval $[T-30, T]$.

For each value of $F$, we initialize the model with a small perturbation of size "eps=0.5" to the fixed point found in **Exc 4.2**.  Answer the following questions:
<ol>
   <li> For what values of $F$ does it look like the fixed point is stable?  How is this reflected in the energy in the system?</li>
   <li> For what values of $F$ does it look like the system settles to periodic motion? How is this reflected in the energy in the system?</li>
   <li> For what values of $F$ does the evolution become "chaotic"? How is this reflected in the energy in the system?
   <li> The classical choice for $F$ is $F=8$.  What kind of behavior is exhibited?
<ol>


```python
# For all i, any n: s(x,n) := x[i+n], circularly.
def s(x,n):
    return np.roll(x,-n)

def animate_lorenz_95(F=0.8,T=0):
    # Initial conditions: perturbations
    eps=.5
    m=10
    x0 = ones(m)
    x0 = x0 * F
    x0[0] += eps
    
    def dxdt(x,t):
        return (s(x,1)-s(x,-2))*s(x,-1) - x + F
    
    tt = linspace(0, T, int(T/.1) + 1)
    xx = odeint(lambda x,t: dxdt(x,t), x0, tt)
    energy =  .5 * np.sum(xx*2, axis=1)
    xx = np.concatenate([xx,
                         np.reshape(xx[:,0], [len(xx[:,0]), 1])],axis=1)
    

    # Plot multiple
    fig = plt.figure(figsize=(16,6))
    ax1 = fig.add_axes([.08, .095,  .4, .89])
    ax2 = fig.add_axes([.525, .095, .4, .89])

    Lag = 4
    colors = plt.cm.cubehelix(0.1+0.6*linspace(0,1,Lag))
    for k in range(Lag,0,-1):
        ax1.plot(xx[max(0,len(xx)-k)],c=colors[Lag-k])

    ax1.set_ylim(-10,20)
    ax1.set_xlabel(r'Sector $j$', size=30)
    ax1.set_xticks(range(0,12,2))
    ax1.set_ylabel(r'$x^j$', size=30)
    ax1.tick_params(
        labelsize=20)
    ax2.plot(tt[-300:], energy[-300:])
    ax2.set_xlabel('Time T',size=30)
    ax2.yaxis.set_label_position("right")
    ax2.set_ylabel('Total Energy',size=30, rotation=270)
    ax2.yaxis.set_label_coords(1.175,.5)
    ax2.tick_params(
        axis='x',
        labelsize=20)
    ax2.tick_params(
        axis='y',
        labelsize=20,
        right=True,
        labelright=True,
        left=False,
        labelleft=False
    )
    tics = [np.round(i,decimals=1) for i in linspace(np.min(energy[-300:]) - 1, np.max(energy[-300:]) + 1, 5)]
    ax2.set_yticks(tics)
    ax2.set_yticklabels([str(i).zfill(2) for i in tics])
    
    plt.show()
    
interact(animate_lorenz_95,T=(0.2,60.2,0.1),F=(0,12,.2));
```

### The Lorenz (1963) system

The <b>[Lorenz-63 system](https://journals.ametsoc.org/doi/abs/10.1175/1520-0469%281963%29020%3C0130%3ADNF%3E2.0.CO%3B2)</b>, commonly known as the "butterfly attractor", is a simplified mathematical model for atmospheric convection respresenting real physics.  The Lorenz equations are derived from the Oberbeck-Boussinesq approximation to the equations describing fluid circulation in a shallow layer of fluid, heated uniformly from below and cooled uniformly from above - this describes Rayleigh-Bénard convection.

The Lorenz-63 system is given by the 3 coupled ordinary differential equations (ODE):
<h3>
$$\begin{aligned}
\dot{x} & = \sigma(y-x) \\
\dot{y} & = \rho x - y - xz \\
\dot{z} & = -\beta z + xy
\end{aligned}$$
</h3>
where 
<h3>$$\dot{\ast} \triangleq \frac{{\rm d} }{{\rm d} t}$$</h3>

As a test case for DA, the state vector is <span style='font-size:1.25em'>$\mathbf{x} = (x,y,z)$</span>, and the parameters are typically set to


```python
SIGMA = 10.0
BETA  = 8/3
RHO   = 28.0
```

The equations relate the properties of a two-dimensional fluid layer uniformly warmed from below and cooled from above. In particular, the equations describe the rate of change of three quantities with respect to time: x is proportional to the rate of convection, y to the horizontal temperature variation, and z to the vertical temperature variation.  

The dynamics can be written as follows


```python
def dxdt(xyz, t0, sigma=SIGMA, beta=BETA, rho=RHO):
    """Compute the time-derivative of the Lorenz-63 system."""
    x, y, z = xyz
    return array([
        sigma * (y - x),
        x * (rho - z) - y,
        x * y - beta * z
    ])
```

#### Numerical computation of trajectories

Below is code to numerically integrate the differential equations and plot the solutions. This function has arguments that control the parameters of the differential equation <span style='font-size:1.25em'>$(\sigma,\beta,\rho)$</span>.  In the following we will study how small perturbations of a "control" trajectory change over time.

Additional parameters in the code inlcude:
<ul>
    <li> "N", defining the number of perturbations;</li>
    <li> "eps", defining the size of perturbations;</li>
    <li> "T", defining the length of the forward evolution.  **Note**: we will only plot the trajectory along times $[T-10, T]$ for $T>10$</li>
</ul>

**Exc 4.6**: Use the code below to investigate sensititivy to initial conditions.  Answer the following questions:
<ol>
    <li> For small pertubations, eps=0.01, how long does it take to see the nearby initial conditions lose track of the control trajectory?
    <li> Does it appear that there are parameter configurations where there are **stable** fixed points? </li>
    <li> Does it appear that there are parameter configurations where there are **stable** periodic behaviors?  **Hint**: what pattern emerges with $\rho = 350$?  What about $\rho=100.5$?</li>
    <li>When all trajectories are drawn to a limiting behavior, we call this behavior **globally stable**.  Do the periodic solutions appear to be globally or locally stable?  Try multiple values of epsilon, and consider what happens when there is more than one periodic behavior.</li>
</ol>


```python
def animate_lorenz(sigma=SIGMA, beta=BETA, rho=RHO, N=2, eps=0.01, T=0.1):    
    
    # Initial conditions: perturbations around some "proto" state
    seed(1)
    x0_proto = array([-6.1, 1.2, 32.5])
    x0 = x0_proto + eps*randn((N, 3))

    # Compute trajectories
    tt = linspace(0, T, int(T/.01)+1)               # Time instances for trajectory
    d2 = lambda x,t: dxdt(x,t, sigma,beta,rho)      # Define dxdt(x,t) with fixed params.
    xx = array([odeint(d2, x0i, tt) for x0i in x0]) # Integrate
    
    
    # PLOTTING
    ax = plt.figure(figsize=(16,8)).add_subplot(111, projection='3d')
    
    colors = plt.cm.jet(linspace(0,1,N))
    for i in range(N):
        # plot each ensemble member, but only the last 1000 steps of the trajectory
        ax.plot(*(xx[i,-1000:,:].T),'-'  ,c=colors[i])
        ax.scatter3D(*xx[i,-1,:],s=40,c=colors[i])

    ax.axis('off')
    plt.show()

    
w = interactive(animate_lorenz, sigma=(0.,50), rho=(0.,350, .5), beta=(0.,5),
                N=(1,10), eps=(0.01,10.01, .1),T=(0.1,100))

w
```

In the standard configuration of the Lorenz-63 model we see that small perturbations quickly diverge from control trajectories.  But even with perfect knowledge of the state of the atmosphere, our numerical approximations of its evolution would quickly degrade the forecast skill to zero.  As a proof of concept, suppose the "<b>true</b>" atmosphere is equal to the Lorenz-63 system, evolved via the simple <a href="https://en.wikipedia.org/wiki/Euler_method" target="blank"><b>forward Euler method</b></a>, with a time step of 


```python
h_true = 0.0001
```

Suppose we are given the <b>exact</b> state of the atmosphere, defined as


```python
xyz_exact = array([-6.1, 1.2, 32.5])
```

But we must use the <a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method" target="blank"><b>Order 4.0 Runge-Kutta</b></a> scheme,


```python
def l63_rk4_step(xyz, h):
    """ calculate the evolution of Lorenz-63 one step forward via RK-4"""
    
    k_xyz_1 = dxdt(xyz, h, sigma=SIGMA, beta=BETA, rho=RHO)
    k_xyz_2 = dxdt(xyz + k_xyz_1 * (h / 2.0), h, sigma=SIGMA, beta=BETA, rho=RHO)
    k_xyz_3 = dxdt(xyz + k_xyz_2 * (h / 2.0), h, sigma=SIGMA, beta=BETA, rho=RHO)
    k_xyz_4 = dxdt(xyz + k_xyz_3 * h, h, sigma=SIGMA, beta=BETA, rho=RHO)

    xyz_step = xyz + (h / 6.0) * (k_xyz_1 + 2 * k_xyz_2 + 2 * k_xyz_3 + k_xyz_4)

    return xyz_step
```

with a time step of


```python
h_approximate = 0.1
```

to derive an <b>approximate forecast</b>, due to computational constraints.  Note that the discretization error of <b>both</b> the forward Euler (with time step h_true=0.0001) and the Order 4.0 Runge-Kutta (with time step h_approximate=0.1) is on the order $\mathcal{O}\left(10^{-4}\right)$.

**Exc 4.8**: Use the function "<b>dxdt</b>" defined above to code the forward Euler method.  Fill in the missing line in the code below. 


```python
def l63_forward_euler_step(xyz, h):
    """x_step is the one-step-forward state, derived from the initial condition xyz"""
    
    ### Fill in missing line here ###
    
    return xyz_step
```


```python
## Example solution

# show_answer('forward_euler')
```

**Exc 4.10**: Verify that your solution to Exc 4.6 works, using the GUI slider below.  The following code will generate the two paths from the **same** initial condition:
<ol>
    <li>the "true" atmosphere, generated by the forward Euler scheme and time step of 0.0001;</li>
    <li>the approximate atmospher, generated by the Runge-Kutta scheme with a time step of 0.1;</li>
</ol>
where the discretization error of each scheme is on the order of $\mathcal{O}\left(10^{-4}\right)$.



```python
def animate_approximation_divergence(T=0.1):    
    
    ## Compute trajectories
    xyz_approx_step = copy.copy(xyz_exact)
    xyz_true_step = copy.copy(xyz_exact)
    
    # define the number of integration steps
    true_steps = int(T / h_true)
    approx_steps = int(T / h_approximate)
    
    # define storage for the approximate trajectory
    xyz_approx = zeros([approx_steps + 1, 3])
    
    # define storage for the true trajectory, but where we only store the same time steps as the approximate one
    xyz_true = zeros([approx_steps + 1, 3])
    
    for i in range(approx_steps + 1):
        # store the value for each the true and approximate trajectory
        xyz_true[i, :] = xyz_true_step
        xyz_approx[i, :] = xyz_approx_step
            
        # forward propagate the approximate trajectory only at increments of 0.1
        xyz_approx_step = l63_rk4_step(xyz_approx_step, h_approximate)
        
        for j in range(1000):
            # forward propagate the true trajectory at every step of 0.0001
            xyz_true_step = l63_forward_euler_step(xyz_true_step, h_true)
        
            
            
    # PLOTTING
    ax = plt.figure(figsize=(10,5)).add_subplot(111, projection='3d')
    xx = np.dstack([xyz_true, xyz_approx])
    
    colors = plt.cm.jet(linspace(0,1,2))
    for i in range(2):
        # plot each of the trajectories, integrated with separate rules
        ax.plot(*(xx[-4:,:,i].T),'-'  ,c=colors[i])
        ax.scatter3D(*xx[-1,:,i], s=40, c=colors[i])

    ax.set_xlim((-15, 15))
    ax.set_ylim((-25, 25))
    ax.set_zlim((15, 40))
    ax.view_init(30, 120)
    ax.axis('off')
    plt.show()
    

w = interactive(animate_approximation_divergence, T=(0.4,4))
w
```

**Exc 4.12**: What do you notice about this plot?  What does this say about small simulation errors in chaotic systems?

### Next: [Dynamics of ensembles and perturbations](T2 - Dynamics of ensembles and perturbations.ipynb)

bfe9cbd477a0226fabed707caafb531269bc57d5
ipynb
Jupyter Notebook
tutorials/DA and the Dynamics of Ensemble Based Forecasting/T1 - Instability, stability and chaos in dynamical systems.ipynb
brajard/DAPPER
1a513b2f23041b15fb335aeb17906607bf2a5350
MIT
2021-07-31T10:13:11.000Z
2022-01-14T16:52:04.000Z
tutorials/DA and the Dynamics of Ensemble Based Forecasting/T1 - Instability, stability and chaos in dynamical systems.ipynb
franktoffel/dapper
373a27273ea109f349e5edcdcef0cfe0b83b925e
MIT
tutorials/DA and the Dynamics of Ensemble Based Forecasting/T1 - Instability, stability and chaos in dynamical systems.ipynb
franktoffel/dapper
373a27273ea109f349e5edcdcef0cfe0b83b925e
MIT
2020-01-25T16:35:00.000Z
2021-04-08T03:20:48.000Z
Qwen/Qwen-72B
1. YES
2. YES