Solving a system of ODE over different intervals (with different conditions on parameters)

Good Evening dear Matlab Community,

Now that my first part of code works, I have the ambition to make it a little complicated but I can't figure out how to implement my plan in matlab:

in the code below you can see that I'm solving a system of two differential equations. Now my challenge is that I would like to solve it with different conditions on certain parameters over the full interval l. I would like to divide the full interval l into sections i each itself divided into two zones, j.

Each section i starts with a conductive zone (j=1) of say l=3 meter for which:

Aap=lcont
Acp=Aap
hap=17
kap=0.018

And ends with a convection zone (j=2) of l=2 meter for which following is true:

Aap=2*lcont
Acp=0
hap=20
kap=0.02

After this, comes the conductive zone of the next section...

I am confused about how to define indices for the sections and zones and whether I need to include loops within the functions corresponding to my ODE in order to solve my system. I have tried (and failed) by introducing following lines within the script of my ODE functions:

for j=1 
Aap=lcont
Acp=Aap
hap=17
kap=0.018
end

for j=2
Aap=2*lcont
Acp=0
hap=20
kap=0.02
end 

If I define the array J=[1 2], how can I make sure that the same indice is used to solve both ODE for each section?

Ideally I would even like to be able to build-in variation for the parameters of the conductive zone between two different sections (like for section 1 and j=1, hap=17 and for section 2 and j=1, hap=5) but I guess once I understand how to implement it for one case, I'll be able to extend it to further! :)

[l,y]=ode45(@myODE,[0 45],[0.48,30]);
plot(l,y(:,:))

function dy = myODE(l,y)
global u
u=y(1);
global Tp
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end

function dudl = myODE2(l,u,Tp)
lcont=2.855;             
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;

dudl=-(kap*Aap)*(Ppw/(Tp+273.15)-0.19);
end

function dTpdl=myODE1(l,u,Tp)
lcont=2.855;    
hcp0=0.65;                 
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
DHs = 0.1*(u^1.1)*((Tp+275.15)^2)*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3*Tp);
DHevap = DHs+DHvap;

dTpdl=((hcp0+0.955*u)*(140-Tp)*Acp+hap*(100-Tp)*Aap+DHevap*(-kap*Aap)*(Ppw/(Tp+273)-0.19))/(1.4+4.18*u);                                                                                           
end

[Jan]

for j=1 
Aap=lcont
Acp=Aap
hap=17
kap=0.018
end

This is no real loop, because it runs once only. The for line is equivalent to j=1, but j is not used inside. Please learn, what loops are: doc for, doc while.

I do not undestand, what you want to achieve. "build-in variation for the parameters of the conductive zone between two different sections" is not clear. The same for: "solve it with different conditions on certain parameters over the full interval l. I would like to divide the full interval l into sections i each itself divided into two zones, j." What are "zones" and "intervals"?

[Lenilein]

Hi Jan,

You are right to refer me to the documentation, the loop as I expressed it doesn't lead anywhere!

Regarding what I want to achieve, my two differential equations are expressed in function of a length element dl (and not time element). They are solved over the entire length of a paper web.

This paper web is for some parts (from l=0 to l=a, l=a+b to l=2a+b, l=2a+2b to l=3a+2b, etc.) in direct contact with drying cylinders (that I called conductive zones) and rest of web (from l=a to l=a+b, l=2a+b to l=2a+2b, l=3a+2b to l=3a+3b, etc.) is in so-called convective zones.

The values of several parameters are different depending if the web is in a conductive or convective zone.

Here a simplified schematic view of the web:

[Jan]

This seems to mean, that you have to split the integration "time" (the physical meaning does not matter) into intervals like in the 1st example in my answer:

tPool = [0,a,a+b,2*a+b, ...];

Now integrate the function in these intervals and use the final value of one integration as initial value to the next one.

Answer 1

[Jan]

Avoid using globals to provide parameters. See <Answers: Anonymous for params>

You can integerate over different time intervals by using a loop:

