Using Newton's method for 2 equations for finding concentration at different points
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Torsten
2022-11-7
编辑:Torsten
2022-11-7
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 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!