I have numerically integrated the (reduced) homogeneous equation of a damped harmonic oscillator in order to see how the error propagates.

$$\frac{d^2 X}{d\phi^2} + \frac{1}{Q}\frac{dX}{d\phi}+X(\phi) = 0$$ where $\phi = 2 \pi \tau$ and $\tau$ is the reduced time

In the process I have also extracted the logarithm of the energy and plotted it against time, for different values of the quality factor, Q. One can expect the energy to decay slower for larger Q.

$$ E = P^2 + X^2$$ where $P = \dot{X}$

What I found is that for very large values of Q, the energy of the oscillator begins to increase. I checked my code several times and I couldn't see anything wrong with it so I suspected that it purely something to do with numbers so I began playing around with them, but first here is the graph.

enter image description here

What I found is odd; say I'm using a step size of $1 \times 10^{-n}$ . If Q's magnitude matches or exceeds $10^{n}$ the energy begins to increase! What is causing this? What is the connection between Q and the step size?

Also, should the energy be oscillating? (If you zoom in on the red and green lines on the graph you will see that it is oscillating much slower there.)

Here is how I am integrating the equation for those interested:

  1. Let $\dot{X} = P$
  2. Thus original equation becomes $\dot{P} + \frac{1}{Q} P + X = 0$
  3. Now we have two 1st order DEs.
  4. Using discrete steps we can roughly find the next point of $X$ and $P$ by considering that $\dot{X} = \frac{dX}{d\phi} = \frac{X_{n+1} - X_n}{\Delta \phi}$ and approaching the other eqution in the same manner.
  5. We can then rearrange for $X_{n+1}$ and $P_{n+1}$ and create a code that calculates those points, given some basic initial conditions.

You are basically using Forward Euler timestepping, which is not stable. It is much better to use Backward Euler. Even better, if you care about the total energy of the system being conserved, you need to use a symplectic integration method.

  • 1
    $\begingroup$ The wiki article on symplectic integration is unnecessarily technical for an intro to the subject; you could look up the Verlet method for starters. $\endgroup$ Jun 25 '14 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.