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

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)

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

Author

Sylvain Mangiarotti

Examples

Run this code

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