矢量化梯形和辛普森积分规则。Trapezoidal 只是从另一个使用 logspace 而不是 linspace 的项目复制和粘贴，以便它可以利用非均匀网格。
def trap(func,a,b,num):
xlinear=np.linspace(a,b,num)
slicevol=np.diff(xlinear)
output=integrand(xlinear)
output=output[:,:-1]+output[:,1:]
return np.dot(output,slicevol)/2
def simpson(func,a,b,num):
a=float(a)
b=float(b)
h=(b-a)/num
output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1)
output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1)
output+=np.sum(integrand(b),axis=1)
output+=np.sum(integrand(a),axis=1)
return output*h/3
def integrand(rlin):
first=np.arange(1,11)[:,None]
second=np.arange(2,12)[:,None]
return np.vstack((rlin*first,np.power(rlin,second)))
检查梯形和辛普森规则累积相对误差：
b=float(100)
first=np.arange(1,11)*(b**2)/2
second=np.power(b,np.arange(3,13))/np.arange(3,13)
exact=np.vstack((first,second))
for x in range(3):
num=x*100+100
tr=trap(integrand,0,b,num).reshape(2,-1)
si=simpson(integrand,0,b,num).reshape(2,-1)
rel_trap=np.sum(abs((tr-exact)/exact))*100
rel_simp=np.sum(abs((si-exact)/exact))*100
print 'Number of points:',num,'Trap Rel',round(rel_trap,6),'Simp Rel',round(rel_simp,6)
Number of points: 100 Trap Rel 0.4846 Simp Rel 0.000171
Number of points: 200 Trap Rel 0.119944 Simp Rel 1.1e-05
Number of points: 300 Trap Rel 0.053131 Simp Rel 2e-06
时间。请注意，两个梯形规则都使用 200 点，而基于上述收敛性，辛普森一家的时间仅为 100。对不起，我没有同情：
s="""
import numpy as np
from scipy.integrate import trapz
def integrand(rlin):
first=np.arange(1,11)[:,None]
second=np.arange(2,12)[:,None]
return np.vstack((rlin*first,np.power(rlin,second)))
def trap(func,a,b,num):
xlinear=np.linspace(a,b,num)
slicevol=np.diff(xlinear)
output=integrand(xlinear)
output=output[:,:-1]+output[:,1:]
return np.dot(output,slicevol)/2
def simpson(func,a,b,num):
a=float(a)
b=float(b)
h=(b-a)/num
output=4*np.sum(integrand(a+h*np.arange(1,num,2)),axis=1)
output+=2*np.sum(integrand(a+h*np.arange(2,num-1,2)),axis=1)
output+=np.sum(integrand(b),axis=1)
output+=np.sum(integrand(a),axis=1)
return output*h/3
def simpson2(func,a,b,num):
a=float(a)
b=float(b)
h=(b-a)/num
p1=a+h*np.arange(1,num,2)
p2=a+h*np.arange(2,num-1,2)
points=np.hstack((p1,p2,a,b))
mult=np.hstack((np.repeat(4,p1.shape[0]),np.repeat(2,p2.shape[0]),1,1))
return np.dot(integrand(points),mult)*h/3
def x2(x):
return x**2
def x3(x):
return x**3
def x4(x):
return x**4
def x5(x):
return x**5
def x5(x):
return x**5
def x6(x):
return x**6
def x7(x):
return x**7
def x8(x):
return x**8
def x9(x):
return x**9
def x10(x):
return x**10
def x11(x):
return x**11
def xt5(x):
return 5*x
"""
zhenya="""
a=[xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,xt5,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11]
vectorize(quad)(a, 0, 100)
"""
usethedeathstar="""
g=lambda x: np.array([[x,2*x,3*x,4*x,5*x,6*x,7*x,8*x,9*x,10*x],[x**2,x**3,x**4,x**5,x**6,x**7,x**8,x**9,x**10,x**11]])
xv=np.linspace(0,100,200)
trapz(g(xv))
"""
vectrap="""
trap(integrand,0,100,200)
"""
vecsimp="""
simpson(integrand,0,100,100)
"""
vecsimp2="""
simpson2(integrand,0,100,100)
"""
print 'zhenya took',timeit.timeit(zhenya,setup=s,number=100),'seconds.'
print 'usethedeathstar took',timeit.timeit(usethedeathstar,setup=s,number=100),'seconds.'
print 'vectrap took',timeit.timeit(vectrap,setup=s,number=100),'seconds.'
print 'vecsimp took',timeit.timeit(vecsimp,setup=s,number=100),'seconds.'
print 'vecsimp2 took',timeit.timeit(vecsimp2,setup=s,number=100),'seconds.'
结果：
zhenya took 0.0500509738922 seconds.
usethedeathstar took 0.109386920929 seconds.
vectrap took 0.041011095047 seconds.
vecsimp took 0.0376999378204 seconds.
vecsimp2 took 0.0311458110809 seconds.
在时间安排中需要指出的是，真雅的答案应该更准确。我相信一切都是正确的，如果需要更改，请告诉我。
如果您提供您将使用的功能和范围，我可能会为您的系统提供更好的东西。您也有兴趣使用额外的核心/节点吗？