I am trying to solve a system of three differential equations simultaneously

I am attempting to sovle three differental equations simultaneously. The end goal is to plot the trajectory of (R2, R1) on an x,y plot.

Equation 1: dR1/dt = k1(a1*Rw - R1)

Equation 2: dR2/dt = k2(a2*Rw - R2)

Equation 3: (Xw*dRw)/dt = (Rwin-Rw)*u - (X1(dR1))/dt - (X2(dR2))/dt

ideally I would like to generate code to solve these three equations to where I could define the variables as needed but if this helps,

u=0:1
R1=6
R2=5.5
Rw=-0.5
a1=1.001
a2=0.9998
X1=X2 (these are mole fractions)
K1/K2 = 5 (these are rate constants)
Xw = 0.01 to 0.999 with steps of 0.01, 0.1, 0.2, 0.333, 0.5, 0.666, 0.8, 0.999

Atttached is the original paper for the equations and the way I would like to graph them for reference. The three ways I want to graph the data are "Closed" system, "Buffered" system, and Open System.

Thank you!! If this question belongs somwhere else please let me know.

[Sam Chak]

Can you rearrange the equations to this form?

appears to rise linearly. What is rate of change of when expressed in the form of as a function of time?

Answer 1

[Star Strider]

It ‘sort of’ has an analytic solution (if you want to let it run for a while, and then can make sense of that result), however if you want a numeric result from it, provide values for the variables and integrate it numerically —

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
dRw = diff(Rw);
Equation1 = dR1 == k1*(a1*Rw - R1)

Equation1(t) = 

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation2(t) = 

Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2))

Equation3(t) = 

% S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
% R1 = S.R1
% Rw = S.R2
% Rw = S.Rw

[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3)

VF = 

Subs = 

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2})        % Use This Result With The Numeric ODE Integrator Of Your Choice

DifEqns = function_handle with value:
    @(t,Y,k1,k2,X1,X2,u,Rwin,Xw,a1,a2)[k2.*(a2.*Y(3)-Y(1));k1.*(a1.*Y(3)-Y(2));(Rwin.*u-u.*Y(3)+X1.*k1.*Y(2)+X2.*k2.*Y(1)-X1.*a1.*k1.*Y(3)-X2.*a2.*k2.*Y(3))./Xw]

.

[Erin Summerlin-Donofrio]

If I am following corectly, this makes R1, R2 and R1 into Y1, Y2, and Y3?

I've attempted running the ode but obviously something in the way I have configured the code isn't correct. I'd appreciate any help.

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
dRw = diff(Rw);
Equation1 = dR1 == k1*(a1*Rw - R1);
Equation2 = dR2 == k2*(a2*Rw - R2);
Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2));
S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
R1 = S.R1;
R2 = S.R2;
Rw = S.Rw;
[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3);
DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2});        % Use This Result With 
%this is the ODE example I was trying to follow below
% y'=2*t;
%tspan = [0 5];
%y0 = 0;
%[t,y] = ode45(@(t,y) 2*t, tspan, y0);
%My variables defined
k1=5;
k2=1;
%R1=6.0; This is the initial starting R1
%R2=5.5; This is the initial starting R2
X1=0.8;
X2=0.2;
u=[0 1];
Rwin=-14;
Xw=[0 1];
a1=1.001;
a2=0.9988;
t = [0 1000]; 
Y0=0; %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw
%The ode code
[t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2] = ode45(@(t, Y) DifEqns,  t, Y0);

%plotting 
plot(R2,R1);

[Torsten]

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
dRw = diff(Rw);
Equation1 = dR1 == k1*(a1*Rw - R1);
Equation2 = dR2 == k2*(a2*Rw - R2);
Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2));
[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3);
DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2})  ;   
k1=5;
k2=1;
X1=0.8;
X2=0.2;
u=1;
Rwin=-14;
Xw=1;
a1=1.001;
a2=0.9988;
t = [0 1000]; 
Y0=[0,0,0];
[T,Y]=ode15s(@(t,y)DifEqns(t,y,k1,k2,X1,X2,u,Rwin,Xw,a1,a2),t,Y0);
plot(T,Y)

