position :math:`x` in 2D and a parameter :math:`\mu` in 2D too. It therefore
enables a pedagogical illustration of optimal measurement points.
-This function depends on the position :math:`x=(x_1,x_2)\in\Omega=[0.1,0.9]^2`
-in the 2D plane, and on the parameter
+This function :math:`G` depends on the position
+:math:`x=(x_1,x_2)\in\Omega=[0.1,0.9]^2` in the 2D plane, and on the parameter
:math:`\mu=(\mu_1,\mu_2)\in\mathcal{D}=[-1,-0.01]^2` of dimension 2 :
.. math:: G(x;\mu) = \frac{1}{\sqrt{(x_1 - \mu_1)^2 + (x_2 - \mu_2)^2 + 0.1^2}}
20x20 points. It is available in ADAO built-in test models under the name
``TwoDimensionalInverseDistanceCS2010``, together with the spatial and
parametric domain default definition. So here we first build a set of
-simulations of :math:`G`, then we look for the best locations for measurements
-to obtain an DEIM interpolation representation of the fields, by applying the
-DEIM-type decomposition algorithm to it, and then derive some simple
-illustrations. We choose to look for an arbitrary number ``nbmeasures`` of 15
-measurement locations.
+simulations of :math:`G` for various :math:`\mu` parameter values, then we look
+for the best locations for measurements to obtain an DEIM interpolation
+representation of the fields, by applying the DEIM-type decomposition algorithm
+to it, and then derive some simple illustrations. We choose to look for an
+arbitrary number ``nbmeasures`` of 15 measurement locations.
It can be seen that the singular values decrease steadily down to numerical
noise, indicating that around a hundred basis elements are needed to fully
-represent the information contained in the set of :math:`G` simulations.
-Furthermore, the optimal measurement points in the :math:`\Omega_G` domain are
-inhomogeneously distributed, favoring the spatial zone near the
-:math:`(0.1,0.1)` corner in which the :math:`G` function varies more.
+represent the information contained in the set of :math:`G` parametric
+simulations. Furthermore, the optimal measurement points in the
+:math:`\Omega_G` domain are inhomogeneously distributed, favoring the spatial
+zone near the :math:`(0.1,0.1)` corner in which the :math:`G` function varies
+more.
:align: center
Python, 3.6.5, 3.12.6
- Numpy, 1.14.3, 2.1.2
+ Numpy, 1.14.3, 2.1.3
Scipy, 0.19.1, 1.14.1
MatplotLib, 2.2.2, 3.9.2
GnuplotPy, 1.8, 1.8
donc une illustration pédagogique pour représenter les points optimaux de
mesure.
-Cette fonction dépend de la position :math:`x=(x_1,x_2)\in\Omega=[0.1,0.9]^2`
-dans le plan 2D, et du paramètre
+Cette fonction :math:`G` dépend de la position
+:math:`x=(x_1,x_2)\in\Omega=[0.1,0.9]^2` dans le plan 2D, et du paramètre
:math:`\mu=(\mu_1,\mu_2)\in\mathcal{D}=[-1,-0.01]^2` de dimension 2 :
.. math:: G(x;\mu) = \frac{1}{\sqrt{(x_1 - \mu_1)^2 + (x_2 - \mu_2)^2 + 0.1^2}}
La fonction est représenté sur une grille spatiale régulière :math:`\Omega_G`
de taille 20x20 points. Elle est disponible dans les modèles de tests intégrés
pour ADAO sous le nom ``TwoDimensionalInverseDistanceCS2010``. On construit
-donc ici tout d'abord un ensemble de simulations de :math:`G`, puis on cherche
-les meilleurs positions de mesures pour obtenir une représentation par
-interpolation DEIM des champs, en appliquant l'algorithme de décomposition de
-type DEIM, et on en tire ensuite des illustrations simples. On choisit de
-rechercher un nombre arbitraire ``nbmeasures`` de 15 positions de mesures.
+donc ici tout d'abord un ensemble de simulations de :math:`G` pour différents
+paramètres :math:`\mu`, puis on cherche les meilleurs positions de mesures pour
+obtenir une représentation par interpolation DEIM des champs, en appliquant
+l'algorithme de décomposition de type DEIM, et on en tire ensuite des
+illustrations simples. On choisit de rechercher un nombre arbitraire
+``nbmeasures`` de 15 positions de mesures.
On observe ainsi que les valeurs singulières décroissent régulièrement jusqu'au
bruit numérique, indiquant qu'il faut environ une centaine d'éléments de base
pour complètement représenter l'information contenue dans l'ensemble des
-simulations de :math:`G`. Par ailleurs, les points de mesure optimaux dans le
-domaine :math:`\Omega_G` sont répartis de manière inhomogène, privilégiant la
-zone spatiale proche du coin :math:`(0.1,0.1)` dans laquelle la fonction
-:math:`G` varie plus.
+simulations paramétriques de :math:`G`. Par ailleurs, les points de mesure
+optimaux dans le domaine :math:`\Omega_G` sont répartis de manière inhomogène,
+privilégiant la zone spatiale proche du coin :math:`(0.1,0.1)` dans laquelle la
+fonction :math:`G` varie le plus.
:align: center
Python, 3.6.5, 3.12.6
- Numpy, 1.14.3, 2.1.2
+ Numpy, 1.14.3, 2.1.3
Scipy, 0.19.1, 1.14.1
MatplotLib, 2.2.2, 3.9.2
GnuplotPy, 1.8, 1.8
"""
Classe de simulateur ODE d'ordre 1 pour modèles dynamiques :
- dy / dt = F_µ(t, y)
+ dy/dt = F_µ(t, y)
avec y = f(t) et µ les paramètres intrinsèques. t est couramment le temps,
mais il peut être une variable quelconque non temporelle.
autonomous=None,
):
"None default values are mandatory to allow overriding"
- self.set_canonical_description()
- self.set_description(mu, integrator, dt, t0, tf, y0, autonomous)
+ if hasattr(self, "CanonicalDescription"):
+ self.CanonicalDescription()
+ self._description(mu, integrator, dt, t0, tf, y0, autonomous)
# --------------------------------------------------------------------------
- # User defined ODE model
+ # User defined ODE model and canonical description
def ODEModel(self, t, y):
"""
"""
raise NotImplementedError()
- def set_canonical_description(self):
+ def CanonicalDescription(self):
"""
User settings for default or recommended EDO description
- Setters that >>> can <<< be used:
- - self.set_mu
- - self.set_integrator
- - self.set_dt
- - self.set_t0
- - self.set_tf
- - self.set_y0
- - self.set_autonomous
+ Setters/Getters that >>> can <<< be used:
+ - self.Parameters
+ - self.Integrator
+ - self.IntegrationStep
+ - self.ObservationStep
+ - self.InitialTime
+ - self.FinalTime
+ - self.InitialCondition
+ - self.Autonomous
"""
pass
# --------------------------------------------------------------------------
- def set_description(
- self,
- mu=None,
- integrator=None,
- dt=None,
- t0=None,
- tf=None,
- y0=None,
- autonomous=False,
- ):
- "Explicit setting of EDO description"
- self.set_mu(mu)
- self.set_integrator(integrator)
- self.set_dt(dt)
- self.set_t0(t0)
- self.set_tf(tf)
- self.set_y0(y0)
- self.set_autonomous(autonomous)
-
- def set_mu(self, mu=None):
- "Set EDO intrinsic parameters µ"
- if mu is not None:
- self._mu = numpy.ravel(mu)
+ @property
+ def Parameters(self):
return self._mu
- def set_integrator(self, integrator=None):
- "Set integration scheme name"
- if integrator is None:
+ @Parameters.setter
+ def Parameters(self, value):
+ if value is None:
+ pass
+ else:
+ self._mu = numpy.ravel(value)
+
+ # -------
+ @property
+ def Integrator(self):
+ return self._integrator
+
+ @Integrator.setter
+ def Integrator(self, value):
+ if value is None:
pass
- elif not (integrator in self.__integrator_list):
+ elif not (value in self.__integrator_list):
raise ValueError(
- "Wrong value %s for integrator in set_integrator call. \nAvailable integrators are: %s"
- % (integrator, self.__integrator_list)
+ "wrong value %s set for the integrator scheme. \nAvailable integrator scheme are: %s"
+ % (value, self.__integrator_list)
)
else:
- self._integrator = str(integrator)
- return self._integrator
+ self._integrator = str(value)
- def set_dt(self, value=None):
- "Set integration step size dt"
- if value is not None:
- self._dt = max(2.0e-16, abs(float(value)))
+ # -------
+ @property
+ def IntegrationStep(self):
return self._dt
- def set_t0(self, value=None):
- "Set initial integration time"
- if value is not None:
+ @IntegrationStep.setter
+ def IntegrationStep(self, value):
+ if value is None:
+ pass
+ elif float(value) > 0:
+ self._dt = max(2.0e-16, float(value))
+ else:
+ raise ValueError("integration step has to be strictly positive")
+
+ # -------
+ @property
+ def ObservationStep(self):
+ return self._do
+
+ @ObservationStep.setter
+ def ObservationStep(self, value):
+ if value is None:
+ pass
+ elif float(value) > 0:
+ self._do = max(2.0e-16, float(value))
+ if hasattr(self, "_dt") and self._do < self._dt:
+ raise ValueError(
+ "Observation step is inconsistent with integration one"
+ )
+ else:
+ raise ValueError("observation step has to be strictly positive")
+
+ # -------
+ @property
+ def InitialTime(self):
+ return self._t0
+
+ @InitialTime.setter
+ def InitialTime(self, value):
+ if value is None:
+ pass
+ else:
self._t0 = float(value)
if hasattr(self, "_tf") and self._t0 > self._tf:
- raise ValueError("Initial time has to remain less than final time")
- return self._t0
+ raise ValueError("initial time has to remain less than final time")
+
+ # -------
+ @property
+ def FinalTime(self):
+ return self._tf
- def set_tf(self, value=None):
- "Set final integration time"
- if value is not None:
+ @FinalTime.setter
+ def FinalTime(self, value):
+ if value is None:
+ pass
+ else:
self._tf = float(value)
if hasattr(self, "_t0") and self._t0 > self._tf:
- raise ValueError("Initial time has to remain less than final time")
- return self._tf
+ raise ValueError("initial time has to remain less than final time")
- def set_y0(self, value):
- "Set integration initial condition"
- if value is not None:
- self._y0 = numpy.ravel(value)
+ # -------
+ @property
+ def InitialCondition(self):
return self._y0
- def set_autonomous(self, value=None):
+ @InitialCondition.setter
+ def InitialCondition(self, value):
+ if value is None:
+ pass
+ else:
+ self._y0 = numpy.ravel(value)
+
+ # -------
+ @property
+ def Autonomous(self):
+ return self._autonomous
+
+ @Autonomous.setter
+ def Autonomous(self, value):
"Declare the system to be autonomous"
- if value is not None:
+ if value is None:
+ pass
+ else:
self._autonomous = bool(value)
- return self._autonomous
- def set_do(self, value=None):
- "Set observation step size do"
- if value is not None:
- self._do = max(2.0e-16, abs(float(value)))
- return self._do
+ # --------------------------------------------------------------------------
+
+ def _description(
+ self,
+ mu=None,
+ integrator=None,
+ dt=None,
+ t0=None,
+ tf=None,
+ y0=None,
+ autonomous=False,
+ ):
+ "Explicit setting of EDO description at init time"
+ self.Parameters = mu
+ self.Integrator = integrator
+ self.IntegrationStep = dt
+ self.InitialTime = t0
+ self.FinalTime = tf
+ self.InitialCondition = y0
+ self.Autonomous = autonomous
# --------------------------------------------------------------------------
_euler_step = _rk1_step
- def Integration(self, y0=None, t0=None, tf=None, mu=None):
+ def Integration(self, y1=None, t1=None, t2=None, mu=None):
"""
- Exécute l'intégration du modèle entre t0 et tf, en partant de y0,
- via le schéma d'intégration choisi
+ Exécute l'intégration du modèle entre t1 et t2, en partant de y1, via
+ le schéma d'intégration choisi.
"""
- if y0 is not None:
- self.set_y0(y0)
- if t0 is not None:
- self.set_t0(t0)
- if tf is not None:
- self.set_tf(tf)
- if mu is not None:
- self.set_mu(mu)
+ if y1 is None:
+ _ly0 = self._y0
+ else:
+ _ly0 = numpy.ravel(y1)
+ if t1 is None:
+ _lt0 = self._t0
+ else:
+ _lt0 = float(t1)
+ if t2 is None:
+ _ltf = self._tf
+ else:
+ _ltf = float(t2)
+ self.Parameters = mu
if (
(self._mu is None)
or (self._integrator is None)
or (self._dt is None)
- or (self._t0 is None)
- or (self._tf is None)
- or (self._y0 is None)
+ or (_lt0 is None)
+ or (_ltf is None)
+ or (_ly0 is None)
):
raise ValueError(
"Some integration input information are None and not defined\n(%s, %s, %s, %s, %s, %s)"
self._mu,
self._integrator,
self._dt,
- self._t0,
- self._tf,
- self._y0,
+ _lt0,
+ _ltf,
+ _ly0,
)
)
#
ODE = self.ODEModel
- times = numpy.arange(self._t0, self._tf + self._dt / 2, self._dt)
+ times = numpy.arange(_lt0, _ltf + self._dt / 2, self._dt)
if self._integrator == "odeint":
# intégration 'automatique' dans le cas d'un système pouvant être
# problématique avec rk4 ou euler (comme Van Der Pol)
trajectory = odeint(
ODE,
- numpy.array(self._y0, dtype=float),
+ numpy.array(_ly0, dtype=float),
times,
tfirst=True,
)
sol = solve_ivp(
ODE,
- (self._t0, self._tf),
- numpy.array(self._y0, dtype=float),
+ (_lt0, _ltf),
+ numpy.array(_ly0, dtype=float),
t_eval=times,
)
trajectory = sol.y.T
% self._integrator
)
#
- t = self._t0
- y = self._y0
- trajectory = numpy.array([self._y0])
+ t = _lt0
+ y = _ly0
+ trajectory = numpy.array([_ly0])
#
- while t < self._tf - self._dt / 2:
+ while t < _ltf - self._dt / 2:
[t, y] = integration_step(t, y, self._dt, ODE)
trajectory = numpy.concatenate((trajectory, numpy.array([y])), axis=0)
#
return [times, trajectory]
def ForecastedPath(self, y1=None, t1=None, t2=None, mu=None):
- "Trajectoire de t1 à t2, en partant de yn, pour un paramètre donné mu"
+ "Trajectoire de t1 à t2, en partant de y1, pour un paramètre donné mu"
#
_, trajectory_from_t1_to_t2 = self.Integration(y1, t1, t2, mu)
#
return trajectory_from_t1_to_t2[-1, :]
def StateTransition(self, y1=None):
- "État y[n+1] intégré depuis y[n] sur pas constant ou non"
- if self.set_autonomous():
+ "État y[n+1] intégré depuis y[n] sur un pas self.ObservationStep"
+ if self.Autonomous:
if not hasattr(self, "_do") or self._do is None:
raise ValueError(
" StateTransition requires an observation step size to be given"
)
- return self.ForecastedState(y1, t1=0.0, t2=self.set_do())
+ return self.ForecastedState(y1, 0.0, self.ObservationStep, self.Parameters)
else:
raise NotImplementedError(
" StateTransition has to be provided by the user in case of non-autonomous ODE"
)
- def HistoryBoard(self, t_s, i_s, y_s, filename="history_board_of_trajectory.pdf"):
+ def HistoryBoard(
+ self, t_s, y_s, i_s=None, filename="figure_of_trajectory.pdf", suptitle="", title="", xlabel="Time", ylabel="State variables", cmap="gist_gray_r",
+ ):
"""
t_s : série des instants t
- i_s : série des indices i des variables
y_s : série des valeurs 1D des variables du système dynamique, pour
chaque pas de temps, sous forme d'un tableau 2D de type:
SDyn(i,t) = SDyn[i][t] = [SDyn[i] pour chaque t]
+ i_s : série des indices i des variables
"""
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
+ #
+ if i_s is None:
+ i_s = range(y_s.shape[0])
levels = MaxNLocator(nbins=25).tick_values(
numpy.ravel(y_s).min(), numpy.ravel(y_s).max()
)
fig, ax = plt.subplots(figsize=(15, 5))
fig.subplots_adjust(bottom=0.1, left=0.05, right=0.95, top=0.9)
im = plt.contourf(
- t_s, i_s, y_s, levels=levels, cmap=plt.get_cmap("gist_gray_r")
+ t_s, i_s, y_s, levels=levels, cmap=plt.get_cmap(cmap)
)
fig.colorbar(im, ax=ax)
- plt.title("Model trajectory with %i variables" % len(y_s[:, 0]))
- plt.xlabel("Time")
- plt.ylabel("State variables")
+ if len(suptitle) > 0:
+ plt.suptitle(suptitle)
+ if len(title) > 0:
+ plt.title(title)
+ else:
+ plt.title("Model trajectory with %i variables" % len(y_s[:, 0]))
+ plt.xlabel(xlabel)
+ plt.ylabel(ylabel)
if filename is None:
plt.show()
else:
import sys
import unittest
-import math
import numpy
from daCore.BasicObjects import DynamicalSimulator
# ==============================================================================
-class Lorenz1963(DynamicalSimulator):
+class EDS(DynamicalSimulator):
"""
Three-dimensional parametrized nonlinear ODE system depending on µ=(σ,ρ,β):
- ∂x/∂t = σ (y − x)
- ∂y/∂t = ρ x − y − x z
- ∂z/∂t = x y − β z
+ dx/dt = σ (y − x)
+ dy/dt = ρ x − y − x z
+ dz/dt = x y − β z
with t ∈ [0, 40] the time interval, x(t), y(t), z(t) the dependent
variables, and with σ=10, ρ=28, and β=8/3 the commonly used parameter
values. The initial conditions for (x, y, z) at t=0 for the reference
case are (0, 1, 0).
- This is the well known parametrized coupled system of three nonlinear
- ordinary differential equations:
+ This is a parametrized coupled system of three nonlinear ordinary
+ differential equations described in:
Lorenz, E. N. (1963). Deterministic nonperiodic flow. Journal of the
Atmospheric Sciences, 20, 130–141.
doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
"""
- def set_canonical_description(self):
- self.set_mu((10.0, 28.0, 8.0 / 3.0)) # µ = (σ, ρ, β)
- self.set_integrator("rk4")
- self.set_dt(0.01)
- self.set_t0(0.0)
- self.set_tf(40)
- self.set_y0((0.0, 1.0, 0.0))
- self.set_autonomous(True)
+ def CanonicalDescription(self):
+ self.Parameters = (10.0, 28.0, 8.0 / 3.0) # µ = (σ, ρ, β)
+ self.Integrator = "rk4"
+ self.IntegrationStep = 0.01
+ self.InitialTime = 0.0
+ self.FinalTime = 40
+ self.InitialCondition = (0.0, 1.0, 0.0)
+ self.Autonomous = True
return True
def ODEModel(self, t, Y):
- "ODE dY/dt = F(Y,t)"
- sigma, rho, beta = self.set_mu()
+ "ODE dy/dt = F_µ(t, y)"
+ sigma, rho, beta = self.Parameters
x, y, z = map(float, Y)
#
- rx = sigma * (y - x)
- ry = x * rho - y - x * z
- rz = x * y - beta * z
+ dxdt = sigma * (y - x)
+ dydt = x * rho - y - x * z
+ dzdt = x * y - beta * z
#
- return numpy.array([rx, ry, rz])
+ return numpy.array([dxdt, dydt, dzdt])
# ==============================================================================
@classmethod
def setUpClass(cls):
print("\nAUTODIAGNOSTIC\n==============\n")
- print(" " + Lorenz1963().__doc__.strip())
+ print(" " + EDS().__doc__.strip())
def test001(self):
- numpy.random.seed(123456789)
- ODE = Lorenz1963() # Default parameters
+ ODE = EDS() # Default parameters
+ ODE.ObservationStep = 0.2
trajectory = ODE.ForecastedPath()
+ lastvalue = numpy.array([16.48799962, 14.01693428, 40.30448848])
#
print()
- self.assertTrue(trajectory.shape[0] == 1 + int(ODE.set_tf() / ODE.set_dt()))
self.assertTrue(
- abs(
- max(
- trajectory[-1]
- - numpy.array([16.48799962, 14.01693428, 40.30448848])
- )
- )
- <= 1.0e-8,
- msg=" Last value is not equal to the reference one",
+ trajectory.shape[0] == 1 + int(ODE.FinalTime / ODE.IntegrationStep)
+ )
+ erreur = abs(max(trajectory[-1] - lastvalue))
+ self.assertTrue(
+ erreur <= 1.0e-8,
+ msg=" Last value is not equal to the reference one. Error = %.2e"
+ % erreur,
)
print(" Last value is equal to the reference one")
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2008-2024 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library 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 library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+import sys
+import unittest
+import numpy as np
+from daCore.BasicObjects import DynamicalSimulator
+
+
+# ==============================================================================
+class EDS(DynamicalSimulator):
+ """
+ N-dimensional parametrized nonlinear ODE system depending on a F forcing
+ term:
+
+ dy[i]/dt = (y[i+1] − y[i-2]) * y[i-1] − y[i] + F
+
+ with i=1,...,N and periodic conditions coming from circle-like domain
+ for the N variables: y[0]=y[N], y[-1]=y[N-1], y[1]=y[N+1]. With a
+ forcing term F value of 8, the dynamics of this model is chaotic.
+
+ This is a parametrized system of nonlinear ordinary differential equations,
+ first described during a seminar at the European Centre for Medium-Range
+ Weather Forecasts (ECMWF) in the Autumn of 1995 (Seminar on Predictability,
+ 4-8 September 1995), the proceedings of which were published as Lorenz
+ (1996):
+ Lorenz, E.N. Predictability: A problem partly solved. In Proceedings of
+ the Seminar on Predictability, ECMWF, Reading, UK, 9–11 September 1996,
+ Volume 1, pp. 1-18.
+ See:
+ https://www.ecmwf.int/en/elibrary/75462-predictability-problem-partly-solved
+ The system is chaotic and numerically very sensible, even on the order of
+ numerical calculations.
+ """
+
+ def CanonicalDescription(self):
+ N = 40
+ self.Parameters = (N, 8.0) # N, F
+ self.Integrator = "odeint"
+ self.IntegrationStep = 0.05
+ self.InitialTime = 0.0
+ self.FinalTime = 1.0
+ self.InitialCondition = np.arange(N)
+ self.Autonomous = True
+ return True
+
+ def ODEModel(self, t, Y):
+ "ODE dy/dt = F_µ(t, y)"
+ N, F = self.Parameters
+ N = int(N)
+ F = float(F)
+ Y = np.ravel(Y)
+ dydt = np.zeros(np.shape(Y))
+ assert len(Y) == N, "%i =/= %i" % (len(Y), N)
+ #
+ # Boundary case equations (rank 1, 2, N):
+ dydt[0] = (Y[1] - Y[N - 2]) * Y[N - 1] - Y[0]
+ dydt[1] = (Y[2] - Y[N - 1]) * Y[0] - Y[1]
+ dydt[-1] = (Y[0] - Y[-3]) * Y[-2] - Y[-1]
+ # General indices (rank 3 to N-1)
+ dydt[2:-1] = (Y[3:] - Y[:-3]) * Y[1:-2] - Y[2:-1]
+ # Adding forcing
+ dydt = dydt + F
+ #
+ return dydt
+
+
+# ==============================================================================
+class LocalTest(unittest.TestCase):
+ @classmethod
+ def setUpClass(cls):
+ print("\nAUTODIAGNOSTIC\n==============\n")
+ print(" " + EDS().__doc__.strip())
+
+ def test001(self):
+ ODE = EDS() # Default parameters
+ ODE.ObservationStep = 0.2
+ trajectory = ODE.ForecastedPath()
+ lastvalue = np.array(
+ [
+ 4.87526450e00,
+ 1.47956922e00,
+ -3.59900230e00,
+ 1.04256938e00,
+ 3.63094386e00,
+ 1.88075245e01,
+ -1.67896983e01,
+ 9.47857563e00,
+ -4.10761433e01,
+ -7.22094136e00,
+ 1.82999964e01,
+ -9.50258813e00,
+ 1.97970038e01,
+ -9.14157688e-01,
+ -3.23915931e00,
+ 9.08515096e00,
+ 1.90398329e01,
+ 6.63360465e00,
+ 3.62060706e00,
+ 1.40450966e00,
+ -4.23014939e00,
+ 3.12067035e00,
+ -7.45630836e00,
+ 5.29044964e00,
+ -4.95916975e00,
+ 1.49074059e00,
+ -5.34660284e00,
+ 7.88846593e00,
+ 9.82088697e00,
+ 9.01760941e-01,
+ 7.27891147e00,
+ 4.66438119e00,
+ -1.22290917e00,
+ 3.38177733e00,
+ 6.19234382e00,
+ 3.66288450e00,
+ 2.51247077e00,
+ 3.61265491e00,
+ 4.56523506e00,
+ 5.03793331e00,
+ ]
+ )
+ #
+ print()
+ self.assertTrue(
+ trajectory.shape[0] == 1 + int(ODE.FinalTime / ODE.IntegrationStep)
+ )
+ erreur = abs(max(trajectory[-1] - lastvalue))
+ self.assertTrue(
+ erreur <= 1.0e-7,
+ msg=" Last value is not equal to the reference one. Error = %.2e"
+ % erreur,
+ )
+ print(" Last value is equal to the reference one")
+
+ def tearDown(cls):
+ print("\n Tests are finished\n")
+
+
+# ==============================================================================
+if __name__ == "__main__":
+ sys.stderr = sys.stdout
+ unittest.main(verbosity=0)
