LinOSS-IM: Implicit Time Integration

• LinOSS-IM is a high-order implicit time integration method offering unconditional stability and tunable numerical dissipation for stiff problems.
• It employs rational approximation and partial fraction expansion to efficiently recover acceleration and derived quantities with high accuracy.
• Designed for PDEs, structural dynamics, and sequence modeling, it decouples explicit and implicit terms to ensure robust and efficient integration.

Implicit Time Integration (LinOSS-IM) is a high-order unconditionally stable time discretization framework designed for stiff linear and nonlinear initial value problems, particularly in PDEs, structural dynamics, and sequence modeling. LinOSS-IM combines algebraic stability, controllable high-frequency numerical dissipation, efficient recovery of derived quantities (notably acceleration in transient dynamics), and practical implementation via a mix of rational approximation and partial fraction expansion.

1. Foundational Principles and Mathematical Formulation

The LinOSS-IM methodology originates from semi-discretized dynamical systems of the form:

$$M\,\ddot{u}(t) + C\,\dot{u}(t) + K\,u(t) = f_I(u,\dot{u}) + f_E(t)$$

where $M$ (mass), $C$ (damping), and $K$ (stiffness) matrices are typically symmetric and positive definite, and $f_I$, $f_E$ denote nonlinear internal and external forcing. Introducing the augmented state $z=[\dot{u};u]$, recasting in dimensionless time $s=(t{-}t_{n-1})/\Delta t$ yields a first-order ODE:

$\frac{dz}{ds} = A\,z(s) + \overline{f}(z(s))$

