Hello Victor!
I understand that you have implemented Newton’s method to solve for two concentrations, CA1 and CA2, with a dependency on alpha (α).
To plot how CA2 varies as α ranges from 0 to 1, you can wrap your Newton solver in a loop over α and store the final CA2 value from each run. Once the loop finishes, you can simply plot CA2 against α.
I have attached the code for the same below, in which I have wrapped the logic (written by you) in a “for” loop. Also, I have defined the helper functions (fx, J) outside the loop to keep things clean.
function TEST11
clear
alphas = linspace(0,1,100); % 100 alpha values from 0 to 1
CA2_vals = zeros(size(alphas)); % To store CA2 for each alpha
for a = 1:length(alphas)
al = alphas(a); % Current alpha
n=2;
Jac=zeros(n); fv=zeros(n,1); f=zeros(n,1); Cd=zeros(n,1); C=zeros(n,1);
C=[0.5; 0.3];
CA0=1; V=1; Q=50; k1=1; k2=1;
Rel_dx=0.01;
itr_max=100; tol=1e-4;
for k=1:itr_max
Jac=J(C, al, Q, CA0, k1, V, Rel_dx);
fv=fx(C, al, Q, CA0, k1, k2, V);
f=-fv';
d=Jac\f;
C=C+d;
Rerr=d./C;
max_err=max(abs(Rerr));
if max_err<tol; break; end
end
if k==itr_max
fprintf('Newton failed at alpha = %.2f\n', al);
end
CA2_vals(a) = C(2); % Store CA2 for this alpha
end
% Plot CA2 vs alpha after loop
figure
plot(alphas, CA2_vals, 'b-', 'LineWidth', 2);
xlabel('\alpha'); ylabel('CA2');
title('CA2 vs \alpha');
grid on
end
function Jac = J(C, al, Q, CA0, k1, V, Rel_dx)
n = length(C);
fv = fx(C, al, Q, CA0, k1, 1, V); % '1' is a placeholder for k2 which is not used in fx1
Jac = zeros(n);
for j = 1:n
Ctemp = C;
Ctemp(j) = Ctemp(j)*(1+Rel_dx);
fvd = fx(Ctemp, al, Q, CA0, k1, 1, V);
Jac(:, j) = (fvd - fv) / (Ctemp(j) - C(j));
end
end
function fv = fx(C, al, Q, CA0, k1, k2, V)
fv = zeros(1, 2);
fv(1) = Q*CA0 - Q*C(1) - k1*C(1)^2*(al*V);
fv(2) = Q*C(1) - Q*C(2) - k2*C(2)^2*(1 - al)*V;
end
TEST11
