how accurate is the function fsolve?

I have a nonlinear equation to solve and I know (from the textbook) that the exact solution is zero but when I use this method I get 0.9891. Here's how I did it:

function F = myFunction(z) 
x = z(1); 
% Probability of False Alarm [Pfa] 
Pfa=0.5; 
F(1) = exp(-x^2/2)/(x^2/2) - sqrt(2*pi)*Pfa; 
end

Then on command window: >> zg=[1]; z=fsolve(@myFunction,zg) gives: z=0.9891; However if I change Pfa to 0.001, the exact solution is 3 and it gives a much better solution of 2.9939!

Answer 1

[John D'Errico]

I'll be lazy, and write the function as a function handle.

Pfa=0.5;
F = @(x) exp(-x.^2/2)./(x.^2/2) - sqrt(2*pi)*Pfa;
ezplot(F,[.5,1.5])
grid on

fsolve(F,1)
ans =
       0.989138491394299

To me, that solution seems pretty accurate compared to the plot. Lets try a symbolic solution.

vpasolve(exp(-x.^2/2)./(x.^2/2) - sqrt(2*pi)*Pfa == 0,1)
ans =
0.98913855820552302829982713582152

Fsolve uses a tolerance. (read the help docs for fsolve) The default tolerance is 1.e-6, which is roughly where I see some deviation.

Anyway, you claim that when Pfa == 0.5, the exact solution is zero. Of course that is completely wrong. At zero that function has a singularity, essentially 1/0 in the limit, so it must produce a +inf as a result as x approaches 0.

By the way, for PFA = 0.001, the exact solution is not 3.

Pfa = 0.001;
vpasolve(exp(-x.^2/2)./(x.^2/2) - sqrt(2*pi)*Pfa == 0,1)
ans =
2.9958361659324325479490777418211

Anyway, what seems to be the problem? Besides the error in either your textbook, or what you think you understood from that textbook. Or, maybe the function that you wrote is not what was actually in that textbook. I cannot know for sure. But fsolve seems to have worked reasonably well.

[Diogo]

Hi John and thanks for the reply. Well, the problem is that according to my textbook, when Pfa=0.5 the solution is 0, not 0.9891, and this is a huge difference!

Here's how it is on the textbook: For example, if Pfa=0.5, then the threshold is found from:

    0.5 = integral('(1/(sqrt(2*pi)))*exp(-0.5*t.^2)',x,Inf)

I reached the function "F(1) = exp(-x^2/2)/(x^2/2) - sqrt(2*pi)*Pfa;" by developing this integral manually and then input the function into the fsolve. In the above expression, Pfa is 0.5.

Here, x is the unknown variable I'm trying to find. Then it says, "as x=0, Then... " so I just tried to reach the same solution using fsolve but like I said before, didn't quite work out! Unless the textbook is wrong...

This picture shows better what I'm trying to say here:

[Steven Lord]

Rephrasing what John said, you're dividing by x^2/2 when you compute F(1) so F(1) is going to be Inf at x = 0.

x = 0;
firstPartOfF = exp(-x^2/2)/(x^2/2)

The second part of F(1) is a finite constant, so subtracting it from firstPartOfF won't make a difference in this case. Therefore one of the following must be true:

1. The textbook is incorrect.
2. The problem you posted in Answers is not the same as the problem as written in your textbook, and 0 is a solution for the problem in the textbook but not the problem in the Answers post.
3. The equation that you generated to solve the problem (and which you're passing into fsolve) is not a correct solution to the problem in the textbook.

Based on the snippet you copied from your textbook, I think you don't need to use fsolve. Take a look at the erf, erfc, erfinv, and erfcinv functions.

[Star Strider]

The function in your code is incorrect. You cannot integrate the error function analytically. It is what it is (although a series expression for it exists). You have to do the integration:

To illlustrate:

fcn = @(gma,cdprb) integral(@(x) exp(-x.^2/2), gma, Inf)/sqrt(2*pi) - cdprb;

cdprob = [0.5, 0.001];
for k1 = 1:length(cdprob)
    gma(k1) = fsolve(@(gma)fcn(gma,cdprob(k1)), 1);
end
fprintf('p = %.3f,  gamma = %.3f\n', [cdprob', gma']')

p = 0.500,  gamma = 0.000
p = 0.001,  gamma = 3.078

[Diogo]

That's a great solution! I really messed up developing the integration I should have seen that. Thanks for writing and I'm sorry about the other post, it's my first time using mathworks website.

[John D'Errico]

What you have written is simply not the analytical integral of the normal density function. As others have pointed, there is no simple expression for that integral, but there are ways to solve the problem in MATLAB. I'd suggest using the stats toolbox functions normcdf or norminv. If you don't have that TB, then there are still many ways to solve the problem.

[Star Strider]

Diogo — Thank you!

No apologies necessary. It would have been helpful to us if you had posted the image earlier. We would then know what you actually want to do.

The function you want (and that does this integration for you) is the ‘complementary error function’, erfc. The erfc function calculates it for you. It is a core MATLAB function, so no Toolboxes are necessary. I used the anonymous function and integral to demonstrate how to use fsolve to calculate it, since that appears to be the essence of your question.

If you have the opportunity, take a course in partial differential equations. That’s where I first encountered the maths behind the error function, and that it cannot be analytically integrated.