Re: lsode use

[Céline]

The problem is that I need to ask to the user of the script some parameters,
So I need to pass some parameters to the function... And I'm really lost, I
don't know how to do that :-(

I post you all the code and try to be clear on what is what since there are
French words
I know that the structure is not the best etc., but I'm really a beginner...

%Definition of the parameters
Xc=1;             %Humidité des levures pour laquelle la fin de la phase de 
                  %séchage à vitesse constante est observée
p=101325;         %Pression totale (= atmosphérique) [Pa]
ke=0.12;          %Coefficient de transfert externe (déterminé exp) [m/s]
Mw=0.018;         %Masse molaire de l'eau [kg/mol]
Ma=0.02896;       %Masse molaire de l'air [kg/mol]
a=6;              %Aire de la surface extérieure des levures par unité de
masse  
                  %de matière sèche [m²/kg]
Rg=8.31;          %Constante universelle des gaz parfaits [J/mol/K]
eps=0.3;          %Porosité du lit, doit être compris entre [0.3;0.6]
tau=5;            %Tortuosité, doit être comprise entre [2;10]
D=0.000025;       %Coefficient de diffusion de la vapeur d'eau dans le gaz
[m²/s]
Lk=2250000;       %Chaleur latente de vaporisation de l'eau [J/kg]
cpa=1000;         %Chaleur spécifique massique de l'air sec [J/K/kg]
Tsat0=373.15;     %Température de référence pour saturation [K]
psat0=101315;     %Pression de saturation de l'eau à la température T0 [Pa]
Lm=2470000*0.018; %Chaleur latente de vaporisation (molaire)

%Variable parameters
%A lot of terms are in comments to test the code rapidly but this is the
part where I ask a lot of parameters to the user
%X0=input("Entrez la teneur initiale en eau des levures:\n");  %(2.45 valeur
ref)
X0=2.45;
%Xr=input("Entrez la teneur residuelle en eau des levures:\n");%(0.1 valeur
ref)
Xr=0.1;
Xc=1;             %Valeur ref d'humidité critique
%R=input("Entrez le rayon moyen des grains de levure (mm):\n");
%R=R/1000;         %Conversion en m
R=0.0005;
%M=input("Entrez la masse de levure non-sechee (kg):\n");
M=1500;
Ms=M/(1+X0);      %Calcul de la masse de solide sec après séchage complet
%G=input("Entrez le debit d'air 'sec' a l'entree du lit fluidise (L/h):\n");
%G=(G/1000)*1.2;   %Conversion en kg/h (38500 valeur ref)
G=38500;
%ts=input("Entrez le temps de sechage total (s):\n");
ts=50;

%Initial conditions
Y0=0.009;         %Humidité initiale de l'air (valeur ref)
T0=273.15+90;     %A CHANGER - Température initiale de l'air

%Initialization of Y and T
Y=Y0;
T=T0;

Jres=[];         %Creation of the result matrix for the evaporation rate
Yres=[Y];        %Creation of the result matrix for air humidity content

%Calculation of the saturation pressure and of the evaporation rate -
initialization
psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0)));
J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));

%Creation of "time matrix" (temps = time)
step=0.01;
temps=[0:step:ts+step]; %Création de la matrice temps - Vecteur ligne
temps=temps.';          %Transposée de la matrice temps --> Vecteur colonne
i=1;

for t=0:step:ts
     x=@(X0,J) humsolide(X0,t,J);
    [X, ISTATE, MSG]=lsode(x,X0,temps);
    psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0)));
    if X(i)>=Xc
      %Evaporation rate calculation
      J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
    else X(i)<=Xc
      %Parameters calculations
      Ri=R*(((X-Xr)/(Xc-Xr))^(1/3));
      ki=((eps*D)/(tau*((R)^2)))*(1/((1/Ri)-(1/R)));
      %Evaporation rate calculation
      J=a*Mw*(1/((1/ke)+(1/ki)))*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y))); 
  end
  Y=((Ms*J)/G)+Y0;
  T=T0-((Ms*J*Lk)/(cpa*G));
  %Adding the results in the result matrix 
  Yres=[Yres;Y];
  Jres=[Jres;J];
  i+=1;
end

