ode45 for coupled system (2nd order ODE)

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

ode45 for coupled system (2nd order ODE)

Joza
Hi guys. I haven't used ode45 before, so I was hoping to get some help here regarding its use.

I want to solve a coupled system, namely dy2/dt = -sin(y1) and dy1/dt = y2.

ode45 will solve this for me for y1 as a function of t, that's nice. But say I want to obtain (y1, dy1/dt) ... is there a nice way of doing this, or will ode45 only put out (y1, t) ?

Reply | Threaded
Open this post in threaded view
|

RE: ode45 for coupled system (2nd order ODE)

Allen.Windhorn-2
Joza,

> -----Original Message-----
> From: [hidden email] [mailto:help-octave-
>
> Hi guys. I haven't used ode45 before, so I was hoping to get
> some help here regarding its use.
>
> I want to solve a coupled system, namely dy2/dt = -sin(y1) and
> dy1/dt = y2.
>
> ode45 will solve this for me for y1 as a function of t, that's
> nice.  But say I want to obtain (y1, dy1/dt) ... is there a nice
> way of doing this, or will ode45 only put out (y1, t) ?

You have y2, which is dy1/dt.  If you need all the derivatives,
you already have a function to use in ode45 like

fdydt = @(t,y) [-sin(y(2); y(1)];

so your derivatives can be calculated as fdydt(t,y).

I tried to use ode45 for Matlab compatibility, but had some trouble
with convergence, so I have been using lsode instead.  The system
you are showing here is not stiff, so ode45 should work OK, but you
may want to try other methods as you get more complex problems.

Regards,
Allen
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: ode45 for coupled system (2nd order ODE)

Jordi Gutiérrez Hermoso-2
In reply to this post by Joza
On 9 November 2012 09:29, Joza <[hidden email]> wrote:
> Hi guys. I haven't used ode45 before, so I was hoping to get some help here
> regarding its use.
>
> I want to solve a coupled system, namely dy2/dt = -sin(y1) and dy1/dt = y2.
>
> ode45 will solve this for me for y1 as a function of t, that's nice.

Btw, I know Matlab users are taught to use ode45 as the default ODE
solver, but in Octave there is also lsode. Its interface is similar to
ode45, but it can usually handle a wider class of problems.

- Jordi G. H.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

RE: ode45 for coupled system (2nd order ODE)

Joza
In reply to this post by Allen.Windhorn-2
Why would I use -sin(dy1/dt) ?The sin function only applies to y1, not its derivative?
Reply | Threaded
Open this post in threaded view
|

RE: ode45 for coupled system (2nd order ODE)

Joza
Here is my code anyway:

function du = f1(t, y)
        dy = zeros(2,1);
        dy(1) = y(2);
        dy(2) = -sin(y(1));
end

...

[t, y1] = ode45(@f1, [0, 100], [1, 0]);

This solves the system for y1. But I would like to also have dy1/dt .... how can I obtain that with ode45?
Reply | Threaded
Open this post in threaded view
|

Re: ode45 for coupled system (2nd order ODE)

Juan Pablo Carbajal-2
On Fri, Nov 9, 2012 at 6:41 PM, Joza <[hidden email]> wrote:

> Here is my code anyway:
>
> function du = f1(t, y)
>         dy = zeros(2,1);
>         dy(1) = y(2);
>         dy(2) = -sin(y(1));
> end
>
> ...
>
> [t, y1] = ode45(@f1, [0, 100], [1, 0]);
>
> This solves the system for y1. But I would like to also have dy1/dt .... how
> can I obtain that with ode45?
>
>
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/ode45-for-coupled-system-2nd-order-ODE-tp4646324p4646344.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

dy1/dt according to your equation is nothing but y1(,:)

If you need other derivatives, you can evaluate your f1 function.

Reocmendation:
You can write your ode/lsode functions in the following way
function du = f1(y, t)
         dy = zeros(2,size(y,2));
         dy(1,:) = y(2,:);
         dy(2,:) = -sin(y(1,:));
 end

then once you have solved your equation
[t, sol] = lsode(@f1, [0, 100], [1, 0]);

you can run
dsol = f1(sol', t)

to get all derivatives.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

RE: ode45 for coupled system (2nd order ODE)

Allen.Windhorn-2
In reply to this post by Joza
Joza,

> -----Original Message-----
> From: [hidden email] [mailto:help-octave-
>
> Here is my code anyway:
>
> function du = f1(t, y)
           ^^ this should be dy

> dy = zeros(2,1);
> dy(1) = y(2);
> dy(2) = -sin(y(1));
> end
>
> ...
>
> [t, y1] = ode45(@f1, [0, 100], [1, 0]);
>
> This solves the system for y1. But I would like to also have
> dy1/dt .... how can I obtain that with ode45?

Well, dy1/dt IS y2, from the code above.  And you are really
computing y1 AND y2 (look at the value of "y1", it is actually

octave:11> y1(1:20,:)
ans =

   1.00000   0.00000
   0.97534  -0.20289
   0.90040  -0.40248
   0.77762  -0.58699
   0.61078  -0.74686
   0.40903  -0.86857
   0.21155  -0.93531
   0.02071  -0.95863
  -0.16367  -0.94481
  -0.34044  -0.89700
  -0.50747  -0.81692
  -0.66159  -0.70529
  -0.79795  -0.56192
  -0.90883  -0.38565
  -0.98170  -0.17496
  -0.99920   0.03680
  -0.96441   0.24332
  -0.87976   0.44054
  -0.74789   0.62101
  -0.57295   0.77460 ...

where the left column is y1 and the right is y2, which is dy1/dt
for the values of time in t).

If you also needed dy2/dt, after using ode45, you could just do

dy = f1(t, y1)

if f1 would handle a matrix input.  I know there's a way to do
that but I can't think what it is at the moment.

Regards,
Allen
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave