How to solve ODEs with time dependent parameters by ode45 method
4 views (last 30 days)
Show older comments
Deva Narayanan
on 17 Sep 2020
Answered: Star Strider
on 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 Comments
Accepted Answer
Star Strider
on 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 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!