如何在Mathematica 8中并行计算积分

Question

有没有人知道如何使用所有核心计算积分？我需要使用并行化或并行表格，但怎么做呢？

 f[r_] := Sum[(((-1)^n*(2*r - 2*n - 7)!!)/(2^n*n!*(r - 2*n - 1)!))*
 x^(r - 2*n - 1), {n, 0, r/2}]; 


 Nw := Transpose[Table[f[j], {i, 1}, {j, 5, 200, 1}]]; 

 X1 = Integrate[Nw . Transpose[Nw], {x, -1, 1}]; 

 Y1 = Integrate[D[Nw, {x, 2}] . Transpose[D[Nw, {x, 2}]], {x, -1, 1}]; 

 X1//MatrixForm
 Y1//MatrixForm

- George Mills

---

还有在 math.SE 上的相关问题和 superuser 上的问题。 - Simon

Answer 1

我将列表的集成更改为集成列表，以便我可以使用ParallelTable：

X1par=ParallelTable[Integrate[i, {x, -1, 1}], {i, Nw.Transpose[Nw]}];

X1par==X1

(* ===> True *)

Y1par = ParallelTable[Integrate[i,{x,-1,1}],{i,D[Nw,{x,2}].Transpose[D[Nw,{x,2}]]}]

Y1 == Y1par

(* ===> True *)

在我的测试中，使用 {j, 5, 30, 1} 代替 {j, 5, 200, 1} 来限制一些时间，这样的速度大约是在我的四核处理器上快3.4倍。但如果采用以下方式，速度甚至可以更快:

X2par = Parallelize[Integrate[#, {x, -1, 1}] & /@ (Nw.Transpose[Nw])]

X2par == X1par == X1

(* ===> True *)

这样可以快6.8倍，2.3倍归功于Parallelize。在并行执行时，Timing和AbsoluteTiming不是非常可靠。我使用AbsoluteTime在每行前后取差值。

---

编辑

我们不应忘记ParallelMap：

在最粗糙的列表级别（1）：

ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {1}]  

在最深的列表级别（最细粒度的并行化）：

ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {2}]

- Sjoerd C. de Vries

---

能否使用简化后的公式加速运算？f[r_]=FullSimplify[Sum[(((-1)^n*(2r-2n-7)!!)/(2^nn!(r-2*n-1)!))x^(r-2n-1),{n,0,r/2}],r>0&&r[Element]Integers]简化后的公式为 $$ f(r)=-\frac{\sqrt{\pi } (-1)^r 2^{r-3} x^{r-1} ,_2\tilde{F}_1\left(\frac{1-r}{2},1-\frac{r}{2};\frac{7}{2}-r;\frac{1}{x^2}\right)}{\Gamma (r)}. $$ 但当使用以下公式时，却不起作用：f[r_]:=-1/Gammar^r 2^(-3+r) Sqrt[[Pi]] ξ^(-1+r) Hypergeometric2F1Regularized[(1-r)/2,1-r/2,7/2-r,1/ξ^2]; - George Mills

@Sjoerd，你的万人派对即将到来！你买足够的啤酒了吗？ - Dr. belisarius

@George 在你的评论中提到的函数中，Sqrt[[Pi]] 应该是 Sqrt[Pi]，[Xi] 应该是 [Xi]。此外，如果你定义了 Xi，它会加快运行速度。在这种情况下，仅并行化的开销似乎超过了并行化带来的任何收益。不确定为什么会这样。 - Sjoerd C. de Vries

@Sjoerd，那我就自己喝啤酒了 :) 预先恭喜你！ - Dr. belisarius

@mr.wizard 谢谢你推动我前进！哇，我现在感觉如此不同。 - Sjoerd C. de Vries

Answer 2

如果我们先扩展矩阵元素，再帮助整合一下，即使需要花费一点点的精力，也能够完成任务。

在装有Windows和Mathematica 8.0.4的四核笔记本电脑上，对于DIM=200，以下代码运行大约需要13分钟；对于DIM=50，代码运行时间为6秒。

$starttime = AbsoluteTime[]; Quiet[LaunchKernels[]]; 
DIM = 200; 
Print["$Version = ", $Version, "  |||  ", "Number of Kernels : ", Length[Kernels[]]]; 
f[r_] := f[r] = Sum[(((-1)^n*(-(2*n) + 2*r - 7)!!)*x^(-(2*n) + r - 1))/(2^n*n!*(-(2*n) + r - 1)!), {n, 0, r/2}]; 
Nw = Transpose[Table[f[j], {i, 1}, {j, 5, DIM, 1}]]; 
nw2 = Nw . Transpose[Nw]; 
Print["Seconds for expanding Nw.Transpose[Nm] ", Round[First[AbsoluteTiming[nw3 = ParallelMap[Expand, nw2]; ]]]]; 
Print["do the integral once: ", Integrate[x^n, {x, -1, 1}, Assumptions -> n > -1]]; 
Print["the integration can be written as a simple rule: ", intrule = (pol_Plus)?(PolynomialQ[#1, x] & ) :> 
     (Select[pol,  !FreeQ[#1, x] & ] /. x^(n_.) /; n > -1 :> ((-1)^n + 1)/(n + 1)) + 2*(pol /. x -> 0)]; 
Print["Seconds for integrating Nw.Transpose[Nw] : ", Round[First[AbsoluteTiming[X1 = ParallelTable[row /. intrule, {row, nw3}]; ]]]]; 
Print["expanding: ", Round[First[AbsoluteTiming[preY1 = ParallelMap[Expand, D[Nw, {x, 2}] . Transpose[D[Nw, {x, 2}]]]; ]]]]; 
Print["Seconds for integrating : ", Round[First[AbsoluteTiming[Y1 = ParallelTable[py /. intrule, {py, preY1}]; ]]]]; 
Print["X1 = ", (Shallow[#1, {4, 4}] & )[X1]]; 
Print["Y1 = ", (Shallow[#1, {4, 4}] & )[Y1]]; 
Print["seq Y1 : ", Simplify[FindSequenceFunction[Diagonal[Y1], n]]]; 
Print["seq X1 0 : ",Simplify[FindSequenceFunction[Diagonal[X1, 0], n]]]; 
Print["seq X1 2: ",Simplify[FindSequenceFunction[Diagonal[X1, 2], n]]]; 
Print["seq X1 4: ",Simplify[FindSequenceFunction[Diagonal[X1, 4], n]]]; 
Print["overall time needed in seconds: ", Round[AbsoluteTime[] - $starttime]]; 

- Rolf Mertig

---

尊敬的Rolf先生，如果我需要X1的两个积分，该如何修改此代码？X1 = aIntegrate[Nw . Transpose[Nw], {x, -1, 0.3}]+bIntegrate[Nw . Transpose[Nw], {x, 0.3, 1}]，其中a和b是常数。非常感谢您提供这些详细的信息，人们应该向您学习如何解决问题。非常感谢。 - George Mills