Trajectory simulation

CLODE can simulate ODE trajectories using the TrajectorySimulator class.

Example - FitzHugh-Nagumo oscillator

The following example simulates the FitzHugh-Nagumo oscillator using the RK45 integrator.

Python

>>> import clode
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from typing import List
>>> def fitzhugh_nagumo(
...     time: float,
...     variables: List[float],
...     parameters: List[float],
...     derivatives: list[float],
...     aux: list[float],
...     wiener: list[float],
... ) -> None:
...     V: float = variables[0]
...     w: float = variables[1]
... 
...     a: float = parameters[0]
...     b: float = parameters[1]
...     current: float = parameters[2]
...     epsilon: float = parameters[3]
... 
...     dV: float = V - V ** 3 / 3 - w + current
...     dw: float = epsilon * (V + a - b * w)
... 
...     derivatives[0] = dV
...     derivatives[1] = dw
... 
>>> a = 0.7
>>> variables = {"V": 1.0, "w": 0.0}
>>> parameters = {"a": a, "b": 0.8, "current": 0.0, "epsilon": 1.0 / 12.5}
>>> simulator = clode.TrajectorySimulator(
...     rhs_equation=fitzhugh_nagumo,
...     variables=variables,
...     parameters=parameters,
...     t_span=(0, 200),
...     dt=0.02,
... )
>>> ensemble_parameters = {"current": np.arange(0.0, 0.6, 0.1)}
>>> simulator.set_ensemble(parameters=ensemble_parameters)
>>> trajectories = simulator.trajectory()
>>> plt.figure(figsize=(8, 6))
<Figure size 800x600 with 0 Axes>
>>> for index in range(len(trajectories)):
...     label = f"I={ensemble_parameters['current'][index]:.1f}"
...     plt.plot(trajectories[index].x[:, 0], trajectories[index].x[:, 1], label=label)
... 
[<matplotlib.lines.Line2D object at 0x10ff6cdd0>]
[<matplotlib.lines.Line2D object at 0x11561a3f0>]
[<matplotlib.lines.Line2D object at 0x115619010>]
[<matplotlib.lines.Line2D object at 0x115618530>]
[<matplotlib.lines.Line2D object at 0x11561b1d0>]
[<matplotlib.lines.Line2D object at 0x115619e50>]
>>> plt.xlabel("V")
Text(0.5, 0, 'V')
>>> plt.ylabel("w")
Text(0, 0.5, 'w')
>>> plt.legend()
<matplotlib.legend.Legend object at 0x117a10290>
>>> plt.title("FitzHugh-Nagumo phase plane")
Text(0.5, 1.0, 'FitzHugh-Nagumo phase plane')
>>> plt.show()
>>> plt.figure(figsize=(8, 6))
<Figure size 800x600 with 0 Axes>
>>> for index in range(0, len(trajectories), 2):
...     label = f"I={ensemble_parameters['current'][index]}"
...     plt.plot(trajectories[index].t, trajectories[index].x[:, 0], label=label)
... 
[<matplotlib.lines.Line2D object at 0x114e68a10>]
[<matplotlib.lines.Line2D object at 0x114e6aa50>]
[<matplotlib.lines.Line2D object at 0x114e6a2d0>]
>>> plt.xlabel("t")
Text(0.5, 0, 't')
>>> plt.ylabel("V")
Text(0, 0.5, 'V')
>>> plt.legend()
<matplotlib.legend.Legend object at 0x115381490>
>>> plt.title("FitzHugh-Nagumo time series")
Text(0.5, 1.0, 'FitzHugh-Nagumo time series')
>>> plt.show()