Nonlinear Galerkin methods for a system of PDEs with Turing instabilities

Abstract

We address and discuss the application of nonlinear Galerkin methods for the model reduction and numerical solution of partial differential equations (PDE) with Turing instabilities in comparison with standard (linear) Galerkin methods. The model considered is a system of PDEs modelling the pattern formation in vegetation dynamics. In particular, by constructing the approximate inertial manifold on the basis of the spectral decomposition of the solution, we implement the so-called Euler–Galerkin method and we compare its efficiency and accuracy versus the linear Galerkin methods. We compare the efficiency of the methods by (a) the accuracy of the computed bifurcation points, and, (b) by the computation of the Hausdorff distance between the limit sets obtained by the Galerkin methods and the ones obtained with a reference finite difference scheme. The efficiency with respect to the required CPU time is also accessed. For our illustrations we used three different ODE time integrators, from the Matlab ODE suite. Our results indicate that the performance of the Euler–Galerkin method is superior compared to the linear Galerkin method when either explicit or linearly implicit time integration scheme are adopted. For the particular problem considered, we found that the dimension of approximate inertial manifold is strongly affected by the lenght of the spatial domain. Indeeed, we show that the number of modes required to accurately describe the long time Turing pattern forming solutions increases as the domain increases.

Appendix

1.1 Computation of the Hausdorff distance

We present a simple algorithm for parallel calculation of Hausdorff distance. The code was implemented in matlab environment with the use of a processor intel XEON E52630 [from Institute for Research on Combustion (IRC)]. Matlab(2016b) parallelization use 12 cores under the command parfor. The serial implementation of the code takes in CPU time about 86.83 sec, while the simple parallel code (with the use of parfor loop) takes in mean value 8.95 sec (almost 10 times faster as it was expected).

The Hausdorff distance is defined as

$$\begin{aligned} d_{{{\mathrm H}}}(X,Y)=\max \{\,\sup _{{x\in X}}dist(x,Y),\,\sup _{{y\in Y}}dist(y,X)\,\}{\text{, }}\! \end{aligned}$$

