Search code examples
matlabnumerical-integration

Integral within exponential within integral

I am trying to numerically integrate a few different expressions, but they all have an integral within an exponential, surrounded by another integral, i.e.

integral image

I can't seem to get this to solve using the code below:

syms s p t
W = 3.18*10^(-22);
b = 10^(-23);
X = 300;
L = 0.374;
intlim = 0.589;

myfuncirc = @(s,p,t) (-W).*sqrt((L.^2)-(s.^2)).*(((-p.*cos(t))+sqrt(1-((p.^2).*((sin(t)).^2)))-sqrt(L.^2-s.^2)).^(-2));

s_min = 0;
s_max = L;
t_min = 0;
t_max = pi;

integral(@(p)(p.*exp(-(integral2(@(s,t) myfuncirc(s,t,p),s_min,s_max,t_min,t_max))/(b.*X))),0,intlim,'Arrayvalued',true)

I get the error message shown below, but I'm expecting a number between 0 and 1:

Warning: Infinite or Not-a-Number value encountered. In funfun\private\integralCalc>iterateArrayValued at 267 In funfun\private\integralCalc>vadapt at 130 In funfun\private\integralCalc at 75 In integral at 88

I'm also attempting a similar integration in the following form, but also not getting expected answers:

second integral image

pmax = y;
pmin = 0;
ymax = 1;
ymin = @(x) x;
xmax = 1;
xmin = 0;

integral3(@(x,y,p) (exp(-(integral(@(s)myfun(s,p),0,lam,'ArrayValued',true)./(k.*T)))),xmin,xmax,ymin,ymax,pmin,pmax,'Method','iterated')

Warning: Infinite or Not-a-Number value encountered. > In funfun\private\integralCalc>iterateScalarValued at 349 In funfun\private\integralCalc>vadapt at 132 In funfun\private\integralCalc at 75 In funfun\private\integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xiones(size(y)),y),y1i,y2i,opstruct.integralOptions) at 17 In funfun\private\integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xiones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) at 17 In funfun\private\integralCalc>iterateScalarValued at 314 In funfun\private\integralCalc>vadapt at 132 In funfun\private\integralCalc at 75 In funfun\private\integral2Calc>integral2i at 20 In funfun\private\integral2Calc at 7 In integral3>innerintegral at 137 In funfun\private\integralCalc>iterateScalarValued at 314 In funfun\private\integralCalc>vadapt at 132 In funfun\private\integralCalc at 75 In integral3 at 121 Warning: The integration was unsuccessful. > In integral3 at 125

Solution

  • Problem 1 is you define

    myfuncirc = @(s,p,t) ...

    But when it is called it is called as 

    myfuncirc(s,t,p)

    The arguments are not in the same order. You should change the definition of myfuncirc to

    myfuncirc = @(s, t, p) ...

    It's harder to tell for the second one because what you posted doesn't give the results you printed, it errors - several values are undefined such as myfun and k. However, I think the problem is that pmax is not defined as a function handle, but as a value - whatever is in y at the time pmax is assigned to. I suspect you want

    pmax = @(x,y) y;