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{:});
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 =
