How to solve ODEs with time dependent parameters by ode45 method
显示 更早的评论
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 个评论
Star Strider
2020-9-17
What about that solution does not work in your current problem?
Deva Narayanan
2020-9-17
采纳的回答
Star Strider
2020-9-17
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Structural Mechanics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
