Problem of integral combined with fzero

Question

Chao-Zhen Liu 2018-9-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/419005-problem-of-integral-combined-with-fzero

评论： Chao-Zhen Liu 2018-9-17

采纳的回答： Torsten

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-9-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/419005-problem-of-integral-combined-with-fzero#answer_336810

在 MATLAB Online 中打开

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.

16 个评论
显示 14更早的评论隐藏 14更早的评论

Chao-Zhen Liu 2018-9-14

在 MATLAB Online 中打开

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!

Torsten 2018-9-17

在 MATLAB Online 中打开

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.

Torsten 2018-9-17

编辑：Torsten 2018-9-17

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

Correct.

Chao-Zhen Liu 2018-9-17

Hi Torsten,

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

I will try it.

请先登录，再进行评论。

Problem of integral combined with fzero

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

16 个评论
显示 14更早的评论隐藏 14更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Problem of integral combined with fzero

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

16 个评论 显示 14更早的评论隐藏 14更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

16 个评论
显示 14更早的评论隐藏 14更早的评论