"Peter " <pgillies3@gmail.com> wrote in message <iuqpp6$917$1@newscl01ah.mathworks.com>...
> I'm new to using ode functions in matlab. I've got the basic stuff working but have found something strange.
> Here is my core function:
>
> ============
> function dxdt = harm1(tr,xr)
> % Simple harmonic oscillator from Thompson & Stewart
> dxdt = [xr(2); 4*pi^2*sin(xr(1))];
> end
> =============
>
> It is an undamped, unforced pendulum. The result should maintain the same amplitude (???) I have 2 results from 2 initial conditions.
> As I extend the time series (e.g. to 500 secs) I find that one of the results shows a very significant reduction in amplitude, and the other shows a small increase in amplitude. (The reason for running a long time series  apart from just curiosity  is to get a clear fft plot).
> The starting conditions are x=1.5, dx/dt=0; x=1.55, dx/dt=0
> The amplitude of condition 1 at the start is, of course, 1.5  but at 500 secs is 0.93.
> The amplitude of condition 2 at the start is 1.55  and at 500 secs is 1.59
> My phase plot (x vs dx/dt) shows clear gradual decline in amplitude of condition 1 and clear gradual increase of condition 2.
> Is this an ode45 issue? Is there something I need to do with long time series?
> Any help would be appreciated, thanks.
          
What you are seeing is undoubtedly an accumulative effect of errors from the finite difference approximations being made, as well as of round off errors. If you are going to be running this through such long periods as 500 seconds, you should most certainly set the Reltol and Abstol tolerances in the options setting at correspondingly small amounts. The tighter you set them, the smaller the errors from finite differencing. Remember, the ode solvers for all their sophistication, must use finite differences to approximate derivatives.
Of course there is an absolute limit imposed on the accuracy from inevitable round off errors, and there is nothing you can do about that, given the 53bit accuracy of the double precision numbers that are being used.
By the way, I would not call a pendulum a "simple harmonic oscillator". That terminology, in my opinion, is reserved for the case:
dxdt = [xr(2); 4*pi^2*xr(1)];
in which the solution is a pure sine function.
Finally, I question the need for running this particular computation past a single cycle. The theoretical solution is some periodic function and what is there to be gained by going through the computation more than one period, once you have determined what that period is? Only if you introduced, for example friction, would there be a call for many repeated cycles to see how fast the amplitude dies.
Roger Stafford
