0 个投票

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.

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

 采纳的回答

Mike Hosea
Mike Hosea 2012-3-20

2 个投票

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

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Kala
Kala 2012-3-27
编辑：Walter Roberson 2023-6-17
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
Mike Hosea 2012-3-27
It should work.
Kala
Kala 2012-3-31
编辑：Walter Roberson 2023-6-17
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
Mike Hosea 2012-3-31
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.

请先登录，再进行评论。

更多回答（2 个）

Hussein Thary
Hussein Thary 2023-6-17
Ran in:

0 个投票

Ran in:
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);

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Walter Roberson
Walter Roberson 2023-6-17
Ran in:
Ran in:
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

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by