strange solution -Newton Raphson

2 views (last 30 days)
Naceur Khraief
Naceur Khraief on 8 Jan 2019
Commented: Naceur Khraief on 9 Jan 2019
I am tried to resolve this numerical equation, I found a solution out of the intervall (2400, 3002). The solution is "3001398". help please.
xguess = 2400; % initial guess
emax = 1e-7; % Maximum error
imax = 1000; % maximum number of iterations
maxstep = 3002; % maximum x-step size
minslope = .01; % minimum value for the slope (prevents divide by zero)
F1=@(x) 2400 + (1- (sqrt(pi)*erf(x))./exp(-x.^2)/2)./(exp(-(x-2774)/228).^2)/2/sqrt(2*pi)*223;
syms x
dF1 = matlabFunction(diff(F1(x),x));
lim = @(x,xlim) max(-xlim,min(xlim,x)); % limiter function
i=1;
x=xguess;
while (i<imax) && (abs(10))>=emax % needs abs(deltax) (can be + or -)
yguess = F1(x); % Get the function value for the x guess
slope = dF1(x); % Get the function slope for x-guess
if(abs(slope)<minslope) % Don't allow slope to go to zero (trap divide by zero condition)
slope = minslope*sign(slope);
end
yerr = yguess; % Since the desired value is zero, the function value represents the error value
deltax = -yerr/slope; % calculate the x-step
xstep = lim(deltax,maxstep); % limit the x-step (deltax) to +/- maxstep
x = x + xstep; % apply the (limited) x-step to the x-guess value and repeat the iteration
i=i+1; % update the increment counter
end

Answers (1)

John D'Errico
John D'Errico on 8 Jan 2019
Edited: John D'Errico on 8 Jan 2019
It always makes sense to test your function. Does it make sense?
Your starting value is 2400.
F1=@(x) 2400 + (1- (sqrt(pi)*erf(x))./exp(-x.^2)/2)./(exp(-(x-2774)/228).^2)/2/sqrt(2*pi)*223;
F1(2400)
ans =
-Inf
What are the odds you will get anything meaningful out of this iteration? (Zero.)
So, does your function make any sense at all in terms of numerical methods? In terms of floating point arithmetic using doubles? Not really. At leasty not if you think the function will do somethign meaningful for x on the order of 2400.
For example, I see terms like erf(x), exp(-x^2/2), etc. really? Do you have a clue what the value of erf(2400) is? Can you distinguish that number from 1? Even the symbolic toolbox cannot do so. So then I tried computing erf(2400) using my own HPF toolbox, in 100000 decimal digits of precision. It still returned exactly 1.
erf(hpf(2400,100000)) - 1
ans =
0
Or, how about exp(-x^2/2) when x is on the order of 2400?
x = hpf(2400,25);
exp(-x^2/2)
ans =
7.800431631420969539097625e-1250769
Do I need to go further? Spend some time and think about if your function makes numerical sense.
  1 Comment
Naceur Khraief
Naceur Khraief on 9 Jan 2019
The original form of this function is as follow:
F1(x)= 2400 + (1-F(x))/F'(x)
where F(x) is the distribution function of truncated normal in the rang [2292, 3332] defined as follow:
F(x)= (sqrt(pi)*erf(x)/ exp(-x^2)/2)
F'(x) is the truncated normal probability density function (derivative of F(x))
the mean is mu=2774.5 and standard deviation 228.83.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!