Newton Raphson implementation, not converging, maybe error in implementation

27 views (last 30 days)
VIGNESH BALAJI
VIGNESH BALAJI on 7 Feb 2024
Edited: Matt J on 9 Feb 2024
Ran in:
I have implemented Newton-Raphson method for my non-linear function with one variable "aa" and others are constant.
I read that Newton-Raphson method does not depend much on the initial guess. My implementation is not converging and it is not getting me a solution. I know the solution for this problem based on other method as 0.0142. I would like to know the mistakes in implementation, please take a look at it and let me know your suggestions.
% Predefined parameters
dd = 0.1069;
vv = 0.7889;
L = 1;
% Define the function and its derivative
func = @(aa) ( (-(sqrt(L^2-vv^2)/2)) + (aa*sinh((dd/2)/aa)) );
func_derivative = @(aa) (sinh(dd/(2*aa)) - ((-dd/(2*aa))*(cosh( (dd/(2*aa))/aa ))) );
% Initial guess
aa_guess = 0.01; % You can choose any initial guess
% Tolerance and maximum iterations
tol = 1e-6;
max_iter = 1000;
% Newton-Raphson method
iter = 0;
while true
iter = iter + 1;
% Evaluate the function and its derivative at the current guess
f = func(aa_guess);
%f_prime = func_derivative(aa_guess);
f_prime = (func(aa_guess+1e-6)-f)*1e6
% Newton-Raphson update rule
aa_new = aa_guess - f / f_prime;
% Display iteration results
fprintf('Iteration %d: aa = %f, f(aa) = %f, f''(aa) = %f\n', iter, aa_guess, f, f_prime);
% Check for convergence
if abs(aa_new - aa_guess) < tol || iter >= max_iter
break;
end
% Update the guess
aa_guess = aa_new;
end
f_prime = -455.1300
Iteration 1: aa = 0.010000, f(aa) = 0.740505, f'(aa) = -455.130017
f_prime = -178.3815
Iteration 2: aa = 0.011627, f(aa) = 0.269331, f'(aa) = -178.381458
f_prime = -89.7550
Iteration 3: aa = 0.013137, f(aa) = 0.076755, f'(aa) = -89.754985
f_prime = -64.3456
Iteration 4: aa = 0.013992, f(aa) = 0.011643, f'(aa) = -64.345556
f_prime = -60.2285
Iteration 5: aa = 0.014173, f(aa) = 0.000376, f'(aa) = -60.228463
f_prime = -60.0928
Iteration 6: aa = 0.014179, f(aa) = 0.000000, f'(aa) = -60.092768
if iter >= max_iter
disp('Newton-Raphson method did not converge within maximum iterations.');
else
fprintf('Root found at aa = %f after %d iterations.\n', aa_guess, iter);
end
Root found at aa = 0.014179 after 6 iterations.
  4 Comments
Show 2 older comments
Dyuman Joshi
Dyuman Joshi on 7 Feb 2024
"I am unable to understand why you multiplied it with 10^-6 after taking the difference ?"
That is not multiplication with 10^-6, but division by 10^-6, which results in multiplication by 10^6. Note the sign of the exponent.
Basically that is f' = (f(x+del) - f(x))/del
John D'Errico
John D'Errico on 7 Feb 2024
Edited: John D'Errico on 7 Feb 2024
Ran in:
@VIGNESH BALAJI
Are there methods that depend less on the initial guess? YES, and sometimes, no. And this is big reason why you don't want to write your own solvers, but also why you need to learn about the solvers, and how to use them properly.
For example, fzero is a better choice of solver almost always than writing a Newton method of your own. It will converge rapidly when the problem is well behaved. But it also has ways of avoiding problems.
dd = 0.1069;
vv = 0.7889;
L = 1;
% Define the function and its derivative
func = @(aa) ( (-(sqrt(L^2-vv^2)/2)) + (aa*sinh((dd/2)/aa)) );
[xval,fval,exitflag] = fzero(func,0.1)
xval = 0.0142
fval = -5.5511e-17
exitflag = 1
Even better is if you will provide a bracket around the root, as that will insure it finds a proper solution.
[xval,fval,exitflag] = fzero(func,[1e-2,0.1])
xval = 0.0142
fval = -5.5511e-17
exitflag = 1
USE EXISTING TOOLS WHENEVER POSSIBLE. Unless of course, you can do a better job than the professional who wrote that tool.
Finally, you should note that had you used fzero here, it would have saved you perhaps hours of time writing code, and trying then to figure out where you went wrong. And it took me one line of code to solve the problem. I never even needed to compute the derivative, which you apparently got wrong.