[Erin Summerlin-Donofrio]

So from the code output I see that the three Y values give the correct shape of the curve I'm looking for if I plot Y1 as x and Y2 as y -- Thank you!!

How could I modify this to where the initial starting point for Y1= 5.5 and for Y2 = 6.0 instead of starting at zero? The curve is for Y2, Y1 (which was initially R2, R1).

[Star Strider]

First, you included the arguments in the output of the ode45 call that should be inputs to the function, so I changed that. (Correct idea, wrong implementation.)

Second, ‘u’ can be 0 (although that give a 0 output for everything), however ‘Xw’ cannot, because if it is, everything is either Inf (with u~=0) or NaN (with u=0) because of the divide-by-zero.

Third, because there are 3 differential equations, there have to be 3 initial conditions. I changed ‘Y0’ to provide them.

With both ‘u’ and ‘Xw’ set to be equal to 1 (or anything other than 0), you get acceptable results —

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
dRw = diff(Rw);
Equation1 = dR1 == k1*(a1*Rw - R1);
Equation2 = dR2 == k2*(a2*Rw - R2);
Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2));
 % S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
 % R1 = S.R1;
 % R2 = S.R2;
 % Rw = S.Rw;
[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3)

VF = 

Subs = 

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2});        % Use This Result With 
%this is the ODE example I was trying to follow below
% y'=2*t;
%tspan = [0 5];
%y0 = 0;
%[t,y] = ode45(@(t,y) 2*t, tspan, y0);
%My variables defined
k1=5;
k2=1;
%R1=6.0; This is the initial starting R1
%R2=5.5; This is the initial starting R2
X1=0.8;
X2=0.2;
u=[0 1];
u = u(2);
Rwin=-14;
Xw=[0 1];
Xw = Xw(2);
a1=1.001;
a2=0.9988;
t = [0 1000]; 
Y0 = zeros(3,1); %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw
%The ode code
[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2),  t, Y0);
R2 = Y(:,1);
R1 = Y(:,2);
Rw =  Y(:,3);

% SizeY = size(Y)
% NrNaN = nnz(isnan(Y))

%plotting 
semilogx(t, Y);
grid
xlabel('Time')
ylabel('Amplitude')
legend(string(Subs))

figure
plot(R2,R1)
grid
xlabel('R_2')
ylabel('R_1')

Everything now works!

.

[Erin Summerlin-Donofrio]

Thank you!!!

I do have a couple of questions if you don't mind.

Is this syntax choosing the second value for u? also for Xw in the above code. My apologies, for Xw I meant very close to zero, not zero. As in the above question I'd like to vary u and Xw (from 0.005 to 0.9999, for example).

u=[0 1];
u = u(2);

And this is used to call R1 as Y1, etc?

R2 = Y(:,1);
R1 = Y(:,2);
Rw =  Y(:,3);

Thank you!

[Star Strider]

My pleasure!

The order of the differential equations is: [R2 R1 Rw], so the initial conditions must match that order. That means:

Y0 = [5.5; 6.0; 0];

To vary ‘u’ and ‘Xw’ it will be necessary to run the code separately for each combination of them. Since each are (1x2) vectors, running the code 4 times, once with each combination. You can do this in two nested loops. If there are more than 2 values for each, the code will have to be run for each combination. Again, nested loops.

for k1 = 1:numel(u)

for k2 = 1:numel(Xw)

[t{k1,k2}, Y{k1,k2}] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, X1, X2, u(k1), Rwin, Xw(k2), a1, a2), t, Y0);

end

end

This creates cell arrays for ‘t’ and ‘Y’ and saves them for each run.

