6
$\begingroup$

I am working with the simple example of an oscillator:

$$(1) \; \; \ddot{u} + u = 0, \; \; u(0) = u_0$$

I know that Forward Euler does not preserve an invariant of the above system: $$(2) \; \; \dot{u}^2 + u^2 = 1$$

I can get the eigenvalues of a matrix that approximately solves for the solution of $(1)$ using Forward Euler if I were to write it as a first order system of equations. However, I don't understand what I can deduce about the ability of Forward Euler to preserve $(2)$, given these eigenvalues and the stability properties of Forward Euler?

I would like to be able to understand this in general, for any scheme, not necessarily FE.

So for instance, I know that the Trapezoidal Rule scheme, unlike Forward Euler, is good at preserving $(2)$, but its difficult for me to see what is qualitatively different about the eigenvalues of the matrix for the problem based on the Trapezoidal Rule...

The Trapezoidal Rule is L-stable (that is, $Re(G(k\lambda)) < 0$ for $\dot{u} = \lambda u$ with $Re(\lambda) < 0$), but so is Backwards Euler, and it does not preserve $(2)$. What's the pattern here?

$\endgroup$
6
$\begingroup$

In a certain sense, @Geoff Oxberry is correct in saying that stability and preservation of quadratic invariants are not directly related. For instance, there exist explicit methods that will preserve energy for your problem (and they are certainly not $A$-stable).

However, in another sense there is a relation between the two.

Well-posedness

First, note that your problem as described is not well posed, and the invariant you give does not follow from the problem description. I believe you should have instead

$$\ddot{u} + u = 0$$

with initial conditions

$$u(0) = u_0$$ and $$\dot{u}(0) = v_0$$ where $u_0^2 + v_0^2 = 1$. Then the problem is well-posed and the invariant condition you have written follows.

Conservation of energy

Now, if you rewrite this second-order equation as a first-order system $$\begin{pmatrix} \dot{u} \\ \dot{v} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} u \\ v\end{pmatrix},$$

you find that the eigenvalues are $\pm i$. If we write the system above as $\dot{w} = A w$, the exact solution is given by $$w(t) = \exp(At)w(0)$$ and the reason you have energy conservation is that $|e^{\lambda t}|=1$ for both eigenvalues of $A$.

Discrete conservation

After discretizing with a one-step method, you will have instead an equation of the form $$w_{n+1} = R(\Delta t A)w_n,$$

where $\Delta t$ is the numerical time step and $R(z)$ is the stability function of the method (generally a polynomial or rational function). Since $A$ is normal, to have discrete conservation you need that $$|R(\pm i \Delta t)|=1.$$.

The connection

One way to interpret the condition just stated is to say that the eigenvalues of the ODE system, scaled by the numerical timestep, must lie on the boundary of the absolute stability region of the numerical integrator. That's because the stability region is defined precisely as the set of complex numbers $z$ for which $|R(z)| \le 1$.

Thus, with just the right time step, even the 4th-order Runge-Kutta method can be used to preserve energy for this problem, since its stability region boundary intersects the imaginary axis. Of course nobody recommends this, because it only works for one special value of $\Delta t$.

The trapezoidal rule conserves energy for any time step, because its region of absolute stability is exactly the left half-plane; i.e., the entire imaginary axis is on its stability region boundary, so the condition above is fulfilled for any $\Delta t$.

Meanwhile, it's impossible to conserve energy with forward or backward Euler, since those methods have a stability region whose boundary only touches the imaginary axis at the origin.

$\endgroup$
  • $\begingroup$ This turned out to be exactly the right answer. $\endgroup$ – user89 Nov 2 '14 at 18:39
5
$\begingroup$

Stability and preservation of invariants are unrelated. Hairer does a nice job of proving how various methods preserve different kinds of invariants in his book Geometric Integration. Chapter III is the most relevant to your question, where he demonstrates, for example, that RK methods and linear multistep methods preserve linear first invariants.

$\endgroup$
  • $\begingroup$ Thanks! I added a bit to my post just as you answered. Basically, I described the effect of a couple of different schemes and couldn't see any pattern...this explains it! $\endgroup$ – user89 Oct 31 '14 at 5:47
  • $\begingroup$ Could you confirm my understanding (based on a quick read-through of Hairer's book, searching for "eigenvalue" and "Trapezoidal rule") that similar to how the L-stability or A-stability of the scheme have no bearing on their ability to preserve invariants, neither do the eigenvalues of the matrix of the scheme applied to the problem of interest? $\endgroup$ – user89 Oct 31 '14 at 5:57
  • $\begingroup$ From reading through the proofs, eigenvalues have nothing to do with preservation of invariants. $\endgroup$ – Geoff Oxberry Oct 31 '14 at 6:02
  • $\begingroup$ The discussion in these comments is not correct -- the question is about a quadratic invariant, and that is certainly related to eigenvalues, at least if the Jacobian is normal. See my answer. $\endgroup$ – David Ketcheson Nov 2 '14 at 4:51
  • 1
    $\begingroup$ @user89: Depends on whether you're talking about the eigenvalues of the Jacobian matrix of the continuous ODE right-hand side or the eigenvalues of the Jacobian matrix of the nonlinear system you solve at each time step of an implicit method for solving ODEs numerically. The former doesn't depend on the time step; the latter does. $\endgroup$ – Geoff Oxberry Nov 2 '14 at 18:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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