The simulation model

To construct a Lorenz-63 simulation model for pypfilt, we need to create a Model subclass that defines the state vector structure and the state update rule. Because this model is a system of ordinary differential equations, we can derive from the OdeModel class, which provides a convenient wrapper around scipy.integrate.solve_ivp().

Here, we define the state vector to contain six floating-point values:

\[x_t = [\sigma, \rho, \beta, x, y, z]^T\]

• The state vector structure is defined by the field_types() method, which returns a list of (name, type) pairs.

• The right-hand side is defined by the d_dt() method, which keeps the parameters \(\sigma\), \(\rho\), and \(\beta\) fixed, and calculates the rate of change for \(x(t)\), \(y(t)\), and \(z(t)\).

class Lorenz63(OdeModel):
    def field_types(self, ctx):
        r"""
        Define the state vector :math:`[\sigma, \rho, \beta, x, y, z]^T`.
        """
        return [
            ('sigma', float),
            ('rho', float),
            ('beta', float),
            ('x', float),
            ('y', float),
            ('z', float),
        ]

    def d_dt(self, time, xt, ctx, is_forecast):
        """
        The right-hand side of the ODE system.

        :param time: The current time.
        :param xt: The particle state vectors.
        :param ctx: The simulation context.
        :param is_forecast: True if this is a forecasting simulation.
        """
        rates = np.zeros(xt.shape, xt.dtype)
        rates['x'] = xt['sigma'] * (xt['y'] - xt['x'])
        rates['y'] = xt['x'] * (xt['rho'] - xt['z']) - xt['y']
        rates['z'] = xt['x'] * xt['y'] - xt['beta'] * xt['z']
        return rates