ODEs across a series of rows

Question

SuppoCheng 2016-9-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/302361-odes-across-a-series-of-rows

评论： Star Strider 2016-9-8

Hey, so I'm trying to integrate across a series of rows. Currently, my code spits out an array (Cu_int{1,1}) that has a conditional variable among it. This one graphs perfectly.

However, subsequent arrays (Cu_int{1,2} and (Cu_int{1,3})), do not. Cu_int{1,2} is dependent on the second column of Cu_int{1,1} as the initial value of C_in, that is: Cu_int{1,1}(1,2) = Cu_in for (k1 = 2), which gives a value of Cu_int{1,2}(1,2), which is = Cu_in for (k1 = 3), giving Cu_int{1,3}(1,2).

Thus, Cu_int{1,1}(2,2) to calculate for Cu_int{1,2}(2,2) to calculate for Cu_int{1,3}(2,2).

Sorry for the poor explanation, but I can't seem to wrap my head around getting this loop to work.

Thanks!

clearvars;
close all;
Funderin = [1.7153, 1.052, 0.7085];
Funderout = [1.052, 0.7085, 0.4813];
Fsep = [0.6633, 0.3435, 0.2272];
MassUnders = [14.73, 9.918, 6.738];
StartTime = 0
TimeStep = 0.1
InjectionTime = 3
FinalTime = 50
%tv = linspace(0, 3);    %time vector
tv = StartTime:TimeStep:InjectionTime;
Cu_in = 0.9717  %tracer flow rate kgTr/s
Cui = 0;        %initial concentration
Cu_int = [];
rhs = @(t,Cu,Cu_in,Fin,Fout,Fsp,MassU) (Fin* Cu_in - Fout*Cu - Fsp*Cu)/MassU;   %creates "call" function for ode
for k1 = 1:1
      %Start first half of the integration (inlets)
      %Cu_int{k1} = [];    %this creates the array?
      tspan = tv;
      [t,Cu] = ode45( @(t,Cu) rhs(t,Cu,Cu_in,Funderin(k1),Funderout(k1),Fsep(k1),MassUnders(k1)), tspan, Cui);
      %Start second half of the integration (outlets)
      Cu_int{k1} = [t, Cu];
      Cu_in = 0;
      Cui = Cu(end);  %Cuin uses conc from last 
      %tspan = tv + t(end);
      tspan = (InjectionTime + TimeStep):TimeStep:FinalTime
      [t,Cu] = ode45( @(t,Cu) rhs(t,Cu,Cu_in,Funderin(k1),Funderout(k1),Fsep(k1),MassUnders(k1)), tspan, Cui);
      Cui = 0;
      Cu_in = Cu(end);                                    % ‘Cu_in’ Assignment ?
      Cu_int{k1} = [Cu_int{k1}; t, Cu];
end
for k1 = 2:length(Funderin)
      tspan = StartTime:TimeStep:FinalTime;
      [t, Cu] = ode45( @(t,Cu) rhs(t,Cu,Cu_in,Funderin(k1),Funderout(k1),Fsep(k1),MassUnders(k1)), tspan, Cui);
      for p = 1:(size(Cu_int{1,1},1))
      Cu_in = Cu_int{1,k1-1}(p,2);  
      end
      Cu_int{k1} = [t, Cu]; %this creates the array
end
figure(1)
plot(Cu_int{1}(:,1),Cu_int{1}(:,2));
hold on
for k1 = 2:length(Funderin)
    plot(Cu_int{k1}(:,1),Cu_int{k1}(:,2));
end
hold off
xlabel('Time'); ylabel('Cu');

3 个评论
显示 1更早的评论隐藏 1更早的评论

SuppoCheng 2016-9-8

编辑：SuppoCheng 2016-9-8

It's a tracer study in a tanks-in-series model. What it assumes is, the concentration of tracer (or dye) in the first tank instantly disperses/mixes with the volume of the tank (in this case, MassUnders), setting it as the initial condition for the next tank. Since there is dye being added for a certain period of time, the concentration goes up.

As soon as dye stops being added (peak of the first curve), the concentration will start to diminish, but will continue feeding to the next tank. At the same time, the second tank will disperse the dye before passing it onto the next tank, and so on. This is why the peak concentration of each tank is 'delayed' as shown below:

What I've realised is, I can apply the ODE to each subsequent tank as soon as the first point in the first tank is solved, which is why I'm trying a heap of nested for-loops now.

