So it’s holiday time and I find some times to write, my objective today would be to present and idea regarding how to use the first and the second Dahlquist barrier theorem to prove some theoretical limit of symplectic integrator.
In this first section we will focus our attention to the numerical method to solve the wave equation in one dimension with Cauchy initial condition, what we mean when saying this is that we are searching a function at least such that:
where , and are given and is not necessarily the derivative of . In particular we will call homogeneous wave initial value problem (IVP) the previous equation with , and we will begin by studying it.
Before dealing directly with the wave equation is interesting to introduce briefly the characteristic method. This is a resolution technique for the first order parabolic partial differential equation (PDE), but we will see that can be applied as well to the wave equation. Let’s consider the advocation equation:
we fix the point and define the function . If we define and we can compute:
this tell’s us that the solution to the ordinary differential equation (ODE) is equivalent to solving the advocation equation. If we know that is constant and therefore which is equivalent to ask and therefore:
Now if instead considering the advocation equation alone we consider the following IVP:
we have the solution to the avocation and is given by the characteristic line equation:
this tell’s as well that the solution to the avocation equation is constant along a line, that we will call characteristic line and has equation , for any .
Now let’s go back to the homogeneous wave IVP, and we can observe that the one dimensional wave equation can be decomposed as follow:
and therefore we know that the solution of the wave equation is as well solution of both forward and backward advocation equation. Therefore the characteristic line of the wave equation are both and . Now our will be to compute to prove that is constantly null:
and if we substitute in the last equation we get that . Now we can integrate to obtain a nice expression for :
where and are two function. Now let’s observe that since we have that and is as well possible to write down:
therefore combining what we have just seen we get the following result:
this last equation give as an analytical solution to the wave equation and is known as the d’Alembert formula. Analogously a formula for the inhomogeneous wave IVP can be obtained, provided that is at least :
Let us focus our attention to the d’Alembert formula, when and are “nice” function the analytical solution can be computed easily by hand but when the function aren’t nice we can use a wide number of numerical integration technique and quadrature formula to compute the integral in d’Alembert formula. For lack of a better name we will call here this class of solvers d’Alembert solvers.
We will here show the result of implementing a d’Alembert solvers that computes integral in the d’Alembert formula using equi-spaced node and the trapezoid quadrature formula, but is important to notice that far better result can be obtained using higher order Newton-Cotes methods, Gaussian quadrature etc.
We will here try to explain a different approach, known with the name of time stepping. Let’s define an equi-spaced spatial mesh , the time stepping method consist in approximating numerically using a stiffness matrix , ie:
where is a discretisation of the according to the spatial approximation chosen. In this report we will only treat the finite difference as a spatial approximation method, more information regarding finite difference can be found in the second part of this report. Now once we observe that still depends upon time, ie is easy to see that the IVP equation becomes:
which is am typical numerical ODE problem of second order. One might think about splitting the second order ODE system into a system of two first order ODE to be solved coupled using method such as Euler, Crank-Nicolson, Runghe-Kutta etc etc. This certainly is a valid idea that we will explore in the third part of our report.
Our aim here will be to address the problem of solving the IVP using second order ODE integration technique. In particular this section focuses on the Leapfrog integration. Let’s introduce to begin with an equi-spaced time mesh of time step and adopt the following notation , to define an other mesh of step such that . Now we are ready to define the Leapfrog integration scheme:
If we locate the Leapfrog integration on the mesh by observing that since then , and that the following Taylor expansions holds:
we get the following equation to describe the Leapfrog integration:
Our aim is here to study the stability of the Leapfrog integration using a the theory developed in B1. Our aim will be to see Leapfrog integration when solving the second order differential equation:
as way of solving the system of differential equations:
where . We know that equation that characterize the leapfrog method, and we can see the second equation as being the one solving the second equation of the first order formulation while the third equation solves the first equation of the first order formulation.
In particular is important to notice that while the second equation of the first order formulation dosen’t present particular problems the first one can be reformulated to obtain component by component a set of differential equation of the form .
The idea behind this that if we suppose our stiffness matrix is diagonalizable we can rewrite as:
where is the digitalization of by the matrix , and this means that we can study instead of system of equation that describe the first order formulation the following formulation:
now we rewrite the Leapfrog scheme equations as a linear multi step method (more information can be found in B1):
supposing that the eigenvalues of are negative we get that both the linear multi step method respect the roots criteria, since they have characteristic polynomials:
and both polynomials have roots in the unit disk.
Therefore the Leapfrog scheme is 0-stable but what we are interested to check now is if it is as well absolutely stable, as defined in B1 which is different from being -stable, as defined by Hairer-Lubich and as well in B1.
To study the Leapfrog scheme we will notice that the first equation of that describe the Leapfrog scheme can be obtained by applying the explicit Euler scheme to the problem:
and under suitable hypothesis on , which we are not going to present in detail here, we need for the solution not to explode with that:
Last but not least we would like to observe that if as stiffness the finite difference matrix in 1D (more information can be found in the second part of this report) we know that its eigenvalues are:
is clear that is in the negative portion of the complex plane and the second requirement mentioned in the definition of -stability is tight if which is equivalent to ask:
which turns out to be the Courant–Friedrichs–Lewy condition.
In this section our aim is to present a wide class of methods, among which we could find as well the Leapfrog integration scheme. The class of methods here presented is known as sympletic integrator and have the peculiar property of conserving some sort of numerical energy.
If we image that describes the motion of particles in space, we introduce the Hamiltonian of this system and Hamilton’s equations:
now we would like to integrate numerically the above mentioned system of equation, this means that the numerical solution of the above system will be associated to an perturbed Hamiltonian, which is conserved since solution to the Hamilton equations conserves .
We will here present the same approach as in B2, which is base upon the idea that we can split the Hamiltonian as follow:
now we introduce the vector and observe that Hamilton equations can be rewritten as:
where are the Poisson brackets. We now introduce a notation useful in future:
rewriting those equation as we can easily find a solution to the ODE using the exponential matrix, since which is a property that comes directly from the assumption that :
we now wont to approximate the above mentioned operator as follow:
if the problem is quick to solve and painless , for the coefficients we are searching for are and . In particular in B3 solutions until are presented, while in B2 the Baker-Campbell-Hausdorff formula is used to a compute and for any symmetric symplectic integrator.
If we perform a Taylor expansion of the exponentials involved in the previous equations we get that:
we end up with the following equations:
now if we substitute this recursively in the previous equations we get:
which in Lagrangian coordinates becomes:
and those equations defines recursively a symplectic integrator. For example if we have as coefficents we have symplectic Euler integration which is defined by the following equations:
Instead if we choose as coefficients we end up with the Verlet integration scheme, which is defined by the following equations:
Last not but least we find out the Leapfrog integration if we choose .
We will here call symplectic integrator any integrator built with the technique presented, that can be found as well in B2 and B3.
Using the same argument proposed in the previous section, if we rewrite a symplectic scheme as a linear multi step scheme is possible to use the Dahlquist’s second barrier we can find out that no symplectic method can be -stable or -stable.
In particular we now know that any symplectic integrator as a multi step linear method, therefore we can use some of the instruments developed for the stability of multi step linear method with regard to symplectic integrator. The instruments we are speaking about are the first and the second Dahlquist Barrier:
First Dahlquist Barrier There are no multi step linear method, consisting of step which have a order of convergence if is odd and of if is even.
Second Dahlquist Barrier There are no explicit multi step linear method which are -stable or -stable. Further more there is no multi step linear method which is -stable and has order of convergence greater then .
In particular we’d like to formulate the following theorem, that will be based upon the first and the second Dahlquist barrier.
Theorem There are no symplectic integrator obtained by iteration of the Negri technique that have a order of convergence if is odd and of if is even. Further more there is no symplectic integrator that is -stable and has order of convergence greater then .
Proof We will show here just a sketch of the proof that is based on what we have seen previously:
- We reduce a symplectic integrator to a system of equations representing a multi step linear method, as we have done for the Leapfrog integration.
- We are justified to study equation by equation the stability of the symplectic integrator, as we proven in the Leapfrog Integration paragraph.
- We apply the second and the first Dahlquist barrier results to each equation representing a multi step linear method to prove the thesis above mentioned.
B1. Quarteroni, Alfio, Riccardo Sacco, and Fausto Saleri. Matematica numerica. Springer Science & Business Media, 2008.
B2. Yoshida, Haruo. “Construction of higher order symplectic integrators.” Physics letters A 150.5-7 (1990): 262-268.
B3. Negri, F. “Lie algebras and canonical integration.” Department of Physics, University of Maryland, prepffnt (1988).