Second Order Time Integration for Stiff Linear System Avoiding any Explicit Step

Question

I have a linear system

$$ \dot x(t) = Ax(t), \quad x(0)=x_0 \tag{*} $$

with $A$ being Hurwitz (i.e. the solutions may oscillate but will eventually tend to zero) but really stiff. $A$ might be large, but I can do explicit solves.

I want to integrate it in time. Since it is coupled with a flow solver, stability is my primary concern. However, since my flow is integrated with order 2, I would fancy a 2nd order scheme for $(*)$ too.

I found that the implicit trapezoidal rule gives accurate and stable solutions. But the oscillations in the initial phase will likely destroy my simulation.

That's why I went for BDF2. Here, however, I face the problem of the initialization (I need a second order approximation of the value at the first time step computed with a different scheme). The canonical choice of Heun's method led to extreme overshoots in the approximation. Things slightly improve when I use an Implicit Euler step for the prediction. But the correction step is still explicit and produces the overshoot.

In both cases, a smaller time step improves things -- but I don't want to go to small.

---

So my question is: Is there a method that gives my a second order increment of

$$x(0)=x_0\to x_1\approx x(0+h)$$

that completely avoids the explicit application of $A$?

---

Below are plots of my numerical tests. I have cropped the y-axis to better see what is happening.

Trapezoidal rule:

BDF2 with one step Heun for initialization:

For comparison: Implicit Euler:

---

EDIT-2: some more plots

Actually, in the above plots, only the output of the system was shown, namely $y=Cx$, where $x$ is the state and $C$ is a constant matrix.

Plot of $x$ computed with the Implicit Euler scheme:

Plot of $y=Cx$ computed with the trapezoidal rule

Close up for the trapezoidal rule on the initial phase -- large time step:

Close up for the trapezoidal rule on the initial phase -- small time step:

EDIT: some words on the system (matrices)...

The system derives from a semi-discrete approximation of the incompressible Navier-Stokes equations that is coupled to a controller. In theory, it reads

\begin{align} \dot v &= \Pi \tilde A(v) + \Pi f (v) + B\hat C \hat v \\ \dot{\hat v} &= \hat A \hat v + \hat B C v \end{align}

with a nonlinearity $f$ and a projector $\Pi=I-J^T(JJ^T)J$ for some matrix $J$.

In the practical realization, I cannot use the projector $\Pi$ but it's realization through solving saddle point systems. E.g., one step of the implicit trapezoidal rule will read like $$ \begin{bmatrix} v_{k+1} \\ \hat v_{k+1} \\ \sim \end{bmatrix} = \begin{bmatrix} I - \frac{h}{2} \tilde A & - \frac{h}{2}B\hat C & -J^T \\ - \frac{h}{2}\hat BC & I - \frac{h}{2}\hat A & 0 \\ -J & 0 & 0 \end{bmatrix}^{-1} \begin{bmatrix} I + \frac{h}{2} \tilde A & \frac{h}{2}B\hat C & 0 \\ \frac{h}{2}\hat BC & I + \frac{h}{2}\hat A & 0 \\ 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} v_{k} \\ \hat v_{k} \\ 0 \end{bmatrix} $$

Here, $\tilde A$, $J$ are large but sparse matrix of size, say, $10^5$. The matrix $\hat A$ is small (of size ~100).

Answer 1

See comments for the original discussion.

The lack of L-stability of the trapezoidal rule in the first step is the source of your problem. A simpler and famous toy problem that shows the point is the Curtiss-Hirschfelder equation $$y'=-2000(y-\cos(t)) \\y(0)=0$$

Integrating this with Backward Euler and the Trapezoidal rule gives the following result (I wrote the easy routine myself for illustration purposes, without using odeint or other libraries)

$T=0.5$" />

$T=10.0$" />

obtained with the following Python snippet

import numpy as np
import matplotlib.pyplot as plt
    
# Trapz rule is NOT L-stable! (While Backward Euler is)
# Test equation: y'(t)=-2000*(y-cos(t)), y(0)=0
C = -2000

def Fun(t,x):
    return C*(x-np.cos(t))
    
tf=10.0
t0=0
ts=100
k=(tf-t0)/ts
t=np.linspace(t0,tf,ts+1)
y=np.zeros((ts+1)) #trapz
yE=np.zeros((ts+1)) #Backward Euler
y0=0
y[0]=y0
yE[0]=y0

for i in range (0,ts):
    y[i+1] = (y[i]+0.5*k*(Fun(t[i],y[i]) -C*np.cos(t[i+1]) ) ) /(1.0-0.5*k*C)
    yE[i+1] = (yE[i] + k*C*(-np.cos(t[i+1])) )/(1-k*C)

plt.plot(t,y,marker='o',label='Trapezoidal')
plt.plot(t,yE,marker='d',label='Backward Euler')
plt.title("Integration up to T="+str(tf))
plt.legend()
plt.show()

As pointed out in the reference from @cfdlab: https://doi.org/10.1016/j.jcp.2019.04.070 Indeed we can read:

BDF2 should be started with a second-order L-stable scheme such as SDIRK.

which is exactly your problem. The good point is that BDF2 is L-stable, so as long as you take an L-stable method with order 2 as first step, you'll have a second order scheme which won't have those wild oscillations at the beginning

Another way could be to choose as first step the exact solution which would require to compute $e^{h A}x_0$ with some suitable routine for the action of the matrix exponential.