How to solve a system of two paired differential equations

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!

Answer 1

[Torsten]

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

[Lenilein]

Thank you very much Torsten for helping out! I'll analyse and test the code as soon as I get time.

Answer 2

[Star Strider]

‘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.

[Star Strider]

@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]

I will give you other opportunities to share your knowledge in near future Star Stride ;)

Answer 3

[Lenilein]

Good Morning,

Thanks for giving an answer to my questions and pointing out ODE EventLocation, I think that is what I was looking for!

If I didn't send you my code before it is because I was truely hoping that you would not need to dig into the equations to give me some hint.

Here come myODE1 and myODE2 equations:

function dTpdl=myODE1(l,Tp)
dTpdl=zeros(1,1);

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 drying Acp=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
Tp=ones(1);     %for giving size of matrix and faster solution finding
u=ones(1);

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)*Gf*Acp*(-kap*Aap*Mw*0.001/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273))/((v/60)*Gf*Acp*(cf+cw*u)))); 

function dudl=myODE2(l,Tp)
dudl=zeros(1,1);

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 drying Acp=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 
Tp=ones(1,1);   
u=ones(1);      

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*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273)))

[Torsten]

By setting Tp=ones(1) and u=ones(1) in your ODE functions, you always work with Tp=1 and u=1 in your equations. I think this is not what you intend to do.

Furthermore, I don't see where

dTpdl=fct(dudl,Tp) 

because you don't use dudl to calculate dTpdl.

[Lenilein]

Thank you Torsten for answering!

I'm new to Matlab and still not sure about when the value of a variable is overwritten: I definitely don't want Tp and u to be constants in my equations! Following your comment, I removed those lines. My thought at first was that it would be a good idea to define the size of my array but I obviously used a wrong way.

Regarding your second question I noticed typos in my code, after correction

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))

dudl=-(0.001*(kap*Aap*Mw)/((v/60)*R*Gf*Acp))*((Ppw/(Tp+273))-(Paw/(Ta+273)));

so I can replace part of this expression of dTpdl with dudl as follows:

dTpdl=((hcp0+955*u)*0.001*(Tc-Tp)*Acp+hap*0.001*(Ta-Tp)*Aap+(v/60)*DHevap*dudl)/((v/60)*Gf*Acp*(cf+cw*u))

That's why I wrote in first place that dTpdl=fct(dudl,Tp)

[John D'Errico]

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.

[Lenilein]

my bad John, I realized this is not the best practice shortly before you wrote it!