Newton Raphson implementation, not converging, maybe error in implementation
Afficher commentaires plus anciens
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
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
4 commentaires
VIGNESH BALAJI
le 7 Fév 2024
Dyuman Joshi
le 7 Fév 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
le 7 Fév 2024
Modifié(e) : John D'Errico
le 7 Fév 2024
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)
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])
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.
Réponse acceptée
John D'Errico
le 7 Fév 2024
Modifié(e) : John D'Errico
le 7 Fév 2024
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)
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)
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 commentaires
VIGNESH BALAJI
le 7 Fév 2024
Torsten
le 7 Fév 2024
If I am correct fsolve (True region method), will also suffer from this incorrect search region problem just like fzero and newton-raphson. Am I correct ?
Yes.
I have explained my problem in detail. I would really like your suggestion to solve this problem for all possible range of configurations for my nonlinear function.
Find two points xl and xr for a with f(xl)*f(xr) < 0 and call fzero with the interval [xl xr] as initial guess. To find these points, you could successively scan the range of values for "a" that are physically senseful.
VIGNESH BALAJI
le 7 Fév 2024
dd = 0.1069;
vv = 0.7889;
L = 1;
func = @(a) ( (-(sqrt(L^2-vv^2)/2)) + (a.*sinh((dd/2)./a)) );
amin = 1e-5; % minimum lower limit for a
amax = 1; % maximum upper limit for a
deltaa = 1e-3; % Scan size for the interval [amin amax]
a = amin:deltaa:amax;
fa = func(a);
index = find(diff(fa>0),1); % find first sign change position
a = fzero(func,[a(index) a(index+1)])
John D'Errico
le 7 Fév 2024
Modifié(e) : John D'Errico
le 7 Fév 2024
fzero is a more robust solver (as long as there is only one variable) than fsolve. With more than 1 variable, you have no choice as fzero would not apply. If you have a bracket around the solution where there is a change of sign, then fzero will always generate a solution.
If you have no bracket around the root, then fzero uses a fairly simple scheme to try to find a bracket (but it is not terribly intelligent in this respect.) Once it has one, then it can just drop into the default algorithm.
Is there some magic way you can be assured of finding a bracket? NO. I can always generate a function where a bracket can never be found, even though technically, one exists. There is no mathematical trick that will insure a bracket will be found. And one thing I observe with your function, if you go too far out, then you will see numerical problems.
Can you just sample the function at some set of points, and look to see if there is a pair that form a brcket? Yes. Of course. This is a simple scheme. As long as your function is reasonably smooth, you do not even need the points to be terribly close together. And if it is fairly smooth, very often the simple scheme fzero uses internally is sufficient to find the bracket on its own.
You said that sometimes, you had issues with fzero, on some problems. Is there any way to ALWAYS insure a solution will be found for ANY problem? OF COURSE NOT! For some problems, numerical issues may interfere. You may be experiencing underflows or overflows. You may be able to scale things to avoid that, but we have not been shown which which cases had a problem. Possibly, for some choices of parameters, no solution exists at all.
VIGNESH BALAJI
le 7 Fév 2024
VIGNESH BALAJI
le 7 Fév 2024
@John D'Errico Okay, I get it you mean fzero (Bisection method) is better than fsolve (True Region method).
The main way to solve my problem is to find a right Bracket. Can you elaborate your statement if I go too far in my function, then I will see numerical problems.
I just got an idea, my variation in the 2 constant parameters of the function are based on real world scenarios (have a length constraint, thereby they have a minimum and a maximum range). In that case, I can plot this function with the minimum (dd = 0.01, vv=0.98) and maximum range(dd=0.98, vv = 0) and also I can get some points in the middle. This can help me to plot and see how the function varies, a visual representation is always good to understand it. I would really like to see the existence of valleys (local minimas) and flat points (I didn't understand the concept of flat points earlier, can you please explain it again).
I made a small code for this visualisation
close all
clear all
clc
x1 = 0;
y1 = 1.0;
x2 = 0.01;
y2 = 1.98;
%% plot sinh
dd = abs(x1-x2);
vv = abs(y1-y2);
L = 1.0;
func = @(a) ( (-(sqrt(L^2-vv^2)/2)) + (a.*sinh((dd/2)./a)) );
figure();
[X,Y] = fplot(func);
idx = find( diff( sign( Y(1,:) ) ) ) + 1;
plot(X,Y);
yline(0, '--','Color', 'r');
ylim([-0.5 1.5]) % ADDED
grid
When I looked at the plots for some of the range of values for dd (horizontal distance) and vv (vertical distance). I was able to understand the funcction better and I was unable to see the valley or flat point in these plots. Please do help me spot them.
- Normal case dd = 0.5, vv = 0
2. Maximum range dd = 0.98, vv = 0
3. This is the plot where the problem started with incorrect search region dd = 0.1069 , vv = 0.788932 .
In this plot, I earlier had a start search region of 0.2, I guess this is in the flat line and when it was made to 0.01 fzero was able to find the solution as it can see the direction of the curve in that region. I guess the problem with flat line is it does not give the direction of the curve, Am I correct ?
4. This dataset d = 0.01, v = 0.98 has a problem in zero crossing as it does not cross zero. I am not sure whether any method can find this corner case
In this case, fzero did not find a solution whereas fsolve found a solution. I guess the reason could be fzero looks for exact zero crossing and fsolve is looking for a region close to zero. Based on this I am not sure why you still call fzero as robust when compared to fsolve ? Please let me know your suggestion on this case ?
This is the code i used to arrive at a solution for this corner case
dd = 0.01; %abs(x1-x2);
vv = 0.98; %abs(y1-y2);
L = 1.0;
x0 = 0.1;
func = @(a) ( (-(sqrt(L^2-vv^2)/2)) + (a.*sinh((dd/2)./a)) );
a = abs(fsolve(func, x0));
Answer a = 9.322285907071188e-04
Thanks a lot for your answers, I believe I am very close to having a complete understanding. Please do answer me the questions I have asked above about flat points, fsolve and my approach to visualise and select a good enough bracket (as it is not always possible to get a right bracket).
I am looking forward for your answers :) :)
First, that you can always choose any start point. WRONG.
In general, that is true, but I believe that Newton-Raphson is globally convergent for strictly convex functions, which in this case it is (assuming, aa>0, dd>0). You would have to threshold the iteration updates to some lower bound lb where fun(lb)>0,
aa_new=max(lb, aa_guess - f / f_prime;)
but this should converge from any starting aa>0.
It might be a problem to get a suitable lb or to formulate more complicated updates for aa:
% 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.1;
%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
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
VIGNESH BALAJI
le 9 Fév 2024
VIGNESH BALAJI
le 9 Fév 2024
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
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
Plus de réponses (1)
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
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
Catégories
En savoir plus sur Solver-Based Nonlinear Optimization dans Centre d'aide et File Exchange
le 7 Fév 2024
le 9 Fév 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
