numicano: Numerical Integration of models in ODE of polynomial form

Description Usage Arguments Value Author(s) See Also Examples

View source: R/numicano.R

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
numicano(
  nVar,
  dMax,
  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.

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#############
# 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)')

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

## End(Not run)

GPoM documentation built on Feb. 18, 2020, 5:08 p.m.