Source code for gemseo.problems.ode.van_der_pol

# Copyright 2021 IRT Saint Exupéry,
# 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
# 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
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)


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

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

    from numpy.typing import NDArray

[docs] 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