Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ode45 problems

Subject: ode45 problems

From: Peter

Date: 3 Jul, 2011 22:17:10

Message: 1 of 4

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.

Subject: ode45 problems

From: Roger Stafford

Date: 3 Jul, 2011 23:15:10

Message: 2 of 4

"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 53-bit 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

Subject: ode45 problems

From: Peter

Date: 4 Jul, 2011 03:07:09

Message: 3 of 4

Thanks very much Roger for the full response.

That makes sense, and I will have a play with those options.

The reason I am running many cycles is to compare different initial conditions both in the time and the frequency domain. This might be unrealistic due to the issues you point out.

Subject: ode45 problems

From: Peter

Date: 4 Jul, 2011 05:08:07

Message: 4 of 4

I have used RelTol=1e-6 & AbsTol = 1e-4 - excellent results (over 500secs) including crystal clear spectral plots.

Thanks again.

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us