Learn R Programming

GPoM (version 1.0)

numicano: Numerical Integration of models in ODE of polynomial form

Description

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

Usage

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

Arguments

nVar

The model dimension expected. This parameter will be deduced from the input data (series) if series is a matrix. If series is a vector, the expected dimension nVar should be provided.

dMax

Maximum degree of the polynomial functions allowed in the model (see poLabs).

Istep

Number of integration time steps

onestep

Time step length

KL

Matricial formulation of the model to integrate numerically

PolyTerms

Vectorial formulation of the model (only for models of canonical form)

v0

Initial conditions (a vector which length should correspond to the model dimension nVar)

method

Integration method (See deSolve), by default method = 'ode45'.

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)

Examples

Run this code
# NOT RUN {
#############
# 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(nVar, dMax, K)
# 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)')

# }
# NOT RUN {
#############
# 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, 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(nVar, dMax, reconstr$KL)
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab