Reaction-diffusion systems, invariants and pattern computation

Let $m\ge 1$, let $U\subseteq {\mathbb{R}}^m$ be open, and let $F:U\to {\mathbb{R}}^m$ be locally Lipschitz. Write $F_i$ for its $i$th component. The autonomous reaction equation is $z^{\mathrm{\prime }}(t)=F(z(t))$. The field is quasi-positive on $U\cap [0,\mathrm{\infty }{)}^m$ if, for every $i$ and every $z\in U\cap [0,\mathrm{\infty }{)}^m$ with $z_i=0$, one has $F_i(z)\ge 0$. This condition concerns only the coordinate faces of the nonnegative orthant.

For a rectangle $R={\prod }_{i=1}^m[a_i,b_i]\subset U$ with $a_i<b_i$, the field is inward-pointing on $R$ if $F_i(z)\ge 0$ whenever $z\in R$ and $z_i=a_i$, and $F_i(z)\le 0$ whenever $z\in R$ and $z_i=b_i$. If $a_i=0$ for every $i$, quasi-positivity supplies the lower-face half of this condition, but the upper-face inequalities remain separate hypotheses. Neither boundary condition by itself asserts that a numerical update preserves $R$.

For smoothness-dependent results below, $F\in C^4(U,{\mathbb{R}}^m)$ will be stated explicitly. No smoothness beyond local Lipschitz continuity is built into this definition.

Fix $d\ge 1$, side lengths $L_j>0$, and integers $N_j\ge 3$. Put $h_j=L_j\mathrm{/}{N}_j$, ${\mathbb{T}}_L^d={\prod }_{j=1}^d(\mathbb{R}\mathrm{/}{L}_j\mathbb{Z})$, and $G_h={\prod }_{j=1}^d(\mathbb{Z}\mathrm{/}{N}_j\mathbb{Z})$. Grid nodes are $x_i=((i_j+\frac{1}{2})h_j{)}_{j=1}^d$ modulo the side lengths, grid indices wrap modulo $N_j$, and $P_{h}u=(u(x_i){)}_{i\in G_h}$ denotes sampling. The cell volume is $m_h={\prod }_j{h}_j$. The positive periodic heat stencil is the unit-diffusivity specialization of the periodic operator in Fully discrete θ heat finite-difference scheme:

For $s$ species with diffusivities $d_r\ge 0$, the reaction-diffusion grid generator acts on species $r$ by $-d_r{A}_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$. Each grid node has weight $m_h={\prod }_j{h}_j$; the stencil has zero row and column sums, so each species' weighted sum is conserved by diffusion. Its matrix representation satisfies $L_h=m_h{A}_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$, with symmetric $L_h$, nonpositive off-diagonal entries and $L_h\mathbf{1}=0$.

The sign convention is that $A_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$ is nonnegative and $-d_r{A}_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$ generates diffusion. The exact eigenvalues of the positive stencil are supplied by Periodic heat-stencil eigenvalues and θ amplification factor. This grid is the only spatial discretization used in the periodic error theorem; no mesh-uniform claim for general finite-volume meshes is made there.

Let $F:U\to {\mathbb{R}}^m$ be locally Lipschitz, and let $R={\prod }_{i=1}^m[a_i,b_i]\subset U$ with finite $a_i<b_i$. The rectangle is reaction-inward if, for every $i$, $F_i(z)\ge 0$ on the lower face $z_i=a_i$ and $F_i(z)\le 0$ on the upper face $z_i=b_i$, with all other coordinates in their intervals. It is reaction-invariant if the solution map ${\mathrm{\Phi }}_F^t$ is defined for all $t\ge 0$ on $R$ and ${\mathrm{\Phi }}_F^t(R)\subseteq R$ for every $t\ge 0$.

For a rectangle with $a_i=0$ for all $i$, quasi-positivity from Reaction vector fields and quasi-positivity implies the lower-face condition. That condition tests only points in $R$ and does not imply quasi-positivity on all of $U\cap [0,\mathrm{\infty }{)}^m$. The upper-face condition is independent and must be checked separately. The inward-sign condition is a hypothesis, not by definition an invariance theorem.

Let $F:U\to {\mathbb{R}}^m$ be the reaction vector field of Reaction vector fields and quasi-positivity, and let $R\subset U$ be a state set. A vector $c\in [0,\mathrm{\infty }{)}^m$ is a stoichiometric conservation vector on $R$ if $c^{T}F(z)=0$ for every $z\in R$. It is a stoichiometric dissipation vector on $R$ if $c^{T}F(z)\le 0$ for every $z\in R$. The corresponding pointwise mass is $c^{T}z$.

