Plot how would I take my out put C(2) and plot it vs Al (Alpha) This newtons method works to get Ca2 and CA1, but I want to bang CA2 vs Alpha for alpha 0-1. Do I use a matrix?

function TEST11

clear

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; al=0.11;

Rel_dx=0.01;

itr_max=100; tol=10^-4;

for k=1:itr_max

Jac=J(C);

fv=fx(C);

f=-1.*fv';

d=Jac\f;

C=C+d;

Rerr=d./C;

max_err=max(abs(Rerr));

fprintf('CA1=%f & CA2=%f max error=%f\n',C,max_err)

if max_err<tol; break; end;

end

if k==itr_max

fprintf('Newtons method failed to converge in %d iternations\n', itr_max)

end

if i<itr_max

fprintf('The solution is CA1=%f CA2=%f in %d iterations\n', C,k)

end

function Jac=J(C)

fv=fx(C);

for j=1:n

C(j)=C(j)*(1+Rel_dx);

fvd=fx(C);

for i=1:n

Jac(i,j)=(fvd(i)-fv(i))/C(j)/Rel_dx;

end

C(j)=C(j)/(1+Rel_dx);

end

end

function fv=fx(C)

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

end

Answer 1

[Aashray]

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

●