subs doesnt substitue correctly in MATLAB-CodePudding

The code below is supposed to simulate the Newton Raphson method to solve an equation, but the substitution of the main variable (V) works for the equation, but not for its derivative. Can someone please help me deal with this error? I can't find the reason subs works just one line above but not the following line - Thanks

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;

CodePudding user response：

The problem you have is not that subs(DF,V,oldguess) doesn't give any value, it is that it gives 0 after iterating for a number of times.

Add the line disp([f_V, fd_V]) to the loop. This is what I see:

[ -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.

The very first iteration takes you from your initial guess of 10 to -43951, and each subsequent iteration moves your guess further and further towards minus infinity. At some point you are so far out that the derivative is too small to properly represent even in the vpa data type, and is rounded to 0.

One thing you should always do before attempting to find a root is to determine if a function has roots. This one doesn't, it is pretty close to -10,000 everywhere except close to 0 where it goes towards minus infinity on both sides, and at 0 where the function cannot be evaluated because you divide by 0. You can see this by old-fashioned function analysis, and if you don't learn this in school any more, just do fplot(F) to get an overview of the function's behavior.