python3.6でシンプソンの1/3公式を利用したプログラムの実装がわかりません。

Question

###前提・実現したいこと
python3.6で1/3シンプソンの公式を利用したプログラムを実装中です。
その際、きざみ幅がhの場合の近似値vと、きざみ幅をその半分にした場合の近似値wとを計算し、|v^w| <= 10^(-10)となった場合にプログラムを終了し、その際の計算結果を出力するというものです。

###発生している問題・エラーメッセージ

|v^w| <= 10^(-10)となった時の分割数と近似値を出力してbreakするように実装したつもりなのですが、
実行時に何も表示されなくて困っています。

###該当のソースコード

import math

def f1(x):
    return (x**3-2*x**2+x-1)
def f2(x):
    return (math.cos(x))

N = 10000        # 分割数

a_list = [-2,0]
b_list = [3,math.pi/2]
cnt = 0
for i in range(1,N,1):
    h = (b_list[0]-a_list[0])/i     # きざみ幅
    h2 = h/2
    v = (h/3) * (f1(h*i) + 4*f1(h*(i+1)) + f1(h*(i+2)))
    w = (h2/3) * (f1(h2*i) + 4*f1(h2*(i+1)) + f1(h2*(i+2)))
    if abs(v-w) <= 10**(-10):
        print("(1)\n")
        print("分割数 N = "+str(n)+"\n")
        print("近似値 = "+str(v)+"\n")
        break
for i in range(1,N,10):
    h = (b_list[1]-a_list[1])/i
    h2 = h/2
    v = (h/3) * ((f2(h*i) + 4*f2(h*(i+1)) + f2(h*(i+2))))
    w = (h2/3) * ((f2(h2*i) + 4*f2(h2*(i+1)) + f2(h2*(i+2))))
    if abs(v-w) <= 10**(-10):
        print("(2)\n")
        print("分割数 N = "+str(n)+"\n")
        print("近似値 = "+str(v)+"\n")
        break

###試したこと
for文内で|v^w|の値を出力をしてみましたが、とても10^(-10)になりそうもなく、小数点1桁くらいの大きい値が出力されました。

###補足情報(言語/FW/ツール等のバージョンなど)
なし

Answer 1

可能性１：
分割数10000程度では誤差が10^(-10)にはならないので、もっとNを大きくするか許容誤差を大きくするかしなければならない
（単純にiで割り算するのではなく、ループを繰り返すごとに割る数を2倍するような処理のほうが良いかもしれない）

可能性２：
コンピューターで扱う数値の精度や演算誤差の問題
浮動小数点演算、その問題と制限

Comments

[motyaa]

回答ありがとうございます。
計算方法にミスがありましたが、解決することができました。

Answer 2

判定条件が不明ですね。
積分としての和がありません。

私にはどう計算しているのかよく理解できませんでした。
もう少し詳しく説明していただけると要望にあった回答ができるかもしれません。

とりあえずWikipediaを参照して式を変更したところ、Scipyの積分モジュールと同じ計算結果を得ることができました。

import math
from scipy import integrate
import numpy as np

def f1(x):
    return x**3-2*x**2+x-1
def f2(x):
    return np.cos(x)

N0 = 100000
a_list = [-2,0]
b_list = [3,math.pi/2]

x = np.linspace(a_list[0], b_list[0], N0)
ans0 = integrate.simps(f1(x), x)
print('ans0', ans0)

x = np.linspace(a_list[1], b_list[1], N0)
ans1 = integrate.simps(f2(x), x)
print('ans1', ans1)

acc = -10

N = 100000
cnt = 0
ans0 = -99999
for i in range(2,N,2):
    h = (b_list[0]-a_list[0])/i
    x = np.arange(i+1)*h + a_list[0]
    x0 = x[2:i:2]
    x1 = x[1:i:2]
    v = (h/3) * (f1(x[0]) + 2*np.sum(f1(x0)) + 4*np.sum(f1(x1)) + f1(x[-1]))
    if abs(v-ans0) <= 10**(acc):
        print("(1)")
        print("分割数 N = " + str(i))
        print("近似値 = " + str(v))
        break
    ans0 = v

ans1 = -99999
for i in range(2,N,2):
    h = (b_list[1]-a_list[1])/i
    x = np.arange(i+1)*h + a_list[1]
    x0 = x[2:i:2]
    x1 = x[1:i:2]
    v = (h/3) * (f2(x[0]) + 2*np.sum(f2(x0)) + 4*np.sum(f2(x1)) + f2(x[-1]))
    if abs(v-ans1) <= 10**(acc):
        print("(2)")
        print("分割数 N = " + str(i))
        print("近似値 = " + str(v))
        break
    ans1 = v

Comments

[motyaa]

回答ありがとうございます。
頂いたソースコードで実行したところ、得たい結果を出力することができました。
定積分の結果とプログラムで得られた近似値の誤差も指定範囲内でしたので、解決致しました。
大変ありがとうございした。

Answer 3

つw = (h/3) → (h2/3)

Comments

[motyaa]

ご指摘ありがとうございます。
頂いたバグの部分を修正したのですが、変わらず何も出力されず、if文に入っていないようです。
また、プログラムを更新致しましたのでもし可能でしたらご確認下さい。