• We have seen \(dy(t)/dt \approx y_{t+1} - y_t\). For discrete data points, the differentiation operator is discretized as matrix. For example, for \(3\) data points \(\mathbf{y}=[y_1, y_2, y_3]^T\), the backward difference matrix: \[\mathbf{D}=\left[\begin{array}{ccc} 1\\ -1 & 1\\ & -1 & 1\\ & & -1 \end{array}\right],\:\mathbf{D}\mathbf{y}=\left[\begin{array}{c} y_{1}\\ y_{2}-y_{1}\\ y_{3}-y_{2}\\ -y_{3} \end{array}\right]\]

• Forward difference matrix: \[\mathbf{D}^{T}=\left[\begin{array}{cccc} 1 & -1\\ & 1 & -1\\ & & 1 & -1 \end{array}\right],\:\mathbf{D}^{T}\mathbf{y}=\left[\begin{array}{c} y_{1}-y_{2}\\ y_{2}-y_{3}\\ y_{3}-y_{4} \end{array}\right]\]

• The symmetric matrix \(\mathbf{S}=\mathbf{D}^T \mathbf{D}\) is a second difference matrix: \[\mathbf{S}=\left[\begin{array}{ccc} 2 & -1\\ -1 & 2 & -1\\ & -1 & 2 \end{array}\right],\:\mathbf{S}\mathbf{y}=\left[\begin{array}{c} 2y_{1}-y_{2}\\ -y_{1}+2y_{2}-y_{3}\\ -y_{2}+2y_{3} \end{array}\right]\]

• The backward difference matrix for \(\mathbf{y}\), namely \(\mathbf{Dy}\) corresponds to \(dy/dt\).

• The forward difference matrix for \(\mathbf{y}\), namely \(\mathbf{D}^T\mathbf{y}\) corresponds to \(-dy/dt\), the dual operator of \(dy/dt\).

• The tri-diagonal second difference matrix for \(\mathbf{y}\), namely \(\mathbf{Sy}\) corresponds to second derivative \(- dy^2/dt^2\): \[- [(y_t - y_{t-1}) - (y_{t+1}- y_t)]= y_{t-1} - 2 y_{t} + y_{t+1} \approx \frac{\partial^2 y(t)}{\partial t^2}.\]

• The \(N\times N\) matrix \(\mathbf{S}\) is positive definite. \(\mathbf{S}\) corresponds to \(-d^2 /dt^2\) and \(-d^2 /dt^2 = (-d/dt)(d/dt)=\mathbf{D}^T \mathbf{D}\).

• Simple continuous states: \(s \in \mathbb{R}\).

• Continuous state combined with continuous time.

• Suppose \(y(t,s)\) is a continuous function of time \(t\in [0,\infty)\) and of a spatial variable \(s \in [0,\infty)\).

• We can take derivatives w.r.t. \(t\), \(\partial y(t,s)/\partial t\), and w.r.t. \(s\), \(\partial y(t,s)/\partial s\)

• Diffusion equation (a partial differential equation - PDE): 1st order change of \(t\), \(\partial y(t,s)/\partial t\) (time variation) is compensated by 2nd order change of \(s\), \(\partial^2 y(t,s)/\partial s^2\) (spatial variation): \[\frac{\partial y(t,s)}{\partial t} = \frac{\partial^2 y(t,s)}{\partial^2 s}\]

• Put the 1D spatial domain into a bounded interval and discretize it: \(s\in \{s_1,s_2,\dots, s_N \}\) where \(s_1<s_2\cdots<s_N\) for \(s_i \in \mathbb{R}\).

• \(y(t,s)\) becomes \(\mathbf{y}(t) = [y_1(t),y_2(t),\dots, y_N(t)]^T\) where \(y_i(t) = y(t,s_i)\).

• We know that \(-d^2 /dt^2\) can be discretized as \(N \times N\) tri-diagonal matrix \(\mathbf{S}\). So holding continuous \(t\), the diffusion equation is (roughly speaking) discretized as a system of \(N\) ODEs \[ \frac{\partial \mathbf{y}(t)}{\partial t} = - \mathbf{S} \mathbf{y}(t)\] where the \(i\)-th (\(1<i<N\)) function in \(\mathbf{y}\) is \[ \frac{d y_i(t)}{dt} = y_{i-1}(t) -2y_{i}(t) + y_{i+1}(t).\]

library(ReacTran) # Package for grid generation and discretization 
N = 100; sgrid = setup.grid.1D(x.down=1, x.up=0, N=N); s = sgrid$x.mid; 
head(s)

## [1] 0.005 0.015 0.025 0.035 0.045 0.055

Diffusion = function(t, y, parms){
  tran = tran.1D(C=y, D=1, v=0, dx=sgrid)
  list(dy = tran$dC)
}
yinit = sin(pi*s); times = seq(0, 0.2, by=0.01)
out = ode(y= yinit, times = times, func = Diffusion, parms =NULL)

• 1D diffusive-advective transport equation: \(D \frac{\partial^2 C}{\partial x^2} - v\frac{\partial C}{\partial x}.\)

plot(s, out[1, 2:(N+1)], type = "l", ylab = "y", xlab = "Distance s", col="grey") 
for (i in 2:20){lines(s, out[i, 2:(N+1)])}

