Why does the following simple Harmonic oscillator (without damping) behave as a damped harmonic oscillator?

function resonance

omega = 100; %

b = 0.0; %

A = 0.0; % driving amplitude per unit mass

omega0 = 0; % driving frequency

tBegin = 0; % time begin

tEnd = 80; % time end

x0 = 0.2; % initial position

v0 = 0.8; % initial velocitie

a = omega^2; % calculate a coeficient from resonant frequency

% Use Runge-Kutta 45 integrator to solve the ODE

[t,w] = ode45(@derivatives, [tBegin tEnd], [x0 v0]);

x = w(:,1); % extract positions from first column of w matrix

v = w(:,2); % extract velocities from second column of w matrix

plot(t,x);

title('Damped, Driven Harmonic Oscillator');

ylabel('position (m)');

xlabel('time (s)');

    % Function defining derivatives dx/dt and dv/dt

    % uses the parameters a, b, A, omega0 in main program but changeth them not

    function derivs = derivatives(tf,wf)

        xf = wf(1);            % wf(1) stores x

        vf = wf(2);            % wf(2) stores v

        dxdt = vf;                                     % set dx/dt = velocity

        dvdt = -a * xf - b * vf + A * sin(omega0*tf);  % set dv/dt = acceleration

        derivs = [dxdt; dvdt];  % return the derivatives

    end

end

%

Answer 1

[David Goodmanson]

Hi Sameer,

I believe this is a standard numerical accuracy issue. Here is some shortened code that is basically the same as yours. If you run it as is, you get similar behavior to yours.

[t x] = ode45(@(t,x) fun(t,x), [0 2000], [1 0]);
plot(t,x(:,1))
function dxdt = fun(t,x)
dxdt = [0 1;-5 0]*x;
end

If you replace the first line with the two lines

 options = odeset('reltol',1e-5);
 [t x] = ode45(@(t,x) fun(t,x), [0 2000], [1 0],options);

then with the tighter tolerance the amplitude is basically constant all the way across. Of course it is going to take a bit longer. There is also an abstol you can vary, and it appears that the default values for reltol and abstol are 1e-3 and 1e-6. The odeset function without any arguments gives you a list of what all the optional parameters are.