The semi-implicit Euler method can be applied for stiff ODE systems with the initial condition e.g., \[ \frac{du}{dt} = l(t,u) \quad \text{and}\quad u(t_0) = u_0 \] where \( l(t,u) \) represents the PDE such that the weak formulation of the stationary elliptic problem can be written as \[ B \dot u = L(u). \]
To determine the step size, one Newton iteration per discretization step is needed to solve. This approach is known as semi-implicit or linear-implicit discretization methods and requires a Jacobian calculation of \( L \).
We consider the general problem formulation \[ B \dot u = L(u) \] with constant matrix \( B \). Imposing the implicit Euler method and linearizing the term \( L(u^{k+1}) \) by Taylor expansion of first order yields the so-called linearly implicit Euler method (or as well called semi-implicit, semi-explicit, modified, symplectic method) \[ \begin{aligned} B \frac{u^{k+1}-u^k}{\tau} &= L(u^{k+1}) \\ &= L(u^k) + L’(u_0)(u^{k+1}-u^k). \end{aligned} \] with linearization point \( u_0 \) for the Jacobian and stepsize \( \tau \).
Note: The linearization point \( u_0 \), which is the initial value used to compute the Jacobian denoted by \( L’ \), can be replaced by the current iterate \( u_k \) at each time step. If so, we need to assemble the left-hand side at each iteration even though \( B \) does not depend on \( u \) due to Jacobian computation. Otherwise we can compute the Jacobian at \( u_0 \) and assemble the left-hand side only once, which will speed up the computation in each iteration, however it might converge slower to the solution.
The problem will then be formulated as \[ B \delta u^k = \tau L(u^k) + \tau L’(u_0) \delta u^k \] which leads us after a simple rearrangement to \[ (B-\tau L’(u_0))\delta u^k = \tau L(u^k) \] with the displacement \( \delta u^k = u^{k+1}-u^k\). By denoting \( D:=(B-\tau L’) \), this expression corresponds to one Newton iteration \( u^{k+1} = u^k - \tau D^{-1}L(u^k) \) for the nonlinear system, arising from the implicit Euler discretization. The left-hand side \( D:=B - \tau L’(u_0) \) is assembled as matrix, and the right-hand side \( \tau L(u_k) \) is assembled as vector by the assembler.
Example: In the tutorial for a time-dependent PDE, the problem has the specific form \[ (B-\tau A)\delta u^k = \tau (Au^k-f). \] with stiffness matrix A, mass matrix B and source f. The stiffness matrix does not depend on u, hence the choice of the linearization point in this case is irrelevant. One solution step to solve this system of ODEs could then be performed by a direct inverse operation (as the simplest method among others) \[ \delta u^k = (B-\tau A)^{-1}\tau (Au^k-f) \] resulting in the new iterate \( u^{k+1} = u^k + \delta u^k \).
If the mass matrix \( B \) depends on \( u \), then \( B(u) \) has to be an invertable Nemyckii operator (a diagonal operator, mapping \( u(x) \) to \( B(x,u(x)) u(x) \). Then the linearly implicit Euler scheme is \[ (B(u_0) - \tau L’(u_0) + \tau_k B^{-1}(u_0)B’(u_0)\mathrm{diag}L(u_0)) \delta u_k = \tau_k L(u)+\frac{\tau_k}{\tau_{k-1}}(B(u_0)-B(u_k))\delta u_{k-1}. \] Derivation: see Lubich and Roche (Rosenbrock Methods for Differential-algebraic Systems with Solution-dependent Singular Matrix Multiplying the Derivative, Computing 43:325-342, 1990)
The class SemiImplicitEulerStep defines the weak formulation of a stationary elliptic problem, resulting from time discretization of a parabolic PDE, and provides the implementation of \( B-\tau L’(u_0) \) and right-hand side \( \tau L(u^k) \) for assembly, based on the d1, d2 and b2 caches of the equation. The ParabolicEquation concept defines the interface that is accessed by the SemiLinearizationAtInner, SemiLinearizationAt or SemiLinearization adapters.
The class SemiLinearizationAt, a proxy class for the linearization of semi-linear time stepping schemes, requires a second linearization point \( u_0 \) which is used to compute the Jacobian denoted by \( L’ \) occuring in JU: ?
\[ B \dot u - L’ u = L(u) - L’ u. \]
In case of a subdomain problem and including the inner boundaries functions, SemiLinearizationAtInner is applied as a proxy class where InnerBoundaryCache class is included as well. Otherwise, the rest is similar to class SemiLinearizationAt. More details will be added later.