To then get ‘Y’ for the second value of ‘u’ and the first value of ‘Xw’ the addressing would be:

Y21 = Y{2,1};

If you want them all to have the same row lengths, define:

t = linspace(0, 1000, N);

where ‘N’ can be any reasonable number (500, 1000, 1500 ...). That may make it easier to work with them, and it only requires one reference toto the time vector.

This:

R2 = Y(:,1);

R1 = Y(:,2);

Rw = Y(:,3);

or, keeping with the ongoing example:

R2 = Y21(:,1);

R1 = Y21(:,2);

Rw = Y21(:,3);

selects the varioous values from the ‘Y’ matrix so you can plot them individually. That is easier than having to refer to each column each time you need to use it.

.

Answer 2

[Sam Chak]

Hi @Erin Summerlin-Donofrio

The third differential equation can be rearranged conveniently, allowing all three state equations to be incorporated into the 'odeModel()' function, which can be called by the ode45 solver. The remaining task is to ensure the correct input of parameters. Here is a code snippet for reference. You may need to design the flow signal u (also known as the manipulated variable) to achieve the desired output responses in , , and .

Edit: The code has been updated with the provided parameters and has been slightly modified using a 'for-loop' method to simulate the dynamical system iteratively for various values of the parameter , following similar approach demonstrated in Criss et al. (1987).

%% Settings
tspan   = [0, 10];              % simulation time
R0      = [6.0; 5.5; -0.5];     % initial values of R1, R2, Rw
Xw      = [0.01, 0.1, 0.2, 1/3, 0.5, 2/3, 0.8, 0.999];  % 8 elements in Xw
i       = 1;                    % for placing text position
a       = 1;                    % for placing text position
b       = 60;                   % for placing text position

%% Call ode45 to solve odeModel for different values of Xw
for j   = 1:numel(Xw)           % looping 8 times
    hold on,                    % hold current figure
    [t, R]  = ode45(@(t, R) odeModel(t, R, Xw(j)), tspan, R0);
    plot(R(:,2), R(:,1))        % plot R1 (y-axis) vs. R2 (x-axis)
    Ra      = R(:,2);           % for placing text position
    Rb      = R(:,1);           % for placing text position
    text(Rb(-i*a+b), Ra(-i*a+b), "Xw="+string(Xw(j)), 'FontSize', 8)
    i       = i + 1;            % for placing text position
end
axis([2 6.5 1.5 6.5])
grid on, xlabel('R_{2}(t)'), ylabel('R_{1}(t)')
title('Mineral–Water Oxygen Isotope Exchange')

%% Dynamics of Mineral–Water Oxygen Isotope Exchange
function dRdt = odeModel(~, R, Xw)
    % definitions
    R1  = R(1);
    R2  = R(2);
    Rw  = R(3);
    
    % parameters
    k1  = 5;
    k2  = 1;
    a1  = 1.001;
    a2  = 0.9988;
    Rwin=-0.5;
    X1  = 0.5;
    X2  = 0.5;
    u   = 0;
    
    % state equations
    dR1 = k1*(a1*Rw - R1);
    dR2 = k2*(a2*Rw - R2);
    dRw = ((Rwin - Rw)*u - X1*dR1 - X2*dR2)/Xw;
    
    % derivative vector
    dRdt= [dR1; 
           dR2; 
           dRw];
end

[Erin Summerlin-Donofrio]

So below is the code I modified with the values I want for the variables. That appears to work for any choice of Rw between 0.01 and 1, as well as for u=0 and u=1.

But if I want to graph (R2, R1) how would I modify the R0 line to do so? I tried plot(R2, R1) before end but it didn't work.

