Source code for gemseo.problems.ode.van_der_pol

# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
# Contributors:
#    INITIAL AUTHORS - API and implementation and/or documentation
#        :author: Isabelle Santos
#    OTHER AUTHORS   - MACROSCOPIC CHANGES
r"""The Van der Pol (VDP) problem describing an oscillator with non-linear damping.

Van der Pol, B. & Van Der Mark, J.
Frequency Demultiplication.
Nature 120, 363-364 (1927).

The Van der Pol problem is written as follows:

.. math::

    \frac{d^2 x(t)}{dt^2} -
    \mu (1-x(t)^2) \frac{dx(t)}{dt} + x = 0

where :math:`x(t)` is the position coordinate as a function of time, and
:math:`\mu` is a scalar parameter indicating the stiffness.

This problem can be rewrittent in a 2-dimensional form with only first-order
derivatives. Let :math:`y = \frac{dx}{dt}` and
:math:`s = \begin{pmatrix}x\\y\end{pmatrix}`. Then the Van der Pol problem is:

.. math::

    \frac{ds}{dt} = f(s, t)

with

.. math::

    f : s = \begin{pmatrix} x \\ y \end{pmatrix} \mapsto
    \begin{pmatrix} y \\ \mu (1-x^2) y - x \end{pmatrix}

The jacobian of this function can be expressed analytically:

.. math::

    \mathrm{Jac}\, f = \begin{pmatrix}
        0 & 1 \\
        -2\mu xy - 1 & \mu (1 - x^2)
    \end{pmatrix}

There is no exact solution to the Van der Pol oscillator problem in terms of
known tabulated functions (see Panayotounakos *et al.*, On the Lack of Analytic
Solutions of the Van Der Pol Oscillator. ZAMM 83, nᵒ 9 (1 septembre 2003)).
"""

from __future__ import annotations

from typing import TYPE_CHECKING

from numpy import array
from numpy import zeros

from gemseo.algos.ode.ode_problem import ODEProblem

if TYPE_CHECKING:
    from numpy.typing import NDArray




class VanDerPol(ODEProblem):
    """Representation of an oscillator with non-linear damping."""

    _mu: float
    r"""Stiffness parameter.

    Van der Pol is stiffer with larger values of :math:`\mu`.
    """

    state_vect: NDArray[float]
    r"""State vector :math:`s=(x, \dot{x})` of the system."""

    def __init__(
        self,
        initial_time: float = 0,
        final_time: float = 0.5,
        mu: float = 1.0e3,
        use_jacobian: bool = True,
        state_vector: NDArray[float] = None,
    ) -> None:
        """
        Args:
            mu: The stiffness parameter.
            initial_time: The start of the integration interval.
            final_time: The end of the integration interval.
            use_jacobian: Whether to use the analytical expression of the Jacobian.
                If false, use finite differences to estimate the Jacobian.
            state_vector: The state vector of the system.
        """  # noqa: D205 D212
        self._mu = mu

        if state_vector is None:
            state_vector = zeros(2)
        self.state_vect = state_vector

        initial_state = array([
            2 + self.state_vect[0],
            -2 / 3
            + 10 / (81 * self._mu)
            - 292 / (2187 * self._mu * self._mu)
            + self.state_vect[1],
        ])

        jac = self.__compute_rhs_jacobian if use_jacobian else None
        super().__init__(
            func=self.__compute_rhs,
            jac=jac,
            initial_state=initial_state,
            initial_time=initial_time,
            final_time=final_time,
        )

    def __compute_rhs(self, time: float, state: NDArray[float]) -> NDArray[float]:  # noqa:U100
        """Compute the right-hand side of the ODE.

        Args:
            time: The time at which the function is evaluated.
            state: The state in which the function is evaluated.

        Returns:
             The value of the right-hand side of the ODE.
        """
        return array([state[1], self._mu * state[1] * (1 - state[0] ** 2) - state[0]])

    def __compute_rhs_jacobian(
        self, time: float, state: NDArray[float]
    ) -> NDArray[float]:  # noqa:U100
        """Compute the Jacobian of the right-hand side of the ODE.

        Args:
            time: The time at which the function is evaluated.
            state: The state in which the function is evaluated.

        Returns:
             The Jacobian of the right-hand side of the ODE.
        """
        jac = zeros((2, 2))
        jac[1, 0] = -self._mu * 2 * state[1] * state[0] - 1
        jac[0, 1] = 1
        jac[1, 1] = self._mu * (1 - state[0] * state[0])
        return jac