<jupyter_start><jupyter_text>Dynamics of ensembles and perturbationsIn the following, we will neglect the important issues of model deficiencies and the inherent stochasticity of certain dynamics.  Instead, we will focus on the simplified scenario where we assume:     we can perfectly model and compute the purely deterministic dynamics; and     prediction error originates soley from the uncertainty in initial conditions.Understanding that perturbations rapidly diverge even in a chaotic system as described above, this led to the transition from single-trajectory forecasts to ensemble-based forecasts.   Ensembles are used to "average out" our initialization errors, and to understand the variability and uncertainty in forecasts. The advantages of ensemble based forecasts over single trajectory forecasts historically led to a search for perturbations that are most representative of the error growth in operational forecasts.  Prediction centers have sought to initialize ensemble-based forecasts in a way to best capture the variability induced by the dynamical chaos.  Two major techniques emerged,    "bred vectors", and     forcing singular vectors".    These lead to different formulations of the classical "Lyapunov vectors".  We do not stress here what a "Lyapunov vector" is, rather we will discover their nature experimentally in the following work.  This will lead to a formal definition of one type of Lyapunov vectors by the end of the exercises.  "Breeding" growing modesSuppose we have a smooth, nonlinear dynamical system,$$\begin{align}\dot{\mathbf{x}} = f(\mathbf{x}) & & \mathbf{x} \in \mathbb{R}^n,\end{align}$$and a precise estimate of an initial condition $\mathbf{x}^c_0$, from which we want to make a forecast.  Suppose, also, that there are future observations that we will assimilate. It was suggested by Toth and Kalnay to use the evolution of the initial estimate $\mathbf{x}^c_0$ as a control trajectory, while introducing small perturbations to generate an ensemble,$$\begin{align}\mathbf{x}^i_0 = \mathbf{x}^c_0 + \boldsymbol{\delta}^i_0 \end{align}$$where $\left \rvert \boldsymbol{\delta}^i \right \rvert  = \epsilon \ll 1$.The ensemble is evolved in parallel to the control trajectory.  Between times $t_{k-1}$ and $t_k$, this takes the form$$\begin{align}\widehat{\mathbf{x}}^c_k &= \mathbf{x}_{k-1}^c + \int_{t_{k-1}}^{t_k} f(x) {\rm d}t \\\widehat{\mathbf{x}}^i_k &= \mathbf{x}_{k-1}^i + \int_{t_{k-1}}^{t_k} f(x) {\rm d}t.\end{align}$$ At the point of analyzing new observations we form a new estimate for the control trajectory, taking $\widehat{\mathbf{x}}_k^c$ to $\mathbf{x}_k^c$. the perturbations are rescaled back to their original small size while maintaining their directions.  That is to say,$$\begin{align}\widehat{\boldsymbol{\delta}}_k^i \triangleq \mathbf{x}_k^c - \widehat{\mathbf{x}}_k^i, & &\boldsymbol{\delta}_k^i \triangleq \frac{\epsilon}{\left\rvert \widehat{\boldsymbol{\delta}}^i_k\right\rvert} \widehat{\boldsymbol{\delta}}^i_k.\end{align}$$ "Breeding growing modes" is designed to simulate how the modes of fast growing error are maintained and propagated through the successive use of short range forecasts.  The resulting perturbations are thus meant to represent a perturbation field of the "errors of the day", i.e., uncertainties in the initial condition at the present time that result from the repeated cycle of forecasts and analyses. **Exc 4.14**: Run the code below and use the sliders to examine behavior of successive "breeding" of growing modes. The parameter B stands for the number of breeding cycles.  The parameter eps stands for the re-scaling parameter $\epsilon$ defined above.  The parameter N is the number of perturbations.  The the plots on the left hand side show the evolution of the control trajectory and the perturbed trajectories along each breeding cycle.  The right hand side plots $\pm 1$ times the normalized perturbations,$$ \frac{ \pm\widehat{\boldsymbol{\delta}}_k^i}{\left\rvert \widehat{\boldsymbol{\delta}}^i_k\right\rvert},$$giving the directions of the perturbations, plotted as lines through the unit sphere.  **Answer the following questions**:     For small values of $\epsilon$, what is significant about the long term behavior of the directions of the perturbations?      Does this behavior change with large N, i.e., more directions for the pertubations?      How does this behavior change when $\epsilon$ is increased?  Do the directions of the perturbations depend on N for large $\epsilon$?<jupyter_code>SIGMA = 10.0
BETA  = 8/3
RHO   = 28.0

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
    ])