Sign in to comment.

Sign in to answer this question.

Accepted Answer

John D'Errico
John D'Errico on 7 Feb 2024
Edited: John D'Errico on 7 Feb 2024
Ran in:
There are some significant misunderstandings on your part that I want to clear up, as well as perhaps help you solve your problem.
First, that you can always choose any start point. WRONG.
It is not true that a Newton scheme can start anywhere. Some start points may cause your method to diverge to +/- inf. Some start points may cause your scheme to converge to a non-solution. Some start points may cause numerical problems in the evaluation. So it is not at all true that any point is acceptable!
For example, consider the simple cubic polynomial...
syms x
y = expand((x+1)*(x-1)*(x-2) + 1)
y = 
fplot(y,[-2,3])
yline(0)
Clearly, it has only one real root, near -1.5. If you start the solver in the region of that root, it will converge to the root we want.
But imagine you decide to start the solver out around x==2? Now it will probably get stuck in that valley. A Newton scheme may just oscillate in that valley around 1.5, and never escape.
And suppose we decide to start the solver at one of two flat spots on the curve, where the slope is exactly zero? It will diverge to plus or minus infinity at the first step, due to the resulting divide by zero.
The concept of a basin of attraction is an important one. It is the set of starting points that will converge to a valid solution. As long as you start the solver in a spot that lies in the basin of attraction, AND that does not create a numerical problem, like an underflow or overflow in double precision arithmetic, then yes, you can succeed. But you stated several times that it does not matter where you start. And that is simply wrong.
Now we can look at your specific problem.
dd = 0.1069;
vv = 0.7889;
L = 1;
% Define the function and its derivative
% Use the dotted operators to allow the code to be vectorized
func = @(aa) ( (-(sqrt(L^2-vv^2)/2)) + (aa.*sinh((dd/2)./aa)) );
% PLOT THE FUNCTION!!!!!!!!!!!! ALWAYS PLOT EVERYTHING!!!!!!!!!
fplot(func)
It looks like your function has a spike in it around zero. Zoom in.
fplot(func,[-.1,.1])
ylim([-.2,.2])
yline(0)
grid on
So 0.01 might not be a terrible guess as a start point. A little low, but not bad.
Next, check if your derivative is correct. I can drop the first term, since it is not a function of aa.
syms DD AA
diff(AA*sinh((DD/2)/AA),AA)
ans = 
And you wrote this:
(sinh(dd/(2*aa)) - ((-dd/(2*aa))*(cosh( (dd/(2*aa))/aa ))) );
which appears to be wrong. So fix your derivative. It looks like you divided by aa too many times.
And use fewer parens. Using more parens is not always a good thing, if it makes your code unreadable, to the extent that you cannot read it yourself.
  12 Comments
