numiMultiX: Numerical Integration polynomial ODEs with Multiple eXternal...

In GPoM: Generalized Polynomial Modelling

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

GPoM documentation built on July 9, 2023, 6:23 p.m.