Function argument plug-in puzzle

I made a function intd(). Given all the arguments except x, I am supposed 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 at xx(58) when I attempt the sequence xx=[1:0.1:10] for x. However, when I plug in the value of xx(58), which is 6.7, for x in h(), it gets through and gives me a perfectly reasonable answer. I am puzzled at why it makes a difference in my function, using a value in a sequence or a straight number for the argument.

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);

xx=[1:0.1:10]; h(xx(58))

%%%%%%%%% 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) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

>>xx(58)

ans =

    6.7000

>> h(6.7)

ans =

    0.0500

Answer 1

[the cyclist]

One thing to notice is that if you type

>> format long

you will see that xx(58) is actually equal to 6.699999999999999 and not 6.7. Furthermore, hh(6.699999999999999) causes an error, just like you are seeing.

So, what you are seeing is something of a numerical stability issue, caused by the floating point representation. (One needs to be very careful doing loops over real numbers.) See, for example, this thread on floating point representation: http://www.mathworks.com/matlabcentral/answers/69-why-does-1-2-3-1-3-not-equal-zero.

It looks like this particular case can be fixed if you use

>> xx = (1:100)/10

because then the 58th element will be 6.7.

I did not dig deeply into why such a tiny difference causes the error.

[Amykelly]

The redefined xx did work well. Is that to say that my function could not handle the floating point issue? How could I make my function insensitive to the floating issue?

On the other hand, my ultimate goal is to find a root to function h(). Then I used fzero(h, initial value) and got the exact same error message. I figured maybe running fzero brings the floating issue back. If I couldn't change my own function, is it possible to modify fzero to avoid the isse and do similarly what the redefined xx does? Or would there be any other root-finding functions that innately has no floating issues?

I appreciate your answer!

[the cyclist]

I think the fundamental issue is your function E(t) inside intd(). It looks like it is prone to numerical stability issue, due to the subtractions and ratios.

I don't really have a solution for you. One thing I can suggest is that if you type

 >> dbstop if error

before running your code, then execution will automatically stop when it hits the error and you will go into debug mode. Then you can see what the values are of each of your variables, and maybe diagnose more specifically what is going wrong.

Type

>> dbclear if error

when you are done, to turn off that debugging feature.