Show 10 older comments
VIGNESH BALAJI
VIGNESH BALAJI on 9 Feb 2024
@John D'Errico and @Torsten Both of your suggestions and code snippets helped me to solve my actual problem of finding the roots of the equation when the parameters vary. I final settled with by understanding the need for finding a good search region and toggling between fzero (if I have a zero crossing, sign change) and if not I started using fsolve with a decent initial guess. I will accept this answer and close the thread. Thanks a lot :)
Matt J
Matt J on 9 Feb 2024
Ran in:
It might be a problem to get a suitable lb or to formulate more complicated updates for aa:
Well, bisection faces the same problem. The algorithm has to find an interval where a sign change occurs at the end points. In this case, we know the function goes to Inf at aa=0, so it should be sufficient to try lb od descending magnitude, lb=0.1, 0.01,0.001,... Until a qualified lower bound is found. As shown below, it didn't take long to find one.
% Predefined parameters
dd = 0.1069;
vv = 0.7889;
L = 1;
% Define the function and its derivative
func = @(aa) ( (-(sqrt(L^2-vv^2)/2)) + (aa.*sinh((dd/2)./aa)) );
func_derivative = @(aa) sinh(dd./(2*aa)) - dd*cosh( dd./(2*aa))./(2*aa);
% Initial guess
aa_guess = 1; % You can choose any initial guess
lb = 0.01;
%lb = 0.01;
% Tolerance and maximum iterations
tol = 1e-6;
max_iter = 1000;
% Newton-Raphson method
iter = 0;
while true
iter = iter + 1;
% Evaluate the function and its derivative at the current guess
f = func(aa_guess);
f_prime = func_derivative(aa_guess);
% Newton-Raphson update rule
aa_new = max(lb,aa_guess - f / f_prime);
% Display iteration results
fprintf('Iteration %d: aa = %f, f(aa) = %f, f''(aa) = %f\n', iter, aa_guess, f, f_prime);
% Check for convergence
if (abs(aa_new - aa_guess) < tol && abs(func(aa_new)) < tol) || iter >= max_iter
break;
end
% Update the guess
aa_guess = aa_new;
end
Iteration 1: aa = 1.000000, f(aa) = -0.253785, f'(aa) = -0.000051 Iteration 2: aa = 0.010000, f(aa) = 0.740505, f'(aa) = -455.279643 Iteration 3: aa = 0.011626, f(aa) = 0.269426, f'(aa) = -178.474720 Iteration 4: aa = 0.013136, f(aa) = 0.076826, f'(aa) = -89.802472 Iteration 5: aa = 0.013992, f(aa) = 0.011672, f'(aa) = -64.368193 Iteration 6: aa = 0.014173, f(aa) = 0.000380, f'(aa) = -60.240739 Iteration 7: aa = 0.014179, f(aa) = 0.000000, f'(aa) = -60.103653
if iter >= max_iter
disp('Newton-Raphson method did not converge within maximum iterations.');
else
fprintf('Root found at aa = %f after %d iterations.\n', aa_guess, iter);
end
Root found at aa = 0.014179 after 7 iterations.

Sign in to comment.

More Answers (1)

