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

Question

Erin Summerlin-Donofrio 2024-3-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2096491-i-am-trying-to-solve-a-system-of-three-differential-equations-simultaneously

评论： Star Strider 2024-5-8

criss1987.pdf

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.

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

Sam Chak 2024-3-20

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 2024-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2096491-i-am-trying-to-solve-a-system-of-three-differential-equations-simultaneously#answer_1428041

在 MATLAB Online 中打开

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]

.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Erin Summerlin-Donofrio 2024-3-20

编辑：Walter Roberson 2024-3-21

在 MATLAB Online 中打开

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);

Star Strider 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

编辑：Erin Summerlin-Donofrio 2024-3-20

在 MATLAB Online 中打开

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 2024-3-20

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

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 2024-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2096491-i-am-trying-to-solve-a-system-of-three-differential-equations-simultaneously#answer_1428426

编辑：Sam Chak 2024-3-21

在 MATLAB Online 中打开

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

16 个评论
显示 14更早的评论隐藏 14更早的评论

Erin Summerlin-Donofrio 2024-3-25

在 MATLAB Online 中打开

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;

Erin Summerlin-Donofrio 2024-5-3

@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 2024-5-3

在 MATLAB Online 中打开

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 2024-5-7

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?

Star Strider 2024-5-7

在 MATLAB Online 中打开

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 2024-5-8

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 2024-5-8

As always, my pleasure!

请先登录，再进行评论。

Answer 3

Sam Chak 2024-5-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2096491-i-am-trying-to-solve-a-system-of-three-differential-equations-simultaneously#answer_1452106

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

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

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（2 个）

16 个评论
显示 14更早的评论隐藏 14更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

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

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（2 个）

16 个评论 显示 14更早的评论隐藏 14更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

16 个评论
显示 14更早的评论隐藏 14更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论