NIntegrate stumbling over If statement

I tried to integrate the following function开发者_JAVA技巧:

test[a_] := 
 Module[{dnda, cg, atg, acg, alphag, betag, sig, b1, b2, dndavsg, a01,
 a02, bc5},
 alphag = -1.96;
 betag = -0.813;
 atg = 6.93*^-6;
 acg = 0.000348;
 cg = 2.95*^-13;
 dnda = (cg/a) (a/atg)^alphag;
 If[betag >= 0, 
       dnda = dnda (1 + betag a/atg), 
       dnda = dnda/(1 - betag a/atg)
 ];
 If[a > atg, dnda = dnda Exp[((atg - a)/acg)^3]];
 a01 = 3.5 10^-8;
 a02 = 3 10^-7;
 sig = 0.4;
 b1 = 2.0496 10^-7;
 b2 = 9.6005 10^-11;
 bc5 = 4.;
 dndavsg = (b1/a) Exp[-0.5 (Log[a/a01]/sig)^2] + 
           (b2/a) Exp[-0.5 (Log[a/a02]/sig)^2];
 If[dndavsg >= 0.0001 dnda, dnda = dnda + bc5 dndavsg];
 dnda]

Curiously, NIntegrate stumples upon the If in the function definition:

In[604]:= NIntegrate[test[x]\[Pi] x^2,{x,3.5 10^-8,6. 10^-8}]
Out[604]= 1.95204*10^-23  

In[605]:= NIntegrate[Interpolation[Table[{x,test[x]\[Pi] x^2},{x,3 10^-8,10. 10^-8,.01 10^-8}]][x],{x,3.5 10^-8,6. 10^-8}]
Out[605]= 2.18089*10^-21

I am curious why this is the case? I am aware that Integrate has issues with If statements, but naively I would assume that NIntegrate boils down to calculating tables of numbers from the function definitions. How does this conflict with the If ?

I am aware that by replacing the last If statement in the definition by e.g. a Piecewise statement helps NIntegrate to get the correct result.

 dnda = Piecewise[{{dnda + bc5 dndavsg, dndavsg >= 0.0001 dnda}, {dnda,dndavsg < 0.0001 dnda}}]

Do you know any other ways short of rewriting the function definition to coerce NIntegrate to swallow the If?

Solution

There's nothing wrong with using Ifs, but: As a general rule (precaution), whenever you pass a function to NIntegrate, make sure that the function does not evaluate for symbolic arguments. Define it as:

Clear[test] (* get rid of previous definition *)
test[x_?NumericQ] := ...

Why this happens

The key is in: what happens when you evaluate test for a symbolic argument, e.g. test[x]? Mathematica will need to evaluate something similar to If[x > 0, a, b], i.e. an If where the condition is neither True or False. x > 0 does not evaluate to either True or False because x is a Symbol with no value. The result is that If will evaluate neither its 2nd or 3rd arguments (a or b).

Not too difficult to avoid the If statements:

mytest[a_] := 
 Module[{dnda, cg, atg, acg, alphag, betag, sig, b1, b2, dndavsg, a01,
    a02, bc5}, alphag = -1.96;
  betag = -0.813;
  atg = 6.93*^-6;
  acg = 0.000348;
  cg = 2.95*^-13;
  dnda = (cg/a) (a/atg)^alphag;
  (*If[betag>=0,dnda=dnda (1+betag a/atg),dnda=dnda/(1-betag a/
  atg)]; *)
  dnda = dnda (1 + Sign[betag] betag a/atg);
  (* If[a>atg,dnda=dnda Exp[((atg-a)/acg)^3]]; *)

  dnda = dnda Exp[Min[0, ((atg - a)/acg)^3]];
  a01 = 3.5 10^-8;
  a02 = 3 10^-7;
  sig = 0.4;
  b1 = 2.0496 10^-7;
      b2 = 9.6005 10^-11;
  bc5 = 4.;
  dndavsg = (b1/a) Exp[-0.5 (Log[a/a01]/sig)^2] + (b2/
       a) Exp[-0.5 (Log[a/a02]/sig)^2];
  (*If[dndavsg>=0.0001 dnda,dnda=dnda+bc5 dndavsg];*)

  dnda = dnda + HeavisideTheta[dndavsg - 0.0001 dnda] bc5 dndavsg;
  dnda]

In[11]:= NIntegrate[mytest[x] [Pi] x^2, {x, 3.5 10^-8, 6. 10^-8}]

Out[11]= 2.18111*10^-21

In[12]:= NIntegrate[ Interpolation[ Table[{x, mytest[x] [Pi] x^2}, {x, 3 10^-8, 10. 10^-8, .01 10^-8}]][x], {x, 3.5 10^-8, 6. 10^-8}]

Out[12]= 2.18111*10^-21