Using Newton's method for 2 equations for finding concentration at different points
1 view (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Torsten
on 7 Nov 2022
Edited: Torsten
on 7 Nov 2022
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
sol
0 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!