image(out, grid=s)

- Diffusion with logistic growth: \[\frac{\partial y(t,x)}{\partial t} = D \frac{\partial^2 y(t,x)}{\partial x^2} + r y(t,x) (1-y(t,x))\] for \(0 <x <L\).

• On the boundary: \(\partial y(t,0)/\partial x =0\) and \(\partial y(t, L)/\partial x =0\)

• Initial condition: \(y(0,x)=f(x)\) for some function \(f(x)\).

• Discretize w.r.t. time \[\frac{\partial y(t,x)}{\partial t} = y(t+\Delta t, x) -y(t,x).\]

• Discretize w.r.t. space \[\frac{\partial^2 y(t,x)}{\partial x^2} = \frac{y( t, x-\Delta x) - 2 y(t,x) + y(t, x+\Delta x)}{(\Delta x)^2} \]

- The system of equations:\[\begin{align} y(t+1, 0) &= y(t, 0) + \frac{D \Delta t}{(\Delta x)^2} (-2y(t,0)+2y(t,1)) + r\Delta t y(t,0) (1- y(t,0)) \\ y(t+1, j) &= y(t, j) + \frac{D \Delta t}{(\Delta x)^2} (y(t,j-1)-2y(t,j)+y(t,j+1)) + r \Delta t y(t,j) (1- y(t,j)) \\ y(t+1, N+1) &= y(t, N+1) + \frac{D \Delta t}{(\Delta x)^2} (2y(t,N)-2y(t,N+1)) + r \Delta t y(t,N+1) (1- y(t,N+1)).\end{align}\]

D = 0.001; r=0.5; L=2; N=40; lambda = 0.3 # discretized factor
dx = L/(N+1); dt = lambda*(dx*dx)/D;
c1 = D*dt/(dx*dx); c2 = r*dt; 
yold = rep(0, N+2); ynew = rep(0, N+2);

# initial function f(x)
x = dx*c(0:(N+1));
for(j in 1:(N+2)){
  xx = x[j];
  if (xx<L/4) {
    yold[j] = 0;
  }
  else if (xx<L/2) {
    yold[j] = (xx-L/4)/(L/4);
  }
  else if (xx<3*L/4){
    yold[j] = (3*L/4 - xx)/(L/4);
  }
  else { yold[j] = 0}
}
plot(x, yold, ylab ="y(t,x)", xlab="0<x<2", type="l", col="blue")

tend = 20; t = dt; b=0.01
timeinv = seq(t,tend,by=dt)
for(t in timeinv){
  # 1st equation
  ynew[1] = yold[1] + c1*(-2*yold[1] + 2*yold[2]) + c2*yold[1]*(1-yold[1]);
  # 2nd equation
  ynew[2:(N+1)] = yold[2:(N+1)] + c1*(yold[1:N] - 2*yold[2:(N+1)] + yold[3:(N+2)])
  + c2*yold[2:(N+1)]*(1-yold[2:(N+1)]);
  # 3rd equation
  ynew[N+2] = yold[N+2] + c1*(2*yold[N+1] - 2*yold[N+2]) + c2*yold[N+2]*(1-yold[N+2]);
  # plot lines: morph from blue to red
  b=b+0.03; lines(x, ynew,  col=rgb(b, 0, 1-b, 0.5));
  yold = ynew;
}

• A simple logistic growth model (\(K=1\)) \[\frac{dy(t)}{dt} = r y(t) \left(1 - y(t) \right).\]

• Now consider two communities \(y_1(t)\) and \(y_2(t)\) (Turing’s morphogenesis): \[y_{1}(t)= r(y_{2})\left[y_{1}(1-y_{1})\right]=r(y_{2})f(y_{1})\\ y_{2}(t)= r(y_{1})\left[y_{2}(1-y_{2})\right]=r(y_{1})f(y_{2})\]

• If \(r(y_2)=r(y_1)=r\) constant growth rates for both, then these communities grow independently. Otherwise, there exists competitions or interactions.

• The system is similar to the prey-predator model.

• Modify the tri-diagonal second difference matrix \(\mathbf{S}\) to a birth-death matrix: \[\mathbf{F} = \left[\begin{array}{ccccc} 1 & \epsilon/2\\ & 1-\epsilon & \epsilon/2\\ & \epsilon/2 & 1-\epsilon & \ddots\\ & & \ddots & 1-\epsilon & \epsilon/2\\ & & & \epsilon/2 & 1-\epsilon/2 \end{array}\right]\]

• Simulate the coupled logistic equations induced by the birth-death matrix with \(N=10\), \(\epsilon=0.05\), \(y_5(0)=0.5\): \[\mathbf{y}(t) = f(\mathbf{F}^T\mathbf{y}(t))= \mathbf{F}^T\mathbf{y}(t) ( \mathbf{1} - \mathbf{F}^T\mathbf{y}(t)) \] with \(\mathbf{y}(t)=[y_1(t),\dots,y_{10}(t)]^T\). Simulate another one with \(\epsilon=0.5\).

• Finite differences

• Continuous states and discretization

• Diffusion

• Coupled ODEs