numicano: 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"
)

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)

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

Author

Sylvain Mangiarotti

Examples

Run this code

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

# \donttest{
#############
# 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)
# }



# \donttest{
#############
# 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)')

# }