ode45 Matlab

Rapha · #1 $Old$ November 26th, 2008, 09:41 PM

Hello everyone!

I need to solve a IVP, for example

$\dot{x} = x + 2\dot{y} + y$

$\dot{y} = y - 2\dot{x} + 5x$

x(0) = 0.995
y(0) = 0

(I actually don't know if there is a solution...)

But how to solve it by using ode45? $\dot{x}$ depends on $\dot{y}$ and y ...

By the way,
I know how to solve

$\dot{v}(t) = 10 - 5/2*v^2(t), t\in [0, 10]$

v0 = v(0) = 0

for example:

[ t , v ] = ode45 (@( t , v ) 10-5v^2/2 , [0, 10], v0 ) ;

----
How to solve the IVP with 2 equations?

Best regards,
Rapha

CaptainBlack · #2 $Old$ November 26th, 2008, 10:07 PM

Quote:

Originally Posted by Rapha $View Post$

Hello everyone!

I need to solve a IVP, for example

$\dot{x} = x + 2\dot{y} + y$

$\dot{y} = y - 2\dot{x} + 5x$

x(0) = 0.995
y(0) = 0

(I actually don't know if there is a solution...)

But how to solve it by using ode45? $\dot{x}$ depends on $\dot{y}$ and y ...

By the way,
I know how to solve

$\dot{v}(t) = 10 - 5/2*v^2(t), t\in [0, 10]$

v0 = v(0) = 0

for example:

[ t , v ] = ode45 (@( t , v ) 10-5v^2/2 , [0, 10], v0 ) ;

----
How to solve the IVP with 2 equations?

Best regards,
Rapha

Eliminate the $\dot{y}$ term in the first equation by substituting from the second, and $\dot{x}$ for the second by substituting from the first:

$\dot{x}=\frac{11x+2y}{5}$

$\dot{y}=\frac{-2x-y}{5}$

So now writing this as a first order vector ODE:

$\dot{\bold{x}}=\left[\begin{array}{cc}11/5 & 2/5\\-2/5 & -1/5\end{array}\right] {\bold{x}}$

where:

${\bold{x}}=\left[\begin{array}{c}x\\y\end{array}\right]$

So now can you get ode45 to work?

CB

Rapha · #3 $Old$ November 27th, 2008, 07:42 AM

Hello.

Quote:

So now can you get ode45 to work?

that was very helpful, thank you very much (and yes, i could get ode45 to work).

But actually i kinda screwed it(sorry, i did not exactly realize what my problem was and asked for something different), because my problem is :

$\ddot{x} = x + 2 \dot{y} - 0.9 \frac{x+0.1}{D_1} - 0.1 \frac{x-0.9}{D_2}$

$\ddot{y} = y - 2 \dot{x} - 0.9 \frac{y}{D_1} - 0.1 \frac{y}{D_2}$

with $D_1 = \sqrt{(((x+0.1)^2+y^2)^3)}$

$D_2 = \sqrt{(((x-0.9)^2+y^2)^3)}$

i tried to find the IVP 1. order:

$x_1 = x$

$\dot{x_1} = \dot{x} = x_2$

$y_1 = y$

$y_2 = \dot{y_1} = \dot{y}$

=>
$\dot{x_2} = x_1 + 2 y_2 - 0.9 \frac{x+0.1}{D_1} - 0.1 \frac{x-0.9}{D_2}$

$\dot{y_2} = y_1 - 2 x_2 - 0.9 \frac{y_1}{D_1} - 0.1 \frac{y_1}{D_2}$

I did not realize that there are 4 functions: $y_2(t), y_1(t), x_1(t), x_2(t)$

I still want to solve this by using ode45.

Regards
Rapha

CaptainBlack · #4 $Old$ November 27th, 2008, 09:58 PM

Quote:

Originally Posted by Rapha $View Post$

Hello.

that was very helpful, thank you very much (and yes, i could get ode45 to work).

But actually i kinda screwed it(sorry, i did not exactly realize what my problem was and asked for something different), because my problem is :

$\ddot{x} = x + 2 \dot{y} - 0.9 \frac{x+0.1}{D_1} - 0.1 \frac{x-0.9}{D_2}$

$\ddot{y} = y - 2 \dot{x} - 0.9 \frac{y}{D_1} - 0.1 \frac{y}{D_2}$

with $D_1 = \sqrt{(((x+0.1)^2+y^2)^3)}$

$D_2 = \sqrt{(((x-0.9)^2+y^2)^3)}$

i tried to find the IVP 1. order:

$x_1 = x$

$\dot{x_1} = \dot{x} = x_2$

$y_1 = y$

$y_2 = \dot{y_1} = \dot{y}$

=>
$\dot{x_2} = x_1 + 2 y_2 - 0.9 \frac{x+0.1}{D_1} - 0.1 \frac{x-0.9}{D_2}$

$\dot{y_2} = y_1 - 2 x_2 - 0.9 \frac{y_1}{D_1} - 0.1 \frac{y_1}{D_2}$

I did not realize that there are 4 functions: $y_2(t), y_1(t), x_1(t), x_2(t)$

I still want to solve this by using ode45.

Regards
Rapha

Now we have as state vector:

${\bold{x}}=\left[ \begin{array}{c}x\\y\\ \dot{x}\\ \dot{y} \end{array} \right]$

and derivative:

$\dot{\bold{x}}=\left[ \begin{array}{c} \bold{x}_3\\ \bold{x}_4\\ \bold{x}_1-2\bold{x}_4-0.9(\bold{x}_1+0.1)/D_1-0.1(\bold{x}_1-0.9)/D_2 \\ \bold{x}_2-2\bold{x}_3-0.9\bold{x}_2/D_1-0.1\bold{x}_2/D_2 \end{array} \right]$

where:

$D_1 = \sqrt{(((\bold{x}_1+0.1)^2+\bold{x}_2^2)^3)}$

$D_2 = \sqrt{(((\bold{x}_1-0.9)^2+\bold{x}_2^2)^3)}$

You will need to write a Matlab function to evaluate the derivative and pass that to ode45 (the derivative is now too complicated to be easily passed as an anonymous function).

(You cound simplify the derivative by collecting terms)

CB

Rapha · #5 $Old$ November 29th, 2008, 04:35 AM

Hello CaptainBlack!

Thank you so much for your helpful comment!

Quote:

Originally Posted by CaptainBlack $View Post$

You will need to write a Matlab function to evaluate the derivative and pass that to ode45 (the derivative is now too complicated to be easily passed as an anonymous function).

(You cound simplify the derivative by collecting terms)

I tried it but it still doesn't work.

Code:

function dx = ivp(t, x)

dx =  \dot{x} * x;

clearly represented, I wrote $\dot{x}$ for
$\dot{\bold{x}}=\left[ \begin{array}{c}\bold{x}_3\\\bold{x}_4\\ \bold{x}_1-2\bold{x}_4-0.9(\bold{x}_1+0.1)/D_1-0.1(\bold{x}_1-0.9)/D_2 \\ \bold{x}_2-2\bold{x}_3-0.9\bold{x}_2/D_1-0.1\bold{x}_2/D_2 \end{array} \right]$

then:

[t, x] = ode45(@ivp, [0, 10], [??,??,??,??] )

where $t \in [0, 10]$ and ?? = [x_3(0); x_4(0);, ...(0); ...(0)]

this is stupid, but I don't know the initial values :-(

CaptainBlack · #6 $Old$ November 29th, 2008, 06:59 AM

Quote:

Originally Posted by Rapha $View Post$

Hello CaptainBlack!

Thank you so much for your helpful comment!

I tried it but it still doesn't work.

Code:

function dx = ivp(t, x)
 
dx =  \dot{x} * x;

clearly represented, I wrote $\dot{x}$ for
$\dot{\bold{x}}=\left[ \begin{array}{c} \bold{x}_3\\ \bold{x}_4\\ \bold{x}_1-2\bold{x}_4-0.9(\bold{x}_1+0.1)/D_1-0.1(\bold{x}_1-0.9)/D_2 \\ \bold{x}_2-2\bold{x}_3-0.9\bold{x}_2/D_1-0.1\bold{x}_2/D_2 \end{array} \right]$

then:

[t, x] = ode45(@ivp, [0, 10], [??,??,??,??] )

where $t \in [0, 10]$ and ?? = [x_3(0); x_4(0);, ...(0); ...(0)]

this is stupid, but I don't know the initial values :-(

Somewere in the staement of the problem there should be initial values.

But you can try any initial values to test the code [0;0;1;1] should do.

the following seems to work on FreeMat:

Code:

function dx=deriv(t,x)
   d1=sqrt(((x(1)+0.1)^2+x(2)^2)^3);
   d2=sqrt(((x(1)-0.9)^2+x(2)^2)^3);
   dx=zeros(4,1);
   dx(1)=x(3);
   dx(2)=x(4);
   dx(3)=x(1)-2*x(4)-0.9*(x(1)+0.1)/d1-0.1*(x(1)-0.9)/d2;
   dx(4)=x(2)-2*x(3)-0.9*x(2)/d1-0.1*x(2)/d2;

with calling statement:

Code:

[t,x]=ode45(@deriv ,[0,10],[0;0;1;1]);

(note the absurdty that we have state and derivative as colum vectors, but the output is an array of row vectors, one for each time point, but it does not work if we have row vectors for the state and derivative!!)

CB

Rapha · #7 $Old$ November 29th, 2008, 09:21 PM

Quote:

Originally Posted by CaptainBlack $View Post$

Somewere in the staement of the problem there should be initial values.

Oooops
I forgot to post them :-(

Quote:

Originally Posted by CaptainBlack $View Post$

the following seems to work on FreeMat:

Code:

function dx=deriv(t,x)
   d1=sqrt(((x(1)+0.1)^2+x(2)^2)^3);
   d2=sqrt(((x(1)-0.9)^2+x(2)^2)^3);
   dx=zeros(4,1);
   dx(1)=x(3);
   dx(2)=x(4);
   dx(3)=x(1)-2*x(4)-0.9*(x(1)+0.1)/d1-0.1*(x(1)-0.9)/d2;
   dx(4)=x(2)-2*x(3)-0.9*x(2)/d1-0.1*x(2)/d2;

with calling statement:

Alright, thank you. You are great!

Rapha