Solving a simple harmonic oscillator--but the solution is decaying!

Hello,

When solving the harmonic oscillator equation, I am noticing a slow decay in the absolute value of my solution. I thought I was making sure I initialized the problem correctly, but the problem persists.

More concretely, I am solving for the solution of

(d/dt)^2 f = -w^2*f, where w^2 = k^2 + m^2.

The code defining my function handle for the ode solver is

function output = harmonic(t,f)
k=1;
m=1;
fdot=f(2);

fdotdot=-(k^2 + m^2)*f(1);

output=[fdot fdotdot];
output=output(:);
end

I have an additional condition that the Wronskian of f, f*conj(fdot) - conf(f)*fdot, has to remain at i. I initialized the problem by requiring that

f(0)=1/sqrt(2)*(k^2 + m^2)^.25,

and that

fdot(0) = -i*(k^2 + m^2)^.25 / sqrt(2).

If these two expressions are used to evaluate the Wronskian, it is indeed i. The equation of motion should only rotate the phase of f, so the Wronskian should be preserved. Moreover, the absolute value of f and fdot should be preserved as well.

However, when I run

[t f]=ode45(e,[-5 1000],input);

for example with this code, and plot the absolute value of f, I see that it is going down slowly! Why could this be?

Edit1: I have even taken sample code from websites, such as

http://msg.byu.edu/matlab_chaos.php

and run sho.m, for example. If I extend the time period, the decay in the amplitude of the solution is apparent. Extending it to about 15,000 makes it very obvious. Is this a problem with MATLAB itself? Can it not solve the harmonic oscillator?

Also I've played with the error tolerances, setting both the relative and absolute tolerances to lower than default settings. It didn't help.

[Arun]

In case anyone ever looks at this, the decay had to do with the type of odesolver employed. The solver uses the numerical decay to ensure stability of the solution. If you want to get around that, ode23 can be used, but that leads to other problems.

[Bryan Callaway]

I have observed a similar problem in Simulink. I find that ODE23 solver makes it workse, but 23t seems to solve the problem. For the simple submodel pictured below:

The output in the scope (with solver set to "auto") is as follows:

...With a slight but visible decay. Going through the solvers:

45: comparable to "auto".

23: Rapid decay - reaches zero in about 4 seconds.

113: Rapid decay - linear, should reach zero in about 30 seconds.

15S: Rapid growth; appears to multiply by an exponential function.

23 S: Nearly reaches opposite end of travel, reverses slightly early, then flatlines at a negative value for the remaining 10 seconds. In a sense, no decay, but nothing like the modelled behavior.

23t: expected result, an oscillation with no visible change in behavior over the ten second sampling period.

23tb: comparable to "auto".

[David Goodmanson]

Hi Arun,

this has nothing to do with the effect you are talking about, but w for the usual harmonic oscillator is sqrt(k/m) not sqrt(k^2+m^2)