import csv
import numpy as np
from math import exp
"""
様々な近似法による数値積分の比較を行うスクリプト。
概要:
指数関数 `exp(x)` の定積分をリーマン和、台形公式、シンプソン公式、シンプソン3/8公式、ブール（Bode）公式といった
異なる数値積分法を用いて計算し、解析解と比較する。結果はCSVファイルに書き出される。
詳細説明:
このスクリプトは、与えられた関数 `func(x)` (ここでは `exp(x)`) の区間 `[x, x+h]` における定積分を、
複数の数値積分手法で近似計算する。
各 `x` の値について、数値積分結果と解析解（`integ(x+h) - integ(x)`）を比較し、
その結果を標準出力とCSVファイル `integ_order.csv` に出力する。
計算範囲とステップサイズは `x0`, `nx`, `h` パラメータで設定される。
関連リンク:
:doc:`integ_order_usage`
"""
[ドキュメント]
def func(x):
"""被積分関数 f(x) を定義する。
概要:
被積分関数 `f(x)` を定義する。
詳細説明:
この例では `exp(x)` を返す。
:param x: float 入力値。
:returns: float `exp(x)` の値。
"""
return exp(x)
[ドキュメント]
def integ(x):
"""解析的な原始関数 F(x) を定義する。
概要:
解析的な原始関数 `F(x)` を定義する。
詳細説明:
この例では `exp(x)` を返す。これは `exp(x)` の原始関数である。
:param x: float 入力値。
:returns: float `exp(x)` の値。
"""
return exp(x)
[ドキュメント]
def integ_rieman(func, x, h):
"""リーマン和を用いて数値積分を行う。
概要:
リーマン和を用いて関数 `func` の区間 `[x, x+h]` における数値積分を近似する。
詳細説明:
区間 `[x, x+h]` における左端の点 `x` での関数値 `func(x)` を用いて、
長方形の面積として積分を近似する（左リーマン和）。
:param func: callable 積分対象の関数。
:param x: float 積分区間の開始点。
:param h: float 積分区間の幅。
:returns: float リーマン和による積分の近似値。
"""
return func(x) * h
[ドキュメント]
def integ_trapezoid(func, x, h):
"""台形公式を用いて数値積分を行う。
概要:
台形公式を用いて関数 `func` の区間 `[x, x+h]` における数値積分を近似する。
詳細説明:
区間 `[x, x+h]` の両端 `x` と `x+h` での関数値を用いて、
台形の面積として積分を近似する。
:param func: callable 積分対象の関数。
:param x: float 積分区間の開始点。
:param h: float 積分区間の幅。
:returns: float 台形公式による積分の近似値。
"""
return (func(x) + func(x + h)) * 0.5 * h
[ドキュメント]
def integ_simpson(func, x, h):
"""シンプソン1/3公式を用いて数値積分を行う。
概要:
シンプソン1/3公式を用いて関数 `func` の区間 `[x, x+h]` における数値積分を近似する。
詳細説明:
区間 `[x, x+h]` を2つの部分区間に分け、3つの点（`x`, `x+0.5*h`, `x+h`）での関数値を用いて
二次多項式で近似し、積分を計算する。
:param func: callable 積分対象の関数。
:param x: float 積分区間の開始点。
:param h: float 積分区間の幅。
:returns: float シンプソン1/3公式による積分の近似値。
"""
return (func(x) + 4.0 * func(x + 0.5 * h) + func(x + h)) / 6.0 * h
[ドキュメント]
def integ_simpson38(func, x, h):
"""シンプソン3/8公式を用いて数値積分を行う。
概要:
シンプソン3/8公式を用いて関数 `func` の区間 `[x, x+h]` における数値積分を近似する。
詳細説明:
区間 `[x, x+h]` を3つの部分区間に分け、4つの点（`x`, `x+h/3.0`, `x+2.0*h/3.0`, `x+h`）での関数値を用いて
三次多項式で近似し、積分を計算する。
:param func: callable 積分対象の関数。
:param x: float 積分区間の開始点。
:param h: float 積分区間の幅。
:returns: float シンプソン3/8公式による積分の近似値。
"""
return (
3.0 * func(x)
+ 9.0 * func(x + h / 3.0)
+ 9.0 * func(x + 2.0 * h / 3.0)
+ 3.0 * func(x + h)
) / 24.0 * h
[ドキュメント]
def integ_bode(func, x, h):
"""ブール（Bode）公式を用いて数値積分を行う。
概要:
ブール（Bode）公式を用いて関数 `func` の区間 `[x, x+h]` における数値積分を近似する。
詳細説明:
区間 `[x, x+h]` を4つの部分区間に分け、5つの点（`x`, `x+h/4.0`, `x+h/2.0`, `x+3.0*h/4.0`, `x+h`）での関数値を用いて
四次多項式で近似し、積分を計算する。この公式はシンプソン公式よりも高精度な近似を提供する。
:param func: callable 積分対象の関数。
:param x: float 積分区間の開始点。
:param h: float 積分区間の幅。
:returns: float ブール公式による積分の近似値。
"""
return (
14.0 * func(x)
+ 64.0 * func(x + h / 4.0)
+ 24.0 * func(x + h / 2.0)
+ 64.0 * func(x + 3.0 * h / 4.0)
+ 14.0 * func(x + h)
) / 180.0 * h
# ===================
# parameters
# ===================
outfile = 'integ_order.csv'
x0 = 0.0
nx = 10
h = 0.5
[ドキュメント]
def main():
"""数値積分の比較を実行し、結果を標準出力とCSVファイルに出力する。"""
print("Numerical integration using different approximations")
print("Write to [{}]".format(outfile))
with open(outfile, 'w', newline='') as f:
fout = csv.writer(f, lineterminator='\n')
fout.writerow([
'x', 'f', 'integ(f)[x,x+h]',
'Riemtan',
'Trapezoid',
'Simpson',
'Simpson 3/8',
'Bode'
])
print("")
print("{:^3}:\t{:^10}\t{:^10}\t"
"{:^12}\t{:^12}\t{:^12}\t{:^12}\t{:^12}"
.format(
'x', 'f', 'integ(f)[x,x+h]',
'Riemtan',
'Trapezoid',
'Simpson',
'Simpson 3/8',
'Bode'
))
for ix in range(nx):
x = x0 + ix * h
fx = func(x)
Fx = integ(x + h) - integ(x)
print("{:^3}:\t{:^10.4f}\t{:^10.4f}\t".format(x, fx, Fx), end='')
res = [
integ_rieman(func, x, h),
integ_trapezoid(func, x, h),
integ_simpson(func, x, h),
integ_simpson38(func, x, h),
integ_bode(func, x, h),
]
print("{:^12.6f}\t{:^12.6f}\t{:^12.6f}\t{:^12.6f}\t{:^12.6f}\t".format(
res[0], res[1], res[2], res[3], res[4]
), end='')
fout.writerow([x, fx, Fx] + res)
print("")
if __name__ == "__main__":
main()