Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.
share|improve this question

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.

share|improve this answer
The wiki article on symplectic integration is unnecessarily technical for an intro to the subject; you could look up the Verlet method for starters. – Daniel Shapero Jun 25 '14 at 17:39

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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