4 views (last 30 days)

Show older comments

V=1000; Q=50; Ca0=1; k=1;

for a=0:1:0.1

f=@(x) [Q*Ca0-Q*x(1)-k*x(1)^2*(a*V); Q*x(1)-Q*x(2)-k*x(2)^2*(1-a)*V];

fsolve(f,[0.5,0.5])

end

plot(a,x(2))

Alan Stevens
on 16 Sep 2021

Your first equation is a simple quadratic in x(1); your second is a quadratic in x(2) that depends on x(1), so, assuming you are only interested in the positive roots, these can be solved as follows:

V=1000; Q=50; Ca0=1; k=1;

a = 0:0.01:1;

x1 = zeros(1,numel(a));

x2 = zeros(1,numel(a));

for i=1:numel(a)

% assuming you want positive values of x1 and x2

if a(i) == 0

A = k*V;

x1(i) = Ca0;

x2(i) = (-Q + sqrt(Q^2 + 4*A*Q*x1(i)))/(2*A);

elseif a(i) == 1

A = k*V;

x1(i) = (-Q + sqrt(Q^2 + 4*A*Q*Ca0))/(2*A);

x2(i) = x1(i);

else

A1 = k*a(i)*V;

A2 = k*(1-a(i))*V;

x1(i) = (-Q + sqrt(Q^2 + 4*A1*Q*Ca0))/(2*A1);

x2(i) = (-Q + sqrt(Q^2 + 4*A2*Q*x1(i)))/(2*A2);

end

end

subplot(2,1,1)

plot(a,x1),grid

xlabel('a'),ylabel('x1')

subplot(2,1,2)

plot(a,x2),grid

xlabel('a'),ylabel('x2')

% The two equations can be expressed as:

% k*a*V*x1^2 + Q*x1 - Q*Ca0 = 0

% k*(1-a)*V*x2^2 + Q*x2 - Q*x1 = 0

Alan Stevens
on 16 Sep 2021

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!