Matt J
Matt J on 9 Feb 2024
Edited: Matt J on 9 Feb 2024
Ran in:
A simple transformation of variables z=dd/2/aa makes this problem, and its solution by NR, much easier. In particular, you can then just chose 0 as the lower bound, lb, and any non-negative initial guess z>=0 should work..
% Predefined parameters
dd = 0.1069;
vv = 0.7889;
L = 1;
const= -sqrt(L^2-vv^2)/2;
% Define the function and its derivative
func = @(z) const*z + dd.*sinh(z)/2;
func_derivative = @(z) 1+dd*cosh(z)/2;
% Initial guess
z_guess = 1; % You can choose any initial guess
lb = 0;
%lb = 0.01;
% Tolerance and maximum iterations
tol = 1e-6;
max_iter = 1000;
% Newton-Raphson method
iter = 0;
while true
iter = iter + 1;
% Evaluate the function and its derivative at the current guess
f = func(z_guess);
f_prime = func_derivative(z_guess);
% Newton-Raphson update rule
z_new = max(lb,z_guess - f / f_prime);
% Display iteration results
fprintf('Iteration %d: aa = %f, f(aa) = %f, f''(aa) = %f\n', iter, z_guess, f, f_prime);
% Check for convergence
if (abs(z_new - z_guess) < tol && abs(func(z_new)) < tol) || iter >= max_iter
break;
end
% Update the guess
z_guess = z_new;
end
Iteration 1: aa = 1.000000, f(aa) = -0.244446, f'(aa) = 1.082478 Iteration 2: aa = 1.225821, f(aa) = -0.293440, f'(aa) = 1.098895 Iteration 3: aa = 1.492853, f(aa) = -0.345781, f'(aa) = 1.124926 Iteration 4: aa = 1.800234, f(aa) = -0.395843, f'(aa) = 1.166131 Iteration 5: aa = 2.139684, f(aa) = -0.433511, f'(aa) = 1.230221 Iteration 6: aa = 2.492068, f(aa) = -0.444921, f'(aa) = 1.325216 Iteration 7: aa = 2.827803, f(aa) = -0.418580, f'(aa) = 1.453454 Iteration 8: aa = 3.115793, f(aa) = -0.355863, f'(aa) = 1.603869 Iteration 9: aa = 3.337670, f(aa) = -0.274083, f'(aa) = 1.753351 Iteration 10: aa = 3.493990, f(aa) = -0.194671, f'(aa) = 1.880519 Iteration 11: aa = 3.597510, f(aa) = -0.130451, f'(aa) = 1.976387 Iteration 12: aa = 3.663514, f(aa) = -0.084114, f'(aa) = 2.042911 Iteration 13: aa = 3.704688, f(aa) = -0.052930, f'(aa) = 2.086691 Iteration 14: aa = 3.730053, f(aa) = -0.032807, f'(aa) = 2.114575 Iteration 15: aa = 3.745568, f(aa) = -0.020147, f'(aa) = 2.131982 Iteration 16: aa = 3.755018, f(aa) = -0.012303, f'(aa) = 2.142718 Iteration 17: aa = 3.760759, f(aa) = -0.007487, f'(aa) = 2.149290 Iteration 18: aa = 3.764243, f(aa) = -0.004547, f'(aa) = 2.153297 Iteration 19: aa = 3.766354, f(aa) = -0.002758, f'(aa) = 2.155732 Iteration 20: aa = 3.767634, f(aa) = -0.001671, f'(aa) = 2.157210 Iteration 21: aa = 3.768409, f(aa) = -0.001013, f'(aa) = 2.158106 Iteration 22: aa = 3.768878, f(aa) = -0.000613, f'(aa) = 2.158649 Iteration 23: aa = 3.769162, f(aa) = -0.000371, f'(aa) = 2.158978 Iteration 24: aa = 3.769334, f(aa) = -0.000225, f'(aa) = 2.159177 Iteration 25: aa = 3.769438, f(aa) = -0.000136, f'(aa) = 2.159297 Iteration 26: aa = 3.769501, f(aa) = -0.000082, f'(aa) = 2.159370 Iteration 27: aa = 3.769539, f(aa) = -0.000050, f'(aa) = 2.159414 Iteration 28: aa = 3.769562, f(aa) = -0.000030, f'(aa) = 2.159441 Iteration 29: aa = 3.769576, f(aa) = -0.000018, f'(aa) = 2.159457 Iteration 30: aa = 3.769585, f(aa) = -0.000011, f'(aa) = 2.159467 Iteration 31: aa = 3.769590, f(aa) = -0.000007, f'(aa) = 2.159473 Iteration 32: aa = 3.769593, f(aa) = -0.000004, f'(aa) = 2.159477 Iteration 33: aa = 3.769595, f(aa) = -0.000002, f'(aa) = 2.159479 Iteration 34: aa = 3.769596, f(aa) = -0.000001, f'(aa) = 2.159480
aa= dd/2/z_guess; %undo change of variables
if iter >= max_iter
disp('Newton-Raphson method did not converge within maximum iterations.');
else
fprintf('Root found at aa = %f after %d iterations.\n', aa, iter);
end
Root found at aa = 0.014179 after 34 iterations.

Categories

Find more on Sparse Matrices in Help Center and File Exchange

Tags

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!