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