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

[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