python怎么算积分_Python 中的科学计算模块 SciPy 里的积分函数效率如何？用这个是否比自己用 Fortran 或 C 写积分程序然后用 python 调用要快？...

6

2016-05-11 11:16:18 +08:00 3

效率已经很高了, 但无法满足较为苛刻的计算问题, 毕竟是 python 扩展(套上了 gil 的脚镣)

scipy 的 quad 函数底层是 f77 写的, 它会根据给定的参数决定具体调用哪个 fortran 的 subroutine 积分,

def quad(func, a, b, args=(), full_output=0, epsabs=1.49e-8, epsrel=1.49e-8,

limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50,

limlst=50):

if not isinstance(args, tuple):

args = (args,)

if (weight is None):

retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)

else:

retval = _quad_weight(func,a,b,args,full_output,epsabs,epsrel,limlst,limit,maxp1,weight,wvar,wopts)

首先分有权重跟无权重积分的情形, 这里只看看无权重的情形.无权重时调用_quad 函数,

def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):

......

省略多行处理 Inf 的代码

......

if points is None:

if infbounds == 0:

return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)

else:

return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)

else:

if infbounds != 0:

raise ValueError("Infinity inputs cannot be used with break points.")

else:

nl = len(points)

the_points = numpy.zeros((nl+2,), float)

the_points[:nl] = points

return _quadpack._qagpe(func,a,b,the_points,args,full_output,epsabs,epsrel,limit)

现在可以看到会调用 3 个不同的 fortran 函数 _qagse, _qagie, _qagpe 之一

找到 scipy 源码中的 f77 文件, 定位到 dqagse.f, dqagie.f, dqagpe.f 这几个文件,

发现:

1. dqagse.f 中调用的 subroutine dqk21; dqk21.f 中显示使用的积分方法是 gauss-kronrod rules( https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula)

2. dqagie.f 中调用的 subroutine dqk15i; dqk15i.f 中显示使用的积分方法也是 gauss-kronrod rules

3. dqagpe.f 中也是调用的 subroutine dqk21

由此可知, 这几个都是用的高斯积分方法, 从数值分析可以知道, 高斯方法的确比梯形法, 辛普森, 龙格库塔法有优越性: 对被积函数的拟合精度高, 收敛快. 从 fortran 代码中可以看出, 里面做了大量的错误检查跟精度的处理, 因为要适用于各种不同的平台, 以及各种不同的奇葩函数.

如果具体问题具体分析, 抛开各种检查的步骤以及平台相关的处理, 只是实现一个积分算法比这种现成的库函数更快还是极有可能的, 我用 c 简单写了一个自适应幸普森积分方法:

double

simp(integrand func, /*in*/

double left_end, /*in*/

double right_end, /*in*/

int n /*in*/

)

{

double result = 0.0;

double result0 = func(left_end) + func(right_end);

double h = (right_end - left_end) / n;

double even_sum = 0.0;

double odd_sum = 0.0;

double x;

for (int i = 0; i < n; ++i){

x = left_end + i*h;

if (i%2){

odd_sum += func(x);

}

else{

even_sum += func(x);

}

}

result = h * (result0 + 2.0 * even_sum + 4.0 * odd_sum) / 3.0;

return result;

}

double

adaptive_simp(integrand func, /*in*/

double left_end, /*in*/

double right_end, /*in*/

double tol /*in*/

)

{

int n = 20;

int max_splits = 20;

double result0, result;

result0 = simp(func, left_end, right_end, n);

for (int i = 0; i < max_splits; ++i){

n = n * 2;

result = simp(func, left_end, right_end, n);

if (fabs(result - result0) < tol)

break;

result0 = result;

}

return result;

并不困难, 很小的两个函数.

scipy 用 quadpack.py 对_quadpack.pyd 做了一层包装, 如果从 scipy import 的 quad 的话, 会比直接调用_quadpack.pyd 中相应的函数多一点点 python 的开销, 在 cython 中使用 quad, _quadpack.pyd 中的_qagse, 以及上面 c 实现的 simp(gcc 添加-O3 --fast-math 选项编译)在给定相同的误差容限下对 sinx(使用 c 标准库 sin)在[0.0, pi]区间积分 100W 次

-------------------------------------------

run quad 1000000 times, use 5.559261s

run _qagse 1000000 times, use 4.717242s

run simp 1000000 times, use 3.680429s

result from quad: 2.0

result from _qagse: 2.0

result from simp: 2.0000000001

-------------------------------------------

可以看到, 结果还是可以的.(没用 fortran 试了, 太久没写生疏了, 而且放到 cython 中还是要转一次 c, 比较麻烦)

另外在上面的 simp c 函数中累加部分可以很容易的用 openmp 并行化,曾经我对一个物体分割后积分,最后所有积分累加,这儿的累加也可以在 cython 中释放掉 gil 利用 openmp 并行化, 充分利用多核. 我那个计算问题一个结果大概要 100W 次积分, 还要放到上万次的迭代中(个人 PC 上算). 当现有的库函数无法满足效率要求时, 手写一个是明智的. 曾经还想过利用 opencl 把这个问题放到显卡上去烧, 但发现难以实现, 因为在迭代的时候总是有太多的逻辑判断跟信息交换, 只把可以纯粹的并行的部分放上去呢, 又会大大增加 CPU 跟 GPU 的信息交换, 效率又会极大的下降, 所以还是压榨的 CPU.