下面的代码应该模拟牛顿·拉夫森方法来求解方程,但是替换主变量(V)对方程有效,但对其导数无效。有人能帮我处理一下这个错误吗?我找不到subs工作的原因只是上面一行,而不是下面一行-谢谢
close all;
clear all;
clc;
i = 0;
root = 0;
T = 50;
i = i+1;
R = 0.4615;
b = 0.00169;
a = 1.703;
P = 10000;
syms V
F = (R*T)/(V-b)- a/(V^2)- P;
DF = diff(F,V);
oldguess = 10;
newguess = 10;
realerror = 100;
acceptable_error = 0.0001;
while realerror > acceptable_error
f_V = vpa(subs(F,V,oldguess)); %% f_V receives a value
fd_V= vpa(subs(DF,V,oldguess)); %% fd_V doesnt get any value
newguess = oldguess - f_V/fd_V; %% this will be an error since we are dividing by zero
realerror = ((newguess-oldguess)/newguess)*100;
oldguess = newguess;
end
root = newguess;
您遇到的问题不是subs(DF,V,oldguess)
不给出任何值,而是它在迭代数次后给出0。
添加行disp([f_V, fd_V])
到循环中。这是我看到的:
[ -9997.7091399665843527556156990531, -0.22742201327580833604722207521671]
[ -10000.00052501683105455937390418, -0.000000011945511284499679327573322611372]
[ -10000.000000000027564264394642461, -0.00000000000000000000003292691968007601924354847989367]
[ -10000.000000000000000000000000076, -2.5017500564189169087426021010076e-52]
[ -10000.0, -1.4442073343107664457291836426706e-110]
[ -10000.0, -4.8128331074807003817768165960745e-227]
[ -10000.0, -5.3449459015966837719809956687411e-460]
[ -10000.0, -6.5921690739471393457565679030652e-926]
[ -10000.0, -1.0027631932710798302770360761963e-1857]
[ -10000.0, -2.3202697552555340007454567882706e-3721]
[ -10000.0, -1.2422776383481875816235353595144e-7448]
[ -10000.0, -3.5610579836824354028159535601917e-14903]
[ -10000.0, -2.9261716619964961538294369733188e-29812]
[ -10000.0, -1.9757923974050109050486823460432e-59630]
[ -10000.0, -9.0079160415622390528769344757706e-119267]
[ -10000.0, -1.8723643738280768889550544416625e-238539]
[ -10000.0, -8.0895143138882494056514979405605e-477085]
[ -10000.0, -1.5100335803334613297994112468127e-954175]
[ -10000.0, -5.2615647621928018162324162884447e-1908357]
[ -10000.0, -6.3880977095623306259506020122841e-3816720]
[ -10000.0, -9.4163980840507503900230854098031e-7633446]
[ -10000.0, -2.0460268576440353827658449489008e-15266897]
[ -10000.0, -9.6597162693281748513616574438529e-30533801]
[ -10000.0, -2.1531309821705331994608199009543e-61067607]
[ -10000.0, -1.0697507758377953232431300192808e-122135220]
[ -10000.0, -2.6406262119508412798733809908278e-244270447]
[ -10000.0, 0]
Error using symengine
Division by zero.
第一次迭代使您从最初的猜测10到-43951,随后的每次迭代都使您的猜测越来越接近负无穷。在某些情况下,即使在vpa
数据类型中,导数也太小而无法正确表示,并且四舍五入为0。
在尝试查找根之前,您应该始终做的一件事是确定函数是否有根。这个没有,它在所有地方都很接近- 10000除了在接近0的地方两边都趋向于负无穷,在0处函数不能求值因为除以0。您可以通过老式的函数分析看到这一点,如果您不再在学校里学习这个,只需执行fplot(F)
以获得函数行为的概述。