def animate_bred_vectors(B=0, eps=0.01, N=10):    
    
    # Initial conditions: perturbations around some "proto" state
    sigma=SIGMA 
    beta=BETA 
    rho=RHO
    T=0.05
    
    seed(1)
    x_0 = array([-6.1, 1.2, 32.5])               # define the control
    
    # define the perturbations, randomly generated but of fixed norm epsilon
    perts = randn([N,3])
    perts = array([eps * perts[i] / sqrt(perts[i] @ perts[i]) for i in range(N)])
    delta_x = x_0 + perts
                  
    tt = linspace(0, T, 10)           # Time instances for trajectory
    d2 = lambda x,t: dxdt(x,t, sigma,beta,rho)  # Define dxdt(x,t) with fixed params.    
    
    # for each breeding cycle
    for kk in range(B):
        # Compute trajectories
        x_traj = array([odeint(d2, x_0, tt)])        # integrate the control trajectory
        x_0 = np.squeeze(x_traj[:, -1, :])
        
        delta_x_traj = array([odeint(d2, delta_xi, tt) for delta_xi in delta_x]) # Integrate the perturbations
        perts = delta_x_traj[:, -1, :] - x_0                                     # redefine the perts
        perts = array([eps * perts[i] / sqrt(perts[i] @ perts[i]) for i in range(N)])
        delta_x = x_0 + perts # redefine the initialization of the perturbed trajectories
    
    # PLOTTING
    fig = plt.figure(figsize=(16,8))
    ax1 = plt.subplot(121, projection='3d')
    ax2 = plt.subplot(122, projection='3d')
    
    if B==0:
        ax1.scatter3D(*x_0, s=40, c='k')

    else:
        ax1.plot(*x_traj[0,:,:].T, '-', c='k')
        ax1.scatter3D(*x_traj[0,-1,:].T, '-', s=40, c='k')
            
    colors = plt.cm.jet(linspace(0,1,N))
    for i in range(N):
        # for each breeding cycle
        if B==0:
            # if just the initial conditions, we plot these
            ax1.scatter3D(*delta_x[i,:],s=40,c=colors[i])
            
        else:
            # otherwise, plot the trajectories over a breeding cycle
            ax1.plot(*(delta_x_traj[i,:,:].T),'-'  ,c=colors[i])
            ax1.scatter3D(*delta_x_traj[i,-1,:],s=40,c=colors[i])
                          
        # we plot the normalized perturbations on the unit sphere
        tmp = perts[i,:]/sqrt(perts[i,:] @ perts[i, :])
        p_vect = np.concatenate([np.reshape([0,0,0],[1,3]), np.reshape(tmp,[1,3])], axis=0)
        
        # delta * +1
        ax2.plot(p_vect[:,0], p_vect[:,1], p_vect[:,2],'-'  ,c=colors[i])
        ax2.scatter3D(*tmp[:],s=40,c=colors[i], marker='o')

        # delta * -1
        ax2.plot(-1*p_vect[:,0], -1*p_vect[:,1], -1*p_vect[:,2],'-'  ,c=colors[i])
        ax2.scatter3D(*tmp[:]*(-1),s=40,c=colors[i], marker='o')

    ax1.axis('off')
    ax2.set_xlim((-.9, .9))
    ax2.set_ylim((-.9, .9))
    ax2.set_zlim((-.9, .9))
    plt.show()
    
