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: trapezoidal rule

BDF2 with one step Heun for initialization: enter image description here

For comparison: Implicit Euler: enter image description here

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: enter image description here

Plot of $y=Cx$ computed with the trapezoidal rule enter image description here

Close up for the trapezoidal rule on the initial phase -- large time step: enter image description here

Close up for the trapezoidal rule on the initial phase -- small time step: enter image description here

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

  • $\begingroup$ 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 $\endgroup$
    – VoB
    Commented Dec 16, 2020 at 12:25
  • $\begingroup$ 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. $\endgroup$
    – Jan
    Commented Dec 16, 2020 at 12:59
  • $\begingroup$ Indeed, an exponential integrator does not require the whole $A$ explicitely, but only its action on a vector $v$. $\endgroup$
    – VoB
    Commented Dec 16, 2020 at 13:26
  • 1
    $\begingroup$ 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 $\endgroup$
    – VoB
    Commented Dec 16, 2020 at 14:15
  • 1
    $\begingroup$ This may be relevant to your problem doi.org/10.1016/j.jcp.2019.04.070 $\endgroup$
    – cfdlab
    Commented Dec 18, 2020 at 4:15

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

Integration up to <span class=$T=0.5$" />

Integration up to <span class=$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))
y=np.zeros((ts+1)) #trapz
yE=np.zeros((ts+1)) #Backward Euler

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,yE,marker='d',label='Backward Euler')
plt.title("Integration up to T="+str(tf))

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.

  • $\begingroup$ 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. $\endgroup$
    – Jan
    Commented Dec 18, 2020 at 14:24
  • 1
    $\begingroup$ 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 $\endgroup$
    – VoB
    Commented Dec 18, 2020 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.