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

[winkmal]

Naturally, this cannot work, since you call

myODE(l,y,Aap,Acp,hap,kap)

but your function declaration states

function dy = myODE(l,y)

So you are passing more arguments to the function that it has defined in its declaration. You would need to change it to include those parameters, either explicitly or via varargin, and then handle them correctly inside the function.

Also, I am unsure whether it will work to have both ODEs in separate functions, since it seems their states depend on each other!? Have you tried putting them in one single function? I.e., merge all the code from myODE1 and myODE2 into myODE?