• Differentiation is for an infinitesimal (infinitely small) change.

• Differentiation w.r.t. time: \[ \frac{dy(t)}{dt} = \lim_{\epsilon \rightarrow 0}\frac{y(t)- y(t-\epsilon)}{\epsilon} \]

• Continuous time model. An (autonomous) ODE is often written as \[ \frac{dy}{dt}=f(y), \, f\mbox{ is a sufficiently differentiable map}.\]

• Deterministic dynamics describes how the state \(y\) evolves over time.

• For \(y(t) \in \mathbb{R}\), \(\mathbb{R}\) is called the phase space, \(y(\cdot):[0,T)\rightarrow\mathbb{R}\) is a trajectory (or an orbit). \(f(\cdot)\) is the vector field.

• Some necessary conditions, such as the initial condition of \(y(0)\) or the boundary condition of \(y(T)\) (termination) are needed for the solutions.

• For example, the linear growth model \[\frac{dy}{dt} = y\] has solution \(y(t)=y_0 e^t\) if we know \(y(0)=y_0\).

• An equilibrium point (or called steady state or stationary point) \(y^{*}\): For an ODE \(dy/dt=f(y)\), \(f(y^{*})=0\).

• The point \(y^{*} = 0\) in linear growth model is the equilibrium point.

• Consider a nonlinear function \(f\) in \(dy/dt=f(y)\): \[\frac{dy}{dt} = r y \left(1- \frac{y}{N}\right),\] the logistic equation, where \(r\) is the growth rate and \(N\) is carrying capacity.

• The growth rate can be re-scaled by change of time scaling: \(t=\tilde{t}/r\), then \(dy/d\tilde{t} = y (1- y/N)\).

• Equilibria are found by setting \(y(1-y/N)=0\) when \(y^*=0\) or \(y^*=N\).

• Stable or unstable: An equilibrium \(y^*\) is stable if there exists a neighbourhood \(\mathcal{U}\) of \(y^*\) such that all trajectories starting in \(\mathcal{U}\) converge to \(y^*\). Otherwise, \(y^*\) is unstable.

times = seq(0,100,by=1)
funct = function(t,integ,parms){
  y = integ[1]
  dy =  parms[1]*y*(1-y/parms[2])
  list(dy)}
require(deSolve)
yinit = 1; parms = c(0.1,10) # Also try yinit = 0.01; parms = c(1,10)
logit = ode(yinit,times,funct,parms) 
tail(logit, 2)

##        time        1
## [100,]   99 9.995486
## [101,]  100 9.995916

• \(y^*=0\) is unstable and \(y^*=N\) is stable.

require(rootSolve) 
equil=runsteady(y=yinit, times=c(0,100), func=funct, parms=parms) 
round(equil$y, 3)

## [1] 10

• Determine the stability: local analysis for the ODE \(dy/dt = f(y)\).

