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 件)
参考
カテゴリ
Help Center および File Exchange で Communications Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!