%% Call ode45
tspan    = [0, 10];      % simulation time
R0       = [6.0; 5.5; -0.5];    % initial values of R1, R2, Rw
[t, R]   = ode45(@odeModel, tspan, R0);
% [t, R] = ode45(@(t, R) odeModel(t, R, Xw), tspan, R0);    % for iteration
%% Plot result
plot(t, R), grid on, xlabel('t'), ylabel('{\bf R}(t)')
legend('R_{1}', 'R_{2}', 'R_{w}', 'fontsize', 16)

function dRdt = odeModel(~, R)
    % definitions
    R1  = R(1);
    R2  = R(2);
    Rw  = R(3);
    
    % parameters
    k1  = 5;
    k2  = 1;
    a1  = 1.001;
    a2  = 0.9988;
    Rwin= -0.5;
    X1  = 0.5;
    X2  = 0.5;
    Xw  = 0.1;
    u   = 0;
    
    % state equations
    dR1 = k1*(a1*Rw - R1);
    dR2 = k2*(a2*Rw - R2);
    dRw = ((Rwin - Rw)*u - X1*dR1 - X2*dR2)/Xw;
    
    % derivative vector
    dRdt= [dR1; 
           dR2; 
           dRw];
  
end

[Sam Chak]

Hi @Erin Summerlin-Donofrio

I have made the necessary updates to the code in my previous answer, incorporating the provided parameters. The resulting code will generate a graph similar to the one demonstrated in Criss et al. (1987). It is important to note that if the differential equations are algebraically more complex, the approach demonstrated by Star Strider using 'odeToVectorField()' is the recommended method.

For the remaining task, you may need to determine how to design the flowrate u in order to achieve your desired outcomes, as I am not familiar with the specifics of the isotope exchange system.

[Erin Summerlin-Donofrio]

Absolutely. Thank you editing the code and for providing the graph! I was getting the same answer running the sequence one by one for different Xw but the shape of the curve as Xw approaches 1 didn't match the figure in the paper so I thought I was doing something incorrectly. I think that was a mistake on the author's part. I appreciate you @Star Strider and @Sam Chak both taking the time to help me and for providing the two differnt approaches to the code. I'm trying to learn more syntax and methods for constructing code so offering different approaches to arrive at the same answer is so incredibly helpful.

Thank you so much, again!!

[Erin Summerlin-Donofrio]

So I'm running into an issue somewhere for the correct values of R1 and R2. For Y0 we assign R1 = 6.0 in the first column, second column as R2 = 5.5, and third column Rw = -0.5 as the starting values for Y0. Under plot, R1= (Y:,1) R2= (Y:,2) and Rw= (Y:,3). So that all seems to check out. However when I run the code, it plots R2 as starting at 6 and R1 as starting at 5.5 if you zoom in. Could you help me figure out why? Thanks!

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
dRw = diff(Rw);
Equation1 = dR1 == k1*(a1*Rw - R1)

Equation1(t) = 

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation2(t) = 

Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2))

Equation3(t) = 

%S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
%R1 = S.R1
%R2 = S.R2
%Rw = S.Rw

[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3)

VF = 

Subs = 

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2})        % Use This Result With The Numeric ODE Integrator Of Your Choice

DifEqns = function_handle with value:
    @(t,Y,k1,k2,X1,X2,u,Rwin,Xw,a1,a2)[k2.*(a2.*Y(3)-Y(1));k1.*(a1.*Y(3)-Y(2));(Rwin.*u-u.*Y(3)+X1.*k1.*Y(2)+X2.*k2.*Y(1)-X1.*a1.*k1.*Y(3)-X2.*a2.*k2.*Y(3))./Xw]

k1=5;
k2=1;
X1=0.5;
X2=0.5;
u=0;
Rwin=-0.5;
Xw=1;
a1=1.001;
a2=0.9988;
t = [0 10]; 
Y0=[6.0,5.5,-0.5];
[T,Y]=ode15s(@(t,y)DifEqns(t,y,k1,k2,X1,X2,u,Rwin,Xw,a1,a2),t,Y0);


