numiMultiX {GPoM}

Numerical Integration polynomial ODEs with Multiple eXternal forcing

Description

Function for the numerical integration of Ordinary Differential Equations of polynomial form including single or Multiple external forcing

Usage

numiMultiX(
  nVar,
  dMax,
  Istep = 1000,
  onestep = 1/125,
  KDf,
  extF = extF,
  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. By default, Istep = 1000

onestep

The time step to be used for numerical integration

KDf

The nonautonomous model in its matrix formulation, NA (i.e. not available) values should be provided for forcing variables provided as an external signal

extF

A matrix providing the time vector in the first column, and time series of each forcing in the next ones

v0

The initial conditions. Its length should be in agreement with the dynamical system dimension. Therefore, 0 or NA can be provided for external forcing

method

integration method. By default 'rk4' is used

Value

A list of two variables:

$KDf The nonautonomous 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, numicano, numinoisy

Examples

#############
# Example 1 #
#############
# build a non autonomous model
nVar = 4
dMax = 3
gamma = 0.05
KDf=matrix(0, nrow = d2pMax(nVar = nVar, dMax = dMax), ncol = nVar)
KDf[11,1]  = 1
KDf[2,2]  = 1
KDf[5,2]  = 1
KDf[11,2]  = -gamma
KDf[35,2] = -1
KDf[2,3]  = NA
KDf[2,4]  = NA
visuEq(K = KDf, substit = c('x', 'y', 'u', 'v'))

# build an external forcing
# number of integration time step
Istep <- 500

# time step
smpl <- 1 / 20

# output time vector
tvec <- (0:(Istep-1)) * smpl

# angular frequency (for periodic forcing)
omega = 0.2

# half step time vector (for Runge-Kutta integration)
tvecX <- (0:(Istep*2-2)) * smpl / 2
# generate the forcing (here variables u and v)
extF = cbind(tvecX, -0.1 * cos(tvecX * omega), 0.05 * cos(tvecX * 16/3*omega))

# decimate the data
extFrs <- extF[seq(1,dim(extF)[1],by=50),]
extFrs <- rbind(extFrs,extF[dim(extF)[1],])

# Initial conditions to be used (external variables can be set to 0)
etatInit <- c(-0.616109362 , -0.126882584 , NA, NA)
# model integration
out <- numiMultiX(nVar, dMax, Istep=Istep, onestep=smpl, KDf=KDf,
                       extF,
                       v0=etatInit, method="rk4")
outrs <- numiMultiX(nVar, dMax, Istep=Istep, onestep=smpl, KDf=KDf,
                       extFrs,
                       v0=etatInit, method="rk4")
dev.new
oldpar <- par(no.readonly = TRUE)    
on.exit(par(oldpar))  
par(mfrow = c(2, 2), # 2 x 2 pictures on one plot
    pty = "s")
plot(out$reconstr[,2],out$reconstr[,3],
    xlab = 'x(t)', ylab = 'y(t)', type = 'l', col = 'red')
lines(outrs$reconstr[,2],outrs$reconstr[,3],
    xlab = 'x(t)', ylab = 'y(t)', type = 'l', col = 'green')
plot(out$reconstr[,2],out$reconstr[,4],
    xlab = 'x(t)', ylab = 'u(t)', type = 'l', col = 'red')
plot(out$reconstr[,4],out$reconstr[,5],
    xlab = 'u(t)', ylab = 'v(t)', type = 'l', col = 'red')