I am simulating a damped harmonic oscillator using the RK4 method of numerical integration. I am comparing the simulated results with the analytical ones (for the free evolution case) and obtaining the error as a function of step-size, for different values of the quality factor Q.

$$\ddot{X} + \frac{1}{Q}\dot{X} + X = 0$$ The analytical solution is given by: $$x(\tau) = e^{-\tau / 2Q} \left[c_1 \cos\left(\tau \sqrt{1 - \frac{1}{4Q^2}}\right) + c_2\sin\left(\tau \sqrt{1 - \frac{1}{4Q^2}}\right)\right]$$ where $c_2$ and $c_2$ can be found by using the boundary conditions: $x(0) = 1$ and $\dot{x}(0) = 0$.

For RK4 method, we obtain two first order equations which we can use in the standard RK4 formula: Eqn1: $$\dot{X} = P \Rightarrow \dot{P} = \ddot{X} $$ Substitute into original equation to obtain: $$\Rightarrow \dot{P} ~+ \frac{1}{Q}P + X = 0$$ Which leads to Eqn2: $$\dot{P} = -\left(\frac{P}{Q}+X\right)$$

So these equations allow us to obtain the next point in $X$ and $dX$ for a particular step-size, $h$, using $$y_{n+1} = y_n + h\left(\frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6}\right)$$

Error Investigation

I am interested in the relative error for different step-sizes and Quality factors.

$$Err_{relative} = \frac{\Sigma_i (x_i - X_i)^2}{\Sigma_i x_i^2}$$

We can logically expect the error to increase as the step-size increases, which is what my graphs show. I run this simulation for different values of Q.
I would expect the error to decrease as Q increases, but instead the trajectory is just shifted upwards, i.e: the error is increased but the behaviour with respect to the step size is what we should expect.

enter image description here

It might be hard to see on the graph but those trajectories are for Q = 5 and Q = 95. The trajectory for Q = 5 is lower than the other one. I would have expected it to be the other way round.

Assuming I haven't goofed up my code, are my assumptions correct or is this normal behaviour?

  • $\begingroup$ Sorry, I don't have time for a detailed answer right now. But your intuition is not generally correct (that is, it may be right or wrong, depending on the integrator and the step size chosen). To be brief, your discretization uses a polynomial approximation to the exponential function and the size of the error depends on how well that polynomial approximates the exponential at the eigenvalues of the ODE system. $\endgroup$ – David Ketcheson Nov 27 '14 at 1:51
  • $\begingroup$ See the stability region diagrams in staff.science.uu.nl/~frank011/Classes/numwisk/ch10.pdf (RK4 is one of them). The darkness of the gray-plot within the stability region is the magnitude of the error amplification. The darker the color level for your specific eigen-value times step-size product the less the damping of the truncation error of the method. If the product lies outside the stability region the error explodes and the method is not useable. $\endgroup$ – Tobias Nov 30 '14 at 17:36

Your Answer

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

Browse other questions tagged or ask your own question.