Midpoint method

The Implicit Midpoint method which, for the ODE z = f(z) takes the form 

This indieates that the Implicit Midpoint method has no stability threshold for the stepsize when solving a linear system This remarkable property may inspire some optimism that Implicit Midpoint would allow much larger stepsize in molecular dynamics simulation than does Verlet. 

In aetual practice, this does not work out quite as expected. To understand why not, we need to recognize that the linear stability condition only describes the 

If the Implicit Midpoint method is viewed as the exact solution of a perturbed harmonic oscillator, then the frequency of this oscillator is 

In the diagram below (Fig. 4.2), we plot the ratio against/t and S2. If this ratio is 1, the frequency is accurately represented. The figure shows that the frequency is only well approximated in case hQ is sufficiently small. For hQ 1 the frequency is highly modified. From this we might infer that a molecular dynamics simulation performed using the Implicit Midpoint method is only likely to be accurate in terms of resolution of temporally correlated quantities in case h 2 is relatively small. 

---

A popular implicit discretization is the (implicit) midpoint method [7] which, applied to a system of the type... 

W, g potential functions, k 1, has been discussed in various papers (see, for example, [6, 11, 9, 16, 3]). It has been pointed out that, for step-sizes /j > e = 1/ /k, the midpoint method can become unstable due to resonances [9, 16], i.e., for specific values of k. However, generic instabilities arise if the step-size k is chosen such that is not small [3, 6, 18], For systems with a rotational symmetry this has been shown rigorously in [6j. This effect is generic for highly oscillatory Hamiltonian systems, as argued for in [3] in terms of decoupling transformations and proved for a linear time varying system without symmetry. 

In general, the solution components of the DAE (4) are the correct limits (as K —> oo) of the corresponding slowly varying solution components of the free dynamics only if an additional (conservative) force term is introduced in the constrained system [14, 5]. It turns out [3] that the midpoint method may falsely approximate this correcting force term to zero unless k — 0 e), which leads to a step-size restriction of the same order of magnitude as explicit... 

The purpose of this paper is twofold (i) We summarize possible difficulties with the midpoint method (other than resonance instability, which has been treated extensively elsewhere) by looking at a simple (molecular) model problem, (ii) We investigate the suitability of some energy conserving methods. 

Time-reversible energy conserving methods can be obtained by appropriate modifications to the (time-reversible) midpoint method. Two such modifications are (i) scaling of the force field by a scalar such that total energy... 

Note that this latter method differs from the midpoint method, where one would use r(q +i/2) = (qn+i + qn)/2 instead of (7c) for r +i/2 in (7b). For highly oscillatory systems with k e, this can be a significant difference, because r is discretized directly in (7). An example in 4 below shows that the midpoint method can become unstable while (7) and (6) remain. stable. 

Here a symmetric projection step is used to enforce conservation of energy. Let a(g,p) and b q,p) be two vector-valued functions such that (p a q,p) + U q) b q,p)) is bounded away from zero. Then we propose the following modified midpoint method,... 

This latter modified midpoint method does work well, however, for the long time integration of Hamiltonian systems which are not highly oscillatory. Note that conservation of any other first integral can be enforced in a similar manner. To our knowledge, this method has not been considered in the literature before in the context of Hamiltonian systems, although it is standard among methods for incompressible Navier-Stokes (where its time-reversibility is not an issue, however). 

Note that, in loeal eoordinates. Step 2 is equivalent to integrating the equations (13). Thus, Step 2 can either be performed in loeal or in eartesian coordinates. We consider two different implicit methods for this purpose, namely, the midpoint method and the energy conserving method (6) which, in this example, coineides with the method (7) (because the V term appearing in (6) and (7) for q = qi — q2 is quadratie here). These methods are applied to the formulation in cartesian and in local coordinates and the properties of the resulting propagation maps are discussed next. 

The interesting question is now what happens if the midpoint method is applied to the cartesian formulation (9) with Vi = V2 = 0. The equations are... 

To obtain the unconditional stability of the midpoint method in local coordinates, one would have to consider the decoupling transformation from cartesian to local coordinates for the discrete variables etc. But this transformation, which for the continuous variables is not constant, necessarily is in error which depends on k, not e. The stability properties of the discrete dynamical systems obtained by the midpoint discretization in the different sets of coordinatc.s may therefore be significantly different when it 3> e [3]. 

On the other hand, our computations using (15) indicate that the midpoint method becomes unstable for a > 1. 

Unfortunately, discretization methods with large step sizes applied to such problems tend to miss this additional force term [3]. Furthermore, even if the implicit midpoint method is applied to a formulation in local coordinates, similar problems occur [3]. Since the midpoint scheme and its variants (6) and (7) are basically identical in local coordinates, the same problem can be expected for the energy conserving method (6). To demonstrate this, let us consider the following modified model problem [13] ... 

For the impact assessment, the midpoint method EDIP97 [15] was used, and the following impact categories were included ... 

Boundary value problems in cylindrical and spherical coordinates have an inherent singularity at x = 0. These problems can be tackled using Maple s inbuilt midpoint methods. For example, diffusion of a substrate in an enzyme catalyzed reaction.[6] The governing equation for the dimensionless concentration is... 

Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead... 

Maple identifies the singularity at x = 0 and suggests the midpoint method ... 

A major disadvantage with higher order Runge-Kutta methods is that the derivative must be calculated many times per time step, making these methods more expensive than midpoint methods of comparable order. However, the Runge-Kutta methods are more stable than the multistep methods of the same order and need no data other than the initial condition required by the differential equation itself (self-starting). 

Jll 21 22 — fl22fll2. 10. Show that the implicit midpoint method... 

Fig. 4.2 The frequency ratio 12/i2 for the Implicit Midpoint method. A frequency ratio of 1 indicates that the phase of the harmonic oscillator is accurate...

Implicit Midpoint Preserves First Integrals. Show that the implicit midpoint method exactly preserves quadratic first integrals. That is, if we have a matrix A ... 

Estimates of Tg obtained by the midpoint method for the three heating rates of Fig.3 are, respectively, 7K, IIK and 21K above the midpoint of the Cp vs. T plot. Although it is not shown here. 

The midpoint method is an improvement on the Euler method in that it uses information about the function at a point other than the initial point of the interval. As seen in Figure 6.4, the values of x andy at the midpoint are used to calculate the step across the whole interval, h. 

Figure 6.5. Solution using the midpoint method, step size h = 0.2.

The midpoint method requires two function evaluations in the step,... 

Consequently, the midpoint method increases accuracy by one order. Obviously the increased order of accuracy leads to a better approximation than the Euler solution, as shown in Figure 6.5. 

---