--- title: "GPoM : 1 Conventions" author: "Sylvain Mangiarotti & Mireille Huc" date: "`r Sys.Date()`" output: pdf_document: number_sections: no vignette: > %\VignetteIndexEntry{GPoM : 1 Conventions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} load package: "`r library(GPoM)`" --- The Generalized Global Polynomial Modelling (GPoM) package allows a generic formulation of any Ordinary Differential Equations (ODEs) in polynomial form. The aim of the present vignette `1 Conventions` is to introduce briefly the way to describe a set of polynomial ODEs with GPoM and to show how to perform its numerical integration^[Mangiarotti S, Le Jean F, Chassan M, Drapeau L, Huc M. 2018. GPoM: Generalized Polynomial Modelling. Version 1.1. *Comprehensive R Archive Network*. .]. ## Conventions used to describe a polynomial The polynomial description is based on a convention defined by the function `regOrd` that provides the order of the polynomial terms. This convention depends on the model dimension (that is the number `nVar` of state variables), and on the maximum polynomial degree used for the formulation (defined by the parameter `dMax`). This order can be visualized using the `poLabs` function. For instance, for `nVar = 3` and `dMax = 2`, the convention used to formulate a polynomial will be: ```{r} nVar = 3 dMax = 2 poLabs(nVar = nVar, dMax = dMax) ``` This formulation has `pMax = 10` terms: ```{r, echo=FALSE} nVar = 3 dMax = 2 ``` ```{r} pMax = d2pMax(nVar, dMax) ``` Based on this convention, one single ordinary differential equation (ODE) in polynomial form with `nVar` variables and of maximum polynomial degree `dMax` can be formulated as one single vector using the convention given by `regOrd(nVar, dMax)`. As an example, the equation : $dX_1/dt = 1 + 2 X_1 - 3 X_1X_3 + 4 X_2^2$ has the three variables $X_1$, $X_2$ and $X_3$ (it thus requires at least `nVar = 3`) and is of maximum polynomial degree two due to terms $X_1X_3$ and $X_2^2$ (it thus requires at least `dMax = 2`). Following the convention defined by `poLabs(nVar = 3, dMax = 2)`, it will require the definition of the following vector of parameters: ```{r} param <- c(1, 0, 0, 0, 0, 4, 2, -3, 0, 0) ``` Indeed: ```{r} nVar = 3 dMax = 2 cbind(param, poLabs(nVar, dMax)) ``` The same convention will be used for any other equation. Note that, by default, the notation used for the variables is `X1`, `X2`, etc. However, to facilitate the analysis, alternative notations may also be used using the optional parameter `Xnote`: ```{r} poLabs(3, 2, Xnote = 'y') ``` or for a full choice of the notation: ```{r} poLabs(3, 2, Xnote = c('x','W','y')) ``` ## Definition of a set of polynomial ODE A set of $N$ equations will require the definition of $N$ parameter vectors and will thus be represented by a matrix of `pMax` lines by `nVar` columns. For example, the Rössler system^[O. Rössler, An Equation for Continuous Chaos, Physics Letters, **57A**(5), 1976, 397-398] is defined by a set of three equations $dx/dt = - y - z$ $dy/dt = x + a y$ $dz/dt = b + z (x - c)$. For $(a = 0.52, b = 2, c = 4)$, this system can be decribed by three vectors (one for each equation) ```{r} # parameters a = 0.52 b = 2 c = 4 # equations 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) ``` The model formulation is obtained by concatenating the vectors of the three equations into one single matrix `K` containing all the coefficients of the model: ```{r} K = cbind(Eq1, Eq2, Eq3) ``` The corresponding model equations can be edited in a mathematical form using the function `visuEq()` ```{r} visuEq(K) ``` By default, the notation used in `visuEq` for the variables is `X` with an indicative number, such as `X1`, `X2`, etc. Alternative notation can be used to edit the equations with the `visuEq` function using the optional parameter `substit`. For `substit = 1`, single letters are automatically chosen such as ```{r} visuEq(K, substit = 1) ``` The notation can be defined also manually, such as: ```{r, eval=TRUE} visuEq(K, substit = c("U", "V", "W")) ``` ## Numerical integration The numerical integration of the model defined by matrix `K` can be done using the `numicano` function. It requires the use of the external package `deSolve`. The following parameters are required as input: ```{r} # The initial conditions of the system variables v0 <- c(-0.6, 0.6, 0.4) # the model formulation K (see former section) # the number of integration steps `Istep` nIstep <- 5000 # the time step length `onestep` onestep = 1/50 # the model dimension `nVar` nVar = 3 # the maximum polynomial degree `dMax` dMax = 2 ``` The numerical integration is launched as follows: ```{r} outNumi <- numicano(nVar, dMax, Istep = nIstep, onestep = onestep, KL = K, v0 = v0) ``` The output of the function `numicano` is a list that contains (1) a memory `$KL` of the model parameters ```{r, eval=FALSE} outNumi$KL ``` from which `nVar` and `dMax` (required to reformulate the equations) can be retrieved ```{r, eval = FALSE} # nVar dim(outNumi$K)[2] # dMax pMax <- dim(outNumi$K)[1] p2dMax(nVar, pMaxKnown = pMax) ``` and (2) the simulations `$reconstr`. This matrix has `nVar + 1` columns. The first one is the time, the other ones correspond to the variables of the system $(X1, X2, X3, ...)$ (or $(x,y,z)$ to keep the formulation used previously in the text). Note that all the other input parameters used in `numicano` can be retrieved from the outputs: ```{r, eval=FALSE} # initial conditions head(outNumi$reconstr, 1)[2:(nVar+1)] # time step diff(outNumi$reconstr[1:2,1]) # number of integration time step dim(outNumi$reconstr)[1] ``` The simulated time series can be plotted as follows: ```{r, fig.show='hold', fig.align='center'} plot(outNumi$reconstr[,1], outNumi$reconstr[,2], type='l', main='time series', xlab='t', ylab = 'x(t)') ``` and the plot of the phase portrait as well: ```{r, fig.show='hold', fig.align='center', fig.width=4, fig.height=4} plot(outNumi$reconstr[,2], outNumi$reconstr[,3], type='l', main='phase portrait', xlab='x(t)', ylab = 'y(t)') ``` ## Next steps The aim of the GPoM package is to retreive ODEs from time series using global modelling. Such type of modelling may require a careful data preprocessing. Simple examples of preprocessing will be given in the next vignette `2 Preprocessing`. Examples for applying global modelling to time series will then be presented in vignette `3 Modelling`.