RE: lsode use

[Pirson Céline]

On Wed, Aug 26, 2015 at 11:18 AM, Céline <address@hidden> wrote:
> 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-tp4672270p467
> 2291.html Sent from the Octave - General mailing list archive at 
> Nabble.com.
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/help-octave

Céline, please answer at the bottom.

Forget the issue of asking for parameter values. There are many solutions for 
that and many examples you can look at. Lets focus on your problem of 
integrating the dynamical system.
As I see it, there is conceptual mistake in your for loop. You are integrating 
the system many times over the same time interval.
to integrate the system from 0 to ts, you do NOT need the for loop, just a 
single call to lsode with the vector temps = 0:step:ts

The following script is equivalent to yours (as far as I understand) without 
the parameter bamboozlement. Please try it and understand it.

function xdot = dynsys (x,t)
   if x < 0.5
        J = sin (pi * x);
   else
        J = (1 + 0.8) * exp (- (x - 0.5)^2 / (2 * 0.1^2) ) - 0.8;
   endif
   xdot = J;
endfunction

t = linspace (0,2,100).';

zup = lsode (@dynsys, 0.01, t);
zdn = lsode (@dynsys, 1, t);

plot (t, [zup zdn]);
xlabel ("Time")
ylabel ("X")
axis tight


Note that my function J is continuous and has 1st derivative continuous, so 
lsode doesn't have troubles integrating pass x == 0.5.
lsode can handle some discontinuities (check by editing my script) but you 
should do the integration properly and use an event function in the general 
case.


---------------------------------------------------------------------------------------------------------------------------


At least it's giving me something now,
I've put all the parameters in the humsolide function to make it work...
But X0 is also needed in main...

Now I can't have the values of J and Y for all time, I need to do another loop 
in the main script to calculate it back even if I already calculated it for X 
in the function humsolide...
(and I need to plot also these 2 parameters, so as I was calculating for each 
time J and Y, I added the new value each time at the bottom of the line vector 
being the result matrix, but now it seems that I cannot "extract" Jres and Yres 
from the "humsolide" function)

And if I try to use the "input" function, of course he's asking me the question 
each time he does the derivative :-/
I've tried to use the "global" in front of variables and parameters and putting 
it in the main script but it doesn't work then... (humsolide does not have the 
parameters necessary to calculate)

Do you have an idea to have Jres, Yres and to ask the parameters to the users 
properly ?

function xdot=humsolide(X,t)
%All the parameters without "input"

%Variable parameters - PROBLEM
%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;

%Should also be asked to the user
Y0=0.009;         %Humidité initiale de l'air (valeur ref)
T0=273.15+90;     %A CHANGER - Température initiale de l'air

Y=Y0;
T=T0;

psat=psat0*exp(-(Lm/Rg)*((1/T)-(1/T0)));

if X < Xc
    Ri=R*(((X-Xr)/(Xc-Xr))^(1/3));
    ki=((eps*D)/(tau*((R)^2)))*(1/((1/Ri)-(1/R)));
    %Calcul du taux d'évaporation
    J=a*Mw*(1/((1/ke)+(1/ki)))*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
else
    J=a*Mw*ke*((psat/(Rg*T))-(p/(Rg*T))*(Y/((Mw/Ma)+Y)));
endif
xdot=-J;
Y=((Ms*J)/G)+Y0;
T=T0-((Ms*J*Lk)/(cpa*G));

 Endfunction



And main is now much smaller but now I need to give X0 back to this one also...

%X0=input("Entrez la teneur initiale en eau des levures:\n");  %(2.45 valeur 
ref)
X0=2.45;
%ts=input("Entrez le temps de sechage total (s):\n");
ts=50;
t = linspace (0,ts,100).';

X = lsode (@humsolide, X0, t);

plot (t, X);
xlabel ("Time")
ylabel ("X")
axis tight