Reaction-diffusion systems, invariants and pattern computation — Examples

On an element $[0,h]$, the two affine basis functions are ${\varphi }_1(x)=1-x\mathrm{/}h$ and ${\varphi }_2(x)=x\mathrm{/}h$. Direct integration gives ${\int }_0^h{\varphi }_1^{2}dx={\int }_0^h{\varphi }_2^{2}dx=h\mathrm{/}3,{\int }_0^h{\varphi }_1{\varphi }_{2}dx=h\mathrm{/}6,$ and their derivatives are $-1\mathrm{/}h$ and $1\mathrm{/}h$, so the element stiffness entries are $1\mathrm{/}h,-1\mathrm{/}h,-1\mathrm{/}h,1\mathrm{/}h$. Assembling the two elements yields the displayed matrices: $M=\frac{h}{6}\left(\begin{array}{ccc}2& 1& 0\\ 1& 4& 1\\ 0& 1& 2\end{array}\right),K=\frac{1}{h}\left(\begin{array}{ccc}1& -1& 0\\ -1& 2& -1\\ 0& -1& 1\end{array}\right)\mathrm{.}$ The Neumann condition is natural in this stiffness form and adds no boundary matrix.

Write $M=(h\mathrm{/}6)B$ and $K=C\mathrm{/}h$ for the displayed integer matrices. Direct multiplication verifies $B^{-1}=\frac{1}{12}\left(\begin{array}{ccc}7& -2& 1\\ -2& 4& -2\\ 1& -2& 7\end{array}\right),B^{-1}C=\frac{1}{12}\left(\begin{array}{ccc}9& -12& 3\\ -6& 12& -6\\ 3& -12& 9\end{array}\right)\mathrm{.}$ Thus $M^{-1}K=6B^{-1}C\mathrm{/}{h}^2$, so in particular $(M^{-1}K{)}_{13}=3\mathrm{/}(2h^2)$. With $e_3=(0,0,1{)}^T$, the ODE gives $u^{\mathrm{\prime }}(0)=-\nu M^{-1}K{e}_3$, whose first component is $-3\nu \mathrm{/}(2h^2)<0$.

The linear field $u\mapsto -\nu M^{-1}Ku$ is globally Lipschitz. By [F1], $e^{-t\nu M^{-1}K}{e}_3$ solves the IVP, and [F2] makes it the unique solution. Since $u_1(0)=0$ and $u_1^{\mathrm{\prime }}(0)<0$ by step 2.1, the right difference quotient $u_1(t)\mathrm{/}t$ is negative for all sufficiently small $t>0$. Hence the consistent-mass P1 semidiscrete diffusion is not positivity preserving, despite the nonnegative initial nodal vector.

Direct substitution gives $F(1,3\mathrm{/}2)=(0,0)$. Differentiating the two components gives $J=DF(1,3\mathrm{/}2)=\left(\begin{array}{cc}1\mathrm{/}2& 1\\ -3\mathrm{/}2& -1\end{array}\right),\mathrm{tr}J=-1\mathrm{/}2,\det J=1.$

The coefficient in the dispersion polynomial is ${\beta }_T=d_2{J}_{11}+d_1{J}_{22}=9\mathrm{/}10$ and $d_1{d}_2=1\mathrm{/}5$. Therefore $q(\lambda )=\frac{1}{5}{\lambda }^2-\frac{9}{10}\lambda +1=\frac{1}{5}(\lambda -2)(\lambda -\frac{5}{2})\mathrm{.}$

For the continuous torus mode $k=1$, $\lambda =(2\pi \mathrm{/}L{)}^2=9\mathrm{/}4$. This lies strictly between $2$ and $5\mathrm{/}2$, so step 2.1 gives $q(9\mathrm{/}4)<0$. By [F1], the mode matrix has one positive eigenvalue and one negative eigenvalue. This example proves a linear mode instability only; it does not assert nonlinear pattern existence or stability.

Here $h=L\mathrm{/}3=4\pi \mathrm{/}9$. The grid symbol gives ${\lambda }_0=0$ and ${\lambda }_1={\lambda }_2=\frac{3}{h^2}=\frac{243}{16{\pi }^2}<2,$ where the last inequality follows from $\pi >3$, since $243\mathrm{/}(16{\pi }^2)<243\mathrm{/}144=27\mathrm{/}16<2$.

At ${\lambda }_0=0$, the reaction Jacobian has trace $-1\mathrm{/}2$ and determinant $1$, so it is Hurwitz by the two-dimensional criterion in [F2]. For $\mathrm{\ell }=1,2$, step 1.1 gives $0<{\lambda }_{\mathrm{\ell }}<2$, so both factors in $q({\lambda }_{\mathrm{\ell }})=\frac{1}{5}({\lambda }_{\mathrm{\ell }}-2)({\lambda }_{\mathrm{\ell }}-\frac{5}{2})$ are negative and their product is positive. The mode trace is $-1\mathrm{/}2-(1\mathrm{/}10+2){\lambda }_{\mathrm{\ell }}<0$, so [F2] makes both modes Hurwitz. These are all three grid modes, even though the continuous mode in [F1] is unstable.

On $[0,1]$, $F(0)=F(1)=0$, so the field points inward at both endpoints. Also $\mathrm{\mid }{F}^{\mathrm{\prime }}(u)\mathrm{\mid }=\mathrm{\mid }r(1-2u)\mathrm{\mid }\le r$, giving the coordinate Lipschitz constant $L=r$ and Euler restriction $\tau \le 1\mathrm{/}r$.

For $u\in [0,1]$, set $\alpha =\tau r\in [0,1]$. The Euler value is $u^{+}=u+\alpha u(1-u)\ge 0$, while $1-u^{+}=(1-u)(1-\alpha u)\ge 0$ because both factors are nonnegative. Thus the reaction step preserves $[0,1]$.

If $\alpha =\tau r>1$, choose $u$ with $1\mathrm{/}\alpha <u<1$. Then $1-u^{+}=(1-u)(1-\alpha u)<0$, so the Euler map leaves $[0,1]$. Together with step 2.1 this proves the stated necessary-and-sufficient bound for the reaction map.

The reaction map satisfies the hypothesis of [F3] by step 2.1; applying [F2] and then [F3] proves that the full TPFA Lie update also remains in $[0,1]$. For a grid state $U=(u_i{)}_i$ with positive weights $m_i$, its Euler reaction mass increment is ${\sum }_i{m}_i\tau r{u}_i(1-u_i)$. This is strictly positive for $\tau >0$ exactly when some $u_i\in (0,1)$, and is zero when every $u_i\in \{0,1\}$; for $\tau =0$ it is zero for every state. Thus the mass claim in the Statement holds, and no mass-conservation claim is made.