tPool = [0,1,2,4,8];
y0    = 0;  % Your initial value
tAll  = [];
yAll  = [];
for k = 1:numel(tPool) - 1
    % Your parameter depending on the time step:
    switch k
        case 1
          P = 17;
        case 2
          P = 23;
        case 3
          P = 10992;
        otherwise
          P = 0;
    end
    
    % Call integrator, provide current parameter by anonymous function:
    [t,y] = ode45((t,y) @fcn(t,y,P), [tPool(k:k+1)], y0)
    
    tAll  = cat(1, tAll, t);  % Append new solution to total solution
    yAll  = cat(1, yAll, y);
    
    y0    = y(end, 1);  % Update the initial value
end

If the zones are not predefined by time steps, but by positions, use a while loop until the final time is reached and an event function set by odeset:

t0    = 0;
tEnd  = 8;
y0    = 0;  % Your initial value
tAll  = [];
yAll  = [];
while t0 < tEnd
    % Your parameter depending on the time step:
    if y(1) < 10
      P = 17;
      eventP = 20
    elseif y(1) < 20
      P = 23;
      eventP = 30;
    else   
      P = 10992;
      eventP = 100;
    end
    
    % Create new event function, which stops the integration at the
    % next step:
    opt = odeset('Events', @(t,y) @myEventsFcn(t,y,eventP));
    
    % Call integrator, provide current parameter by anonymous function:
    % tEnd is not reached, because the integration is stopped by the
    % event function except for the last step:
    [t,y] = ode45((t,y) @fcn(t,y,P), [t0, tEnd], y0, opt)
    
    tAll  = cat(1, tAll, t);  % Append new solution to total solution
    yAll  = cat(1, yAll, y);
    
    t0    = t(end);
    y0    = y(end, 1);  % Update the initial value
end

function [position,isterminal,direction] = yourEventFcn(t,y,eventP)
  position   = y(1) - eventP; % The value that we want to be zero
  isterminal = 1;   % Halt integration 
  direction  = 0;   % The zero can be approached from either direction
end

Both examples are just a demonstration and they do not run!!!

[Lenilein]

Thank you Jan, this is very useful!

I'm not sure about the role of tAll and yAll in the first code and I i wonder whether the "@" might be at the wrong place in this expression:

[t,y] = ode45((t,y) @fcn(t,y,P), [tPool(k:k+1)], y0)

I tried to apply your method to my code and I still have some doubts regarding how to integrate the parameters Aap, Acp, hap and kap in the function definition of dudl and dTpdl. Should I add them to the input variables of these functions?

Right now I logically get the error message Too many input arguments. regarding @(l,y)myODE(l,y,Aap,Acp,hap,kap)

I don't exactly understand how to make the link between the three defined functions and the loop.

lcont =2.85;
lconv =1.37;
lPool = [0,lcont,lcont+lconv,2*lcont+lconv,2*lcont+2*lconv,3*lcont+2*lconv];
y0    =[0.48,30];

for k=1:numel(lPool)
    switch k
        case rem(k,2)==1
            Aap=lcont; Acp=Aap; hap=17; kap=0.018;    
        otherwise
            Aap=2*lconv; Acp=0; hap=20; kap=0.02;
    end
    
[l,y]=ode45(@(l,y) myODE(l,y,Aap,Acp,hap,kap),lPool(k:k+1),y0);
y0=y(end,1);

end

plot(l,y(:,:));

function dy = myODE(l,y)
global u
u=y(1);
global Tp
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end

function dudl = myODE2(l,u,Tp)            
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;

dudl=-(kap*Aap)*(Ppw/(Tp+273.15)-0.19);
end

function dTpdl=myODE1(l,u,Tp)
hcp0=0.65;  
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
DHs = 0.1*(u^1.1)*((Tp+275.15)^2)*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3*Tp);
DHevap = DHs+DHvap;

dTpdl=((hcp0+0.955*u)*(140-Tp)*Acp+hap*(100-Tp)*Aap+DHevap*(-kap*Aap)*(Ppw/(Tp+273)-0.19))/(1.4+4.18*u);                                                                                           
end