Problem of integral combined with fzero

Hi everyone,

I have encountered one problem when I want to use fzero to solve a integral equation, the detail is below

    rdata = normrnd(mu,sigma,1,n);
    xbar = mean(rdata);
    sd = std(rdata);
    aLehat = ( max( (xbar-T)*d/Du,(T-xbar)*d/Dl )./d_star )^2 + sd^2/d_star^2;

    delta = (xbar - T)./sd;
    b1 = @(y) sqrt(2./y).*delta;
    if delta >=0
        bL = @(y,c0) sqrt(n).*(Dl./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
        bU = @(y,c0) sqrt(n).*(Du./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
    else
        bL = @(y,c0) sqrt(n).*(Dl./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
        bU = @(y,c0) sqrt(n).*(Du./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
    end
    fun = @(y,c0) (    1./( exp( gammaln(alpha_no) ).*y.^(alpha_no+1) )   ).*exp(-1./y).* (   normcdf( b1(y)+bL(y,c0) ) - normcdf( b1(y)-bU(y,c0) )  );
    c0 = fzero(@(c0) (1-alpha) - integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))),   0.01);

The problem is last line because when I put c0 in the denominator of 2.*(aLehat./ c0) , it could not work. However, if I put the c0 in numerator, it can work but the result is not what I want.

Hope everyone can give me some advice, thanks. Error message is below:

_Error using erfc Input must be real and full.

Error in normcdf>localnormcdf (line 124) p(todo) = 0.5 * erfc(-z ./ sqrt(2));

Error in normcdf (line 46) [varargout{1:max(1,nargout)}] = localnormcdf(uflag,x,varargin{:});

Error in @(y,c0)(1./(exp(gammaln(alpha_no)).*y.^(alpha_no+1))).*exp(-1./y).*(normcdf(b1(y)+bL(y,c0))-normcdf(b1(y)-bU(y,c0)))

Error in @(y)fun(y,c0)

Error in integralCalc/iterateScalarValued (line 314) fx = FUN(t);

Error in integralCalc/vadapt (line 132) [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75) [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88) Q = integralCalc(fun,a,b,opstruct);

Error in @(c0)(1-alpha)-integral(@(y)fun(y,c0),0,(2.*(aLehat./c0)./(n.*aLehat).*((d./Du.*delta).^2+1)))

Error in fzero (line 363) a = x - dx; fa = FunFcn(a,varargin{:});

Answer 1

[Torsten]

Instead of specifying a point (0.01) in the call to "fzero", you should specify a search interval which excludes the point c0=0.

Something like

x0 = [1e-6 1];
x = fzero(fun,x0)

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

It could work now but the solving result is too strange and when I make the lower bound of interval be 1e-5, it does not work and the error message is below:

Error using fzero (line 290) The function values at the interval endpoints must differ in sign.

Can you tell me other advice? Thanks!

[Torsten]

Let f is the function for which "fzero" tries to find a root and [xleft xright] be the search interval you specify.

For "fzero" to work, it must be ensured that

f(xleft)<= 0 and f(xright)>=0

or

f(xleft)>=0 and f(xright)<=0.

Another possible solution for your problem is to immediately return a number different from 0 (e.g. Inf) to "fzero" if "fzero" wants to evaluate your function with an infeasible value for c0.

By the way: Before applying "fzero" to your function, I'd plot it to see whether its shape is reasonable.

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

I realize your answer now , so I try different combination for interval like [1e-6 1] or [0.5 10], but the value is too big or too small respectively.

When I notice this phenomenon, I also try other combination, but it will not solve the equation if interval range is smaller than orignal interval mentioned above.

Does it mean there is some limitations for Matlab to calculate ? Can we have other method to solve it?

Thanks!

[Torsten]

Another possible solution for your problem is to immediately return a number different from 0 (e.g. Inf) to "fzero" if "fzero" wants to evaluate your function with an infeasible value for c0.

Did you try this ?

And did you plot your function ?

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

I am sincerely feel sorry for you! I read your answer carefully but the part you say:

" immediately return a number different from 0 (e.g. Inf) to "fzero" if "fzero" wants to evaluate your function with an infeasible value for c0." and

" Before applying "fzero" to your function, I'd plot it to see whether its shape is reasonable", I did not try it. I am sorry because I do not understand and I omit it...

About the first part can you tell me more detail? How to operate it? It seems like we can change the fzero operation principle from

    f(xleft)<= 0 and f(xright)>=0

to

    f(xleft)<= inf and f(xright)>=inf

but I do not know how to bigger than inf...

About the second part, I also do not how to plot my function? I do not know which function you mention and if this function is

    integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  )))