%Final aim of this script : (in comment to loose less time while testing)
%plot(temps,Jres);
%hold on;
%plot(temps,Yres);
%hold on;
%plot(temps,X);
  
 

and the other one:

function XDOT=humsolide(X0,t,J)
  XDOT=-J;
  endfunction


Thank you a lot for your time and your help








Cöline, 

Just a formal detail. Please answer at the bottom or between the lines 
(as I am doing now). This makes it easy for readers to follow the 
mailing list archives. 


On Wed, Aug 26, 2015 at 10:36 AM, Céline <[hidden email]> wrote: 

> In fact, my J values are different at each time t and are calculated 
> separately at each time. 
> So I need, each time, to give my J value to my ODE to be able to resolve 
> it... 
> 
> I've tried your advices, but I do probably something wrong since it does
> not 
> work at all... 
> 
> Here is the code of the loop (I've all the initialization of parameters
> and 
> variables before this): 
> 
> for t=0:step:ts 
>      x=@(X0,J) humsolide(X0,t,J); 
>     [X, ISTATE, MSG]=lsode(x,X0,temps); 

«  []

What is temps? 


>     psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0))); 
>     if X(i)>=Xc 
>       %evaporation rate calculation 
>       J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y))); 
>     elseif X(i)<=Xc 
>       %parameters calculation 
>       Ri=R*(((X-Xr)/(Xc-Xr))^(1/3)); 
>       ki=((eps*D)/(tau*((R)^2)))*(1/((1/Ri)-(1/R))); 
>       %evaporation rate calculation 
>      
> J=a*Mw*(1/((1/ke)+(1/ki)))*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y))); 
>   end 
>   Y=((Ms*J)/G)+Y0; 
>   T=T0-((Ms*J*Lk)/(cpa*G)); 
>   %Adding the results in "results matrix" 
>   Yres=[Yres;Y]; 
>   Jres=[Jres;J]; 
>   i+=1; 
> end 
> 
> And my function humsolide: 
> 
> function XDOT=humsolide(X0,t,J) 
>   XDOT=-J; 
>   endfunction 
> 
> Do you know what's wrong ? 
> 
> 

«  []
Form what I understand you have a system of the form 

dx/dt = -J(x) 

With J switching functionla form baed on the values of x. Thisi is an 
hybrid dynamical system and you can opt to solve it with lsode or with 
the odepkg 
I did not check how the to forms of J glue to each other (continuous? 
1st derivative continuous?), but assuming J is smooth enough you can 
use lsode in the following way 

function J = Jvalue (x, params) 
    # params is a struct with all your parameter values. Here comes a 
block where you copy them ot variables of a nice name, if you want. 
Otherwise use, e.g. parasm.Mw 
    if x > params.Xc 
        J = # function for x >= xc 
    else # note that you elseif is not necessary and overlapping 
        J = # function for x < xc 
    endif 
endfunction 

function xdot = humsolide (x, t) 
      xdot = Jvalue(x); 
endfunction 

Or you can integrate everything inside humsolide 

function xdot = humsolide (x, t) 
   # params is a struct with all your parameter values. Here comes a 
block where you copy them ot variables of a nice name, if you want. 
Otherwise use, e.g. parasm.Mw 
    if x > params.Xc 
        J = # function for x >= xc 
    else # note that you elseif is not necessary and overlapping 
        J = # function for x < xc 
    endif 
     xdot = J; 
endfunction 

and pass it to lsode as 

func = @(x,t) humsolide (x,t, params) 

The script will then be 

#block where params is defined 

func = @(x, t) humsolide(x,t,params); 
[X, ISTATE, MSG] = lsode(x,X0, 0:step:ts); 

As you see no need for the for loop. 

If your J function is not smooth, then you should ask the question 

Do I know the switching times a priori? 

If the answer is yes, you can integrate using lsode and passing the 
switching times as the argument T_CRIT of lsode (see documentation). 

If the answer is no, then you should try to use odepkg, solver ode45, 
and use the event function. Please refer to the documentation. 

Good luck 






--
View this message in context: 
http://octave.1599824.n4.nabble.com/lsode-solver-problem-tp4672270p4672291.html
Sent from the Octave - General mailing list archive at Nabble.com.