Equations Matlab vs. Scilab

Submitted by vesy on Tue, 06/07/2011 - 19:48

Hello!


I'm new to this forum. I would like to ask you for your help on my following problem Smile I have equation in matlab *.m file, that i would like to run in scilab.
 
MATLAB

function elliptical=equ_mot(t,Y)

if t>=t_start 
  acc=acc_ratio*9.81*sin(t);
else
  acc=0;
end

  elliptical(1,1)=Y(2);
  elliptical(2,1)=(-a*cos(Y(1))*acc+0.5*Y(2)^2*(a^2-b^2)*sin(2*Y(1))-9.81*b*sin(Y(1)))/(a^2*(cos(Y(1)))^2+b^2*sin((Y(1)))^2);

//this is called into primary code with following syntax:
options=odeset('RelTol',1e-6,'AbsTol',1e-6,[1e-6 1e-6])
Y=[0;0];
[t,dif]=ode45('equ_mot',sim_t,Y,options)


I have tried to rewrite equation into scilab as follows:

Y0=[0;0];
Time=20;
sim_t=[0:0.1:Time];

function acc=g(sim_t)
  acc_ratio=0.01
  acc=acc_ratio*9.81*sin(sim_t)
endfunction

function elliptical=f(sim_t,Y)
  elliptical=(-a*cos(Y(1))*g(sim_t)+0.5*Y(2)^2*(a^2-b^2)*sin(2*Y(1))-9.81*b*sin(Y(1)))/(a^2*(cos(Y(1)))^2+b^2*sin((Y(1)))^2);
endfunction

y=ode('rkf',Y0,0,sim_t,1.d-6,1.d-6,f);

scilab returns me following error Sad

>y=ode('rkf',Y0,0,sim_t,1.d-6,1.d-6,f);
                                      !--error 98
 variable returned by scilab argument function is incorrect
at line      17 of exec file called by :
exec("D:/Vesna/DIPLOMSKA (SPH)/SPH/PRIPRAVA ZA POSLAT/NIHALO_en_gib_SCILAB.sce");
in  execstr instruction    called by :
scipad('D:\Vesna\DIPLOM~1\SPH\PRIPRA~1\NIC15D~1.SCE');   

 
I don't know what am I doing wrong Sad I think, that there is something missing, pa I don't know what Sad

thx all for your suggestions
Taxonomy upgrade extras: 
General

sudarson

Wed, 08/17/2011 - 05:08

If you carefully see your code you can get that function " f " is not properly defined. It should be two dimentional as your problem is related with 2nd order differential equation(I guess).

Try this code:

Y0=[0;0];
Time=20;
sim_t=[0:0.1:Time];

function acc=g(sim_t)
acc_ratio=0.01
acc=acc_ratio*9.81*sin(sim_t)
endfunction

function elliptical=f(sim_t,Y)
a=5;b=0.5;
elliptical=[Y(2);(-a*cos(Y(1))*g(sim_t)+0.5*Y(2)^2*(a^2-b^2)*sin(2*Y(1))-9.81*b*sin(Y(1)))/(a^2*(cos(Y(1)))^2+b^2*sin((Y(1)))^2)];
endfunction

y=ode(Y0,0,sim_t,1.d-6,1.d-6,f);
plot2d(sim_t',y')

I did minimal change.

Again you have not given any value of "a" and "b". I have taken some guess value.

Hope this will help you to understand.