with explicit forms for $A$ and $\overline{f}$ dependent on $\Delta t$ and the discretization. LinOSS-IM applies a rational approximation $R(A)=P(A)/Q(A)$ for the matrix exponential $e^{A}$, using a mixed Padé strategy controlled by the spectral radius parameter $\rho_\infty\in[0,1]$ that interpolates between A-stable nondissipative and L-stable highly dissipative schemes (O'Shea et al., 20 Sep 2024).

Through partial fraction decomposition of $Q(A)$, the time-stepping update is reduced to solving $M$ linear systems per time step, each of the form:

$$(r_i^2\,M + r_i\,\Delta t\,C + \Delta t^2\,K)\,x^{(i)}_1 = r_i\,M\,g^{(i)}_1 - \Delta t^2\,K\,g^{(i)}_2 + r_i\,\Delta t^2\,f_{r_i}$$

with associated formulae for the updated velocity and displacement (O'Shea et al., 20 Sep 2024).

2. Stability, Dissipation, and Order of Accuracy

LinOSS-IM offers unconditional stability properties that are tunable via $\rho_\infty$:

• $\rho_\infty=1$: diagonal Padé yields strictly A-stable, zero-numerical-dissipation schemes (order $2M$).
• $0<\rho_\infty<1$: mixed diagonal/subdiagonal Padé provides controllable high-frequency damping (order $2M-1$).
• $\rho_\infty=0$: subdiagonal Padé corresponds to L-stable (maximal dissipation), with order $2M-1$.

This rational approximation preserves accuracy even in the presence of external forcing and physical damping and guarantees that the computed acceleration remains at the full order of accuracy of the primary solution (O'Shea et al., 20 Sep 2024).

For nonlinear problems, LinOSS-IM uses high-order quadrature with force evaluations at Gauss–Lobatto points and Hermite interpolation to maintain the designed accuracy (O'Shea et al., 20 Sep 2024).

3. Efficient Computation of Acceleration and Derived Quantities

A distinctive feature of LinOSS-IM is the direct recovery of acceleration at each step. Specifically:

$$\ddot{u}_n = \rho\,\ddot{u}_{n-1} + \sum_{i=1}^M a_i [r_i\,x^{(i)}_1 - g^{(i)}_1]$$

where all terms are computed as part of the primary update, requiring only vector and scalar operations with no additional solves or inversions (O'Shea et al., 20 Sep 2024). This is especially significant in elastodynamics and wave propagation, where acceleration is often a critical observable. Standard second-order methods can produce polluted acceleration responses, but LinOSS-IM maintains high-quality accuracy across both displacement and acceleration channels (O'Shea et al., 20 Sep 2024).

4. Linearly Stabilized Schemes for Nonlinear PDEs

LinOSS-IM is applicable to stiff nonlinear PDEs by explicit–implicit splitting:

$u_t = F(u) + L\,u$

with nonlinear $F(u)$ treated explicitly and a strong linear stabilizer $p\,L\,u$ handled implicitly. The one-step scheme (SBDF1-type):

$$u^{n+1} = (I - \Delta t\,pL)^{-1} \left[ u^n + \Delta t\left(F(u^n) - pL\,u^n\right) \right]$$

ensures unconditional stability for $\bar{p} \geq \frac{1}{2}$, with higher-order variants (SBDF2) raising the stability threshold to $\bar{p} \geq \frac{3}{4}$ (Chow et al., 2021). These schemes require only one sparse linear solve and one evaluation of $F(\cdot)$ per step. LinOSS-IM decouples explicit (nonlinear) and implicit (linear stabilizer) terms, allowing efficient time integration at large time steps without nonlinear solves typically required by fully implicit methods.

5. Artificial Dissipation and Energy Behavior

LinOSS-IM’s implicit structure naturally produces negative energy production at each step:

$$\|u_{+}\|_M^2 \leq \|u_{0}\|_M^2 - (\Delta t)^2\,\|f(u_{+})\|_M^2$$

guaranteeing full discrete stability and suppressing nonphysical energy growth (Öffner et al., 2016). For explicit schemes, modal filtering or artificial viscosity must be adaptively chosen at each Euler stage to mimic the dissipative effect of implicit integration. This equivalence can be made precise: applying two forward-Euler steps with per-step adaptive filtering yields exactly the same dissipation as one LinOSS-IM update (Öffner et al., 2016). The minimal necessary filter strength is given by:

$$\epsilon = \frac{\Delta t\,\|f(u)\|_M^2}{2\,\sum_{n=0}^p \lambda_n^s\,\tilde u_{+,n}^2\,\|\phi_n\|^2}$$

6. Advanced Implicit Integration and Structural Features

LinOSS-IM incorporates practical high-order features:

• Self-starting: No need for acceleration or prior derivatives at initial time.
• Partial fraction expansion allows for amortized matrix factorization; in “single-multiple-root” formulations, only one matrix factorization per step is required (O'Shea et al., 20 Sep 2024).
• Force sampling via Gauss–Lobatto points matches the rational approximation order, and Hermite interpolation facilitates high-order internal state evaluation without additional ODE solves.

Recent advances include linearly implicit multistep methods (LIMM) that require only one linear system per time step and achieve similar or superior stability compared to BDF, with orders up to five, optimized coefficients, and efficient error estimation and variable-step/variable-order control (Glandon et al., 2020).

7. Applications and Extensions

LinOSS-IM has demonstrated success in diverse areas:

• Stiff nonlinear PDEs: Mean-curvature motion, nonlinear diffusion, with confirmed high-order convergence and stability (Chow et al., 2021).
• Structural dynamics: Wave propagation and vibration in the presence of both linear and nonlinear forces, delivering highly accurate and physically meaningful acceleration responses (O'Shea et al., 20 Sep 2024).
• High-dimensional sequence modeling: By appropriate splitting of implicit and explicit terms, LinOSS-IM is applicable in state-space models for long-range forecasting and learning on long sequences (Rusch et al., 4 Oct 2024). IMEX discretization preserves reversibility and volume in phase space, whereas fully implicit LinOSS-IM introduces controlled dissipation and forgetting.
• Multicomponent reacting flows: Component-splitting implicit integration yields significant performance gains via tailored blockwise flux-vector splitting, scaling linearly with the number of chemical species and enabling robust convergence and accuracy in chemically reacting compressible flow simulations (Zhang et al., 6 Mar 2024).

Table: Illustrative LinOSS-IM Features

Scheme Variant	Dissipation Control	Per-Step Operation Count	Self-Starting	Acceleration Recovery
LinOSS-IM (Padé Diagonal)	None ($\rho_\infty=1$)	$M$ linear solves	Yes	Yes
LinOSS-IM (Padé Mixed)	Tunable ($0<\rho_\infty<1$)	$M$ linear solves	Yes	Yes
SBDF1/SBDF2	Full (by choice of $p$)	1 linear solve	Yes	Indirect
LIMM (order up to 5)	Via multistep coefficients	1 linear solve	Yes	By formula

The table highlights that LinOSS-IM and its related advances achieve both high algebraic order and stability/dissipation tuning with computational efficiency and practical advantages for derived quantities.

---

All statements and formulas are precisely quoted or paraphrased from (O'Shea et al., 20 Sep 2024, Chow et al., 2021, Öffner et al., 2016, Glandon et al., 2020, Zhang et al., 6 Mar 2024, Rusch et al., 4 Oct 2024).