質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

90.98%

  • Python 3.x

    4127questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

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

解決済

回答 3

投稿 編集

  • 評価
  • クリップ 0
  • VIEW 136

motyaa

score 4

前提・実現したいこと

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/ツール等のバージョンなど)

なし

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 3

checkベストアンサー

+1

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

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

とりあえず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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/01/10 01:11

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

    キャンセル

+1

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/01/10 01:12

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

    キャンセル

0

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2018/01/09 23:08

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 90.98%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

  • 解決済

    正の整数nを1から9の数を掛けた数を半角スペース区切りで出力する方法を教えて欲しいです!

    前提・実現したいこと 1~9にnを掛けたものを出力したいのですが上手くいかないのでアドバイスが欲しいです。 n=整数として。 ,がついてきこないように下の方法でやってみました

  • 解決済

    Pythonのthreadingにおける、終了時の処理

    前提・実現したいこと Pythonを使って、動画や画像を表示するGUIを作っています。その中で、みなさまのご支援を頂ながらカメラの画像を取り込んで再生できるところまで来ました。

  • 解決済

    漸化式を求めたいがうまくかけない

    1 3 7 15・・・という数列のn番目の項を求めたいです。(2^(n-1)ずつ値が増える) def get_total(n): if n == 1:

  • 解決済

    if:条件NG時の再計算

    質問事項 質問タイトルがわかり辛くすみません。 pythonの基礎勉強をしています。 下記補足に記したように、whileループの中にif文を2つ作り、各々ifの条件次第で、各々

  • 解決済

    python 文字列を判定する方法

    何らかのstr型が連続して格納されている配列(例['100'],['.'],['abc'])の各要素を判定し、数値ならint型に変換後次要素の判定に移る、数値ではないのであればfa

  • 解決済

    Python 内包表記

    jupyter notebookでとある数列を求めるプログラムを作りました。 for文の中にあるfor文(for j in range(compare.size + 1)...)を

  • 解決済

    python for文

     前提・実現したいこと pythonでプログラム作成をしております。 簡易プログラムは以下です。 for k in all x = str(a) return x if n

  • 解決済

    Python 3.x 辞書のキー値によって変換する場合の高速化

    Pythonにて、辞書(dict({key,value})を使って、list型の全要素をValue値に変換する際の、 高速化が可能かどうかをご教授いただきたいです。 dict1

同じタグがついた質問を見る

  • Python 3.x

    4127questions

    Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。