For a Lebesgue measurable domain $\mathrm{\Omega }\subseteq {\mathbb{R}}^d$ and a measurable field $u:\mathrm{\Omega }\to {\mathbb{R}}^m$ such that ${\int }_{\mathrm{\Omega }}\mathrm{\mid }{c}^{T}u(x)\mathrm{\mid }dx<\mathrm{\infty }$, the continuous total mass is ${\mathcal{M}}_c(u)={\int }_{\mathrm{\Omega }}{c}^{T}u(x)dx$, with Lebesgue integration. In particular, this is defined for bounded measurable fields on finite-volume domains. For any finite index set $I$ with assigned positive weights $m_i$, the discrete total mass of states $U_i\in {\mathbb{R}}^m$ is ${\mathcal{M}}_{c,h}(U)={\sum }_{i\in I}{m}_i{c}^T{U}_i$. The weights and index set are part of the discrete model. Conservation and dissipation always refer to these stated functionals and the stated state set; no claim is made for a different weight vector or for reactions outside $R$.

Let $m\ge 1$, let $D=\mathrm{diag}(d_1,\dots ,d_m)$ with $d_i\ge 0$, and let $F:U\to {\mathbb{R}}^m$ be as in Reaction vector fields and quasi-positivity. A semilinear vector reaction-diffusion system is

on either a flat periodic torus ${\mathbb{T}}_L^d$ or a bounded connected polyhedral domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$ with homogeneous no-flux conditions $d_i{\mathrm{\partial }}_n{u}_i=0$ on each boundary face. Equivalently, ${\mathrm{\partial }}_n{u}_i=0$ is required for components with $d_i>0$; components with $d_i=0$ require no boundary condition. The diffusion matrix is diagonal and independent of $u$; this definition does not include cross-diffusion.

A classical solution is bounded in the rectangle $R$ on $[0,T]$ if $u(x,t)\in R$ at every point and time in its stated solution class, with $R$ as in Inward reaction rectangles. This item fixes the continuous PDE convention; the finite-volume and periodic discrete operators are defined separately.

Let $\mathcal{T}=(\mathcal{K},\mathcal{F})$ be a finite TPFA-admissible cell-face mesh from Finite-volume cell-face meshes and geometric regularity and Orthogonal TPFA diffusion with harmonic face transmissibility. Impose homogeneous no-flux boundary conditions: boundary faces contribute zero flux and therefore no conductance. For adjacent cells $K,L$, set

where the sum is over their common interior faces and is zero when they are not adjacent. Set $c_{KL}=c_{LK}$ and $m_K=\mathrm{\mid }K\mathrm{\mid }$. If more than one face separates a cell pair, summing their conductances gives the same cell-balance flux as retaining the faces separately. Define the symmetric matrix $L_h$ by $(L_h{)}_{KK}={\sum }_{L\mathrm{\ne }K}{c}_{KL}$ and $(L_h{)}_{KL}=-c_{KL}$ for $K\mathrm{\ne }L$.

For species diffusivities $d_r\ge 0$, the no-flux finite-volume diffusion equation is $M{U}_r^{\mathrm{\prime }}+d_r{L}_h{U}_r=0$ for each species $r$, where $M=\mathrm{diag}(\mathrm{\mid }K\mathrm{\mid })$. Its component equation is

The equality is the usual shared-face TPFA balance with outward flux $T_F(U_{K,r}-U_{L,r})$ and no boundary flux. The matrix $L_h$ is symmetric, has zero row sums and nonpositive off-diagonal entries, and obeys $v^{\mathsf{T}}{L}_{h}v={\sum }_{\{K,L\}}{c}_{KL}(v_K-v_L{)}^2\ge 0$, where the sum is over unordered adjacent cell pairs. With $M=\mathrm{diag}(\mathrm{\mid }K\mathrm{\mid })$, the generator $-M^{-1}{L}_h$ is self-adjoint in the weighted pairing $⟨v,w{⟩}_M=v^{\mathsf{T}}Mw$. This definition makes no accuracy claim for arbitrary admissible meshes.

