strange solution -Newton Raphson

1 次查看（过去 30 天）
Naceur Khraief
Naceur Khraief 2019-1-8
评论： Naceur Khraief 2019-1-9
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

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2019-1-8
编辑：John D'Errico 2019-1-8
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 个评论
Naceur Khraief
Naceur Khraief 2019-1-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Error and Exponential Integral Functions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by