Thread Subject:

ODE solver

Hi,

I've been working on a program using ODE's and I've been using ODE45 to get numerical solutions for differential equations. I've been getting incorrect results e.g. when I used
y'' = -y, I get a period that is ~ 2.8 instead of pi. This also happened with various other functions that I've been using.

Is there a solution to this problem?

Thanks!

W dniu 2013-06-28 22:38, William pisze:
> Hi,
>
> I've been working on a program using ODE's and I've been using ODE45 to
> get numerical solutions for differential equations. I've been getting
> incorrect results e.g. when I used
> y'' = -y, I get a period that is ~ 2.8 instead of pi. This also happened
> with various other functions that I've been using.
>
> Is there a solution to this problem?

I'm great wizard from quadruple-precision-in-matlab-land.
My crystal ball told me you have done something wrong[1]. Sadly,
this crystal ball is not HD ready, so I can't read you code:(
Maybe you can post it? :)

> Thanks!

[1]:
% solution of y'' = -y; y(1,:) is y, y(1,:) is velocity (y')

[t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0,10],[0,1]);

% of course ode45 solves first order equation, so we wrote
% y'' = -y as system {y'=u, u'=-y}

plot(t,y(:,1))
grid
%looks nice

fzero(@(x)interp1(t,y(:,1),x),pi)

% ans = 3.14134838301112
% meh, but significantly better than 2.8

options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
[t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0,10],[0,1],options );

root = fzero(@(x)interp1(t,y(:,1),x),pi)
(root-pi)/pi

% root = 3.14159265459964
% ans = 3.21445091288877e-010
% :-)

bartekltg <bartekltg@gmail.com> wrote in message <kqnu6a$k9n$1@node1.news.atman.pl>...
> W dniu 2013-06-28 22:38, William pisze:
> > Hi,
> >
> > I've been working on a program using ODE's and I've been using ODE45 to
> > get numerical solutions for differential equations. I've been getting
> > incorrect results e.g. when I used
> > y'' = -y, I get a period that is ~ 2.8 instead of pi. This also happened
> > with various other functions that I've been using.
> >
> > Is there a solution to this problem?
>
> I'm great wizard from quadruple-precision-in-matlab-land.
> My crystal ball told me you have done something wrong[1]. Sadly,
> this crystal ball is not HD ready, so I can't read you code:(
> Maybe you can post it? :)
>
> > Thanks!
>
> [1]:
> % solution of y'' = -y; y(1,:) is y, y(1,:) is velocity (y')
>
> [t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0,10],[0,1]);
>
> % of course ode45 solves first order equation, so we wrote
> % y'' = -y as system {y'=u, u'=-y}
>
> plot(t,y(:,1))
> grid
> %looks nice
>
> fzero(@(x)interp1(t,y(:,1),x),pi)
>
> % ans = 3.14134838301112
> % meh, but significantly better than 2.8
>
> options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
> [t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0,10],[0,1],options );
>
> root = fzero(@(x)interp1(t,y(:,1),x),pi)
> (root-pi)/pi
>
> % root = 3.14159265459964
> % ans = 3.21445091288877e-010
> % :-)
>
>

So my code looks something like this

[t,phi] = ode45(@diffeq,0:T/100:T,[1,0]);
%T is the period of the diff eq.

then

function phiprime = diffeq(t,phi)
m =1;
g=0;
l=1;
H = 0;
phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];

where i just let phi(1) = y, phi(2) = y'

I've tried this code with g, l, H = 0 which becomes y'' = -y
What exactly does the reltol, abstol do also.

Thanks!

W dniu 2013-06-30 10:32, William pisze:
> bartekltg <bartekltg@gmail.com> wrote in message
> <kqnu6a$k9n$1@node1.news.atman.pl>...
>> W dniu 2013-06-28 22:38, William pisze:
>> > Hi,
>> >
>> > I've been working on a program using ODE's and I've been using ODE45 to
>> > get numerical solutions for differential equations. I've been getting
>> > incorrect results e.g. when I used
>> > y'' = -y, I get a period that is ~ 2.8 instead of pi. This also
>> happened
>> > with various other functions that I've been using.
>> >
>> > Is there a solution to this problem?
>>
>> I'm great wizard from quadruple-precision-in-matlab-land.
>> My crystal ball told me you have done something wrong[1]. Sadly,
>> this crystal ball is not HD ready, so I can't read you code:(
>> Maybe you can post it? :)
>>
>> > Thanks!
>>
>> [1]:
>> % solution of y'' = -y; y(1,:) is y, y(1,:) is velocity (y')
>>
>> [t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0,10],[0,1]);
>>
>> % of course ode45 solves first order equation, so we wrote
>> % y'' = -y as system {y'=u, u'=-y}
>>
>> plot(t,y(:,1))
>> grid
>> %looks nice
>>
>> fzero(@(x)interp1(t,y(:,1),x),pi)
>>
>> % ans = 3.14134838301112
>> % meh, but significantly better than 2.8
>>
>> options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
>> [t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0,10],[0,1],options );
>>
>> root = fzero(@(x)interp1(t,y(:,1),x),pi)
>> (root-pi)/pi
>>
>> % root = 3.14159265459964
>> % ans = 3.21445091288877e-010
>> % :-)
>>
>>
>
> So my code looks something like this
>
> [t,phi] = ode45(@diffeq,0:T/100:T,[1,0]);
> %T is the period of the diff eq.

