Numerical Methods and Invariant Distributions

We have seen that the Liouville equation defines the propagation of the distribution associated to system of differential equations. In a similar way, a numerical method z +i = hiz ) induces an effective distributional propagator. One way to understand this propagator is to interpret the numerical method as being approximately equivalent to the solution of a perturbed differential equation 

The perturbation term may be viewed as the finite truncation of the perturbative expansion obtained from the backward error analysis using the methods of the previous chapter, (r is the order of the numerical method.) Let us assume that the perturbed equation is a realistic model for the numerical solutions, then define the evolving density by 

Under suitable boundedness assumptions on the remainder, we could expect that the distributional error remains bounded by a quantity of order h . Unfortunately, it is not always true that p remains bounded. 

As a simple illustration, consider the harmonic oscillator with Hamiltonian H(q,p) = p /2 - - (ill and the invariant distributions obtained using several numerical methods. Applying Forward Euler to the harmonic oscillator, we have 

As an alternative, consider the Backward Euler method q + = q + hp +i, p +i = p — hq +i. Then it is easily shown that both eigenvalues of the matrix lie in the interior of the unit disk, and hence all solutions of the recurrence relation tend to the origin with increasing n (the origin is an attractive equilibrium point). In this case, all densities evolve toward the Dirac distribution centered at the origin (5[ ](5 p]. The only (distributional) solution of Cjp = 0 is again 3[( ]5[p], which in this case is attractive. 

---