7
$\begingroup$

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:

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

[...]

  1. 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?

$\endgroup$
6
  • $\begingroup$ What's in the Jansen et al. article? Did you read it? What didn't you understand? $\endgroup$
    – Bill Barth
    Sep 22, 2014 at 12:02
  • $\begingroup$ @BillBarth You can have a look here scorec.rpi.edu/~kjansen/genalf.pdf . The article assumes working knowledge of the answer to my question, rather than providing a basic explanation. However, there is a useful clue: "Given the non-linearities present in most interesting engineering systems, it is of great importance to ensure that there is temporal damping for frequencies beyond the chosen resolution level. However, it is equally important that this damping not [affect] the frequencies within the chosen resolution level. This delicate balancing act [...]" $\endgroup$ Sep 24, 2014 at 4:07
  • $\begingroup$ It might be a good idea to limit the discussion to the step-response of a system. The step function has non-zero components for arbitrary large frequencies, and if we wish to find the solution by integrating numerically, we must choose some finite resolution, and, therefore, an upper-bound frequency. $\endgroup$ Sep 24, 2014 at 4:10
  • $\begingroup$ @BillBarth By the way, does this sound like a silly question to you? You should remember that not everybody is an expert in this field. $\endgroup$ Sep 24, 2014 at 5:33
  • $\begingroup$ It doesn't sound like a silly question at all to me, but you gave us nothing to go on. I'm certainly not an expert in ODE integration theory either. $\endgroup$
    – Bill Barth
    Sep 24, 2014 at 11:54

1 Answer 1

6
$\begingroup$

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.

$\endgroup$

Your Answer

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

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