problem to solve an EDO

hi, i'm tring to solve an EDO but i've got this error message : Index exceeds matrix dimensions Error in ==> f_cine at 3 f=abs(y(:,1)-y(:,2)-y(:,3))-sg0 and i don't understand.

Here this is my function describing the system: function dy=f_cine(y,t)

global C1 C2 D1 D2 Q b K n E sig0 ep

f=abs(y(1)-y(2)-y(3))-sg0

if f<0

    dy(4)=0

else

    dy(4)=(f/K)^n

end

dy(2)=C(1)*dy(5)-D1*y(2)*dy(4)

dy(3)=C(2)*dy(5)-D2*y(3)*dy(4)

dy(5)=sign(y(1)-y(2)-y(3))*dy(4) dy(1)=E*(ep-dy(5)) dy=dy' end

thanks for your help

Answer 1

[Colin]

Did you check the dimensions of sg0?

Answer 2

[tony]

yes, it's a reel, defined in my main function.

I put at the beginning y=zeros(5,1) and the error has changed,

??? Attempted to access dy(5); index out of bounds because numel(dy)=4.

Answer 3

[tony]

I'also changed the order like that :

function dy=f_cine(y,t)

global C1 C2 D1 D2 Q b K n E sg0 ep

y=zeros(5,1);

f=abs(y(1)-y(2)-y(3))-sg0

if f<0

    dy(4)=0

else

    dy(4)=(f/K)^n

end

dy(5)=sign(y(1)-y(2)-y(3))*dy(4)

dy(2)=C1*dy(5)-D1*y(2)*dy(4)

dy(3)=C2*dy(5)-D2*y(3)*dy(4)

dy(1)=E*(ep-dy(5))

dy=dy'

end

Now there is no error message, but the solution is fault because all equal zero

Answer 4

[Colin]

next to y = zeros(5,1); you should also initiate dy = zeros(5,1); otherwise Matlab does not know of wich size dy is

How does the RHS of your ODE look like?

Answer 5

[tony]

global C1 C2 D1 D2 Q b K n E sg0 ep

% parametres

E=140000;

sg0=200;

K=800;

n=6;

C1=300000;

C2=25000;

D1=2000;

D2=200;

ep=0,01;

%temps finale

tf=400;

%condition initiale

[t,y]=ode45('f_cine', [0 tf], [0 0 0 0 0]);

y1 = y(:,1)

y2 = y(:,2)

y3 = y(:,3)

y4 = y(:,4)

y5 = y(:,5)

No errors, but the result is obviously wrong

Answer 6

[Colin]

This is, what you coded:

y(5)=sign(y(1)-y(2)-y(3))*dy(4)

dy(2)=C1*dy(5)-D1*y(2)*dy(4)

dy(3)=C2*dy(5)-D2*y(3)*dy(4)

dy(1)=E*(ep-dy(5))

dy=dy'

I guess, what you want to have

dy = zeros(5,1);

dy(1)=E*(ep-y(5))

dy(2)=C1*y(5)-D1*y(2)*y(4)

dy(3)=C2*y(5)-D2*y(3)*y(4)

dy(5)=sign(y(1)-y(2)-y(3))*y(4)

Is that correct? If not, please post the mathematical form of your differential equation.

Answer 7

[tony]

No, my differential equation is

f=abs(sigma-X1-X2)-sig0

if f<0

   dp/dt=0

else

   dp/dt=(f/K)^n

and

dXi/dt=Ci dep/dt-Di*Xi*dp/dt i=1,2

dep/dt=sign(sima-X1-X2)*dp/dt

dsigma/dt=E(de/dt-dep/dt)

de/dt is given and at the end , i would like extract sigma in order to have sigma=f(e)

Answer 8

[Colin]

If I understand this correct, you want to solve an implicit system of differential equations, as for example de/dt is part of your RHS. Maybe you can elliminate that terms by substitution and solve this with ode45.

To solve implict systems of differential equations, have a look to the solver ode15i.

If you have also given de/dt at the final time, you need to solve a boundary value problem, not just an initial value problem. bvp4c for example is a solver, which can solve boundary value problems.

If your method diverges, you should have a closer look to that parts of your differential equations containing sign(). Some methods are using Taylor expansions to integrate the system. If your RHS is not continous differentiable to some order, this might be a problem.

Answer 9

[tony]

yes , i'm using here EDO45 thanks for your help,i'll give a shot a it again