Double integration

12 次查看(过去 30 天)
Kala
Kala 2012-3-20
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
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 个评论
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 个)

Thomas
Thomas 2012-3-20

Hussein Thary
Hussein Thary 2023-6-17
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 个评论
Walter Roberson
Walter Roberson 2023-6-17
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

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Specifying Target for Graphics Output 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by