Other schemes can be used, instead.
Bi = np.linalg.inv(B)
Ri = np.linalg.inv(R)
ua = np.copy(ub)
for iter in range(100):
Dh = JObsOp(ua)
A = Bi + (Dh.T)@Ri@Dh
b = Bi@(ub-ua) + (Dh.T)@Ri@(w-ObsOp(ua))
du = np.linalg.solve(A,b) #solve a linear system
ua = ua + du
if np.linalg.norm(du) <= 1e-4:
break
return ua
```
#### *3.3. Example: Lorenz 63 System*

The Lorenz 63 equations have been utilized as a toy problem in data assimilation studies, capturing some of the interesting mechanisms of weather systems. The three-equation model can be written as

$$\begin{aligned} \frac{\mathbf{dx}}{\mathbf{dt}} &= \sigma(y - \mathbf{x}),\\ \frac{\mathbf{dy}}{\mathbf{dt}} &= \mathbf{x}(\rho - z) - y,\\ \frac{\mathbf{dz}}{\mathbf{dt}} &= xy - \beta z, \end{aligned} \tag{12}$$

where the values of *σ* = 10, *β* = 8/4, *ρ* = 28 are usually used to exhibit a chaotic behavior. If we like to put Equation (12) with the notations introduced in Section 2, we can write **u** = [*x*, *y*, *z*] *<sup>T</sup>* with *n* = 3, and *θ* = [*σ*, *β*, *ρ*] *<sup>T</sup>* with *p* = 3. A Python function describing the dynamics of the Lorenz 63 system is given in Listing 3.

**Listing 3.** A Python function for the Lorenz 63 dynamics. ✞ ☎

```
import numpy as np
def Lorenz63(state,*args): #Lorenz 96 model
sigma = args[0]
beta = args[1]
rho = args[2]
x, y, z = state #Unpack the state vector
f = np.zeros(3) #Derivatives
f[0] = sigma * (y - x)
f[1] = x * (rho - z) - y
f[2] = x * y - beta * z
return f
```
Equation (12) describe the continuous-time evolution of the Lorenz system. In order to obtain the discrete-time mapping *M*(·; ·), a temporal integration scheme has to be applied. In Listing 4, one-step time integration functions are provided in Python using the first-order Euler and the fourth-order Runge–KuKutta schemes. Note that these functions requires a right-hand side function as input, this is mainly the continuous-time model *f*(**u**) (e.g., Listing 3).

✝ ✆

**Listing 4.** Python functions for the time integration using the 1st Euler and the 4th Runge–Kutta schemes. ✞ ☎

```
import numpy as np
def euler(rhs,state,dt,*args):
k1 = rhs(state,*args)
new_state = state + dt*k1
return new_state
def RK4(rhs,state,dt,*args):
k1 = rhs(state,*args)
k2 = rhs(state+k1*dt/2,*args)
k3 = rhs(state+k2*dt/2,*args)
k4 = rhs(state+k3*dt,*args)
new_state = state + (dt/6)*(k1+2*k2+2*k3+k4)
return new_state
```
For twin experiment testing, we suppose a true initial condition of **u***t*(0)=[1, 1, 1] *<sup>T</sup>* and measurements are collected each 0.2 time units for a total time of 2. We suppose that we measure

the full system state (i.e., *h*(**u**) = **u**, *m* = 3, and **H** = **I**3, where **I**<sup>3</sup> is the 3 × 3 identity matrix). Measurements are considered to be contaminated by a white Gaussian noise with a zero mean and a covariance matrix **R** = Diag(*σ*<sup>2</sup> <sup>1</sup> , *<sup>σ</sup>*<sup>2</sup> <sup>2</sup> , *<sup>σ</sup>*<sup>2</sup> <sup>3</sup> ). For simplicity, we let *σ*<sup>1</sup> = *σ*<sup>2</sup> = *σ*<sup>3</sup> = 0.15. For data assimilation testing, we assume that we begin with a perturbed initial condition of **u**(0)=[2, 3, 4] *T*. Then, background state values are computed at *t* = 0.2 by time integration of Equation (12) starting from this false initial condition. Observations at *t* = 0.2 are assimilated to provide the analysis at *t* = 0.2. After that, background state values are computed at *t* = 0.4 by time integration of Equation (12) starting from the analysis at *t* = 0.2, and so on. A sample implementation of the 3DVAR framework is presented in Listing 5, where a fixed background covariance matrix **B** = Diag(0.01, 0.01, 0.01) is assumed. Solution trajectories are presented in Figure 1 for a total time of 10, where observations are only available up to *t* = 2.

**Listing 5.** Implementation of the 3DVAR for the Lorenz 63 system. ✞ ☎

```
import numpy as np
import matplotlib.pyplot as plt
#%% Application: Lorenz 63