How do you now, how big is T?
I think five (approximate) periods will be better.

BTW, tspan other than [tstart, tend] give you more points,
but doesn't increase accuracy.

> then
>
> function phiprime = diffeq(t,phi)
> m =1;
> g=0;
> l=1;
> H = 0;
> phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];
>
> where i just let phi(1) = y, phi(2) = y'

> I've tried this code with g, l, H = 0 which becomes y'' = -y

And it works!

%l=0,h=0,g=0
T=10;
[t,phi] = ode45(@diffeq,0:T/100:T,[1,0]);
plot(t,phi(:,1))
root1 = fzero(@(x)interp1(t,phi(:,1),x),1*pi/2);
root2 = fzero(@(x)interp1(t,phi(:,1),x),3*pi/2);
[root1,root2]-[pi/2,3*pi/2]
%ans = -0.000118529835060111 -0.000454050586505161

Could you paste the code that generates wrong answer?

> What exactly does the reltol, abstol do also.

help odeset
;p

ode45 (and all matlab solvers) is 'adaptative'.
In every iteration compute estimated error, and,
if it is too big, decrease step size (Δt).

Absolute tolerance tell us, how big error
can be, and relative tolerance - how big compedred
to value.
error <= max(RelTol*|y| , AbsTol).

bartekltg <bartekltg@gmail.com> wrote in message <kqp8ic$864$1@node2.news.atman.pl>...
> W dniu 2013-06-30 10:32, William pisze:
> > bartekltg <bartekltg@gmail.com> wrote in message
> > <kqnu6a$k9n$1@node1.news.atman.pl>...
> >> W dniu 2013-06-28 22:38, William pisze:
> >> > Hi,
> >> >
> >> > I've been working on a program using ODE's and I've been using ODE45 to
> >> > get numerical solutions for differential equations. I've been getting
> >> > incorrect results e.g. when I used
> >> > y'' = -y, I get a period that is ~ 2.8 instead of pi. This also
> >> happened
> >> > with various other functions that I've been using.
> >> >
> >> > Is there a solution to this problem?
> >>
> >> I'm great wizard from quadruple-precision-in-matlab-land.
> >> My crystal ball told me you have done something wrong[1]. Sadly,
> >> this crystal ball is not HD ready, so I can't read you code:(
> >> Maybe you can post it? :)
> >>
> >> > Thanks!
> >>
> >> [1]:
> >> % solution of y'' = -y; y(1,:) is y, y(1,:) is velocity (y')
> >>
> >> [t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0,10],[0,1]);
> >>
> >> % of course ode45 solves first order equation, so we wrote
> >> % y'' = -y as system {y'=u, u'=-y}
> >>
> >> plot(t,y(:,1))
> >> grid
> >> %looks nice
> >>
> >> fzero(@(x)interp1(t,y(:,1),x),pi)
> >>
> >> % ans = 3.14134838301112
> >> % meh, but significantly better than 2.8
> >>
> >> options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
> >> [t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0,10],[0,1],options );
> >>
> >> root = fzero(@(x)interp1(t,y(:,1),x),pi)
> >> (root-pi)/pi
> >>
> >> % root = 3.14159265459964
> >> % ans = 3.21445091288877e-010
> >> % :-)
> >>
> >>
> >
> > So my code looks something like this
> >
> > [t,phi] = ode45(@diffeq,0:T/100:T,[1,0]);
> > %T is the period of the diff eq.
>
> How do you now, how big is T?
> I think five (approximate) periods will be better.
>
> BTW, tspan other than [tstart, tend] give you more points,
> but doesn't increase accuracy.
>
> > then
> >
> > function phiprime = diffeq(t,phi)
> > m =1;
> > g=0;
> > l=1;
> > H = 0;
> > phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];
> >
> > where i just let phi(1) = y, phi(2) = y'
>
> > I've tried this code with g, l, H = 0 which becomes y'' = -y
>
> And it works!
>
> %l=0,h=0,g=0
> T=10;
> [t,phi] = ode45(@diffeq,0:T/100:T,[1,0]);
> plot(t,phi(:,1))
> root1 = fzero(@(x)interp1(t,phi(:,1),x),1*pi/2);
> root2 = fzero(@(x)interp1(t,phi(:,1),x),3*pi/2);
> [root1,root2]-[pi/2,3*pi/2]
> %ans = -0.000118529835060111 -0.000454050586505161
>
> Could you paste the code that generates wrong answer?
>
> > What exactly does the reltol, abstol do also.
>
> help odeset
> ;p
>
> ode45 (and all matlab solvers) is 'adaptative'.
> In every iteration compute estimated error, and,
> if it is too big, decrease step size (?t).
>
> Absolute tolerance tell us, how big error
> can be, and relative tolerance - how big compedred
> to value.
> error <= max(RelTol*|y| , AbsTol).
>
>
>
>
>
>

