试图让 Mathematica 逼近一个积分

Question

我试图让 Mathematica 逼近一个积分，它是各种参数的函数。我不需要它非常精确 - 答案将是一个分数，5 位数会很好，但我会满足于少至 2。
问题是主积分中隐藏了一个符号积分，我不能在它上面使用 NIntegrate 因为它是符号的。

    F[x_, c_] := (1 - (1 - x)^c)^c;
    a[n_, c_, x_] := F[a[n - 1, c, x], c];
    a[0, c_, x_] = x;

    MyIntegral[n_,c_] := 
      NIntegrate[Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}],{x,0,1}]

当n 大于 2 且 c 大于 3 左右时，Mathematica 开始挂起（通常因为 n 和 c 都高一点）。

是否有任何技巧可以重写此表达式以便更容易地对其进行评估？我在外部NIntegrate 上使用了不同的WorkingPrecision 和AccuracyGoal 和PrecisionGoal 选项，但是这些都没有帮助内部积分，这就是问题所在。事实上，对于n 和c 的较高值，我什至无法让 Mathematica 展开内导数，即

Expand[D[a[4,6,y],y]] 

挂起。

我正在使用 Mathematica 8 for Students。

如果有人对我如何让 M. 近似此有任何提示，我将不胜感激。

Answer 1

由于您只想要一个数字输出（或者这就是您将得到的），您可以使用 NIntegrate 将符号积分转换为数字积分，如下所示：

Clear[a,myIntegral]
a[n_Integer?Positive, c_Integer?Positive, x_] := 
  a[n, c, x] = (1 - (1 - a[n - 1, c, x])^c)^c;
a[0, c_Integer, x_] = x;

myIntegral[n_, c_] := 
 NIntegrate[D[a[n, c, y], y]*y/(1 - a[n, c, x]), {x, 0, 1}, {y, x, 1},
   WorkingPrecision -> 200, PrecisionGoal -> 5]

这比以符号方式执行集成要快得多。这是一个比较：

尤达：

myIntegral[2,2]//Timing
Out[1]= {0.088441, 0.647376595...}

myIntegral[5,2]//Timing
Out[2]= {1.10486, 0.587502888...}

rcollyer：

MyIntegral[2,2]//Timing
Out[3]= {1.0029, 0.647376}

MyIntegral[5,2]//Timing 
Out[4]= {27.1697, 0.587503006...}
(* Obtained with WorkingPrecision->500, PrecisionGoal->5, MaxRecursion->20 *)

Jand 的功能与 rcollyer 的功能类似。当然，当您增加n 时，您必须将您的WorkingPrecision 增加得比这更高，如you've experienced in your previous question。既然你说你只需要大约 5 位的精度，我已经明确地将 PrecisionGoal 设置为 5。你可以根据需要进行更改。

Answer 2

为了对 cme​​ts 进行编码，我会尝试以下方法。首先，为了消除关于变量n 的无限递归，我会将您的函数重写为

F[x_, c_] := (1 - (1-x)^c)^c;
(* see note below *)
a[n_Integer?Positive, c_, x_] := F[a[n - 1, c, x], c];  
a[0, c_, x_] = x;

那样n==0实际上将成为一个停止点。 ?Positive 表单是PatternTest，可用于将附加条件应用于参数。我怀疑问题是NIntegrate 正在为x 的每个值重新评估内部Integrate，所以我会取消评估，比如

MyIntegral[n_,c_] := 
  With[{ int = Integrate[(D[a[n,c,y],y]*y)/(1-a[n,c,x]),{y,x,1}] },
    NIntegrate[int,{x,0,1}]
  ]

其中With 是几个专门用于创建局部常量的范围构造之一。

您的 cmets 表明内部积分需要很长时间，您是否尝试过简化被积函数，因为它是 a 的导数乘以 a 的函数？对我来说，这似乎是链式规则扩展的结果。

注意：根据 Yoda 在 cmets 中的建议，您可以向 a 添加缓存或记忆机制。将其定义更改为

d:a[n_Integer?Positive, c_, x_] := d = F[a[n - 1, c, x], c];

这里的诀窍是，在d:a[ ... ] 中，d 是一个命名模式，在 d = F[...] 中再次使用它来缓存 a 的值以用于那些特定的参数值。