Using Newton's method for 2 equations for finding concentration at different points

Hi all, I'm currently trying to use Newton's method on two equations (I've attatched them) for finding the concentration of the lung and liver at three different R values 0.6, 0.8, and 0.9 from their intial guesses R = 0.6: Clung = 360.6 and Cliver = 284.6, R = 0.8: Clung = 312.25 and Cliver = 236.25, and R = 0.9: Clung = 215.5 and Cliver = 139.5. What I have so far is shown below, I'm hitting a wall and it wont allow me to create my J without giving an error if I use .*, etc. If I don't use .* in the J it gives me an error on dimensions instead. Any help would be greatly appreciated. Thank you

clearvars

clc

close all

%% Newton's Method to Systems

%% I have 2 equations, 3 sets of initials

maxiter = 100;

tol = 1e-6;

k = 0;

err(1) = 1;

%xold = [0.5; 0.25];

R_Old = [0.6; 0.8; 0.9];

C_Lung_Old = [360.6; 312.25; 215.5];

C_Liver_Old = [284.6; 236.25; 139.5];

Mat_Old = [R_Old,C_Lung_Old,C_Liver_Old];

while k<maxiter && err(k+1)>tol

k = k+1;

f1_1 = -780.*(1-R_Old(1))-R_Old(1).*(0.5.*C_Liver_Old(1)+1.5.*C_Lung_Old(1)) + (8.75.*C_Lung_Old(1)./(2.1 + C_Lung_Old(1))) + 2.*C_Lung_Old(1) ;

f1_2 = -780.*(1-R_Old(2))-R_Old(2).*(0.5.*C_Liver_Old(2)+1.5.*C_Lung_Old(2)) + (8.75.*C_Lung_Old(2)./(2.1 + C_Lung_Old(2))) + 2.*C_Lung_Old(2) ;

f1_3 = -780.*(1-R_Old(3))-R_Old(3).*(0.5.*C_Liver_Old(3)+1.5.*C_Lung_Old(3)) + (8.75.*C_Lung_Old(3)./(2.1 + C_Lung_Old(3))) + 2.*C_Lung_Old(3) ;

f2_1 = -0.5.*C_Lung_Old(1) + 0.322*(118*C_Liver_Old(1)./(7.0 +C_Liver_Old(1))) + 0.5*C_Liver_Old(1) ;

f2_2 = -0.5.*C_Lung_Old(2) + 0.322*(118*C_Liver_Old(2)./(7.0 +C_Liver_Old(2))) + 0.5*C_Liver_Old(2) ;

f2_3 = -0.5.*C_Lung_Old(3) + 0.322*(118*C_Liver_Old(3)./(7.0 +C_Liver_Old(3))) + 0.5*C_Liver_Old(3) ;

nF = [f1_1,f1_2,f1_3;f2_1,f2_2,f2_3];

J = [(C_Liver_Old+3*C_Lung_Old-1560)/2 (20*(3*R_Old-4)*C_Lung_Old*(5*C_Lung_Old+21)+1323*R_Old-5439)/(2*(10*C_Lung_Old+21).^2) 2*R_Old R_Old/2; ...

0 1/2 -((125*C_Liver_Old^2+1750*C_Liver_Old+72618)/(250*(C_Liver_Old+7)^2))];

xx = J\nF;

xnew = xx+Mat_Old;

f1 = -3*xnew(1).^2-xnew(2).^2 +10;

f2 = -xnew(1).^2-xnew(2)+1;

errtemp = [abs(f1) abs(f2)];

err(k+1) = max(errtemp);

Mat_Old = xnew;

end

semilogy(err)

return

Answer 1

[Torsten]

Why don't you use "fsolve" ? Don't you have a license for the optimization toolbox ?

I don't understand J, f1 and f2 in your code.

maxiter = 100;
tol = 1e-6;
R_Old = [0.6; 0.8; 0.9];
C_Lung_Old = [360.6; 312.25; 215.5];
C_Liver_Old = [284.6; 236.25; 139.5];
sol = [];
for i = 1:3
    f = @(x) [-780.*(1-R_Old(i))-R_Old(i)*(0.5*x(2)+1.5*x(1)) + (8.75*x(1)/(2.1 + x(1))) + 2*x(1) ; ...
              -0.5*x(1) + 0.322*(118*x(2)/(7.0 +x(2))) + 0.5*x(2)] ;
    Mat_Old = [C_Lung_Old(i);C_Liver_Old(i)];
    err(1) = 1;
    k = 0;
    while k<maxiter && err(k+1)>tol   
        k = k+1;
        fv = f(Mat_Old);
        df = [(f(Mat_Old+1e-6*[1;0])-fv)/1e-6,(f(Mat_Old+1e-6*[0;1])-fv)/1e-6];
        xx = df\fv;  
        xnew = Mat_Old-xx;   
        err(k+1) = norm(abs(fv));
        Mat_Old = xnew;
    end
    sol = [sol,Mat_Old];
end

k = 3

k = 3

k = 4

sol

sol = 2×3

351.3324 294.6213 185.6450 277.2120 220.9628 114.0475