(Numerically) Integrating bump functions --- on R^1 vs R^2 and R^n

Hi everyone,

I'm still relatively new to MATLAB so any comments would be greatly appreciated! Thanks in advance!

This is a question on numerical integration of bump functions and mollifiers (i.e. http://en.wikipedia.org/wiki/Mollifier ). In particular, I'm interested in computing the very integral that's displayed in the aforementioned Wikipedia article, under the section "Concrete example".

I have no problem integrating this integral on R^1 using the following code.

psi_m1 = @(x) (  x.^2   < 1) .* exp(- 1 ./ ( 1 -  x.^2  ) ) ...
    + (  x.^2  >= 1) .* 0; 

integral(psi_m1, -inf, inf)

From this, I get the answer 0.4440.

Now, suppose I consider the analogous computation on R^2.

psi_m2 = @(x,y) (  x.^2 + y.^2   < 1) .* exp(- 1 ./ ( 1 - ( x.^2 + y.^2 ) ) ) + (  x.^2 + y.^2  >= 1 ) .* 0;  

psi_n = @(x,y) (1 / 2)^(-2) * psi_m2(2 * x, 2 * y); 

integral2(psi_n, -inf, inf, -inf, inf)

I get all sorts of errors of the form:

Warning: Infinite or Not-a-Number value encountered. 
> In funfun/private/integralCalc>iterateScalarValued at 349
  In funfun/private/integralCalc>vadapt at 133
  In funfun/private/integralCalc at 104
  In funfun/private/integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) at 18
  In funfun/private/integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) at 18
  In funfun/private/integralCalc>iterateScalarValued at 314
  In funfun/private/integralCalc>vadapt at 133
  In funfun/private/integralCalc at 104
  In funfun/private/integral2Calc>integral2i at 21
  In funfun/private/integral2Calc at 8
  In integral2 at 107 
Warning: The integration was unsuccessful. 
> In integral2 at 110

Note that I explicitly write out the Euclidean R^2 norm function here as I get all sorts of errors and problems when I attempt to use the norm() function along with integral() and integral2(). But let's brush that aside for a moment.

Question: Anybody have any comments or insights on this problem? As an aside --- this same computation works fine with Mathematica using the NIntegrate[] function there both on R^1, R^2, ...

Thanks!

Answer 1

[Mike Hosea]

Avoid masking the integrand in this manner if possible. It is numerically toxic. INTEGRAL2 is written to accept functional limits. Just as with the 1D problem you should write

>> fx = @(x)exp(-1./(1 - x.^2));
>> integral(fx,-1,1)
ans =
   0.443993816168071

with the 2D problem you should write

>> fxy = @(x,y)exp(-1./abs(1 - (x.^2 + y.^2)));
>> integral2(fxy,-1,1,@(x)-sqrt(1 - x.^2),@(x)sqrt(1 - x.^2))
ans =
   0.466512390886714

>> fxy2 = @(x,y)4*fxy(2*x,2*y);
>> integral2(fxy2,-0.5,0.5,@(x)-sqrt(0.25 - x.^2),@(x)sqrt(0.25 - x.^2))
ans =
   0.466512390886714

Or in polar coordinates:

>> fxy2p = @(r,theta)fxy2(r.*cos(theta),r.*sin(theta)).*r;
>> integral2(fxy2p,0,0.5,0,2*pi)
ans =
   0.466512392735157

[the cyclist]

I didn't see the forest for the trees with this one!

Answer 2

[the cyclist]

During the integration of your function, MATLAB evaluates your function at the values

xc = -0.263261412747471;
yc = -0.425378012941634;

At these "critical" values, your function returns a NaN, because:

(  x.^2 + y.^2   < 1)

is false (i.e. zero numerically), and

exp(- 1 ./ ( 1 - ( x.^2 + y.^2 ) ) )

evaluates as "Inf" (i.e. infinite). Your use of that first condition to suppress the out-of-range value leads to a NaN, and ultimately the inability to calculate the function. Such is the nature of a numerical calculator.

There are probably several ways around this. One simple but somewhat inelegant method would be to use

min(realmax,exp(- 1 ./ ( 1 - ( x.^2 + y.^2 ) ) ))

which will force that expression to be the largest expressible floating point number, so that it will get suppressed when multiplied by zero.