http://www.mathworks.com/matlabcentral/newsreader/view_thread/310048
MATLAB Central Newsreader  ode45 problems
Feed for thread: ode45 problems
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Sun, 03 Jul 2011 22:17:10 +0000
ode45 problems
http://www.mathworks.com/matlabcentral/newsreader/view_thread/310048#843929
Peter
I'm new to using ode functions in matlab. I've got the basic stuff working but have found something strange.<br>
<br>
Here is my core function:<br>
<br>
============<br>
function dxdt = harm1(tr,xr)<br>
% Simple harmonic oscillator from Thompson & Stewart<br>
dxdt = [xr(2); 4*pi^2*sin(xr(1))];<br>
end<br>
=============<br>
<br>
It is an undamped, unforced pendulum. The result should maintain the same amplitude (???) I have 2 results from 2 initial conditions.<br>
<br>
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).<br>
<br>
The starting conditions are x=1.5, dx/dt=0; x=1.55, dx/dt=0<br>
<br>
The amplitude of condition 1 at the start is, of course, 1.5  but at 500 secs is 0.93.<br>
The amplitude of condition 2 at the start is 1.55  and at 500 secs is 1.59<br>
<br>
My phase plot (x vs dx/dt) shows clear gradual decline in amplitude of condition 1 and clear gradual increase of condition 2.<br>
<br>
Is this an ode45 issue? Is there something I need to do with long time series?<br>
<br>
Any help would be appreciated, thanks.

Sun, 03 Jul 2011 23:15:10 +0000
Re: ode45 problems
http://www.mathworks.com/matlabcentral/newsreader/view_thread/310048#843932
Roger Stafford
"Peter " <pgillies3@gmail.com> wrote in message <iuqpp6$917$1@newscl01ah.mathworks.com>...<br>
> I'm new to using ode functions in matlab. I've got the basic stuff working but have found something strange.<br>
> Here is my core function:<br>
> <br>
> ============<br>
> function dxdt = harm1(tr,xr)<br>
> % Simple harmonic oscillator from Thompson & Stewart<br>
> dxdt = [xr(2); 4*pi^2*sin(xr(1))];<br>
> end<br>
> =============<br>
> <br>
> It is an undamped, unforced pendulum. The result should maintain the same amplitude (???) I have 2 results from 2 initial conditions.<br>
> 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).<br>
> The starting conditions are x=1.5, dx/dt=0; x=1.55, dx/dt=0<br>
> The amplitude of condition 1 at the start is, of course, 1.5  but at 500 secs is 0.93.<br>
> The amplitude of condition 2 at the start is 1.55  and at 500 secs is 1.59<br>
> My phase plot (x vs dx/dt) shows clear gradual decline in amplitude of condition 1 and clear gradual increase of condition 2.<br>
> Is this an ode45 issue? Is there something I need to do with long time series?<br>
> Any help would be appreciated, thanks.<br>
          <br>
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.<br>
<br>
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.<br>
<br>
By the way, I would not call a pendulum a "simple harmonic oscillator". That terminology, in my opinion, is reserved for the case:<br>
<br>
dxdt = [xr(2); 4*pi^2*xr(1)];<br>
<br>
in which the solution is a pure sine function.<br>
<br>
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.<br>
<br>
Roger Stafford

Mon, 04 Jul 2011 03:07:09 +0000
Re: ode45 problems
http://www.mathworks.com/matlabcentral/newsreader/view_thread/310048#843939
Peter
Thanks very much Roger for the full response.<br>
<br>
That makes sense, and I will have a play with those options.<br>
<br>
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.

Mon, 04 Jul 2011 05:08:07 +0000
Re: ode45 problems
http://www.mathworks.com/matlabcentral/newsreader/view_thread/310048#843944
Peter
I have used RelTol=1e6 & AbsTol = 1e4  excellent results (over 500secs) including crystal clear spectral plots.<br>
<br>
Thanks again.