Function floating point issue

I made a function intd(). Given all the arguments except x, I first want to plot the function intd() in terms of x. As is shown below, I use a function handle to turn intd() into univariate h().

Suppose (r,s,q,nu,al)=(5, 2, 0.2104, 60, 0.05).Then h() gives me the below error message when I plugged in 6.699999999999999 for x. By comparisons, when I let x be 6.7, h() returns a perfectly reasonable answer. So it looks like my function has some floating point issue. This issue particularly gets in the way when I need to find a root to function h(). I used fzero(h, initial value) to find the root and got the exact same error message.

How could I make my function insensitive to the floating issue? If I had to keep my function, would there be a way to modify fzero so that it could avoid the floating point issue?

function prob=intd(x,r,s,q,nu,al)

a=q.*(1+q.^2).^(-1./2);
b=(1+q.^2).^(-1./2);
E=@(t) fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)- 
a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) ;
up=x.^2./(b.^2)./r;
down=x.^2./r;
prob=fcdf(x.^2./r,r,nu)+quadv(E,down,up)-1+al;
end
h=@(x) intd(x,r,s,q,nu,al);

h(6.699999999999999)
h(6.7)

%________________ Error Message_______________

Error using betainc Inputs must be real, full, and double or single.

Error in fcdf (line 58) p(kk) = betainc(xx, v1(kk)/2, v2(kk)/2,'lower');

Error in intd>@(t)fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) (line 4) E=@(t) fcdf((r-s)./s.*(x.*(r.*t-x.^2).^(1/2)-a.*b.*r.*t).^2./(a.^2.*r.*t-x.^2).^2,s,r-s).*fpdf(t,r,nu) ;

Error in quadv (line 61) y{j} = feval(f, x(j), varargin{:}); %#ok<AGROW>

Error in intd (line 7) prob=fcdf(x.^2./r,r,nu)+quadv(E,down,up)-1+al;

Error in @(x)intd(x,5,2,0.2104,60,0.05) 
%______________________________________________________

>>  h(6.7)

ans =

    0.0500

Answer 1

[Sean de Wolski]

This error doesn't look like a floating point issue but rather a real issue. I.e. at some point it looks like your data became complex.

Run:

dbstop if error

This will stop with the debugger when the error occurs, you will be able to look at the stack information and identify which of your variables has an imaginary component.

[Amykelly]

The issue does come directly from the imaginary part of a value which I am puzzled because no part in my function intd() could possibly generate a complex number. All the numbers I plugged in at fcdf() are real and reasonable.

[Sean de Wolski]

This is why I encourage you to use dbstop if error. It will allow you to see exactly what everything is at the time of the error.

Answer 2

[Image Analyst]

You call anonymous function h, which calls intd(). You pass in only x to h() but intd() requires r,s,q,nu,al in addition. Why are you not passing in those to your anonymous function h()?

[Amykelly]

Because all the other arguments are my parameters. I am only interested in the function with respect to x.