Let $R={\prod }_{i=1}^m[a_i,b_i]$ be the reaction-inward rectangle of Inward reaction rectangles. Assume $F$ from Reaction vector fields and quasi-positivity is continuously differentiable on an open neighborhood of $R$, and its derivative is bounded on some closed neighborhood of $R$. For any $z_0\in R$, the IVP $z^{\mathrm{\prime }}=F(z)$ has a unique solution for every $t\ge 0$, and $z(t)\in R$. If $c$ is a conservation vector on $R$, then $c^{T}z(t)=c^T{z}_0$; if it is a dissipation vector, then $c^{T}z(t)$ is nonincreasing. When every $a_i=0$, quasi-positivity implies the lower-face condition on $R$; upper-face inwardness remains a separate hypothesis.

Let $R={\prod }_{i=1}^m[a_i,b_i]$ be reaction-inward for $F$ from Reaction vector fields and quasi-positivity. Assume for each $i$ there is a coordinate Lipschitz constant $L_i\ge 0$ such that $\mathrm{\mid }{F}_i(z)-F_i(w)\mathrm{\mid }\le L_i\mathrm{\mid }{z}_i-w_i\mathrm{\mid }$ whenever $z,w\in R$ differ only in coordinate $i$. Define

with ${\tau }_R=+\mathrm{\infty }$ if every $L_i=0$. For $z\in R$ and $0\le \tau \le {\tau }_R$, the forward-Euler reaction step $z^{+}=z+\tau F(z)$ lies in $R$. For every conservation vector $c$ on $R$, $c^T{z}^{+}=c^{T}z$; for every dissipation vector, $c^T{z}^{+}\le c^{T}z$.

Let $\mathrm{\Omega }={\prod }_{j=1}^d(0,L_j)$, where $L_j>0$, and let integers $N_j\ge 2$ define $h_j=L_j\mathrm{/}{N}_j$. The finite cell set is $I_h={\prod }_{j=1}^d\{0,\dots ,N_j-1\}$, cell centers are $x_i=((i_j+\frac{1}{2})h_j{)}_{j=1}^d$, and every cell has volume $m_h={\prod }_j{h}_j$. This is the uniform orthogonal Cartesian mesh; adjacent center segments are perpendicular to their shared faces. Take the scalar cellwise coefficient to be ${\lambda }_K=1$ and give boundary faces homogeneous no-flux data, so they contribute no flux.

For a grid vector $v=(v_i{)}_{i\in I_h}$, define the positive cell-centered stencil

where an unavailable neighbor at a boundary is assigned the value of the boundary cell itself. Indeed a face normal to direction $j$ has area $\mathrm{\mid }F\mathrm{\mid }=m_h\mathrm{/}{h}_j$ and adjacent centers are distance $h_j$ apart, so its TPFA transmissibility is $\mathrm{\mid }F\mathrm{\mid }\mathrm{/}{h}_j=m_h\mathrm{/}{h}_j^2$; dividing the cell balance by $m_h$ gives the displayed coefficient. Boundary faces contribute zero flux. Thus $M_h=m_{h}I$ and $L_h^N=m_h{A}_h^N$ are the mass and stiffness matrices of the no-flux specialization in No-flux TPFA diffusion generator. For species diffusivities $d_r\ge 0$, the finite-volume diffusion equation is $U_r^{\mathrm{\prime }}=-d_r{A}_h^N{U}_r$.

Define the doubled periodic grid with $2N_j$ cells per direction on side lengths $2L_j$, using the same widths $h_j$, as in Periodic rectangular reaction-diffusion grid. The reflection map $E_h$ extends $U$ evenly across every boundary face: in one coordinate, its extended index is $i$ for $0\le i<N_j$ and $2N_j-1-i$ for $N_j\le i<2N_j$; apply this coordinatewise. This finite reflection convention is used in the error theorem below.

Let $u\in C^1([0,T];C^2(\bar{\mathrm{\Omega }};{\mathbb{R}}^m))$ solve $u_t=D\mathrm{\Delta }u+F(u)$, with constant diagonal $D=\mathrm{diag}(d_i)$, and let its range lie in $R$. Let $c$ be a conservation or dissipation vector on $R$. One may take (i) a box $\mathrm{\Omega }={\prod }_{j=1}^d[0,L_j]$ with periodic identification or homogeneous Neumann data, or (ii) a bounded domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$, $d\ge 2$, with a specified finite piecewise $C^1$ boundary presentation as in Specified finite piecewise C1 boundary presentations, assuming ${\mathrm{A}\mathrm{C}}_{\omega }$ in this second branch and homogeneous Neumann data on every face outside the edge set. For either branch,

Hence the total mass is constant when $c^{T}F=0$ on $R$ and nonincreasing when $c^{T}F\le 0$ on $R$. The first branch is choice-free; ${\mathrm{A}\mathrm{C}}_{\omega }$ is used only for the stated piecewise $C^1$ divergence-theorem branch.

