Skip to main content

Energy conservation in a finite difference simulation of a simple harmonic oscillator

One key concept in physics is the conservation of energy---in an isolated system, the total energy cannot increase or decrease, it simply remains constant. Thus, the simple lossless harmonic oscillator model defined by the equation

$$ M\frac{d^2 u}{dt^2} = -K u $$

has an associated energy $H$ defined as

$$ H(t) = \frac{M}{2} \left( \frac{du}{dt} \right)^2 + \frac{K}{2} u^2, \qquad \frac{d}{dt} H =0,$$

which is conserved.

Such a concept can also be applied in a numerical simulation, where the existence of a numerical energy-like quantity with certain properties can become a powerful tool in guaranteeing the stability of the scheme.

Discrete operators

A finite difference scheme for the lumped oscillator equation above has been previously discussed here. Before jumping into the analysis of numerical energy, it is worth introducing some useful notation.

Given a discrete variable $u^n$ defined at integer multiples $n$ of a time step $k$, the following discrete operators can be defined:

$$ \delta_{tt} u^n = \frac{u^{n+1} - 2u^n + u^{n-1}}{k^2}, \qquad \delta_{t+} u^n = \frac{u^{n+1} - u^n}{k}, \qquad \delta_{t-} u^n = \frac{u^{n} - u^{n-1}}{k}. $$

These are just a shorthand for the second and first order difference operators. With this notation, the finite difference scheme for the SHO can be written in a much more compact way as

$$ M\delta_{tt} u^n = - K u^n. $$

Discrete energy

It is possible to show that the above scheme as an associated discrete energy quantity defined as

$$ H^{n+1/2} = \frac{M}{2} \left( \delta_{t+} u^n \right)^2 + \frac{K}{2} u^n \, u^{n+1} $$

which is conserved, such that

$$ \delta_{t-} H^{n+1/2} = 0. $$

Notice that the discrete energy is defined at interleaved steps to the displacement $u^n$. The conservation of energy in the discrete case holds assuming infinite precision arithmetic, but breaks down when finite precision representation of numbers (such as floating point) is used. Variations, however, are small, and confined to the last digits of the values, as will be seen shortly.

Numerical simulation

The calculation of the discrete energy can be easily implemented in the code given here. All it is necessary to do is to define a function that calculates the new value of $H^{n+1/2}$ at every time step. In order to simplify the implementation, it is possible to scale the value by $M$ and introduce the parameter $\omega_0^2 = K/M.$

In [1]:
function H_calc(u, u1, k, ω0)
    v = (u-u1)/k
    H = 0.5 * v^2 + 0.5 * ω0^2 * u * u1
    return H
H_calc (generic function with 1 method)

The finite difference code for the SHO with energy is the following.

In [2]:
function sho(f)

# Initial conditions

SR=44100      # Sampling rate
T=1           # Duration of the simulation (s)

# Calculate the time step
# Calculate the number of iterations

# Calculate ω0. This will be useful later

# Check the stability condition
if k>(2/ω0)
    println("Stability condition violated!")

# Initialise and allocate displacement
# Allocate energy vector
H = zeros(NF-1)
H[1] = H_calc(u1, u0, k, ω0)

# Start the Main Loop
for nn=3:NF
  # Update the recursion
  u = 2*u1 - ω0^2*k^2*u1 - u0

  # Store the new value for the position
  # Store the new value for the energy
  H[nn-1] = H_calc(u, u1, k, ω0)
  # Reset the known values of the recursion
  u0, u1 = u1, u;


  # At the end of the loop, return output
  return out, H
sho (generic function with 1 method)

Energy conservation

It is possible now to run a simulation and check if the discrete energy is conserved...

In [3]:
out, H = sho(5000)

using Gadfly

NF = length(H)
t = collect(0:NF-1)/NF
plot(x=t, y=H, Geom.line,
Guide.xlabel("time (s)"), Guide.ylabel("energy (J/kg)"), 
Guide.title("FD simulation - SHO"))
time (s)-1.5-1.0- (J/kg)FD simulation - SHO

As can be seen, the energy curve is not "dead flat", but presents some variations. By normalising the values, however, it is possible to see that such discrepancies are present only in the last few digits of the floating point representation of the numbers (remember that machine epsilon in double precision is roughly $2.22 \cdot 10^{-16}$).

In [4]:
Hnorm = eps(1.0)*(H - H[1])/eps(H[1])

plot(x=t, y=Hnorm, Geom.line,
Guide.xlabel("time (s)"), Guide.ylabel("energy (J/kg)"), 
Guide.title("FD simulation - SHO"))
time (s)-1.5-1.0-