积分测试数值积分

Question

我试图检查收敛使用中点方法作为数值方法。我可以得到一张数据表，但我不确定如何绘制图表或添加日志比例。这是我的代码到目前为止。任何帮助深表感谢。

from math import exp,pi 
import numpy as np 
import matplotlib.pyplot as plt 

def midpoint(f,a,b,n): 
    h=float((b-a)/n) 
    result=0 
    for i in range(n): 
     result+=f((a+ h/2.0)+ i*h) 
    result*= h 
    return result 

#g= lambda y: (y**3)/(exp(y)-1) 
g=lambda y: ((exp(-y)*y**3)/(1-exp(-y))) 
a=0 
b=1000 
print( 'n    midpoint') 
for i in range(1,20): 
    n=2**i 
    m=midpoint(g,a,b,n) 
    print('%7d %.16f' %(n,m)) 

"""r=(pi**4/15) 
plt.plot(n,m,label="numberical") 
plt.plot(n,r,label="analytical") 
plt.xlabel('n') 
plt.ylabel('intergral') 
plt.legend() 
plt.show()""" 

"""plt.semilogx(n,m) 
    plt.show()""" 

Answer 1

我假设你想绘制近似而不是错误与间隔。要做出的关键更改为：

1）n和m应该是数组或列表; 2）由于确切的结果是一个对应于水平线的常数，所以应该从(xmin, r)到(xmax, r); 3）使用ax.set_xscale('log')获取日志缩放x轴。

完整的代码如下。我做了一些更改，尝试使用numpy来提高可读性。

from math import exp,pi 
import numpy as np 
import matplotlib.pyplot as plt 

def midpoint(f,a,b,n): 
    h = float((b - a)/n) 
    result = 0 
    for i in range(n): 
     result += f((a + h/2.0) + i * h) 
    result *= h 
    return result 

g = lambda y: ((exp(-y) * y ** 3)/(1 - exp(-y))) 
a = 0 
b = 1000 
n = 2 ** np.arange(1, 20) 
f = np.vectorize(midpoint) 
m = f(g, a, b, n) 
r = (pi**4/15) 

plt.plot(n, m, label="numberical") 
ax = plt.gca() 
ax.set_xscale('log') 
ax.set_ylim(auto=False) 
xmin, xmax = plt.xlim() 
plt.plot([xmin, xmax], [r, r], label="analytical") 
plt.xlabel('n') 
plt.ylabel('intergral') 
plt.legend() 
plt.show()