"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
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Reply via email to