[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: lsode use
From: |
Juan Pablo Carbajal |
Subject: |
Re: lsode use |
Date: |
Wed, 26 Aug 2015 13:46:01 +0200 |
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/
- Re: lsode use, (continued)
- Re: lsode use, Doug Stewart, 2015/08/25
- Re: lsode use, Céline, 2015/08/25
- Re: lsode use, Olaf Till, 2015/08/25
- Re: lsode use, Céline, 2015/08/26
- Re: lsode use, Juan Pablo Carbajal, 2015/08/26
- Re: lsode use, Céline, 2015/08/26
- Re: lsode use, Juan Pablo Carbajal, 2015/08/26
- Re: lsode use, Céline, 2015/08/26
- Re: lsode use, Juan Pablo Carbajal, 2015/08/26
- RE: lsode use, Pirson Céline, 2015/08/26
- Re: lsode use,
Juan Pablo Carbajal <=
- Re: lsode use, Céline, 2015/08/26
- Re: lsode use, Juan Pablo Carbajal, 2015/08/26