plot(T,Y,LineWidth=1.5)
R1 = Y(:,1);
R2 = Y(:,2);
Rw =  Y(:,3);
% SizeY = size(Y)
% NrNaN = nnz(isnan(Y))
%plotting 


xlabel('Time')
ylabel('Amplitude')
legend(string(Subs))


box on
ax = gca;
ax.LineWidth = 1.5;
ax.FontSize = 20; 

[Star Strider]

‘Could you help me figure out why?’

The order of the differential equations are their order in the ‘Subs’ result, that being in order. So the initial conditions have to be [5.5 6.0 -0.5] to match that ordering.

[Erin Summerlin-Donofrio]

AAH okay thank you!

[Star Strider]

As always, my pleasure!

[Erin Summerlin-Donofrio]

@Star Strider

syms R1(t) R2(t) Rw(t) k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

dRw = diff(Rw);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw*dRw) == (Rwin-Rw)*u - (X1*(dR1)) - (X2*(dR2))

% S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10, R2(0)==R20, Rw(0)==Rw0);

% R1 = S.R1

% Rw = S.R2

% Rw = S.Rw

[VF,Subs] = odeToVectorField(Equation1, Equation2, Equation3)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2})

Y0 = zeros(3,1); %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw

%The ode code

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, X1, X2, u, Rwin, Xw, a1, a2), t, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

Rw = Y(:,3);

How would I modify this code to treat Rw as a constant, where Rw=Rwin? Modify equation 3 somehow? I've attempted it myslef by taking out dRw/dt as a derivative term, but I keep getting the "Unable to convert the initial value problem to an equivalent dynamical system. Either the differential equations cannot be solved for the highest derivatives or inappropriate initial conditions were specified" error.

Thank you!

[Star Strider]

The problem with making as a constant is that ‘Equation3’ is no longer a differential equation. The system is also not a differential algebraic equation because of the remaining derivative terms in ‘Equation3’. The only way to make to be a constant is to eliminate ‘Equation3’ and that means elliminating and . I also experimented with using dsolve on ‘Equation3’ and then using subs with those solutions in ‘Equation1’ and ‘Equation2’, however that also failed. UInless you reformulate ‘Equation1’ and ‘Equation2’ to include and , it may not be possible to integrate your equations and get a reasonable result.

Torsten may have some insight into this, however I do not. (I experimented with this for a few hours.)

[Erin Summerlin-Donofrio]

The resulting equation shown in the paper I'm trying to follow has eliminated X1 and X2. I will try and reformlate equation three without derivatives. Thank you for trying for so long - I very much appreciate it!

[Star Strider]

If you’re eliminating and , you don’t need ‘Equation3’ (unless there’s something I’m failing to understand, since ‘Equation3’ is then a simple algebraic equation that then appears to be irrelevant to the rest of the system, unless you want to solve it for and then substitute it for in the other equations).

In that event, and completely neglecting ‘Equation3’, this should work —

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation1(t) = 

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation2(t) = 

Equation3 = (Xw) == (Rwin-Rw)*u

Equation3 = 

% S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
% R1 = S.R1
% Rw = S.R2
% Rw = S.Rw
[VF,Subs] = odeToVectorField(Equation1, Equation2)

VF = 

Subs = 

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        

DifEqns = function_handle with value:
    @(t,Y,k1,k2,u,Rw,Rwin,Xw,a1,a2)[k2.*(Rw.*a2-Y(1));k1.*(Rw.*a1-Y(2))]

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

Y0 = zeros(2,1); %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw
%The ode code
cv = num2cell(rand(1,8));
[k1, k2, u, Rw, Rwin, Xw, a1, a2] = cv{:}

k1 = 0.5811

k2 = 0.5237

u = 0.7918

Rw = 0.7721

Rwin = 0.0410

Xw = 0.9582

a1 = 0.2903

a2 = 0.4877

tspan = [0 15];
[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2),  tspan, Y0);
R2 = Y(:,1);
R1 = Y(:,2);

