NIntegrate绊倒如果声明

我试着以下功能整合：NIntegrate绊倒如果声明

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]

奇怪的是，NIntegrate stumples于如果函数定义：

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

我很好奇，为什么是这样的话？我知道Integrate与If语句存在问题，但天真地说，我会假设NIntegrate归结为从函数定义中计算数字表。这与If的冲突如何？

我知道，通过将定义中的最后一个If语句替换为分段声明有助于NIntegrate获得正确的结果。

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

你知道短重写函数定义来胁迫NIntegrate吞下如果其他方法吗？

来源

2011-08-12 Markus Roellig

解决方案

没有什么错误使用If S，但：作为一般规则（预防），只要你通过函数NIntegrate，确保该函数不评估符号参数。其定义为：

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

为什么发生这种情况

的关键是：发生了什么，当你对一个象征性的说法，例如评估test？ Mathematica需要评估类似于If[x > 0, a, b]的东西，即If，其中条件既不是True也不是False。 x > 0不计算为True或False，因为x是Symbol没有值。结果是If将不评估其第2或第3个参数（a或b）。

来源

2011-08-12 12:14:21 Szabolcs

谢谢，那是关键。该_？NumericQ测试一直在咬我的背后... –

@Markus，是的，缺少'_？NumericQ'是一个很常见的问题，有没有简单的方法来解决，比如'NIntegrate'功能，使'_ ？NumericQ'不必要。新用户也不容易发现/学习他们需要在类似情况下证明其功能不符合象征性评估。这可能是MathGroup上最常见的投诉。 – Szabolcs

不是太困难，以避免如果语句：

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]

在[11]：= NIntegrate [mytest的[X] [PI]×^ 2，{X，3.5 10^-8， 6. 10^-8}]

缺货[11] = 2.18111 * 10^-21

在[12]：= NIntegrate [插值[表[{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 }]

出[12] = 2.18111×10^-21

来源

2011-08-12 11:09:48

我同意。在大多数情况下，如果不是所有情况下都可以避免。我仍然想知道为什么NIntegrate恨他们这么多...... –

如果'betag <0'，那么'dnda'等于'dnda'除以'（1-betag a/atg）'，所以我认为dnda = dnda（1 + Sign @ betag a/atg）'应该是类似于'dnda = dnda（1 + Signa betag a/atg）^Signβ'的东西。 – Heike

@喜欢是的，那看起来更好。仍然问题仍然存在，为什么问题发生在所有人。 –

NIntegrate绊倒如果声明

回答

相关问题