What is numerical damping in the context of time-dependent FEM solvers?

Question

Comsol Multiphysics (a popular FEM package) includes two time-stepping algorithms (IDA aka BDF, and Generalized-alpha), described in their documentation as follows (quoted here under Fair Use; emphasis mine):

• IDA was created at the Lawrence Livermore National Laboratory (Ref. 4) and is a modernized implementation of the DAE solver DASPK, which uses variable-order variable-step-size backward differentiation formulas (BDF).

• Generalized-α is an implicit, second-order accurate method with a parameter α or ρ∞ (0 ≤ ρ∞ ≤ 1) to control the damping of high frequencies. With ρ∞ = 1, the method has no numerical damping. For linear problems this corresponds to the midpoint rule. ρ∞ = 0 gives the maximal numerical damping; for linear problems the highest frequency is then annihilated in one step. The method was first developed for second-order equations in structural mechanics and later extended to first-order systems (Ref. 7).

[...]

BDF methods have been used for a long time and are known for their stability. However, they can have severe damping effects, especially the lower-order methods. Backward Euler severely damps any high frequencies. Even if you are expecting a solution with sharp gradients, you might get a very smooth solution due to the damping in backward Euler.

[...]

Compared to BDF (with maximum order two), generalized-α causes much less damping and is thereby more accurate. For the same reason it is also less stable.

[...]

References:

4. A.C. Hindmarsh, P.N. Brown, K.E. Grant, S.L. Lee, R. Serban, D.E. Shumaker, and C.S. Woodward, “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers,” ACM T. Math. Software, vol. 31, p. 363, 2005.

[...]

7. K.E. Jansen, C.H. Whiting, G.M. Hulbert, “A Generalized-α Method for Integrating the Filtered Navier–Stokes Equations with a Stabilized Finite Element Method,” Comput. Methods Appl. Mech. Engrg., vol. 190, pp. 305–319, 2000.

This makes good sense, but I'd like to understand the essence of this numerical damping phenomenon. Is there a minimal example somewhere (1D is best, if possible) that shows numerical damping clearly, and allows to see where it comes from?

Answer 1

It is quite straightforward to demonstrate this for implicit Euler (aka backward Euler) on a scalar example. Consider the initial value problem

$$\dot{y}(t) = i \alpha y(t), \ y(0) = 1$$

with solution

$$y(t) = \exp(i \alpha t) = \cos(\alpha t) + i \sin(\alpha t)$$

The solution is a harmonic oscillation with amplitude 1 and this amplitude does not change over time.

Now apply implicit Euler to the ODE, which gives

$$y_{n+1} = y_{n} + h i \alpha y_{n+1}$$

with $h$ being your time step or, after rearranging,

$$y_{n+1} = \frac{y_n}{1 - h i \alpha} = \frac{1}{\alpha^2 h^2+1} \left(1 + i \alpha h \right) y_n$$

Because I assume $y(0) = 1$, that means

$$y_{n+1} = \left( \frac{1}{\alpha^2 h^2 + 1} \right)^n \left( 1 + i \alpha h \right)^n$$

While the first term goes to one in the limit $h \to 0$, we have

$$\frac{1}{\alpha^2 h^2 + 1} < 1$$

for any finite value of $h$. Of course, any actual numerical solution will involve a finite $h$. Therefore, we get

$$\left( \frac{1}{\alpha^2 h^2 + 1} \right)^n \to 0$$

as $n \to \infty$, meaning that in contrast to what should happen, the amplitude of the numerical solution will decrease over time. This effect is referred to as numerical damping because the decrease in amplitude is purely an artifact of the numerical method.

Numerical damping in implicit Euler is very strong. If you implement the scalar problem above, you will see that unless you use very small time steps $h$, after a few oscillations your numerical solution will have almost decayed to zero. Methods like BDF-2 also show numerical damping but much less pronounced than implicit Euler. I'd assume this is why Comsol recommends them in the text you posted.