figure
plot(t, Y)
grid
legend(string(Subs), 'Location','best')

.

[Erin Summerlin-Donofrio]

Yes this is what I was trying to achieve! I think I was confused about what to do for with the derivatives at the end of equation 3 (dR1 and dR2). I have copied the full code above and used hold on to plot curves with different Rw and Rwin values.

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw) == (Rwin-Rw)*u

[VF,Subs] = odeToVectorField(Equation1, Equation2)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

Y0 = [5.5,6.0];

k1=5;

k2=1;

u=1000;

Rw=-0.5;

Rwin=-0.5;

Xw=1;

a1=1;

a2=1;

tspan = [0 15];

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2), tspan, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

figure

plot(R2,R1);

box on

hold on

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw) == (Rwin-Rw)*u

[VF,Subs] = odeToVectorField(Equation1, Equation2)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

Y0 = [5.5,6.0];

k1=5;

k2=1;

u=1000;

Rw=5;

Rwin=5;

Xw=1;

a1=1;

a2=1;

tspan = [0 15];

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2), tspan, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

plot(R2,R1);

hold on

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw) == (Rwin-Rw)*u

[VF,Subs] = odeToVectorField(Equation1, Equation2)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

Y0 = [5.5,6.0];

k1=5;

k2=1;

u=1000;

Rw=-5;

Rwin=-5;

Xw=1;

a1=1;

a2=1;

tspan = [0 15];

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2), tspan, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

plot(R2,R1);

hold on

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw) == (Rwin-Rw)*u

[VF,Subs] = odeToVectorField(Equation1, Equation2)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

Y0 = [5.5,6.0];

k1=5;

k2=1;

u=1000;

Rw=-15;

Rwin=-15;

Xw=1;

a1=1;

a2=1;

tspan = [0 15];

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2), tspan, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

plot(R2,R1);

hold on

syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

dR1 = diff(R1);

dR2 = diff(R2);

Equation1 = dR1 == k1*(a1*Rw - R1)

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation3 = (Xw) == (Rwin-Rw)*u

% S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10, R2(0)==R20, Rw(0)==Rw0);

% R1 = S.R1

% Rw = S.R2

% Rw = S.Rw

[VF,Subs] = odeToVectorField(Equation1, Equation2)

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y

%Y0 = zeros(2,1); %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw

Y0 = [5.5,6.0];

%The ode code

% cv = num2cell(rand(1,8));

% [k1, k2, u, Rw, Rwin, Xw, a1, a2] = cv{:}

k1=5;

k2=1;

u=1000;

Rw=2.5;

Rwin=2.5;

Xw=1;

a1=1;

a2=1;

tspan = [0 15];

[t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2), tspan, Y0);

R2 = Y(:,1);

R1 = Y(:,2);

plot(R2,R1);

How could I draw a line between points of equal time at each curve? For example at t=15, the end points of the curves would be connected by a line. But what about for t=0.25, or t=1? I have been doing this with previous calculations by running the code again and choosing a single t for each curve and then plotting a line through the two x,y for each chosen value of t. Is there a more efficent way to accomplish this?

[Torsten]

Is there a more efficent way to accomplish this?

No.

[Star Strider]

It took a while to figure out how to get this to work optimally. I also made your code a bit more efficient.

Try this —

% syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% dR1 = diff(R1);
% dR2 = diff(R2);
% Equation1 = dR1 == k1*(a1*Rw - R1)
% Equation2 = dR2 == k2*(a2*Rw - R2)
% Equation3 = (Xw) == (Rwin-Rw)*u
% [VF,Subs] = odeToVectorField(Equation1, Equation2)
% DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        
% clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% Y0 = [5.5,6.0];
% k1=5;
% k2=1;
% u=1000;
% Rw=-0.5;
% Rwin=-0.5;
% Xw=1;
% a1=1;
% a2=1;
% tspan = [0 15];
% [t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2),  tspan, Y0);
% R2 = Y(:,1);
% R1 = Y(:,2);
% figure
% plot(R2,R1);
% box on 
% hold on 
% syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% dR1 = diff(R1);
% dR2 = diff(R2);
% Equation1 = dR1 == k1*(a1*Rw - R1)
% Equation2 = dR2 == k2*(a2*Rw - R2)
% Equation3 = (Xw) == (Rwin-Rw)*u
% [VF,Subs] = odeToVectorField(Equation1, Equation2)
% DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        
% clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% Y0 = [5.5,6.0];
% k1=5;
% k2=1;
% u=1000;
% Rw=5;
% Rwin=5;
% Xw=1;
% a1=1;
% a2=1;
% tspan = [0 15];
% [t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2),  tspan, Y0);
% R2 = Y(:,1);
% R1 = Y(:,2);
% plot(R2,R1);
% hold on 
% syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% dR1 = diff(R1);
% dR2 = diff(R2);
% Equation1 = dR1 == k1*(a1*Rw - R1)
% Equation2 = dR2 == k2*(a2*Rw - R2)
% Equation3 = (Xw) == (Rwin-Rw)*u
% [VF,Subs] = odeToVectorField(Equation1, Equation2)
% DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        
% clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% Y0 = [5.5,6.0];
% k1=5;
% k2=1;
% u=1000;
% Rw=-5;
% Rwin=-5;
% Xw=1;
% a1=1;
% a2=1;
% tspan = [0 15];
% [t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2),  tspan, Y0);
% R2 = Y(:,1);
% R1 = Y(:,2);
% plot(R2,R1);
% hold on 
% syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% dR1 = diff(R1);
% dR2 = diff(R2);
% Equation1 = dR1 == k1*(a1*Rw - R1)
% Equation2 = dR2 == k2*(a2*Rw - R2)
% Equation3 = (Xw) == (Rwin-Rw)*u
% [VF,Subs] = odeToVectorField(Equation1, Equation2)
% DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        
% clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
% Y0 = [5.5,6.0];
% k1=5;
% k2=1;
% u=1000;
% Rw=-15;
% Rwin=-15;
% Xw=1;
% a1=1;
% a2=1;
% tspan = [0 15];
% [t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2),  tspan, Y0);
% R2 = Y(:,1);
% R1 = Y(:,2);
% plot(R2,R1);
% hold on 
syms R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
dR1 = diff(R1);
dR2 = diff(R2);
Equation1 = dR1 == k1*(a1*Rw - R1)

Equation1(t) = 

Equation2 = dR2 == k2*(a2*Rw - R2)

Equation2(t) = 

Equation3 = (Xw) == (Rwin-Rw)*u

Equation3 = 

% S = dsolve(Equation1, Equation2, Equation3, R1(0)==R10,  R2(0)==R20,  Rw(0)==Rw0);
% R1 = S.R1
% Rw = S.R2
% Rw = S.Rw
[VF,Subs] = odeToVectorField(Equation1, Equation2)

VF = 

Subs = 

DifEqns = matlabFunction(VF, 'Vars',{t, Y, k1, k2, u, Rw, Rwin, Xw, a1, a2})        

DifEqns = function_handle with value:
    @(t,Y,k1,k2,u,Rw,Rwin,Xw,a1,a2)[k2.*(Rw.*a2-Y(1));k1.*(Rw.*a1-Y(2))]

clear R1(t) R2(t) Rw k1 k2 X1 X2 u Rwin Xw a1 a2 R10 R20 Rw0 Y
%Y0 = zeros(2,1); %Not sure what to do here as we now have Y1, Y2, and Y3 in place of R1, R2, Rw
Y0 = [5.5,6.0];
%The ode code
% cv = num2cell(rand(1,8));
% [k1, k2, u, Rw, Rwin, Xw, a1, a2] = cv{:}
k1=5;
k2=1;
u=1000;
Rw=2.5;
Rwin=2.5;
Xw=1;
a1=1;
a2=1;
tspan = [0 15];

