Follow along with the video below to see how to install our site as a web app on your home screen.
Note: This feature may not be available in some browsers.
>> lsode_options()
Options for LSODE include:
keyword value
------- -----
absolute tolerance 1.49012e-08
relative tolerance 1.49012e-08
integration method stiff
initial step size -1
maximum order -1
maximum step size -1
minimum step size 0
step limit 100000
printf("System of Stiff Ordinary Differential Equations\n");
printf("by hasnaoui\n");
function YDOT = F(Y, T)
EXI=0.10;
C0=EXI-5.0/2.0;
C1=EXI/2.0-1.0/4.0;
C2=-0.50;
C3=-0.50;
C4=EXI-3.0/2.0;
GAM=5.0/3.00;
BETA=EXI/2.0+3.0/4.0;
GU=5.7D-3;
EP=0.10;
ALM=1.0;
XM=1.00;
PM=1.0;
PP=PM*EP+3*XM/ALM**2;
SM=ALM*PP*GU**0.5;
RM=PP/EP;
DEL0=(1-EP**2*(SM**2/2.0+2*(2-BETA)+GU*RM))**0.5;
QQ=ALM*DEL0*(PP-PM*EP)/2.0/GU**0.5;
TAU=PP/PM/EP-1.0;
%VARIABLE DI'INTEGRATION
YDOT(1) = 1.0;
%RESISTIVITE MAGNETIQUE
YDOT(10) = -2.0*Y(1)*exp(-Y(1)**2);
%
%VITESSE TOROIDALE
YDOT(2)=(Y(4)*Y(3)*Y(2)+ ...
TAU/(1.0+TAU)*(Y(9)*Y(6)/Y(10)-C1/BETA*Y(7)*Y(8)) ...
-1.0/(1.0+TAU)*Y(3)*Y(10)) ...
/(2.0*Y(3)*(Y(1)*Y(4)+Y(5)));
%DENSITE VOLUMIQUE
YDOT(3)=((EXI-1.0)*Y(3)*Y(4)*EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5)) ...
-Y(1)*Y(3)* ...
(C2*SM**2*EP**2*Y(3)*Y(4)**2-Y(3)*DEL0**2*Y(2)**2 ...
+Y(3)/(1.0+EP**2*Y(1)**2)**1.50+EP**2*C0*Y(11) ...
+GU*QQ**2*EP**2*Y(7)*(C1*Y(7)-Y(1)*Y(6)/Y(10)) ...
+GU/BETA*(Y(9)-Y(1)/BETA*Y(8))* ...
(-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) ...
-Y(3)* ...
(C3*SM**2*EP**2*Y(3)*Y(4)*Y(5)-Y(1)*Y(3)/ ...
(1.0+EP**2*Y(1)**2)**1.50 ...
-GU*QQ**2*Y(7)*Y(6)/Y(10)-GU/BETA**2/EP**2*Y(8)* ...
(-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) ...
-1.0/BETA*Y(3)*Y(9)**(-1.0-1.0/BETA)*Y(8)* ...
Y(3)*(1.0+Y(1)**2*EP**2)) ...
/(Y(3)*(SM**2*EP**2*(Y(1)*Y(4)+Y(5))**2-Y(9)**(-1.0/BETA)* ...
(1.0+EP**2*Y(1)**2)));
%PRESSION
YDOT(11)=-1.0/BETA*Y(3)*Y(9)**(-1.0-1.0/BETA)*Y(8)+ ...
Y(9)**(-1.0/BETA)*YDOT(3);
%EQUILIBRE RADIALE
YDOT(4)=(C2*SM**2*EP**2*Y(3)*Y(4)**2 ...
-Y(3)*DEL0**2*Y(2)**2 ...
+Y(3)/(1.0+(EP*Y(1))**2)**1.50 ...
+EP**2*C0*Y(11) ...
+GU*QQ**2*EP**2*Y(7)*(C1*Y(7)-Y(1)*Y(6)/Y(10)) ...
+GU/BETA*(Y(9)-Y(1)*Y(8)/BETA)* ...
(-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5))) ...
-EP**2*Y(1)*YDOT(11)) ...
/(EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5)));
%EQUILIBRE VERTICALE
YDOT(5)=(C3*EP**2*SM**2*Y(3)*Y(4)*Y(5)- ...
Y(1)*Y(3)/(1.0+(EP*Y(1))**2)**1.50 ...
-GU*QQ**2*Y(7)*Y(6)/Y(10)- ...
GU*Y(8)/BETA**2/EP**2* ...
(-RM*EP**2/Y(10))*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5))) ...
-YDOT(11)) ...
/(EP**2*SM**2*Y(3)*(Y(1)*Y(4)+Y(5)));
%INDUCTION MAGNETIQUE
YDOT(6)= (EP**2*Y(1)*((C1-1.0)*Y(6)+C1*Y(7)*YDOT(10)) ...
-EP**2*(C1*Y(7)*Y(10)-Y(1)*Y(6))*(C1-1.50) ...
-XM*RM*DEL0/QQ/SM*(1.50/BETA*Y(8)*Y(2)+Y(9)*YDOT(2)) ...
+XM*RM*EP**2*BETA*Y(7)*Y(4) ...
+XM*RM*EP**2*(Y(5)+Y(1)*Y(4))* ...
(Y(6)/Y(10)-Y(7)*YDOT(3)/Y(3))) ...
/(1.0+EP**2*Y(1)**2);
YDOT(7)=Y(6)/Y(10);
%FLUX MAGNETIQUE
YDOT(8)=(EP**2*(BETA*(2.0-BETA)*Y(9)+(2*BETA-3.0)*Y(1)*Y(8)) ...
-RM*EP**2/Y(10)*(BETA*Y(4)*Y(9)-Y(8)*(Y(1)*Y(4)+Y(5)))) ...
/(1.0+(EP*Y(1))**2);
YDOT(9)=Y(8);
endfunction
Y0(1) = 0.0001;
Y0(2) = 1.0;
Y0(3) = 1.0;
Y0(4) = 1.0;
Y0(5) = 0.0;
Y0(6) = -1.;
Y0(7) = 0.0;
Y0(8) = 0.0;
Y0(9) = 1.0;
Y0(10) = 1.0;
Y0(11) = 1.0;
T = linspace(0, 3, 300);
Y = lsode ("F", Y0, T);
for i = 1: columns(Y)
figure(i);
plot(T, Y(:,i));
tit=title ("hasnaoui", "fontsize", 16);
l=xlabel("T", "fontsize", 16);
ylabel(sprintf("y(%d)", i), "fontsize", 16, "rotation", 0);
endfor