How to solve a system of two paired differential equations
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Lenilein
el 2 de En. de 2019
Comentada: Lenilein
el 3 de En. de 2019
Good Evening,
After several hours of unsuccessful wandering through the jungle of matlab posts, I decided to call for help regarding how to solve a system of two coupled non linear differential equations of first order.
To summarize the problem I have two variable Tp and u which both depend on l. My ODE system contains the two following equations:
dTpdl=fct(dudl,Tp)
dudl=fct(u,Tp)
(where dTpdl=dTp/dl and dudl=du/dl)
This is how I proceeded:
in a first code (Editor tab) I defined the function dTpdl=myODE1(l,Tp) and wrote the corresponding detailed equation
in the second code I defined the function dudl=myODE2(l,Tp) and wrote the corresponding detailed equation
In a third tab I called ODE45 Matlab to solve this system:
[l,u]=ode45('myODE2',[0 8],0.5);
[l,Tp]=ode45('myODE1',[0 8],30);
With a plot:
plot(l,u,'c','LineWidth',1.5);hold on;
plot(l,Tp,'k','LineWidth',1.5);
When I run the code of my third tab I do obtain a plot with both variables in function of l. However I noticed that if I comment the first line [l,u]=... it doesn't change anything to my plot! => which made me think that matlab probably doesn't use the solution of myODE2 for solving myODE1. Any idea about what I am missing out?
Furthermore, I would like the solutions u to be in the interval [0,1] but I couldn't figure out how to restrict it.
Many many thanks in advance for helping out a struggling beginner!
0 comentarios
Respuesta aceptada
Torsten
el 3 de En. de 2019
function main
[l,y]=ode45(@myODE,[0 8],[0.5,30]);
plot(l,y(:,1))
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
Más respuestas (2)
Star Strider
el 2 de En. de 2019
‘Any idea about what I am missing out?’
We would have to see ‘myODE1’ and ‘myODE2’ to provide meaningful responses.
With respect to your plot question, MATLAB will plot all columns of the second output of your ode45 call, so you do not need to plot each column separately. If you are calculating both in each function, then there will be no change in the output or plotted values by not evaluating one of the functions.
‘Furthermore, I would like the solutions u to be in the interval [0,1] but I couldn't figure out how to restrict it’
If you want to stop the solver at those points, probably the easiest way is use an ‘event function’, and define that point as an ‘event’. See the documentation on ODE Event Location (link) for a thorough discussion.
2 comentarios
Star Strider
el 3 de En. de 2019
@Lenilein — Unfortunately for me, because of the 15-hour delay in responding to my Answer, anything I have to add at this point is now irrelevant.
It was 02:00 in my time zone when it posted.
Lenilein
el 3 de En. de 2019
Editada: Lenilein
el 3 de En. de 2019
4 comentarios
John D'Errico
el 3 de En. de 2019
Editada: John D'Errico
el 3 de En. de 2019
By the way, please stop adding answers just to respond to a question. Use comments.
And, the immediate answer is exactly as Torsten said. If you replace both u and Tp with ones, then ODE45 never uses the values that are passd in for those variables.
In terms of exceeding 1, you first need to verify that your code is written properly. It is entirely reasonable that this is not a problem in the first place. So first verify the code is working. Only then should you consider secondary issues.
Ver también
Categorías
Más información sobre General Applications en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!