Rwv = [-0.5; -5; -15; 2.5];
Rwinv = [-0.5; -5; -15; 2.5];

for k = 1:numel(Rwv)
    [t, Y] = ode45(@(t, Y) DifEqns(t, Y, k1, k2, u, Rwv(k), Rwinv(k), Xw, a1, a2),  tspan, Y0);
    tc{k} = t;
    R2c{k} = Y(:,1);
    R1c{k} = Y(:,2);
end

tq = [0:0.5:2 15];

figure
hold on
for k = 1:numel(tc)
    hp2(k) = plot(R2c{k},R1c{k}, 'DisplayName',"k = "+k);
    R2R1 = [R2c{k} R1c{k}];
    R2R1i = interp1(tc{k}, R2R1, tq);
    plot(R2R1i(:,1), R2R1i(:,2), 's')
    R2R1c{k} = R2R1i;
end



R2R1m = cell2mat(R2R1c);
for k = 1:size(R2R1m,1)
    R2R1v{k} = reshape(R2R1m(k,:), 2, []);
end
for k  = 1:size(R2R1c{1},1)
    plot(R2R1v{k}(1,:), R2R1v{k}(2,:))
end
grid on
legend(hp2, 'Location','best')
% R2R1c{3}
text(R2R1c{3}(:,1), R2R1c{3}(:,2), compose('t = %.1f', tq), 'Horiz','left', 'Vert','middle', 'Rotation',45)

figure
hold on
for k = 1:numel(tc)
    hp3(k) = plot3(R2c{k}, R1c{k}, tc{k}, 'DisplayName',"k = "+k);
    plot3(R2R1c{k}(:,1), R2R1c{k}(:,2), tq, 's')
end
grid on
xlabel('R_2')
ylabel('R_1')
zlabel('Time')
legend(hp3, 'Location','best')
view(-27,30)

I use the 3D plot to understand what your integrations were doing. I’m leaving it in for now. You can delete it, since it is not important for the code.

.

[Erin Summerlin-Donofrio]

THANK YOU!!

Being able to plot it in 3D is actually really useful for what I'm trying to model and visualize!

Once again I SO appreciate all of your time and effort!!

[Star Strider]

As always, my pleasure!

Answer 3

[Sam Chak]

Hi @Erin Summerlin-Donofrio

This is merely a test. From the perspective of a first-order system, the variable u serves as an inverse of the time constant. Consequently, when the influx of water u is very large, the time constant of system is very small, and the state 'Rw' rapidly converges to the desired reference 'Rwin', nearly instantaneously, as if 'Rw' were constant irrespective of its initial condition.

Here is a demo:.

tspan   = [0 10];
R0      = [4; 3; 2];
[t, R]  = ode45(@(t, R) odefun(t, R, 1, 1, 1, 1, 1, 1000, 1, 1, 1),  tspan, R0);
plot(t, R), grid on, xlabel('t'), legend('R_1', 'R_2', 'R_w', 'fontsize', 14)

%% Buffer System
function dRdt = odefun(t, R, a1, a2, k1, k2, Rwin, u, X1, X2, Xw)
    dRdt    = zeros(3, 1);
    Rw      = R(3);
    dRdt(1) =  - k1*(R(1) - a1*Rw);  % k1 exponential rate response of R1(t) --> a1*Rw if Rw is bounded
    dRdt(2) =  - k2*(R(2) - a2*Rw);  % k2 exponential rate response of R2(t) --> a2*Rw if Rw is bounded
    dRdt(3) = (-  u*(R(3) -  Rwin) - X2*(- k2*(R(2) - a2*Rw)) - X1*(- k1*(R(1) - a1*Rw)))/Xw;
end