= Solution
Choose a fixed time increment $h=\Delta t>0$, let $t_n=nh$, and approximate $\mathbf c(t_n)$ by $\mathbf c_n$. The <Backward Euler method> evaluates the transport term at the new time:
$$
\frac{\mathbf c_{n+1}-\mathbf c_n}{h}=A\mathbf c_{n+1}-\mathbf b.
$$
Rearrange to the <linear system>
$$
\boxed{B\mathbf c_{n+1}=\mathbf c_n-h\mathbf b,\qquad
B=I-hA=\begin{pmatrix}
1+2h&-h&0&-h\\
-2h&1+2h&0&0\\
0&-h&1+2h&-h\\
0&0&-2h&1+2h
\end{pmatrix}.}
$$
Although $A$ is a <singular matrix>, $B$ is an <invertible matrix> for every $h>0$. In each row, its diagonal entry is $1+2h$ and the sum of the <absolute values> of its off-diagonal entries is $2h$; hence it is a <strictly diagonally dominant matrix>. Its <Gershgorin discs> have centre $1+2h$ and radius $2h$, excluding zero. In particular, the zero <eigenvalue> of $A$ becomes an <eigenvalue> $1$ of $B$.
The algorithm uses <Reusing LU factors for multiple right-hand sides>:
* Form $B$ once. Compute a <Pivoted LU decomposition> $PB=L_B U_B$, storing the <permutation matrix> $P$, the unit <lower triangular matrix> $L_B$, and the <upper triangular matrix> $U_B$.
* Set $\mathbf c_0=\mathbf c(0)$.
* At each step form $r_n=\mathbf c_n-h\mathbf b$. Solve $L_B y=Pr_n$ by <forward substitution>, then solve $U_B\mathbf c_{n+1}=y$ by <back substitution>.
* Advance $n$ and repeat until the desired final time. Choose the fixed increment so the final time is an integer multiple of $h$, or refactor $B$ if a shorter final step is used.
Both triangular solves use the factors of \b[$I-hA$]. The <LU factorization> of the singular transport <matrix> $A$ from part (a) is not a solver for this new <linear system>. No <matrix inverse> needs to be formed. Since $A$ and $h$ are fixed, the same <LU factorization> serves every step; only the right-hand side changes. A varying $\mathbf b(t)$ would give $r_n=\mathbf c_n-h\mathbf b(t_{n+1})$ without changing these factors.
As a check, multiplying the update by $\mathbf1^T$ gives
$$
\mathbf1^T\mathbf c_{n+1}=\mathbf1^T\mathbf c_n-h\mathbf1^T\mathbf b,
$$
exactly reproducing the inventory balance at the sampled times, apart from <roundoff error>. The <Backward Euler method> has temporal <global error> $O(h)$ for smooth solutions, so its <A-stability> does not remove the need to select $h$ for accuracy.
Back to article page