Let $u_r\in C^4({\mathbb{T}}_L^d)$, and let $P_h{u}_r$ denote its samples on the periodic rectangular grid of Periodic rectangular reaction-diffusion grid. For every grid node $i$,

For $D=\mathrm{diag}(d_1,\dots ,d_m)$, define the sampled spatial residual $r_h=D(P_h\mathrm{\Delta }u-{\mathrm{\Delta }}_h{P}_{h}u)$, where ${\mathrm{\Delta }}_h=-A_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$ acts componentwise. If the indicated fourth derivatives are uniformly bounded for $0\le t\le T$, then $\mathrm{\parallel }{r}_h(t){\mathrm{\parallel }}_{\mathrm{\infty }}\le C_{\mathrm{s}\mathrm{p}}{h}_{\max }^{}$ with $C_{\mathrm{s}\mathrm{p}}=\frac{d_{\max }}{12}{\sum }_j\mathrm{\parallel }{\mathrm{\partial }}_{jjjj}u{\mathrm{\parallel }}_{L^{\mathrm{\infty }}}$ componentwise, uniformly in the mesh.

Let $u:[0,T]\times {\mathbb{T}}_L^d\to R$ have bounded spatial derivatives through order two, with $R$ a compact rectangle. Let $F\in C^2$ on a neighborhood of $R$, and let ${\mathcal{A}}_h=-D{A}_h^{\mathrm{p}\mathrm{e}\mathrm{r}}$ be the periodic-grid diffusion generator from Periodic rectangular reaction-diffusion grid. Define $F_h(P_{h}u)=P_h(F(u))$ pointwise. Then, uniformly over all periodic grids with $h_j=L_j\mathrm{/}{N}_j$ and $N_j\ge 3$, ${\sup }_{\begin{array}{c}{N}_j\ge 3(1\le j\le d),0\le t\le T\end{array}}{\parallel {\mathcal{A}}_h{F}_h(P_{h}u(t))-DF(P_{h}u(t)){\mathcal{A}}_h{P}_{h}u(t)\parallel }_{\mathrm{\infty }}\le C_{\mathrm{c}\mathrm{o}\mathrm{m}}<\mathrm{\infty },$ where $C_{\mathrm{c}\mathrm{o}\mathrm{m}}$ depends only on the dimension, $D$, the suprema of $D^{2}u$, $DF$, $D^{2}F$, and the compact state rectangle, not on the mesh widths. This is the finite-grid commutator used in the first-order split-error proof; it is not a claim for arbitrary meshes or unbounded operators.

Let $F\in C^1$ and let $y_{\ast }\in {\mathbb{R}}^2$ be a homogeneous equilibrium, $F(y_{\ast })=0$. Set $J=DF(y_{\ast })=(J_{ij})$, assume $\mathrm{tr}J<0$ and $\det J>0$, and let $D=\mathrm{diag}(d_1,d_2)$ with $d_1,d_2>0$. For any spatial eigenfunction of $-\mathrm{\Delta }$ with eigenvalue $\lambda \ge 0$, the linearized reaction-diffusion mode has matrix $J-\lambda D$. For the continuous rectangular torus, $\lambda ={\lambda }_k={\sum }_j(2\pi k_j\mathrm{/}{L}_j{)}^2$ with $k\in {\mathbb{Z}}^d$. For the periodic grid of Periodic rectangular reaction-diffusion grid, the mode indexed by $\mathrm{\ell }$ has $\lambda ={\lambda }_{h,\mathrm{\ell }}=4{\sum }_j{h}_j^{-2}{\sin }^2(\pi {\mathrm{\ell }}_j\mathrm{/}{N}_j)$, with $0\le {\mathrm{\ell }}_j<N_j$.

Define

For either domain or grid and each actual mode, $J-\lambda D$ has a positive-real-part eigenvalue exactly when $q(\lambda )<0$; it is Hurwitz exactly when $q(\lambda )>0$, and has a zero eigenvalue when $q(\lambda )=0$. Thus the continuous or discrete linear instability criterion tests the actual corresponding spectrum. If all nonnegative $\lambda$ are regarded as continuously available, then a diffusion-driven unstable band exists exactly when ${\beta }_T>0$ and ${\beta }_T^2>4d_1{d}_2\det J$; its endpoints are

On a bounded torus or a finite grid, the band criterion alone is not enough: at least one actual Laplacian eigenvalue must lie strictly between ${\lambda }_{-}$ and ${\lambda }_{+}$.