How to solve ODEs with time dependent parameters by ode45 method
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Deva Narayanan
le 17 Sep 2020
Réponse apportée : Star Strider
le 17 Sep 2020
i have found a matlab code in https://in.mathworks.com/ . it works correctly using constant value of Ta=100; Tc=150, v=6.25; . but i want to change these values with respect to l as given in excel file.Is it possible to change this matlab code with respect to l ?
function main
tspan=0:1:8;
IC = [0.5,30];
[l,y]=ode45(@myODE,tspan,IC);
disp([l,y]);
plot(l,real(y(:,1)),'b',l,real (y(:,2)),'r');
title('Solution of cylinder');
xlabel('Time l');
ylabel('Solution y');
legend(' u ','Tp ')
end
function dy = myODE(l,y)
u=y(1);
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end
function dTpdl=myODE1(l,u,Tp)
w=4.23; %width of the paper web in m
lcont=3.01; %length of a contact heat transfer section
Mw=18; %Molar weight of water in g/mol
R=8.314; %Ideal gas Constant in kJ/(kmol*K)
v=6.25; %machine speed in m/min
Gf=0.12; %bone-dry-basis weight of sheet in kg/m2; provided information on drying section: Base weight (g/m2) = 120
cf=1.423; %specific heat of fiber in kJ/(kg.K)
cw=4.1868; %specific heat of water in in kJ/(kg.K)
hcp0=650; %heat transfer parameter between cylinder and sheet in W/(m2.K)
hap=17.7; %heat transfer parameter between air and sheet in W/(m2.K) taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
kap=0.0176; %Transfer parameter between air and sheet in m/s taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
Acp=w*lcont; %contact areas between cylinder and the sheet; in contact dryingAcp=mf/Gf, for convective drying, Acp=0
Aap=Acp ; %contact areas between the sheet and air; in contact drying Aap=mf/Gf, for convective drying, Aap=2mf/Gf
Ta=100; %Air temperature in °C => Temperature of pocket ventilation air [constant]
Tc=150; %Cylinder temperature in °C [constant]
Xa=0.4; %Humidity of air
Psata = 133.322*exp(18.3036-3816.44/(Ta+227.03));
Psatp = 133.322*exp(18.3036-3816.44/(Tp+227.03));
Phi= 1-exp(-(47.58*8^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi ;
Paw= Psata*Xa ; %where Xa is the relative humidity of air
DHs = (R/Mw)*0.10085*(u^1.0585)*(Tp+275.15)^2*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3237*Tp);
DHevap = DHs+DHvap;
dTpdl=((hcp0+955*u)*0.001*(Tc-Tp)*Acp+hap*0.001*(Ta-Tp)*Aap+(v/60)*DHevap*(-kap*Aap*Mw*0.001/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273))))/((v/60)*Gf*Acp*(cf+cw*u));
end
function dudl = myODE2(l,u,Tp)
w=4.23; %width of the paper web in m
lcont=3.01; %length of a contact heat transfer section
Mw=18; %Molar weight of water in g/mol
R=8.314; %Ideal gas Constant in kJ/(kmol*K)
v=6.25; %machine speed in m/min
Gf=0.12; %bone-dry-basis weight of sheet in kg/m2; provided information on drying section: Base weight (g/m2) = 120
cf=1.423; %specific heat of fiber in kJ/(kg.K)
cw=4.1868; %specific heat of water in in kJ/(kg.K)
hcp0=650; %heat transfer parameter between cylinder and sheet in W/(m2.K)
hap=17.7; %heat transfer parameter between air and sheet in W/(m2.K) taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
kap=0.0176; %Transfer parameter between air and sheet in m/s taken for lcont1 from model input data// The temperature and the humidity of air change are far smaller than the above interval in the actual drying process, so we can set hap, kap as a constant to reduce the complexity of the model.
Acp=w*lcont; %contact areas between cylinder and the sheet; in contact dryingAcp=mf/Gf, for convective drying, Acp=0
Aap=Acp ; %contact areas between the sheet and air; in contact drying Aap=mf/Gf, for convective drying, Aap=2mf/Gf
Ta=100; %Air temperature in °C => Temperature of pocket ventilation air [constant]
Tc=150; %Cylinder temperature in °C [constant]
Xa=0.4; %Humidity of air
Psata = 133.322*exp(18.3036-3816.44/(Ta+227.03));
Psatp = 133.322*exp(18.3036-3816.44/(Tp+227.03));
Phi= 1-exp(-(47.58*8^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi ;
Paw= Psata*Xa ; %where Xa is the relative humidity of air
DHs = (R/Mw)*0.10085*(u^1.0585)*(Tp+275.15)^2*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3237*Tp);
DHevap = DHs+DHvap;
dudl=-(0.001*(kap*Aap*Mw)/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273)));
end
2 commentaires
Réponse acceptée
Star Strider
le 17 Sep 2020
I do not see any interp1 calls in the code you posted.
To use a function in a calculation, the function needs to be evaluated.
In my previous code:
interptransm = @(t) interp1(T1.t, T1.transm, t, 'linear','extrap'); % Interpolate
interprecov = @(t) interp1(T1.t, T1.recov, t, 'linear','extrap'); % Interpolate
%so here fun = @(t,y)[S'; I'; R'];
fun = @(t,y)[-interptransm(t)*y(1)*y(2); (interptransm(t)*y(1)*y(2))-(interprecov(t)*y(2)); interprecov(t)*y(2)];
↑ ↑ ↑ ↑ ← EVALUATE THEM HERE
It is necessary to provide the current ‘t’ value passed to the ODE function to each interpolation function in order to get the interpolated result returned. Using the function handle without evaluating it will throw the error you reported.
0 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!