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!