数値積分のためのカバーゲンステスト

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

誤差vs間隔ではなく近似をプロットすると仮定します。 makeの主な変更点は次のとおりです。

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