Re:Re: ODE system solving

[address@hidden]

thanks for fast response!

i try to solve a polymer melt flow (fiber spinning system) with octave 2.1.42

in the file you can see the original system equation, written in the form 
y'=f(z,y)

the other functions (hc,Re,Rz, ecc )are described in the m-file reported here:

function SysO=f(y,z)
g=9810;                                 % mm2/s
rho=9e-07;                              % kg/mm3
Cp=1500;                                % J/Kg k
Ka=3e-5;                                % W/mm k
R1=0.53;                                % mm
R2=0.02;                                % mm
w0=57;                                  % mm/s
Va=1000;                                % mm/s
Ta=25+273;                              % k
Etap=0.434*exp(508.3/y(3));             % Pa/s
Nup=Etap/rho;                           % mm2/s
dNup=-Etap*508.3/(y(3)^2*rho);          
a=0.41;
b=0.61;
Nua=10.2+0.131*(y(3)-273);              % mm2/s
Rz=R1*sqrt(w0/y(1));                    % mm
Re=2*y(1)*Rz/Nua;
Cd=a*Re^b;
tau=(1/2)*Cd*rho*y(1)^2;
hc=0.21*Re^(1/3)*Ka/Rz*(1+(8*Va/Rz))^(1/6);
SysO(1)=y(2); 
SysO(2)=y(1)*y(2)/(3*Nup)-(g/(3*Nup))+(2*tau/(3*Nup*rho*Rz))+((y(2)^2)/y(1)*3*Nup)-(1/Nup)*dNup*y(2)*(-2*hc*(y(3)-Ta)/(Rz*rho*Cp*y(1)));
SysO(3)=-2*hc*(y(3)-Ta)/(Rz*rho*Cp*y(1));
endfunction

then i type the initial condition for w,w' and T (y(1),y(2) and y(3) ), called 
Y0:

Y0=[57 1 520];

initializatin for z:

z=linspace(0,1000,2000)';


and resolution using LSODE

y=lsode("SySO",Y0,z);

i don't use any option (i don't know what options are, i'm a beginner)

then i have to plot w and T versus z:

plot(z,y(1,:),z,y(3,:));   

or something similar...   ;)

________________________________________________________________________
Cerchi un laboratorio fotografico aperto 24 ore su 24?
Stampa le tue foto digitali su Kataweb e le ricevi a domicilio in 48 ore.
http://www.kataweb.it/foto

doc1.doc
Description: Binary data