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

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).

• Your stiff system admits the exact solution $e^A x_0$, being $e^A$ the matrix exponential. You might use the analytical solution in place of $x_1$, i.e. $x_1 = e^A x_0$, but this of course depends on the matrix $A$ as the matrix exponential is not a trivial task, and its accuracy strongly depends on the matrix. Could you provide us your $A$? I think that some ad-hoc technique like exponential integrators could solve the problem
– VoB
Dec 16, 2020 at 12:25
• An exponential integrator might be a remedy. However, there are two more caveats to my problem. (1) the linear system is but a part of an implicit-explicit integration. (2) in my application, I don't have the $A$ explicitly available but only applications of $A$ or $(I-sA)^{-1}$. I will add some explanation of the system to the body. If you still think that it might be worth a try, I can provide the matrices at least of the linear part.
– Jan
Dec 16, 2020 at 12:59
• Indeed, an exponential integrator does not require the whole $A$ explicitely, but only its action on a vector $v$.
– VoB
Dec 16, 2020 at 13:26
• Or you may just compute $e^A x_0$ and then use BFD-2. The point is that to compute $e^A x_0$ you can use a tailored highly accurate routine for the matrix exponential, so that the first step is exact up to machine precision @Jan
– VoB
Dec 16, 2020 at 14:15
• This may be relevant to your problem doi.org/10.1016/j.jcp.2019.04.070 Dec 18, 2020 at 4:15

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=y0
yE=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.

• Would you mind, providing a 2nd plot with a finer time step. I'm curious about the initial oscillations that are also seen in the backward Euler scheme.
– Jan
Dec 18, 2020 at 14:24
• There was an error in my BE :-) Now it's fine, thanks for pointing this out. I've added also the Python snippet with whom I got the plots @Jan
– VoB
Dec 18, 2020 at 18:31