2017-11-17 51 views
0

リーマン和を使って積分を近似し、Pythonでmatplotlibを使ってグラフを作成するプログラムを作成しました。 x軸の上下に等しい面積を持つ関数の場合、結果の領域はゼロになりますが、代わりにプログラムから出力される数値は非常に小さいです。Python Riemann Sumは、正の領域と負の領域でゼロを返しません

次のコードは、奇関数f(x)= x^3を-1から1にグラフ化しているため、領域はゼロでなければなりません。代わりに私のコードは1.68065561477562 e^-15に近似しています。

この原因は何ですか? delta_x、x、またはyの丸め誤差ですか?私は値をゼロに丸めることができることを知っていますが、これを解決する別の問題または方法があるかどうかは疑問です。

私はdelta_xにDecimal.decimalクラスを使用しようとしましたが、ちょっとだけ小さい数字があります。

ザ・Pythonコード:

import matplotlib.pyplot as plt 
import numpy as np 

# Approximates and graphs integral using Riemann Sum 


# example function: f(x) = x^3 
def f_x(x): 
    return x**3 

# integration range from a to b with n rectangles 
a, b, n = -1, 1, 1000 

# calculate delta x, list of x-values, list of y-values, and approximate area under curve 
delta_x = (b - a)/n 

x = np.arange(a, b+delta_x, delta_x) 

y = [f_x(i) for i in x] 

area = sum(y) * delta_x 

# graph using matplotlib 
fig = plt.figure() 
ax = fig.add_subplot(111) 
ax.plot(x, y) 
ax.bar(x, y, delta_x, alpha=.5) 
plt.title('a={}, b={}, n={}'.format(a, b, n)) 
plt.xlabel('A = {}'.format(area)) 
plt.show() 

答えて

1

あなたが計算しているものは、元の意味でリーマン積分ではないことに注意する必要があります。間隔をnビンに分割していますが、n+1ビン(ここではn = 1000,len(x) == 1001)を合計してください。だから、結果はあなたが期待しているものに近いかもしれませんが、確かにそこに着くのは良い方法ではありません。

Riemann sumを使用すると、間隔がnビンに分割され、それらのビンの合計が合計されます(n)。左のリーマン和、右のリーマン和、または可能であれば中間点を計算するかどうかを選択できます。中間点合計がすでに良い結果を与える一方で

import numpy as np 

def f_x(x): 
    return x**3 

# integration range from a to b with n rectangles 
a, b, n = -1, 1, 1000 

delta_x = (b - a)/float(n) 

x_left = np.arange(a, b, delta_x) 
x_right = np.arange(a+delta_x, b+delta_x, delta_x) 
x_mid = np.arange(a+delta_x/2., b+delta_x/2., delta_x) 

print len(x_left), len(x_right), len(x_mid)   ### 1000 1000 1000 


area_left = f_x(x_left).sum() * delta_x 
area_right = f_x(x_right).sum() * delta_x 
area_mid = f_x(x_mid).sum() * delta_x 

print area_left  # -0.002 
print area_right # 0.002 
print area_mid  # 1.81898940355e-15 

、対称関数のためには、

print 0.5*(area_right+area_left) # 1.76204537072e-15 

これをしてもnを選択し、左右の和の平均を取るために良いアイデアのようになり

numpy.arangeはそれだけでエラーを出力します。より良い選択が

print area_left  # -0.002 
print area_right # 0.002 
print area_mid  # 8.52651282912e-17 
print 0.5*(area_right+area_left) # 5.68121938382e-17 

5.68121938382e-17を得numpy.linspace

x_left = np.linspace(a, b-delta_x, n) 
x_right = np.linspace(a+delta_x, b, n) 
x_mid = np.linspace(a+delta_x/2., b-delta_x/2., n) 

を使用されるだろうと、それは完全に0でない理由はかなり0に近い理由は確かにfloating point inaccuraciesです。

その有名な例がこの単に操作はリーマン積分として1E-17の順序の同じエラーを導入することを示すことである代わりに、0の5.5e-17もたらす 0.1 + 0.2 - 0.3あろう。

+0

奇数の部分区間しか使わない場合は十分にうまくいくでしょうし、左のリーマン和を '-1 + delta_x'にすることは' -1-delta_x/2 'から' 1 + delta_x/2 'になり、対称的なままになります。 – eugenhu

0

あなたのコードは、私には有効であると思われます。 がf_xによって返され、arangeの第2引数にdelta_xを追加することによって、1000の部分区間しかないという近似誤差が考慮されます。残念ながら、ここでは丸め誤差が発生していると思います。

1

はい、これは浮動小数点の不正確さによるものです。

np.linspace()は、np.arange()が非整数ステップサイズを使用するときに丸められる方法に応じてエンドポイントを含む場合と含まない場合があるため、より適切な場合があります。 numpy.arange()ドキュメントから

非整数のステップを使用して、例えば0.1として、結果はしばしば一貫性がないであろう。これらの場合には、linspaceを使用する方がよいでしょう。

+0

これは原則的に正しいですが、質問のコードにはここで説明する精度よりも深刻な欠陥があります。私の答えを見てください。 – ImportanceOfBeingErnest

関連する問題