• Consider a perturbation around \(y^*\) \[y(t) = y^* + \epsilon X(t),\,\,\, \epsilon>0\] where \(\epsilon\) is sufficient small. Taylor expansion gives \[\begin{align} \frac{dy}{dt} &= \frac{d(y^* + \epsilon X(t)}{dt} = \epsilon \frac{d(X(t))}{dt} \\ &= f (y^* + \epsilon X(t)) \approx \epsilon \frac{\partial f(y^*)}{\partial y} X(t) + \cdots.\end{align}\]

• Previous linearized ODE \[\frac{dX(t)}{dt} = f'(y^*) X(t).\] It is known that \(f'(y^*)<0\) this ODE is convergent (stable) while \(f'(y^*)>0\) diverges (unstable).

• For the logistic growth \(f(y)=y(1-y/N)\), \(f'(y)=1-2y/N\).

• When \(y^*=0\), the linearized ODE is \(dX/dt = X\) so \(y^*=0\) is unstable.

• When \(y^*=N\), the linearized ODE is \(dX/dt = -X\) so \(y^*=N\) is stable.

• A system of \(N\) first-order ODEs: \[\frac{dy_{1}}{dt} =f_{1}(y_{1},\dots,y_{N}),\\ \vdots \\ \frac{dy_{N}}{dt} =f_{N}(y_{1},\dots,y_{N}).\]
• A compact form \[\frac{d\mathbf{y}}{dt}=\mathbf{f}(\mathbf{y})\] where \(\mathbf{y}=(y_{1},\dots,y_{N})^{T}\) and \(\mathbf{f}(\mathbf{y})=(f_{1}(\mathbf{y}),\dots,f_{N}(\mathbf{y}))^{T}\).

• The vector \(\mathbf{y}(t)\) is called state vector. If \(\mathbf{y}(t) \in \mathbb{R}^N\), then \(\mathbb{R}^N\) is the state space (or phase space) of the system.

• Similarly, the simplest ODE system \(d\mathbf{y}/dt=\mathbf{A}\mathbf{y}\) has the solution \[ \mathbf{y}(t) = e^{\mathbf{A}t}\mathbf{y}_0, \,\,\mbox{when }\mathbf{y}(t_0)=\mathbf{y}_0.\]

• For a general nonlinear vector function \(\mathbf{f}:\mathbb{R}^N \rightarrow \mathbb{R}^N\), the Taylor expansion linerization is also valid: \[\begin{align} \frac{d\mathbf{y}}{dt} &= \frac{d(\mathbf{y}^* + \mathbf{\epsilon}(t))}{dt} = \frac{d(\mathbf{\epsilon}(t))}{dt} \\ &= f (\mathbf{y}^* + \mathbf{\epsilon}) \approx \nabla f (\mathbf{y})|_{\mathbf{y}=\mathbf{y}^*} \times \mathbf{\epsilon} + \cdots.\end{align}\] where \(\nabla f(\mathbf{y})\) is the Jacobian matrix \[\nabla f (\mathbf{y})= \left[ \partial f (\mathbf{y})/\partial y_1, \dots, \partial f (\mathbf{y})/\partial y_N \right]^T = \left( \frac{\partial f_i}{\partial y_j} \right)_{ij}.\]

• Characterise a speed of divergence or convergence of initially nearby phase orbits.

• For the system \(d\mathbf{y}/dt=\mathbf{A}(t)\mathbf{y}\), given a perturbation \(\epsilon_0\), \(\epsilon_0\) evolves (gets large) according to \[\frac{d\mathbf{\epsilon}}{dt}=\mathbf{A}(t)\mathbf{\epsilon}(t).\] The solution at a time instant \(t_0 + \Delta t\) is \[\mathbf{\epsilon}(t_0 + \Delta t) = \exp \left( \int_{t_0}^{t_0 + \Delta t} \mathbf{A}(s)ds \right) \mathbf{\epsilon}(t_0) = \mathbf{M}(t_0,\Delta t)\mathbf{\epsilon}(t_0).\]

• Value assignment operators: <- and =. Traditionally (in R community), the operator <- are said to be used anywhere, whereas the operator = was only allowed at the top level or as one of the subexpressions in a braced list of expressions. However,…

xy <- list(x<-"x", y="y"); x # try y, you will receive an error.

## [1] "x"

xy$y

## [1] "y"

• Hint for beginners: Don’t make use of Global Variables!

• Notice that for any matrix \(\mathbf{A}\) we have \(\mathbf{A}\mathbf{y} = \lambda \mathbf{y}\) where \(\lambda\) is the eigenvalue of \(\mathbf{A}\). \(\lambda\) also satisfies the ODEs: \[\frac{de^{\lambda t}\mathbf{y}_0}{dt}=\mathbf{A}e^{\lambda t}\mathbf{y}_0 \mbox{ implies } \lambda e^{\lambda t}\mathbf{y}_0 = \mathbf{A}e^{\lambda t}\mathbf{y}_0.\] where \(\mathbf{y}_0\) now is the eigenvector.

• Suppose \(N \times N\) matrix \(\mathbf{A}\) has \(N\) independent eigenvectors. Then \(\mathbf{A}\) has \(N\) different eigenvalues. The eigenvectors \(\mathbf{x}_1,\dots, \mathbf{x}_N\) are a basis to express any \(\mathbf{y}_0\): \[\mathbf{y}_0 = c_1 \mathbf{x}_1 + \cdots + c_N \mathbf{x}_N.\]

• The solution has a superposition form \[\mathbf{y}(t) = c_1 e^{\lambda_1 t} \mathbf{x}_1 + \cdots + c_N e^{\lambda_1 t} \mathbf{x}_N.\]

• Note that eigenvalues can be complex \(\lambda = a + \mbox{i}b\).

• Then the real part of \(\lambda\) counts for growth or decay while the imaginary part counts for oscillation or rotation.

• Consider \[\frac{d\mathbf{y}}{dt}=\left[\begin{array}{cc} -2 & 1\\ -1 & -2 \end{array}\right]\mathbf{y},\:\:\mathbf{y}_{0}=\left[\begin{array}{c} 6\\ 2 \end{array}\right]\]

• Find the eigenvalues satisfying \[\det (\mathbf{A} - \lambda \mathbf{I}) = \lambda^2 + 4\lambda + 5 = 0.\] The roots are \(\lambda = (-4 \pm \sqrt{16-20})/2=-2 \pm \mbox{i}\). The eigenvectors are \((1, \pm \mbox{i})\).

• Subsitute into \(\mathbf{y}_0\), there are \(c_1 = 3-\mbox{i}\) and \(c_2 = 3+\mbox{i}\). The solution is \[\mathbf{y}(t) = (3-\mbox{i}) e^{(-2 + \mbox{i})t} \left[\begin{array}{c} 1\\ \mbox{i} \end{array}\right] + (3+\mbox{i}) e^{(-2 - \mbox{i})t}\left[\begin{array}{c} 1\\ -\mbox{i} \end{array}\right] \\ = e^{-2t}\left[\begin{array}{c} 6\cos t+2\sin t\\ 2\cos t-6\sin t \end{array}\right].\]

• For \(N=2\), the representative example of a conservative motion is traveling around a circle. That is, \(\|\mathbf{y}\|=\mathbf{y}^T \mathbf{y}\) stays constant. Thus \[\frac{d}{dt}(\mathbf{y}^T \mathbf{y}) = \mathbf{y}^T (\mathbf{A}^T + \mathbf{A}) \mathbf{y} = 0.\] It implies \(\mathbf{A}^T = -\mathbf{A}\) antisymmetric.

• For stale linear system, the solution \(\mathbf{y}(t)\) goes to zero which means the real part of the eigenvalue is negative \(\mbox{Re}\lambda < 0\).

• Previous example: \(\mbox{Re}\lambda = -2\), stable.

• Consider a second oder ODE: \[\frac{d^2 y}{dt^2} + B\frac{dy}{dt} + Cy = 0\]

• It has exactly the same solution as \[\frac{d\mathbf{y}}{dt}=\left[\begin{array}{c} dy/dt\\ d^{2}y/dt^{2} \end{array}\right]=\left[\begin{array}{cc} 0 & 1\\ -C & -B \end{array}\right]\left[\begin{array}{c} y\\ dy/dt \end{array}\right]=\mathbf{A}\mathbf{y}\]

• The 2nd order ODE: \(d^2y/dt^2 - y + y^3 =0\) can be written as \[\frac{dy}{dt} =x, \,\, \frac{dx}{dt} =y- y^3.\]

• If we introduce a function \(H(x,y)= x^2/2 + y^4/4 - y^2/2\), then we have \[\frac{dy}{dt} =\partial_x H(x,y),\,\, \frac{dx}{dt} = - \partial_y H(x,y).\] The function \(H:\mathbb{R}^2 \rightarrow \mathbb{R}\) is called the Hamiltonian. The above system is called a Hamiltonian system.

• Hamiltonian system: Any solution is constant along a given level set \(\{ H(x,y) = \mbox{constant} \}\).

• Proof: \(d H(x,y) = (dx/dt)\partial_x H + (dy/dt)\partial_y H =0\).

• Van der Pol equation (a self-excited system): \[\frac{d^{2}y}{dt^{2}}-\mu (1-y^{2}){dy \over dt}+y=0,\] with two interesting limits \(\mu \rightarrow 0\) and \(\mu =1/\sqrt{\epsilon} \rightarrow \infty\). It can be rewritten as \[\frac{d y_1}{dt} = y_2,\,\, \frac{d y_2}{dt} = \mu (1-y_1^2) y_2 - y_1.\]

• When \(\mu\) is small (fast motion), we rescale the time \(\tilde{t} = t/\sqrt{\epsilon}\) and \(\tilde{y}_2 = y_2/\sqrt{\epsilon}\) then \[ \frac{d \tilde{y}_2}{d\tilde{t}} = (1-y_1^2) \tilde{y}_2 - y_1, \,\, \frac{d y_1}{d\tilde{t}} = \epsilon \tilde{y}_2.\]

vdpol = function (t, y, mu) {
  list(c( y[2], mu * (1 - y[1]^2) * y[2] - y[1] ))
}
yini = c(y1 = 2, y2 = 0)
slow = ode(y = yini, func = vdpol, times = 0:3000, parms = 1000)
fast = ode(y = yini, func = vdpol, times = seq(0, 30, 0.01), parms = 1)
head(cbind(slow, fast))

##      time       y1            y2 time       y1          y2
## [1,]    0 2.000000  0.0000000000 0.00 2.000000  0.00000000
## [2,]    1 1.999333 -0.0006670373 0.01 1.999901 -0.01970323
## [3,]    2 1.998666 -0.0006674088 0.02 1.999608 -0.03882242
## [4,]    3 1.997998 -0.0006677807 0.03 1.999127 -0.05737271
## [5,]    4 1.997330 -0.0006681535 0.04 1.998462 -0.07537148
## [6,]    5 1.996662 -0.0006685269 0.05 1.997621 -0.09283422

plot(slow, fast, xlim=c(0,100))

equils=runsteady(y=yini, times=c(0,3000), func=vdpol, parms=1000) 
equilf=runsteady(y=yini, times=c(0,0.01), func=vdpol, parms=1) 
round(equils$y, 3)

##     y1     y2 
## -1.504  0.001

round(equilf$y, 3)

##     y1     y2 
##  2.000 -0.027

• Any set of autonomous ODEs can be reduced to a form, solved in respect of the highest derivatives.

• For the system \[\frac{dy_{1}}{dt} =f_{1}(y_{1},\dots,y_{N}), \dots\,, \frac{dy_{N}}{dt} =f_{N}(y_{1},\dots,y_{N})\] the form is \[\frac{dx_{1}}{dt} =x_{2},\, \frac{dx_{2}}{dt} =x_{3},\, \dots,\, \frac{dx_{N}}{dt} =F(x_{1},\dots,x_{N}).\]

• The general result is given by Takens’ theorem (delay embedding theorem).

• We’ve know a system of 1st order ODEs \(\mathbf{x}(t)\) is equal to a single high order ODE \(y(t)\).

• Suppose that we can measure \(y(t)\) and also we are not acquainted with the whole system of \(\mathbf{x}(t)\). Can we obtain the necessary information about the phase space of \(\mathbf{x}(t)\)?

• Answer: embedding, to show how the relevant variables describing the system that evolves in time. (i.e. we have a scalar series, we can construct a vector space whose axes represent all the (relevant) variables.)

• Lorenz system is a simple approximation of a general system. (It based upon a Fourier series approximation of the solution of 3 PDEs, namely Navier-Stokes equations, heat equation.)

• This approximation can be written as a three-dimensional system of ODEs \[\begin{align} \frac{dy_1}{dt} =& - c_1 y_1 + y_2 y_3\\ \frac{dy_2}{dt} =& - c_2 (y_2 - y_3)\\ \frac{dy_3}{dt} =& - - y_1y_2 + c_3 y_2 - y_3 \end{align}\] with parameters \(c_1 = 8/3, c_2 =10, c_3=28\).

chaos = function (t, y, parms) {
  dy1 = -8/3*y[1] + y[2]*y[3]
  dy2 = -10*(y[2]- y[3])
  dy3 = -y[1]*y[2] + 28*y[2] - y[3]
  list(c(dy1, dy2, dy3))}
yini = c(1,1,1); yini2 = yini + c(1e-6, 0, 0)
times = seq(0,100,0.01)
out  = ode(y = yini, times = times, func = chaos, parms = 0)
out2 = ode(y = yini2, times = times, func = chaos, parms = 0)
tail(cbind(out,out2), 3)

##            time        1          2         3   time        1        2
##  [9999,]  99.98 18.97496  -8.977816 -14.15177  99.98 29.63896 7.330461
## [10000,]  99.99 19.79639  -9.501998 -14.80335  99.99 29.14110 7.012169
## [10001,] 100.00 20.73260 -10.036203 -15.40840 100.00 28.63608 6.712828
##                 3
##  [9999,] 4.060228
## [10000,] 3.920508
## [10001,] 3.820715

• In nonlinear systems, a perturbation will either dampen or expand.

• A chaotic system: sensitive dependence on initial conditions. In other words, two very nearby trajectories will diverge exponentially over time.

plot(out, out2, xlim=c(0,30), mfrow=c(1,3))

library(scatterplot3d)
scatterplot3d(out[,2],out[,3],out[,4], type="l")

• Suppose we only observe \(y_3\).

• Use the delay embedding to construct a 8-dim vector \(\mathbf{x}= [y_3(t), y_3(t-1),\dots, y_3(t-8)]^T\).

x = embed(out[,4],8); head(x,3)

##          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]
## [1,] 3.109160 2.738552 2.400161 2.088545 1.798314 1.523999 1.259918
## [2,] 3.517577 3.109160 2.738552 2.400161 2.088545 1.798314 1.523999
## [3,] 3.969623 3.517577 3.109160 2.738552 2.400161 2.088545 1.798314
##          [,8]
## [1,] 1.000000
## [2,] 1.259918
## [3,] 1.523999

par(mfrow=c(1,2))
scatterplot3d(x[,1],x[,4],x[,8], type="l", angle = 70)
plot(x[,1],x[,8], type="l")

• Linear, non-linear and linearized non-linear ODEs

• System of ODEs

• Solutions, stable analysis

• Various features: harmonic, fast-slow, chaotic

• Embedding