Here is the original ODE:

This system is a bit of a different case because of the two 'outlet' flow rates. But the C_outlets are assumed to be the same. F_in, F_out, F_sep and Mass are tank specific parameters (these are expressed in the initial arrays), here Fsep dictates Mass in each tank, but that is solved in another piece of code I've got.

Hopefully this explains it much better. This is a chemical engineering topic and I've tried to explain it as best as possible.

SuppoCheng 2016-9-8

I've added a schematic hoping that it will visually describe the process better:

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2016-9-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/302361-odes-across-a-series-of-rows#answer_234204

在 MATLAB Online 中打开

Lose the loop! You have a system of ordinary differential equations, so you have to solve them as that system. This is easy with ode45.

This is as close as I can get to your desired solution:

Funderin = [1.7153, 1, 0.5];
Funderout = [1, 0.5, 0.2];
Fsep = [0.9406, 0.2, 0.1];
MassUnders = [12.0069, 8, 5];
tspan = linspace(0, 50, 250);
t = 0;
Cu_in = 0.9717   %tracer flow
Cui = 0;        %initial concentration
rhs = @(t,Cu,Cu_in,Fin,Fout,Fsp,MassU) [(Fin(1)*Cu(1)/MassU(1) - Fout(1)*Cu(1) - Fsp(1)*Cu(1))/MassU(1); (Fin(2)*Cu(1) - Fout(2)*Cu(2) - Fsp(2)*Cu(2))/MassU(2); (Fin(3)*Cu(2) - Fout(3)*Cu(3) - Fsp(3)*Cu(3))/MassU(3)];
[t,Cu] = ode45( @(t,Cu) rhs(t,Cu,Cu_in,Funderin,Funderout,Fsep,MassUnders), tspan, [Cu_in; 0; 0]);
figure(1)
plot(t, Cu)
grid
xlabel('Time (\mus)')
ylabel('[C_u] (n\itM\rm)')
legend('Compartment #1', 'Compartment #2', 'Compartment #3', 'Location','NE')

It’s missing the initial concentration ‘step’, since I modeled it as a ‘delta’ function instead of a step for convenience. I will defer to you to modeling the step and getting the concentrations to agree with those in your other simulation. Consider adding an initial compartment to model the initial step.

The plot:

This should get you started.

I’ve done compartmental kinetics, specifically with respect to drug and hormone concentration simulations in pharmacokinetic models, but not exactly what you’re doing.

2 个评论
显示无隐藏无

SuppoCheng 2016-9-8

Unfortunately, I still need to loop it. 3 compartments was chosen because it was a small number to deal with. I'm dealing with a system that's supposed to have an 'N' number of compartments, upwards of 300 depending on the mixing behaviour within the system. But this is an absolutely brilliant start.

Are you referring to compartment modelling? My supervisor coined our research project as a 'compartment model', however, I'm under the impression that tanks-in-series is more appropriate. Compartment modelling is very specific in terms of pharmaceutical applications.

Star Strider 2016-9-8

Thank you.

You can do this with a loop, however the loop must be within the ‘rhs’ function. You would have to create a function file to do that, since an anonymous function would only allow the sort of system I used here. A system of 300 ODEs would create an enormous anonymous function.

Each iteration of the loop (one iteration for each ‘tank’ ODE) would correspond to a row of my anonymous function. An easy way to do this would be to use a single anonymous function and then call it in each loop iteration with new constants and the concentrations from the previous calculation. I would set up a different ODE function for the first compartment (since it has to incorporate the tracer injection dynamics), and then use an if block to use it in the first iteration and skip over it in subsequent iterations. This would set up the system of ODEs you need to use to simulate the dynamics of your connected tanks. Your initial conditions vector would have to match the number of compartments you’re modeling, so use the zeros function for all but the first tank.

Yours could easily become a ‘stiff’ system, so using ode15s or one of the other ‘stiff’ solvers would likely work better than ode45.

Your project reminds me of compartmental modeling, although in compartmental modeling the equations are usually set up differently because the dynamics are different.

请先登录，再进行评论。

ODEs across a series of rows

3 个评论
显示 1更早的评论隐藏 1更早的评论

回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

ODEs across a series of rows

3 个评论 显示 1更早的评论隐藏 1更早的评论

回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

2 个评论
显示无隐藏无