数值积分的收敛性检验

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()