double integral infinite limits

This is easy to modify if you want c and/or d to be a function handle, like for QUAD2D. Also, it should be easy to add tolerances or whatever options you would like.

function y = mydblquad(fun,a,b,c,d)
% Double integral using QUADGK.
innerintegral = @(x)arrayfun( ...
    @(xi,y1i,y2i)quadgk(@(y)fun(xi*ones(size(y)),y),y1i,y2i), ...
    x,c*ones(size(x)),d*ones(size(x)));
y = quadgk(innerintegral,a,b);  
>> mydblquad(@(x,y)1./(1+x.^4+y.^4),-inf,inf,-inf,inf)
  ans =
     5.824747385361142

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

Gautam 2012-8-10

编辑：Gautam 2012-8-10

This was a very helpful post as I was wondering how I would achieve a double integration where the limits could be Inf. Following on Mike's post, I just wanted to confirm that the following modification to Mike's code to handle the case where the upper limit of the inner integral (the integration variable being y) is a function of x:

function y = mydblquad(fun,a,b,c,d) innerintegral = @(x) arrayfun(... @(xi,yli,y2i)quadgk(@(y)fun(xi*ones(size(y)),y),yli,y2i),... x,c*ones(size(x)),feval(d,x)); y = quadgk(innerintegral,a,b);

I tried the above and it seems to work fine. Any comments would be appreciated.

- Gautam

Mike Hosea 2012-8-10

One can indeed extend that method to use arbitrary lower and upper limit functions on the inner integral. I should note that the answer I gave above is somewhat old now. If you have the 2012a release of MATLAB you can use INTEGRAL2 "out of the box".

请先登录，再进行评论。