Energy Analysis

Beginning from the compact representation of the difference scheme given in (3.18), one may derive a discrete conserved energy similar to that obtained for the continuous problem in Section 3.1.2. Multiplying (3.18) by a discrete approximation to the velocity, $\delta_{t\cdot}u$ gives

$$\delta_{t\cdot}u\delta_{tt}u +\omega_{0}^2(\delta_{t\cdot}u)u = 0$$ (3.34)

Employing the product identities given in Section 2.4.2 in the last chapter, one may immediately arrive at

$$\delta_{t+}\left(\frac{1}{2}\left(\delta_{t-}u\right)^2+\frac{1}{2}ue_{t-}u\right)=0$$ (3.35)

or

$$\delta_{t+}{\mathfrak{h}} = 0\qquad{\mathfrak{h}} = {\mathfrak{t}}+{\mathfrak{v}}$$ (3.36)

with

$${\mathfrak{t}} = \frac{1}{2}\left(\delta_{t-}u\right)^2\qquad{\mathfrak{v}} = \frac{\omega_{0}^2}{2}ue_{t-}u$$ (3.37)

${\mathfrak{t}}$ and ${\mathfrak{v}}$ are clearly approximations to the kinetic and potential energies, ${\mathfrak{T}}$ and ${\mathfrak{V}}$ respectively, for the continuous problem, as defined in (3.10).

Equation (3.36) implies that the discrete time series ${\mathfrak{h}}^{n}$, which approximates the total energy ${\mathfrak{H}}$ for the continuous time problem, is constant. Thus the scheme (3.18) is energy conserving, in a discrete sense. In other words,

$${\mathfrak{h}}^{n} = {\rm constant}={\mathfrak{h}}^{0}$$ (3.38)

Though one might think that the existence of a discrete conserved energy would immediately imply good behaviour of the associated finite difference scheme, just as in the continuous time case, it is important to note that the discrete approximation to potential energy, given above in (3.37), is not necessarily positive, and thus neither is ${\mathfrak{h}}$. The determination of positivity conditions on ${\mathfrak{h}}$ leads directly to a numerical stability condition. Expanding out the operator notation, one may write the function ${\mathfrak{h}}^{n}$ at time step $n$ as

$${\mathfrak{h}}^{n} = \frac{1}{2k^2}\left((u^{n})^2+(u^{n-1})^2\right)+\left(\frac{\omega_{0}^2}{2}-\frac{1}{k^2}\right)u^{n}u^{n-1}$$ (3.39)

which is no more than a quadratic form in the variables $u^{n}$ and $u^{n-1}$. Applying the results given in §2.4.3, one may show that the quadratic form above is positive definite when

$$k< \frac{2}{\omega_{0}}$$ (3.40)

which is exactly the condition (3.25) arrived at through frequency domain analysis. Under this condition, one further has that

$$\vert u^{n}\vert\leq \sqrt{\frac{2k^2{\mathfrak{h}}^n}{1-(1-\omega_{0}^2 k^2/2)^2}}=\sqrt{\frac{2k^2{\mathfrak{h}}^0}{1-(1-\omega_{0}^2 k^2/2)^2}}$$ (3.41)

Just as in the case of frequency domain analysis, the solution may be bounded in terms of the initial conditions (or the initial energy) and the system parameters. In fact, the bound above is identical to that given in (3.29). See Problem 3.11.

It is important to make some comments about the nature of this numerical energy conservation property. In infinite precision arithmetic, the quantity ${\mathfrak{h}}^{n}$ remains exactly constant, according to (3.38). That is, the discrete kinetic energy ${\mathfrak{t}}^{n}$ and potential energy ${\mathfrak{v}}^{n}$ sum to the same constant, at any time step $n$. See Figure 3.2(a). But in a finite precision computer implementation, there will be fluctuations in the total energy at the level of the least significant bit; such a fluctuation is shown, in the case of floating point arithmetic, in Figure 3.2(b). In this case, the fluctuations are on the order of $10^{-15}$ of the total energy, and the bit-quantization of these variations is clearly visible. Thus energy is conserved to machine accuracy.

Figure: (a) Variation of the discrete potential energy ${\mathfrak v}$ (solid grey line), discrete kinetic energy ${\mathfrak t}$ (dotted grey line) and total discrete energy ${\mathfrak h}$ (solid black line), plotted against time step $n$, for the output of scheme (3.18). In this case, the values $\omega _{0}=1$ and $k=1/10$ were used, and scheme (3.18) was initialized with the values $u^{0}=0.2$ and $u^{1} = -0.4$. (b) Variation of the error in energy, defined, at time step $n$, as ${\mathfrak h}_{e}^{n} = \left({\mathfrak h}^{n}-{\mathfrak h}^{0}\right)/{\mathfrak h}^{0}$, plotted as black points. Multiples of single bit variation are plotted as grey lines.

The analysis of conservative properties of finite difference schemes is younger than frequency domain analysis, and, as a result, much less well-known. It is, however, far more powerful than the body of frequency domain analysis techniques which form the workhorse of scheme analysis. The energy method [102] is the name given to a class of analysis techniques which do not employ frequency domain concepts. The study of finite difference schemes which possess conservation properties (of energy, and other invariants) experienced a good deal of growth in the 1980s [101,188] and 1990s [139,244,91], and is related to more modern symplectic integration techniques [189,199].