Consider the fundamental model problem for simple harmonic oscillations, $$ \begin{equation} u'' + \omega^2u = 0,\quad u(0)=I,\ u'(0)=0, \tag{362} \end{equation} $$ where \( \omega \) is the frequency of the oscillations (the exact solution is \( u(t)=I\cos\omega t \)). This model can equivalently be formulated as two first-order equations, $$ \begin{align} v' &= -\omega^2 u, \tag{363} \\ u' &= v\tp \tag{364} \end{align} $$ The popular Euler-Cromer scheme for this \( 2\times 2 \) system of ODEs applies an explicit forward difference in (363) and a backward difference in (364): $$ \begin{align} \frac{v^{n+1}- v^n}{\Delta t} &=- \omega^2u^{n}, \tag{365}\\ \frac{u^{n+1} - u^n}{\Delta t} &= v^{n+1}\tp \tag{366} \end{align} $$ For a time domain \( [0,T] \), we have introduced a mesh with points \( 0=t_0 < t_1 < \cdots < t_n=T \). The most common case is a mesh with uniform spacing \( \Delta t \): \( t_n=n\Delta t \). Then \( v^n \) is an approximation to \( v(t) \) at mesh point \( t_n \), and \( u^n \) is an approximation to \( u(t) \) at the same point. Note that the backward difference in (368) leads to an explicit updating formula for \( u^{n+1} \) since \( v^{n+1} \) is already computed: $$ \begin{align} v^{n+1} &= v^n -\Delta t \omega^2u^{n}, \tag{367}\\ u^{n+1} &= u^n + \Delta t v^{n+1}\tp \tag{368} \end{align} $$
The Euler-Cromer scheme is equivalent with the standard second-order accurate scheme for (362): $$ \begin{equation} u^{n+1} = 2u^n - u^{n-1} - \Delta t^2\omega^2 u^n,\ n=1,2,\ldots, \tag{369} \end{equation} $$ but for the first time step, the method for (362) leads to $$ \begin{equation} u^1= u^0 - \frac{1}{2}\Delta t^2\omega^2u^0, \tag{370} \end{equation} $$ while Euler-Cromer gives $$ \begin{equation} u^1= u^0 - \Delta t^2\omega^2u^0, \tag{371} \end{equation} $$ which can be interpreted as a first-order, backward difference approximation of \( u'(0)=0 \) combined with (369). At later time steps, however, the alternating use of forward and backward differences in (367)-(368) leads to a method with error \( \Oof{\Delta t^2} \).
(hpl 13: Do the equations in different sequence, first vel, then pos.)
The fact that the forward and backward differences used the Euler-Cromer method yield a second-order accurate method is not obvious from intuition. A much more intuitive discretization employs solely centered differences and leads to a scheme that is equivalent to the Euler-Cromer scheme. It is in fact fully equivalent to the second-order scheme for (362), also for the first time step. This alternative scheme is based on using a staggered (or alternating) mesh in time.
In a staggered mesh, the unknowns are sought at different points in the mesh. Specifically, \( u \) is sought at integer time points \( t_n \) and \( v \) is sought at \( t_{n+1/2} \) between two \( u \) points. The unknowns are then \( u^1, v^{3/2}, u^2, v^{5/2} \), and so on. We typically use the notation \( u^n \) and \( v^{n+\half} \) for the two unknown mesh functions.
On a staggered mesh it is natural to use centered difference approximations, expressed in operator notation as $$ \begin{align} \lbrack D_t u &= v\rbrack^{n+\half}, \tag{372}\\ \lbrack D_t v &= -\omega u\rbrack^{n+1} \tp \tag{373} \end{align} $$ Writing out the formulas gives $$ \begin{align} u^{n+1} &= u^{n} + \Delta t v^{n+\half}, \tag{374}\\ v^{n+\frac{3}{2}} &= v^{n+\half} -\Delta t \omega^2u^{n+1} \tag{375} \tp \end{align} $$ Of esthetic reasons one often writes these equations at the previous time level to replace the \( \frac{3}{2} \) by \( \half \) in the left-most term in the last equation, $$ \begin{align} u^{n} &= u^{n-1} + \Delta t v^{n-\half}, \tag{376}\\ v^{n+\half} &= v^{n-\half} -\Delta t \omega^2u^{n} \tag{377} \tp \end{align} $$ Such a rewrite is only mathematical cosmetics. The important thing is that this centered scheme has exactly the same computational structure as the forward-backward scheme. We shall use the names forward-backward Euler-Cromer and staggered Euler-Cromer to distinguish the two schemes.
We can eliminate the \( v \) values and get back the centered scheme based on the second-order differential equation, so all these three schemes are equivalent. However, they differ somewhat in the treatment of the initial conditions.
Suppose we have \( u(0)=I \) and \( u'(0)=v(0)=0 \) as mathematical initial conditions. This means \( u^0=I \) and $$ v(0)\approx \half(v^{-\half} + v^{\half}) = 0, \quad\Rightarrow\quad v^{-\half} =- v^\half\tp$$ Using the discretized equation (377) for \( n=0 \) yields $$ v^\half = v^{-\half} -\Delta t\omega^2 I,$$ and eliminating \( v^{-\half} =- v^{\half} \) results in \( v^\half = -\half\Delta t\omega^2I \) and $$ u^1 = u^0 - \half\Delta t^2\omega^2 I,$$ which is exactly the same equation for \( u^1 \) as we had in the centered scheme based on the second-order differential equation (and hence corresponds to a centered difference approximation of the initial condition for \( u'(0) \)). The conclusion is that a staggered mesh is fully equivalent with that scheme, while the forward-backward version gives a slight deviation in the computation of \( u^1 \).
We can redo the derivation of the initial conditions when \( u'(0)=V \): $$ v(0)\approx \half(v^{-\half} + v^{\half}) = V, \quad\Rightarrow\quad v^{-\half} = 2V - v^\half\tp$$ Using this \( v^{-\half} \) in $$ v^\half = v^{-\half} -\Delta t\omega^2 I,$$ then gives \( v^\half = V - \half\Delta t\omega^2 I \). The general initial conditions are therefore $$ \begin{align} u^0 &= I, \tag{378}\\ v^\half &= V - \half\Delta t\omega^2I \tag{379}\tp \end{align} $$
The algorithm goes like this:
Translating the schemes (376)
and (377) to computer code
faces the problem of how to store and access \( v^{n+\half} \),
since arrays only allow integer indices with base 0.
We must then introduce a convention: \( v^{1+\half} \) is stored
in v[n] while \( v^{1-\half} \) is stored in v[n-1].
We can then write the algorithm in Python as
def solver(I, w, dt, T):
dt = float(dt)
Nt = int(round(T/dt))
u = zeros(Nt+1)
v = zeros(Nt+1)
t = linspace(0, Nt*dt, Nt+1) # mesh for u
t_v = t + dt/2 # mesh for v
u[0] = I
v[0] = 0 - 0.5*dt*w**2*u[0]
for n in range(1, Nt+1):
u[n] = u[n-1] + dt*v[n-1]
v[n] = v[n-1] - dt*w**2*u[n]
return u, t, v, t_v
Note that the return \( u \) and \( v \) together with the mesh points such
that the complete mesh function for \( u \) is described by u and t,
while v and t_v represents the mesh function for \( v \).
Some prefer to see a closer relationship between the code and
the mathematics for the quantities with half-integer indices.
For example, we would like to replace the updating equation for
v[n] by
v[n+half] = v[n-half] - dt*w**2*u[n]
This is easy to do if we could be sure that n+half means n and
n-half means n-1. A possible solution is to define half as a
special object such that an integer plus half results in the integer,
while an integer minus half equals the integer minus 1.
A simple Python class may realize the half object:
class HalfInt:
def __radd__(self, other):
return other
def __rsub__(self, other):
return other - 1
half = HalfInt()
The __radd__ function is invoked for all expressions n+half
("right add" with self as half and other as n). Similarly,
the __rsub__ function is invoked for n-half and results in n-1.
Using the half object, we can implement the algorithms in an even
more readable way:
def solver(I, w, dt, T):
"""
Solve u'=v, v' = - w**2*u for t in (0,T], u(0)=I and v(0)=0,
by a central finite difference method with time step dt.
"""
dt = float(dt)
Nt = int(round(T/dt))
u = zeros(Nt+1)
v = zeros(Nt+1)
t = linspace(0, Nt*dt, Nt+1) # mesh for u
t_v = t + dt/2 # mesh for v
u[0] = I
v[0+half] = 0 - 0.5*dt*w**2*u[0]
for n in range(1, Nt+1):
print n, n+half, n-half
u[n] = u[n-1] + dt*v[n-half]
v[n+half] = v[n-half] - dt*w**2*u[n]
return u, t, v, t_v
def test_staggered():
I = 1.2; w = 2.0; T = 5; dt = 2/w
u, t, v, t_v = solver(I, w, dt, T)
from vib_undamped import solver as solver2
u2, t2 = solver2(I, w, dt, T)
error = abs(u - u2).max()
tol = 1E-14
assert error < tol
if __name__ == '__main__':
test_staggered()
Verification of this code is easy as we can just compare the computed
u with the u produced by the solver function in
vib_undamped.py (which solves \( u''+\omega^2u=0 \) directly). The
values should coincide to machine precision since the two numerical
methods are mathematically equivalent. We refer to the file
vib_undamped_staggered.py
for the details of a nose test that checks this property.
The more general model for vibration problems, $$ \begin{equation} mu'' + f(u') + s(u) = F(t),\quad u(0)=I,\ u'(0)=V,\ t\in (0,T], \tag{380} \end{equation} $$ can be rewritten as a first-order ODE system $$ \begin{align} u' &= v, \tag{381} \\ v' &= m^{-1}\left(F(t) - f(v) - s(u)\right)\tp \tag{382} \end{align} $$ It is natural to introduce a staggered mesh (see the section The Euler-Cromer scheme on a staggered mesh) and seek \( u \) at mesh points \( t_n \) (the numerical value is denoted by \( u^n \)) and \( v \) between mesh points at \( t_{n+1/2} \) (the numerical value is denoted by \( v^{n+\half} \)). A centered difference approximation to (381)-(382) can then be written in operator notation as $$ \begin{align} \lbrack D_t u &= v\rbrack^{n-\half}, \tag{383} \\ \lbrack D_tv &= m^{-1}\left(F(t) - f(v) - s(u)\right)\rbrack^n\tp \tag{384} \end{align} $$ Written out, $$ \begin{align} \frac{u^n - u^{n-1}}{\Delta t} &= v^{n-\half}, \tag{385} \\ \frac{v^{n+\half} - v^{n-\half}}{\Delta t} &= m^{-1}\left(F^n - f(v^n) - s(u^n)\right)\tp \tag{386} \end{align} $$ With linear damping, \( f(v)=bv \), we can use an arithmetic mean for \( f(v^n) \): \( f(v^n)\approx = \half(f(v^{n-\half}) + f(v^{n+\half})) \). The system (385)-(386) can then be solved with respect to the unknowns \( u^n \) and \( v^{n+\half} \): $$ \begin{align} u^n & = u^{n-1} + {\Delta t}v^{n-\half}, \tag{387} \\ v^{n+\half} &= \left(1 + \frac{b}{2m}\Delta t\right)^{-1}\left( v^{n-\half} + {\Delta t} m^{-1}\left(F^n - {\half}f(v^{n-\half}) - s(u^n)\right)\right)\tp \tag{388} \end{align} $$
In case of quadratic damping, \( f(v)=b|v|v \), we can use a geometric mean: \( f(v^n)\approx b|v^{n-\half}|v^{n+\half} \). Inserting this approximation in (385)-(386) and solving for the unknowns \( u^n \) and \( v^{n+\half} \) results in $$ \begin{align} u^n & = u^{n-1} + {\Delta t}v^{n-\half}, \tag{389} \\ v^{n+\half} &= (1 + \frac{b}{m}|v^{n-\half}|\Delta t)^{-1}\left( v^{n-\half} + {\Delta t} m^{-1}\left(F^n - s(u^n)\right)\right)\tp \tag{390} \end{align} $$
The initial conditions are derived at the end of the section The Euler-Cromer scheme on a staggered mesh: $$ \begin{align} u^0 &= I, \tag{391}\\ v^\half &= V - \half\Delta t\omega^2I \tag{392}\tp \end{align} $$
We consider the generalized model with quadratic damping, expressed as a system of two first-order equations as in the section A staggered Euler-Cromer scheme for a generalized model: $$ \begin{align*} u^{\prime} &= v,\\ v' &= \frac{1}{m}\left( F(t) - \beta |v|v - s(u)\right)\tp \end{align*} $$ However, contrary to what is done in the section A staggered Euler-Cromer scheme for a generalized model, we want to apply the idea of a forward-backward discretization: \( u \) is marched forward by a one-sided Forward Euler scheme applied to the first equation, and thereafter \( v \) can marched forward by a Backward Euler scheme in the second equation, see in the section The Euler-Cromer method. Express the idea in operator notation and write out the scheme. Unfortunately, the backward difference for the \( v \) equation creates a nonlinearity \( |v^{n+1}|v^{n} \). To linearize this nonlinearity, use the known value \( v^n \) inside the absolute value factor, i.e., \( |v^{n+1}|v^{n}\approx |v^n|v^{n+1} \). Show that the resulting scheme is equivalent to the one in the section A staggered Euler-Cromer scheme for a generalized model for some time level \( n\geq 1 \).
What we learn from this exercise is that the first-order differences
and the linearization trick play together in "the right way" such that
the scheme is as good as when we (in the section A staggered Euler-Cromer scheme for a generalized model)
carefully apply centered differences and a geometric mean on a
staggered mesh to achieve second-order accuracy.
There is a
difference in the handling of the initial conditions, though, as
explained at the end of the section The Euler-Cromer method.
Filename: vib_gen_bwdamping.