Re: lsode use

[Juan Pablo Carbajal]

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-tp4672270p4672291.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.