w = interactive(animate_bred_vectors,B=(0,175,25), eps=(0.01,1.61, .2), N=(1, 161, 20))
w<jupyter_output><empty_output><jupyter_text>**Exc 4.16**: If the "breeding" of perturbations is meant to represent the unstable growth of initial perturbations, what can we learn from their growth rates?  In the following code, fill in the missing lines to define a function that will compute the log-growth rate of the perturbation $\boldsymbol{\delta^i}_k$, i.e., the log growth relative to the length of time in the breeding interval.  This function should return$$\begin{align}\frac{1}{T}\log \left( \frac{\left\rvert \widehat{\boldsymbol{\delta}}^i_k\right \rvert}{\left\rvert \boldsymbol{\delta}_{k-1}^i \right\rvert}\right) \equiv \frac{1}{T}\log \left( \frac{\left\rvert \widehat{\boldsymbol{\delta}}^i_k\right \rvert}{\epsilon}\right),\end{align}$$for a single perturbation.<jupyter_code>def log_growth(x_control_k, x_pert_k, T, eps):
    """function returns array of the log growth values for a single perturbation"""
    nrm = sqrt( (x_pert_k - x_control_k) @ (x_pert_k - x_control_k).T )

    log_growth_rate = (1.0 / T) * log(nrm / eps)
    
    ### Fill in missing line(s) here ###

    
    return log_growth_rate
## Example solution

# show_answer('log_growth')<jupyter_output><empty_output><jupyter_text>**Exc 4.18**: Test your answer to **Exc 4.16**.  Using the code and slider below, investigate the distributions of the log-growth rates of the bred vectors as a function of the number of breeding cycles.  Answer the following questions:     What is the long term behaviour of this distribution?     The leading Lyapunov exponent of the Lorenz-63 system is $\approx 0.9050$, what do you notice about the mean of the log-growth rates over a long number of breeding cycles?     Consider the behavior of small perturbations from Exc 4.14. Can you conjecture what the perturbations are converging to?     What does this suggest about the "representative perturbations" for the error growth in chaotic systems?<jupyter_code>def animate_bred_growth_rates(B=1000):    
    
    # Initial conditions: perturbations around some "proto" state
    sigma=SIGMA 
    beta=BETA 
    rho=RHO
    eps=0.01
    N=1
    T=0.01
    
    seed(1)
    x_0 = array([-6.1, 1.2, 32.5])               # define the control
    
    # define the perturbation, randomly generated but of fixed norm epsilon
    perts = randn(3)
    perts = eps * perts / sqrt(perts @ perts)
    delta_x = x_0 + perts
                  
    tt = linspace(0, T, 20)           # Time instances for trajectory
    d2 = lambda x,t: dxdt(x,t, sigma,beta,rho)  # Define dxdt(x,t) with fixed params.    
    grwt = np.zeros(B)
    
    
    # for each breeding cycle
    for kk in range(B):
        # Compute trajectories
        x_traj = array([odeint(d2, x_0, tt)])        # integrate the control trajectory
        x_0 = np.squeeze(x_traj[:, -1, :])
        
        delta_x_traj = array([odeint(d2, delta_x, tt)]) # Integrate the perturbation
        
        # compute the log growth from the code defined earlier
        grwt[kk] = log_growth(x_0, delta_x_traj[0, -1, :], T, eps)
        
        # redefine perts
        perts = delta_x_traj[0, -1, :] - x_0
        perts = eps * perts / sqrt(perts @ perts.T) 
        delta_x = x_0 + perts
    
    # PLOTTING
    fig = plt.figure(figsize=(16,8))
    ax = plt.subplot(111)
    ax.hist(grwt, bins=linspace(-20,20,4001), density=True)
    ax.set_xlim((-12, 12))
    ax.set_ylim((0, 0.8))
    ax.text(4, 0.6, 'Mean log-growth rate=' + str(np.round(mean(grwt),decimals=4)).zfill(4), size=20)
    
    plt.show()
    
w = interactive(animate_bred_growth_rates,B=(1000,30000, 5800))
w<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11973078285511866, 'fraction_numerical': 0.023025150549061284, 'mean_word_length': 3.2711800302571863, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '81257', 'license': 'MIT', 'n_tokens_mistral': 3752, 'n_tokens_neox': 3488, 'n_words': 1204, 'path': 'tutorials/DA and the Dynamics of Ensemble Based Forecasting/T2 - Dynamics of ensembles and perturbations.ipynb', 'repo_name': 'brajard/DAPPER', 'size': '11292'}