[Scilab-users] 'ode' returns constant values for the variables

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

[Scilab-users] 'ode' returns constant values for the variables

Hello,
I try to solve the sisytema with the help of "ode".
But I can not find where the error is; in the solution it returns constant
values for the variables.
I think the error is that I have not correctly declared the dependence of
the variable variables "t".
The expected solution should allow the graph to be reproduced:

<http://mailinglists.scilab.org/file/t497622/Solution_expected_fromODE.jpg>

clear;
funcprot(0);

//Draghilev's Method

timer();
function Sys=F(x)
    Sys=[];
    x0=x(1);x1=x(2);
    Sys(1)=x1.^2+x0.^2 -8;
    Sys(2)=sin(x0.*x1).*sin(exp(x0.*x1));
endfunction
function Sys=F0(x)
    Sys=[];
    x00=x(1);x01=x(2);
    Sys(1)=x01.^2+x00.^2-8;
    Sys(2)=sin(x00.*x01).*sin(exp(x00.*x01));
endfunction
function Sys=Mp(v)
    [Sys]=F(v(1:2))-v(3)*F0(ics(1:2));
endfunction
//deff('Sys=Mp(v)','[Sys]=F(v(1:2))-v(3).*F0(ics([1:2]))')
function Sys=NF(v)
    Sys=[];
    [Sys]=numderivative(Mp,v);
endfunction
//
//Z=NF
//
function  D=ODEs(ZZ)
A=ZZ(:,1:2);b=ZZ(:,$);
A1=A;A1(:,1)=b;
A2=A;A2(:,2)=b;
D=-det(A);
d(1)=det(A1) ; d(2)=det(A2); d(3)=D;
[D]=-d;
endfunction

ics=[2;1]//[2.784;-0.5];[2;-0.5];
yint1=fsolve(ics,F)
ics=[yint1];
disp(F(yint1),"x->f(x)=0",ics,"x f(x)=0")

ics=[ics;1];
disp(ics)
function [dxdt]=odes(t,x)
[dxdt]=ODEs(NF(x));
endfunction
N=200;
smin=0.0;
smax=15;
h=0.000001;

step=smax/N;
    t=smin:step:smax;
    t0=0;
    atol=h;
    LL= ode("fix",ics,t0,t,atol,odes)
disp(timer(),"tiempo de ejecucion C");

clf();
subplot(1,2,1);
a=get("current_axes")//get the handle of the newly created axes
a.axes_visible="on"; // makes the axes visible
a.font_size=3; //set the tics label font size
a.x_location="middle"; //set the x axis position
plot(t,LL(1,1:$),"-cya",t,LL(2,1:$),"-g",t,LL(3,1:$),"-r")

Gracias



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Rafael Guerra Rafael Guerra
Reply | Threaded
Open this post in threaded view
|

Re: 'ode' returns constant values for the variables

Hi Hermes,

The right side of your differential equation does not seem to depend on 't' as it should.

This does not seem to be a Scilab problem but a math problem and you should in that case point, to those interested, towards a reference describing the differential equations to be solved.

Regards,
Rafael

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Hermes Hermes
Reply | Threaded
Open this post in threaded view
|

Re: 'ode' returns constant values for the variables

In reply to this post by Hermes
"Stiffness occurs when some components of the solution decay much more
rapidly than others."
*Stretched energy function Sys(2)=sin(x0.*x1).*sin(exp(x0.*x1)*

function Sys=F(x)
    Sys=[];
    x0=x(1);x1=x(2);
    Sys(1)=x1*x1+x0*x0-8;
    Sys(2)=sin(x0.*x1).*sin(exp(x0.*x1));
endfunction

function Sys=Mp(v,xx)
    [Sys]=F(v(1:2))-v(3)*F(xx);
endfunction

function Sys=NF(v)
    Sys=[];
    [Sys]=numderivative(list(Mp,ics),v);
endfunction

function  D=ODEs(ZZ)
A=ZZ(:,1:2);b=ZZ(:,$);
A1=A;A1(:,1)=b;
A2=A;A2(:,2)=b;
D=-det(A);//main determinant
d(1)=det(A1) ; d(2)=det(A2); d(3)=D;
[D]=-d;
endfunction
   
ics=[2.784;-0.5];
v0=[ics;1];

function [dxdt]=odes(t,x)
[dxdt]=ODEs(NF(x));
endfunction
N=200;
smin=0.0;
smax=15;
h=0.000001;

step=smax/N;
    t=smin:step:smax;
    t0=0;
    atol=h;
    LL= ode("stiff",v0,t0,t,atol,odes)
disp(timer(),"execution time:");

clf();
subplot(1,2,1);
a=get("current_axes")//get the handle of the newly created axes
a.axes_visible="on"; // makes the axes visible
a.font_size=3; //set the tics label font size
a.x_location="middle"; //set the x axis position
plot(t,LL(1,1:$),"-cya",t,LL(2,1:$),"-g",t,LL(3,1:$),"-r")

subplot(1,2,2);

s= LL(3,:);
ym= s(1:$-1).*s(2:$);
z= find(ym <= 0);
t0 = t(z) - s(z).*(t(z+1)-t(z))./(s(z+1)-s(z));
y01 = interpln([t;LL(1,:)],t0);
y02 = interpln([t;LL(2,:)],t0);

param3d(LL(1,:), LL(2,:), 0 * LL(1,:), alpha=35, theta=45); // circle in the
xy plane
e = gce();
e.foreground =17;
param3d(LL(1,:), LL(2,:), LL(3,:), alpha=35, theta=45) // parametric
function in 3d-space
param3d(y01, y02, 0*y01, alpha=35, theta=45); // zero-crossings
e = gce();
e.line_mode = "off";
e.mark_mode = "on"
e.mark_style = 3;
e.mark_foreground =5 ;

<http://mailinglists.scilab.org/file/t497622/find_roots_syst4Scilab.jpg>

Gracias



--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users