# Energy of damped harmonic oscillator begins to increase with very large Q in numerical integration

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.

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.

How are you integrating the equation? It is notoriously difficult, if not impossible, to numerically integrate Newtonian equations of motion and conserve energy. The need to perform operations sequentially rather than simultaneously introduces errors. Some techniques exist to reduce the problem. This is also happening in your low Q examples, but it's masked by the much larger effect of the decay that you are simulating. The fact that you see oscillations indicates (to me, at least) that your integration technique is not ideal for this problem. In other words, all of these things are artifacts of the discrete approximation to a continuous phenomenon. There was a post very recently on this topic, with some additional info; unfortunately, I haven't found it.

• I am integrating the equation by letting the $\dot{X} = P$ and thus getting two first order DEs. Then I am simply using discrete steps to find the next point: $\dot{X} = (X_{n+1} - X_n) / \Delta t$. I do the same for the P equation. My code operates in a loop which calculates values for P and X for about 10000 iterations. Would you like me to post my code in my question? – turnip Jun 24 '14 at 13:30
• No, that's fine. It will not conserve energy. For starters, see this Wikipedia page and this SE question. . You might try the Runge-Kutta method, which may help and is not difficult to implement. – garyp Jun 24 '14 at 13:37
• My plan is to use the Runge-Kutta 4 method later. (My task at the moment is to use this method). Can you elaborate on the energy conservation bit? – turnip Jun 24 '14 at 13:43
• I'm not an expert, and my applications are not critical. I simply reduce the step size until the results don't change over the time of interest, and energy is reasonably conserved, or I look into other integration methods. Steps that are too small might lead to round-off / floating point errors. I do usually use RK4 for starters. – garyp Jun 24 '14 at 14:26

I think you're seeing two types of errors here.

For low values of $Q$, there is a lot of damping, and thus the values of $\phi$ are low (and decreasing).

I'm not sure what the absolute values of $X$ are, but it might be that they become sufficiently small to make machine errors relevant. For higher values of $Q$, higher values of $X$ arise, and integration errors become an issue.

Do you see the same relationship between the error and the time step for all values of $Q$? For example, if you decrease the time step, does the error for low $Q$ persist?