Re: lsode use

[Juan Pablo Carbajal]

On Wed, Aug 26, 2015 at 1:34 PM, Pirson Céline <address@hidden> wrote:
> 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

Céline, it seems you do not abide to the mailing list rules. I wont
continue helping you.
Check the Octave manual
http://www.gnu.org/software/octave/doc/interpreter/