feat: add Chua system example

docs/gallery_chaos.rst

@@ -205,6 +205,28 @@ Three Dimensional Systems
     [Petras2011]_. The complete setup (with parameters) can be found in
     :download:`examples/gallery/chen.py <../examples/gallery/chen.py>`.
 
+.. card:: Fractional-order Chua System
+    :class-title: sd-text-center
+
+    .. image:: _static/gallery-chua-light.svg
+        :class: only-light
+        :width: 65%
+        :align: center
+        :alt: Fractional-order Chua system phase diagram
+
+    .. image:: _static/gallery-chua-dark.svg
+        :class: only-dark
+        :width: 65%
+        :align: center
+        :alt: Fractional-order Chua system phase diagram
+
+    +++
+
+    This example uses the Caputo derivative and the
+    :class:`~pycaputo.fode.gallery.Chua` class to reproduce
+    Figure 5.6 from [Petras2011]_. The complete setup (with parameters) can be
+    found in :download:`examples/gallery/chua.py <../examples/gallery/chua.py>`.
+
 .. card:: Fractional-order Cellular Neural Network (3 cells) System
     :class-title: sd-text-center
 

examples/gallery/chua.py

@@ -0,0 +1,37 @@
+# SPDX-FileCopyrightText: 2024 Alexandru Fikl <alexfikl@gmail.com>
+# SPDX-License-Identifier: MIT
+
+from __future__ import annotations
+
+import numpy as np
+
+from pycaputo import fracevolve, fracplot
+from pycaputo.fode.gallery import Chua
+
+# setup system (parameters from Figure 5.6 from [Petras2011])
+alpha = (0.93, 0.99, 0.92)
+y0 = np.array([0.2, -0.1, 0.1])
+func = Chua(alpha=10.725, beta=10.593, gamma=0.268, m=(-1.1726, -0.7872))
+
+print(f"alpha {alpha} y0 {y0} parameters {func}")
+
+# setup controller
+from pycaputo.controller import make_fixed_controller
+
+dt = 5.0e-3
+control = make_fixed_controller(dt, tstart=0.0, tfinal=100.0)
+print(control)
+
+# setup stepper
+from pycaputo.fode import caputo
+
+stepper = caputo.PECE(
+    derivative_order=alpha,
+    control=control,
+    source=func.source,
+    y0=(y0,),
+    corrector_iterations=1,
+)
+
+solution = fracevolve(stepper, dtinit=dt)
+fracplot(solution, "gallery-chua")

examples/integrate-and-fire-pif.py

@@ -38,8 +38,7 @@
 logger.info("tspike %.8e tstart %.8e, tfinal %.8e", tspikes[0], tstart, tfinal)
 if tspikes.size < 2:
     raise ValueError(
-        "This example expects at least two spikes. "
-        "Try increasing 'alpha' or 'tfinal'."
+        "This example expects at least two spikes. Try increasing 'alpha' or 'tfinal'."
     )
 
 dtinit = 1.0e-1

scripts/generate-doc-figures.py

@@ -42,6 +42,9 @@
     "gallery/chen.py": {
         "gallery-chen.svg",
     },
+    "gallery/chua.py": {
+        "gallery-chua.svg",
+    },
     "gallery/cnn.py": {
         "gallery-cnn.svg",
     },

src/pycaputo/fode/gallery.py

@@ -157,6 +157,90 @@ def source_jac(self, t: float, y: Array) -> Array:
 # }}}
 
 
+# {{{ Chua
+
+
+@dataclass(frozen=True)
+class Chua(Function):
+    r"""Implements the right-hand side of the Chua system (see Equation 5.13
+    from [Petras2011]_).
+
+    .. math::
+
+        \begin{aligned}
+        D^\alpha[x](t) & =
+            \alpha (y - x + f(x)), \\
+        D^\alpha[y](t) & =
+            x - y + z, \\
+        D^\alpha[z](t) & =
+            -\beta y - \gamma z
+        \end{aligned}
+    """
+
+    alpha: float
+    """Non-dimensional parameter in the Chua system."""
+    beta: float
+    """Non-dimensional parameter in the Chua system."""
+    gamma: float
+    """Non-dimensional parameter in the Chua system."""
+    m: tuple[float, float]
+    """Parameters used to define describing the diode (Equation 5.5)."""
+
+    def source(self, t: float, y: Array) -> Array:
+        m0, m1 = self.m
+        fx = m1 * y[0] + 0.5 * (m0 - m1) * (abs(y[0] + 1) - abs(y[0] - 1))
+
+        return np.array([
+            self.alpha * (y[1] - y[0] - fx),
+            y[0] - y[1] + y[2],
+            -self.beta * y[1] - self.gamma * y[2],
+        ])
+
+    def source_jac(self, t: float, y: Array) -> Array:
+        m0, m1 = self.m
+        df = m0 if abs(y[0]) < 1.0 else m1
+
+        return np.array([
+            [-self.alpha * (df + 1.0), self.alpha, 0.0],
+            [1.0, -1.0, 1.0],
+            [0.0, -self.beta, -self.gamma],
+        ])
+
+    @classmethod
+    def from_dim(
+        cls,
+        C1: float,
+        C2: float,
+        R2: float,
+        RL: float,
+        L1: float,
+        Ga: float,
+        Gb: float,
+    ) -> Chua:
+        """Construct the non-dimensional system from the standard dimensional
+        parameters of the Chua system.
+
+        :arg C1: capacitance (in nano Farad) of first capacitor.
+        :arg C2: capacitance (in nano Farad) of second capacitor.
+        :arg R2: resistance (in Ohm).
+        :arg RL: resistance (in Ohm).
+        :arg L1: inductance (in mili Henry).
+        :arg Ga: slope of nonlinear function (in mili Ampere per Volt).
+        :arg Gb: slope of nonlinear function (in mili Ampere per Volt).
+        """
+        G = 1.0 / R2
+        alpha = C2 / C1
+        beta = C2 / (L1 * G**2)
+        gamma = C2 * R2 / (L1 * G)
+        m0 = Ga / G
+        m1 = Gb / G
+
+        return Chua(alpha=alpha, beta=beta, gamma=gamma, m=(m0, m1))
+
+
+# }}}
+
+
 # {{{ Cellular Neural Network
 
 

src/pycaputo/quadrature/riemann_liouville.py

@@ -747,7 +747,7 @@ def _quad_rl_yuan_agrawal(
 ) -> Array:
     if not callable(f):
         raise TypeError(
-            f"{type(m).__name__!r} requires a callable: " f"f is a {type(f).__name__!r}"
+            f"{type(m).__name__!r} requires a callable: f is a {type(f).__name__!r}"
         )
 
     x = p.x