How to solve ODEs with time dependent parameters by ode45 method

4 views (last 30 days)
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
Deva Narayanan
Deva Narayanan on 17 Sep 2020
i followed your steps in my current problem. but i got following error: "Undefined function 'mtimes' for input arguments of type 'function_handle'". i think , i have a problem to solve the system of odes with time dependent parameters using "interp1" function.

Sign in to comment.

Accepted Answer

Star Strider
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.

More Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!