Positive Conservative Production Destruction Integrators — Examples

The rates satisfy $p_{21}=d_{12}=1000y_1$ and $p_{12}=d_{21}=y_2$, so [F1] gives $y_1^{\mathrm{\prime }}=-1000y_1+y_2$ and $y_2^{\mathrm{\prime }}=1000y_1-y_2$. In the state order $(1,2)$ the coefficient matrix is $Q=\left(\begin{array}{cc}-1000& 1\\ 1000& -1\end{array}\right)$; its off-diagonal entries are nonnegative and both column sums are zero, so it is a finite column generator under [F4].

With $h=1$ and $y^n=(1,1)$, [F2] gives $\left(\begin{array}{cc}1001& -1\\ -1000& 2\end{array}\right)y^{n+1}=\left(\begin{array}{c}1\\ 1\end{array}\right)$. The determinant is $1001\cdot 2-1000=1002$, and direct elimination gives $y_1^{n+1}=3\mathrm{/}1002=1\mathrm{/}334$ and $y_2^{n+1}=2001\mathrm{/}1002=667\mathrm{/}334$. Both coordinates are positive and their sum is $2$, as also predicted by [F3].

For total mass $M=y_1+y_2$ and $s=1001$, the continuous first coordinate obeys $y_1^{\mathrm{\prime }}=-s(y_1-y_1^{\ast })$ with equilibrium $y_1^{\ast }=M\mathrm{/}1001$. The one-step formula in step 1.2 satisfies $y_1^{n+1}-y_1^{\ast }=(y_1^n-y_1^{\ast })\mathrm{/}(1+hs)$; for $h=1$ this positive damping factor is $1\mathrm{/}1002$ although $hs=1001$. This verifies the large-rate damping for the specified example without inferring a general stability theorem.

One forward event consumes mass $2$ from $A$ and produces mass $1$ in each of $B,C$; one reverse event consumes mass $1$ from each of $B,C$ and produces mass $2$ in $A$. Thus the mass transfer rates are $p_{21}=p_{31}=d_{12}=d_{13}=r=y_1$ and $p_{12}=p_{13}=d_{21}=d_{31}=s=y_2{y}_3$, with all other rates zero. These paired rates satisfy [F1], and $y_1+y_2+y_3=2n_A+n_B+n_C$ is the molecular-weight mass.

The transfer factors are $a_{21}=a_{31}=1$, $a_{12}=y_3$, and $a_{13}=y_2$, with all other factors zero; hence the rates have the form required by [F2], even when $y_2$ or $y_3$ is zero. At $y^n=(2,0,0)$, the matrix in [F3] is $I-B(y^n)=\left(\begin{array}{ccc}3& 0& 0\\ -1& 1& 0\\ -1& 0& 1\end{array}\right)$. Solving it against $(2,0,0{)}^T$ gives $y^{n+1}=(2\mathrm{/}3,2\mathrm{/}3,2\mathrm{/}3{)}^T$.

The updated coordinates are nonnegative, and $2\mathrm{/}3+2\mathrm{/}3+2\mathrm{/}3=2=2+0+0$. This is the weighted number-to-mass conservation verified in step 1.1; it also agrees with the general zero-compatible theorem [F4].

The rates are $p_{21}=d_{12}=2y_X$ and $p_{12}=d_{21}=3y_Y$, so [F1] holds. At $y^n=(1,0)$, the standard choice has ${\sigma }_Y=0$, so [F2] returns inadmissible-input before evaluating a rate. The displayed standard-step equation itself would contain $p_{12}(y^n)y_Y^{n+1}\mathrm{/}{\sigma }_Y=0\cdot y_Y^{n+1}\mathrm{/}0$, which is undefined; assigning $0\mathrm{/}0=0$ is not an algebraic evaluation of that expression.

The transfer factors $a_{21}=2$ and $a_{12}=3$ are nonnegative and $C^2$ on $U={\mathbb{R}}^2$, and $y^n=(1,0)\in K$ with $h=1>0$, so the hypotheses of [F5] hold. With $h=1$, [F4] gives $\left(\begin{array}{cc}3& -3\\ -2& 4\end{array}\right)y^{n+1}=\left(\begin{array}{c}1\\ 0\end{array}\right)$, whose determinant is $6$ and whose solution is $(y_X^{n+1},y_Y^{n+1})=(2\mathrm{/}3,1\mathrm{/}3)$.

The computed coordinates are nonnegative and sum to $1$, equal to the initial total. This agrees with [F5] and exhibits a well-defined first-order zero-state extension; it makes no second-order claim for a zero MPRK stage.

Set $m=1$, $p_{11}=d_{11}=0$, $h=1$, $y^n=0$, and ${\sigma }_1=1$. These data satisfy [F1] and the admissible input conditions of [F2], and the system $M{y}^{+}=y^n$ is $y^{+}=0$.

For the candidate $\widehat{y}=-\epsilon$, the exact residual is $r=M\widehat{y}-b=-\epsilon$, so $\mathrm{\mid }r\mathrm{\mid }=\epsilon$ can be made smaller than any prescribed positive tolerance while $\widehat{y}<0$.

Here the exact solve is $y^{\ast }=0$. Whenever the stored-operand and successful-arithmetic hypotheses of [F3] hold, its radius satisfies $\rho \ge \epsilon$, so its sufficient margin $\widehat{y}\ge \rho$ fails. For every real $\epsilon >0$, step 2.1 already gives a negative candidate with exact residual magnitude $\epsilon$, refuting the statement independently of floating-point evaluation.

The rates satisfy $p_{21}=d_{12}=y_1\ge 0$ and all other pair rates vanish, so [F1] gives $y_1^{\mathrm{\prime }}=-y_1$ and $y_2^{\mathrm{\prime }}=y_1$ with conserved sum $y_1+y_2$.

Forward Euler gives $y^1=y^0+hF(y^0)=(1,0)+(3\mathrm{/}2)(-1,1)=(-1\mathrm{/}2,3\mathrm{/}2)$, so the first concentration is negative even though the sum remains $1$.

The system matrix is lower triangular with eigenvalues $-1$ and $0$. For $\lambda =-1$, [F3] gives absolute stability because $0<h\rho =3\mathrm{/}2<2$; for $\lambda =0$, [F2] gives multiplier $R_E(0)=1$. Thus $\mathrm{\mid }1+h\lambda \mathrm{\mid }\le 1$ for both eigenvalues even though step 2.1 has a negative component.

The vector field is linear and smooth, so the compact-tube hypotheses of [F2] hold locally and forward Euler has first-order local consistency. Steps 2.1 and 3.1 show that neither this accuracy order nor scalar absolute stability implies positivity preservation.