Double integration

How to perform double integration of exp(-ax-by)*(x^m)*(y^n)/(cx+dy)where x & y lies between 0 and infinity, a,b,m,n,c,d are positive real numbers.

Thank you.

Answer 1

[Mike Hosea]

If you have the 2012a release, just use INTEGRAL2:

>> a = 2; b = 3; m = 2; n = 3; c = 2; d = 3;
>> f = @(x,y)exp(-a*x-b*y).*(x.^m).*(y.^n)./(c*x+d*y)
f = 
    @(x,y)exp(-a*x-b*y).*(x.^m).*(y.^n)./(c*x+d*y)
>> integral2(f,0,inf,0,inf)
ans =
   0.0030864

Be sure to use .*, .^, and ./ to do elementwise operations instead of matrix ops.

If you don't have 2012a yet, you can use the solution I presented here: http://www.mathworks.com/matlabcentral/answers/14514-double-integral-infinite-limits

[Kala]

Thank you Mr.Mike Hosea for your help. Now I'm facing one more problem. Here is my code...

for l = 0:k-1 %k value is given in the beginning
    for m = 0:r %r value is also given in the beginning
        A1 = (nchoosek(k-1,l)*nchoosek(r,m)*(-1)^(l+m));
        f = @(x,y)exp(-a*x-b*y).*(x.^m).*(y.^n)./(c*x+d*y);
        I = integral2(f,0,inf,0,inf);
        A2 = A1 * I;
    end
end

Is this possible ?? Can we call Integral function inside the loops ?? And a and b values are calculated previously in the same program and our integral has to use the same values of a and b. please help.

[Mike Hosea]

It should work.

[Kala]

Yes sir.. It has worked. Thanks for your help. Now I have one more doubt.Here is my code..

part1=((xij)^(n1*k+alpha))*((yij)^(n2*r+beta))*k/(gamma(n2*r+beta)*gamma(n1*k+alpha));
I2=0;
for l=0:k-1
    for m=0:r
        part2=nchoosek(k-1,l)*nchoosek(r,m)*(-1)^(l+m);
        I=mydblquad(@(v,u) fun2(v,u,xij, yij, n1, n2, k, r, alpha, beta,l,m),0,inf,0,inf);
        I2=I2+part2*I;
    end
end
Rpbcap=I2*part1;
%Baye's estimator for Series
part3=((xij)^(n1*k+alpha))*((yij)^(n2*r+beta))*r/(gamma(n2*r+beta)*gamma(n1*k+alpha));
I3=0;
for m=0:r-1
    part4=nchoosek(r-1,m)*(-1)^(m);
    I1=mydblquad(@(v,u) fun3(v,u,xij, yij, n1, n2, k, r, alpha, beta,m),0,inf,0,inf);
    I3=I3+(part4*I1);
end
Rsbcap=I3*part3;
where
function fun2 = fun2(v,u,xij, yij, n1, n2, k, r, alpha, beta,l,m )
fun2=exp(-xij*u-yij*v).*u.^(n1*k+alpha).*v.^(n2*r+beta-1)./((l+1)*u+m*v);
end
and
function fun3 = fun3(v,u,xij, yij, n1, n2, k, r, alpha, beta, m)
fun3 = exp(-xij*u-yij*v).*u.^(n1*k+alpha-1).*v.^(n2*r+beta)./(k*u+(m+1)*v);
end

Here if I give n1 & n2 values >= 15 your mydblquad is showing NaN.. Is it out of range ? Please help..

[Mike Hosea]

I don't have all your input values to confirm it, but I think your integrand function is returning NaNs for some of the resulting input values. NaNs will be generated when you have overflow and try to calculate inf/inf or inf - inf.

Answer 2

[Thomas]

Use dblquad

http://www.mathworks.com/help/techdoc/ref/dblquad.html

doc dblquad

[Mike Hosea]

DBLQUAD won't handle infinite limits.

Answer 3

[Hussein Thary]

integral

alfa=1;fai=2*pi;theata=2*pi;L=5;k=10;w=2e-2

w = 0.0200

f=@(r, fai)exp(-alfa*L./2).*exp(-j*k*r*theata*cos(fai)).*exp((-r^2./w^2)-(j*fai))

f = function_handle with value:
    @(r,fai)exp(-alfa*L./2).*exp(-j*k*r*theata*cos(fai)).*exp((-r^2./w^2)-(j*fai))

s=integral2(f,0,100,0,2*pi)

Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

Error in solution>@(r,fai)exp(-alfa*L./2).*exp(-j*k*r*theata*cos(fai)).*exp((-r^2./w^2)-(j*fai)) (line 2)
f=@(r, fai)exp(-alfa*L./2).*exp(-j*k*r*theata*cos(fai)).*exp((-r^2./w^2)-(j*fai))

Error in integral2Calc>integral2t/tensor (line 237)
Z1 = FUN(X(VTSTIDX),Y(VTSTIDX)); NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

[Walter Roberson]

alfa=1;fai=2*pi;theata=2*pi;L=5;k=10;w=2e-2

w = 0.0200

f=@(r, fai)exp(-alfa.*L./2).*exp(-j.*k.*r.*theata.*cos(fai)).*exp((-r.^2./w.^2)-(j.*fai));
s=integral2(f,0,100,0,2*pi)

s = 0.0000 - 0.0027i