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.
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?