where $dist(x,Y)=\inf _{{y\in Y}}d(x,y)$. The serial code takes both vector spaces X, Y (one from FD and the other from LGM or EGM and for each $x\in X$, calculates the dist(x, Y), then determines the supremum (and repeats similar, by interchange the X and Y spaces).

We implement the algorithm 1 in a parallel code, partitioning the X space into $N_{cores}$ samples where $N_{cores}$ is the number of cores (i.e.$N_{cores}=12$). Then, for each sample ($X_{sample}$), calculates parallel the $dist(x,Y),x \in X_{sample}$. The partition of X space is implemented according to algorithm 1. The reported CPU times are calculated for matrices (the X and Y vector spaces) of size dimension $5001\times 123$ double.

1.2 The finite difference scheme

Let $N \in \mathbb {N}$, we discretize the domain [0, L] into N equal interval i.e. $h=L/N$ and $x_i=i\cdot h, i=0,1,2\ldots ,N$. Also, we define $B_i^n: = B\left( {{x_i},{t^n}} \right)$ to be the value of biomass B, in place $x_i$ and time $t^n$ (similar we define $W_i^n: = W\left( {{x_i},{t^n}} \right)$ and $T_i^n: = B\left( {{x_i},{t^n}} \right)$ ). The FD scheme for the laplacian operator is

$$\begin{aligned} {\Delta _{{h}}}u\left( {{x_i},{t^n}} \right) :&= \frac{{u\left( {{x_i} + {h_1},{t^n}} \right) - 2u\left( {{x_i},{t^n}} \right) - u\left( {{x_i} - {h_1},{t^n}} \right) }}{{h^2}} \\&=\frac{{u_{i + 1}^n - 2u_i^n + u_{i - 1}^n}}{{h^2}} \end{aligned}$$

where u represents both B, W, T. The mathematical model (3) for each node $x_i$ take the form:

$$\begin{aligned} \left\{ \begin{array}{l} \vartheta _t B_{i}^{n} = {c_B}\left( {B_{i + 1}^n - 2B_i^n + B_{i - 1}^n} \right) + c{\left( {B_i^n} \right) ^2} W_i^n - \left( {d + sT_i^n} \right) B_i^n\\ \vartheta _t W_i^n= {c_W}\left( {W_{i + 1}^n - 2W_i^n + W_{i - 1}^n} \right) + p - r{{\left( {B_i^n} \right) }^2}W_i^n - lW_i^n\\ \vartheta _t T_i^n= {q\left( {d + sT_i^n} \right) {B_i^{n}} - \left( {k + wp} \right) T_i^n} \end{array} \right. \end{aligned}$$

(34)

where the constants are $c_B=\frac{D_B}{h^2}$, $c_W=\frac{D_W}{h^2}$ and $i=1,2,\ldots ,N-1$

For the boundaries points, we implement central differences. For example in the left boundary $x=0,(i=0)$, (34) takes the form

$$\begin{aligned} \left\{ \begin{array}{l} \vartheta _t B_{0}^{n} = {c_B}\left( {2B_{1}^n - 2B_0^n} \right) + c{\left( {B_0^n} \right) ^2} W_0^n - \left( {d + sT_0^n} \right) B_0^n\\ \vartheta _t W_0^n= {c_W}\left( {2W_{1}^n - 2W_0^n} \right) + p - r{{\left( {B_0^n} \right) }^2}W_0^n - lW_0^n\\ \vartheta _t T_0^n= {q\left( {d + sT_0^n} \right) {B_0^{n}} - \left( {k + wp} \right) T_0^n} \end{array} \right. \end{aligned}$$

(35)

and similar for the right boundary $x=L,(i=N)$

$$\begin{aligned} \left\{ \begin{array}{l} \vartheta _t B_{N}^{n} = {c_B}\left( {2B_{N-1}^n - 2B_{N}^n} \right) + c{\left( {B_{N}^n} \right) ^2} W_{N}^n - \left( {d + sT_{N}^n} \right) B_{N}^n\\ \vartheta _t W_{N}^n= {c_W}\left( {2W_{N-1}^n - 2W_{N}^n} \right) + p - r{{\left( {B_{N}^n} \right) }^2}W_{N}^n - lW_{N}^n \\ \vartheta _t T_{N}^n= {q\left( {d + sT_{N}^n} \right) {B_{N}^{n}} - \left( {k + wp} \right) T_{N}^n} \end{array} \right. \end{aligned}$$

(36)

Equation (34–36) define an ordinary differential system of the form

$$\begin{aligned} \frac{d\mathbf f }{dt}=\mathbf f (t,\mathbf u ) \end{aligned}$$

(37)

where $\mathbf u ={[{B_0},{B_1},\ldots {B_{N }},{W_0},{W_1},\ldots {W_{N}},{T_0},{T_1},\ldots {T_{N}}]^T}$. The Jacobian of $\mathbf f$ (the RHS) defines a square block matrix of dimension $3N+3$.

$$\begin{aligned} J = \left( {\begin{array}{*{20}{c}} {{A_{11}}}&{}\quad {{A_{12}}}&{}\quad {{A_{13}}}\\ {{A_{21}}}&{}\quad {{A_{22}}}&{}\quad {{A_{23}}}\\ {{A_{31}}}&{}\quad {{A_{32}}}&{}\quad {{A_{33}}} \end{array}} \right) \end{aligned}$$

Each diagonal block $i.e. A_{i,i}$ has a tridiagonal form and the others matrices ${A_{ij}},i \ne j$ are diagonal. For example: $A_{11}$=

$$\begin{aligned} \begin{bmatrix} {- 2{c_B} + 2{c_B}{B_0}{W_0} - \left( {d + s{T_0}} \right) }&\quad {2c_B}&\quad 0&\quad \dots&\quad 0 \\ {c_B}&\quad { - 2{c_B} + 2{C_B}{B_1}{W_1} - \left( {d + s{T_1}} \right) }&\quad c_B&\quad \dots&\quad 0 \\ \vdots&\vdots&\vdots&\ddots&\vdots \\ 0&\quad 0&\quad 0&\quad \dots&\quad - 2{c_B} + 2{c_B}{B_{N }}{W_{N}} - {d + s{T_{N}}} \end{bmatrix} \end{aligned}$$

${A_{12}} = diag\left( {cB_\iota ^2} \right)$ and ${A_{13}} = diag\left( {s{B_i}} \right) ,i = 0,1,\ldots ,N$