I do not know how to use plot function to plot it because I am not sure how to set up the plot range...

I am not so smart and not so similar to Matlab. I have to be honest and I waste your enthusiasm, sorry! Hope you can give me another chance to get some advice from you and hope you can give me more detail about the two parts what you want me to do.

I am a Asian and it is time to sleep in my hometown, so maybe I could not reply to you in time. Please forgive me. These question is important for me now and sincerely hope you can help me again.

Thanks!

[Chao-Zhen Liu]

Hi Torsten,

Can you help me?

Thanks!

[Chao-Zhen Liu]

Does anyone can give me some advice?

Thanks!

[Torsten]

Try to plot your function and see whether the plot looks reasonable:

      rdata = normrnd(mu,sigma,1,n);
      xbar = mean(rdata);
      sd = std(rdata);
      aLehat = ( max( (xbar-T)*d/Du,(T-xbar)*d/Dl )./d_star )^2 + sd^2/d_star^2;
      delta = (xbar - T)./sd;
      b1 = @(y) sqrt(2./y).*delta;
      if delta >=0
          bL = @(y,c0) sqrt(n).*(Dl./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
          bU = @(y,c0) sqrt(n).*(Du./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
      else
          bL = @(y,c0) sqrt(n).*(Dl./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
          bU = @(y,c0) sqrt(n).*(Du./d).*sqrt((2.*(aLehat./c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))./y-1);
      end
      fun = @(y,c0) (    1./( exp( gammaln(alpha_no) ).*y.^(alpha_no+1) )   ).*exp(-1./y).* (   normcdf( b1(y)+bL(y,c0) ) - normcdf( b1(y)-bU(y,c0) )  );
      plotfun = @(c0)(1-alpha) - integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))),   0.01); 
      c = linspace(1e-5,1,1000)
      for i=1:numel(c)
        fc = plotfun(c(i));
      end
      plot(c,fc)

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

I have try it but the plot is like this:

It means that this function does not cross the x axis so we could not find a root? I am not sure and hope you can tell me more,

Thanks!

[Torsten]

I don't see a function. But if it is there and does not cross the x-axis, then yes: there is no root for c in between 1e-5 and 1.

Best wishes

Torsten.

[Torsten]

Replace the line

plotfun = @(c0)(1-alpha) - integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  ))),   0.01);

by

plotfun = @(c0)(1-alpha) - integral(@(y) fun(y,c0), 0,(2.*(aLehat./ c0)./(n.*aLehat).*(  (d./Du.*delta).^2+1  )));

Best wishes

Torsten.

[Chao-Zhen Liu]

Hi Torsten,

Actually, when I try to plot, I have do the same adjustment like the last comment you say, but the plot is what you see above...

Does any other direction to solve? Like Another possible solution for your problem is to immediately return a number different from 0 (e.g. Inf) to "fzero" if "fzero" wants to evaluate your function with an infeasible value for c0. , about this possible solution, how can I operate it?

Thanks!

[Torsten]

Further, replace

       fc = plotfun(c(i));

by

       fc(i) = plotfun(c(i));

Does it work now ?

[Chao-Zhen Liu]

Hi Torsten,

You are right! It can work now like this: [0.01 10]

[0.01 0.05]

So, it means that around [7 8], there is a solution?

[Torsten]

So, it means that around [7 8], there is a solution?

Correct.

[Chao-Zhen Liu]

Hi Torsten,

Thanks! I think it is something wrong for my code...

I will try it.