So I tried your code out on my matlab and for some reason, i got [root1,root2]-[pi/2,3*pi/2] gave me -0.1647 -0.4933 as opposed to what you got.
So.. maybe this is something wrong with my matlab in general?

W dniu 2013-06-30 20:08, William pisze:


>
> So I tried your code out on my matlab and for some reason, i got
> [root1,root2]-[pi/2,3*pi/2] gave me -0.1647 -0.4933 as opposed to what
> you got.
> So.. maybe this is something wrong with my matlab in general?


Are you sure _everything_ was exactly the same?

"
function phiprime = diffeq(t,phi)
m =1;
g=0;
l=0; % !!!
H = 0;
phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];
"
in file "diffeq.m" ?



Try to clear workspace and work in new, empty folder.

Try also this:


T=10;
[t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0:T/100:T],[1,0]);
roots = [pi/2, 3*pi/2,5*pi/2];

nroots = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
errors = abs(roots - nroots)

options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
[t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0:T/100:T],[1,0],options );

nroots2 = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
errors2 = abs(roots - nroots2)

What answers do you get?

bartekltg <bartekltg@gmail.com> wrote in message <kqq1kj$25h$1@node2.news.atman.pl>...
> W dniu 2013-06-30 20:08, William pisze:
>
>
> >
> > So I tried your code out on my matlab and for some reason, i got
> > [root1,root2]-[pi/2,3*pi/2] gave me -0.1647 -0.4933 as opposed to what
> > you got.
> > So.. maybe this is something wrong with my matlab in general?
>
>
> Are you sure _everything_ was exactly the same?
>
> "
> function phiprime = diffeq(t,phi)
> m =1;
> g=0;
> l=0; % !!!
> H = 0;
> phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];
> "
> in file "diffeq.m" ?
>
>
>
> Try to clear workspace and work in new, empty folder.
>
> Try also this:
>
>
> T=10;
> [t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0:T/100:T],[1,0]);
> roots = [pi/2, 3*pi/2,5*pi/2];
>
> nroots = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
> errors = abs(roots - nroots)
>
> options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
> [t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0:T/100:T],[1,0],options );
>
> nroots2 = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
> errors2 = abs(roots - nroots2)
>
> What answers do you get?
>
>
>
>
>
>

My output looks like this

   1.0e-03 *

    0.1185 0.4541 0.8734

   1.0e-04 *

    0.1434 0.1362 0.0330

I restarted my matlab and it didn't change the output.

W dniu 2013-07-01 09:49, William pisze:
> bartekltg <bartekltg@gmail.com> wrote in message
> <kqq1kj$25h$1@node2.news.atman.pl>...
>> W dniu 2013-06-30 20:08, William pisze:
>>
>>
>> >
>> > So I tried your code out on my matlab and for some reason, i got
>> > [root1,root2]-[pi/2,3*pi/2] gave me -0.1647 -0.4933 as opposed to
>> what
>> > you got.
>> > So.. maybe this is something wrong with my matlab in general?
>>
>>
>> Are you sure _everything_ was exactly the same?
>>
>> "
>> function phiprime = diffeq(t,phi)
>> m =1;
>> g=0;
>> l=0; % !!!
>> H = 0;
>> phiprime=[phi(2); -3*H.*phi(2)-m^2*phi(1) - g*phi(1).^2 - l*phi(1).^3];
>> "
>> in file "diffeq.m" ?
>>
>>
>>
>> Try to clear workspace and work in new, empty folder.
>>
>> Try also this:
>>
>>
>> T=10;
>> [t,y] = ode45 ( @(t,Y) [Y(2);-Y(1)],[0:T/100:T],[1,0]);
>> roots = [pi/2, 3*pi/2,5*pi/2];
>>
>> nroots = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
>> errors = abs(roots - nroots)
>>
>> options = odeset ('RelTol',1e-10,'AbsTol',1e-10); % quite useful:)
>> [t,y] = ode45 ( @(t,x) [x(2);-x(1)],[0:T/100:T],[1,0],options );
>>
>> nroots2 = arrayfun(@(r)fzero(@(x)interp1(t,y(:,1),x),r),roots );
>> errors2 = abs(roots - nroots2)
>>
>> What answers do you get?
>>
>>
>>
>>
>>
>>
>
> My output looks like this
>
> 1.0e-03 *
>
> 0.1185 0.4541 0.8734
>
> 1.0e-04 *
>
> 0.1434 0.1362 0.0330
>
> I restarted my matlab and it didn't change the output.

This is OK.
On my computer:

errors =

   1.0e-003 *

     0.1185 0.4541 0.8734

errors2 =

   1.0e-004 *

     0.1434 0.1362 0.0330

(Try 'format short g' sometimes it looks better)


So your matlab work fine. I still think the problem may by in
your diffeq.m file.