numicano: Numerical Integration of models in ODE of polynomial form

In GPoM: Generalized Polynomial Modelling

Numerical Integration of models in ODE of polynomial form

Description

Function for the numerical integration of Ordinary Differential Equations of polynomial form.

Usage

numicano(
  nVar,
  dMax,
  dMin = 0,
  Istep = 1000,
  onestep = 1/125,
  KL = NULL,
  PolyTerms = NULL,
  v0 = NULL,
  method = "rk4"
)

Arguments

nVar	Number of variables considered in the polynomial formulation.
dMax	Maximum degree of the polynomial formulation.
dMin	The minimum negative degree of the polynomial formulation (0 by default).
Istep	The number of integration time steps
onestep	Time step length
KL	Matrix formulation of the model to integrate numerically
PolyTerms	Vectorial formulation of the model (only for models of canonical form)
v0	The initial conditions (a vector which length should correspond to the model dimension nVar)
method	The integration method (See package deSolve), by default method = 'rk4'.

Value

A list of two variables:

$KL The model in its matrix formulation

$reconstr The integrated trajectory (first column is the time, next columns are the model variables)

Author(s)

Sylvain Mangiarotti

See Also

derivODE2, numinoisy

Examples

#############
# Example 1 #
#############
# For a model of general form (here the rossler model)
# model dimension:
nVar = 3
# maximal polynomial degree
dMax = 2
# Number of parameter number (by default)
pMax <- d2pMax(nVar, dMax)
# convention used for the model formulation
poLabs(nVar, dMax)
# Definition of the Model Function
a = 0.520
b = 2
c = 4
Eq1 <- c(0,-1, 0,-1, 0, 0, 0, 0, 0, 0)
Eq2 <- c(0, 0, 0, a, 0, 0, 1, 0, 0, 0)
Eq3 <- c(b,-c, 0, 0, 0, 0, 0, 1, 0, 0)
K <- cbind(Eq1, Eq2, Eq3)
# Edition of the equations
visuEq(K, nVar, dMax)
# initial conditions
v0 <- c(-0.6, 0.6, 0.4)
# model integration
reconstr <- numicano(nVar, dMax, Istep=1000, onestep=1/50, KL=K,
                      v0=v0, method="ode45")
# Plot of the simulated time series obtained
dev.new()
plot(reconstr$reconstr[,2], reconstr$reconstr[,3], type='l',
      main='phase portrait', xlab='x(t)', ylab = 'y(t)')


#############
# Example 2 #
#############
# For a model of canonical form
# model dimension:
nVar = 4
# maximal polynomial degree
dMax = 3
# Number of parameter number (by default)
pMax <- d2pMax(nVar, dMax)
# Definition of the Model Function
PolyTerms <- c(281000, 0, 0, 0, -2275, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
               861, 0, 0, 0, -878300, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
# terms used in the model
poLabs(nVar, dMax, findIt=(PolyTerms!=0))
# initial conditions
v0 <- c(0.54, 3.76, -90, -5200)
# model integration
reconstr <- numicano(nVar, dMax, Istep=500, onestep=1/250, PolyTerms=PolyTerms,
                     v0=v0, method="ode45")
# Plot of the simulated time series obtained
plot(reconstr$reconstr[,2], reconstr$reconstr[,3], type='l',
     main='phase portrait', xlab='x', ylab = 'dx/dt')
# Edition of the equations
visuEq(reconstr$KL, nVar, dMax)





#############
# Example 3 #
#############
# For a model of general form (here the rossler model)
# model dimension:
nVar = 3
# maximal polynomial degree
dMax = 2
dMin = -1
# Number of parameter number (by default)
pMax <- regOrd(nVar, dMax, dMin)[2]
# convention used for the model formulation
poLabs(nVar, dMax, dMin)
# Definition of the Model Function
a = 0.520
b = 2
c = 4
Eq1 <- c(0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0)
Eq2 <- c(0, 0, 0, a, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0)
Eq3 <- c(b,-c, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0)
K <- cbind(Eq1, Eq2, Eq3)
# Edition of the equations
#visuEq(K, nVar, dMax)
# initial conditions
v0 <- c(-0.6, 0.6, 0.4)
# model integration
reconstr <- numicano(nVar, dMax, dMin, Istep=1000, onestep=1/50, KL=K,
                      v0=v0, method="ode45")
# Plot of the simulated time series obtained
dev.new()
plot(reconstr$reconstr[,2], reconstr$reconstr[,3], type='l',
      main='phase portrait', xlab='x(t)', ylab = 'y(t)')

GPoM documentation built on July 9, 2023, 6:23 p.m.