import numpy as np
from math import exp, sqrt, sin, cos, pi
"""
概要: ガウス・ルジャンドル法による数値積分を行うモジュール。
詳細説明: このモジュールは、ガウス・ルジャンドルの求積公式を用いて特定の関数の定積分を数値的に計算します。
様々な積分点数（2点、3点、4点、5点）に対応しており、定義された関数 `func(x)` とその解析解 `F(x)` を用いて結果を比較します。
関連リンク: :doc:`integ_gauss_legendre_usage`
"""
# define function to be integrated
[ドキュメント]
def func(x):
"""
概要: 積分対象の関数を定義する。
詳細説明: この関数は、ガウス・ルジャンドル法で数値積分される対象の関数 `f(x)` を表します。
現在、`exp(x)` が設定されていますが、コメントアウトされた行を有効にすることで他の関数に変更可能です。
:param x: 独立変数 (float)
:returns: xにおける関数の値 (float)
"""
return exp(x)
# return sqrt(1.0 - x*x)
# define analytical integration of f(x)
[ドキュメント]
def F(x):
"""
概要: 積分対象の関数の解析的な不定積分を定義する。
詳細説明: この関数は、`func(x)` の解析的な不定積分 `F(x)` を表します。
数値積分の結果と比較するための真値を提供するために使用されます。
現在、`exp(x)` の解析解が設定されていますが、コメントアウトされた行を有効にすることで他の関数に変更可能です。
:param x: 独立変数 (float)
:returns: xにおける解析的な不定積分の値 (float)
"""
return exp(x)
# return pi / 2.0
#===================
# parameters
#===================
xmin, xmax = -1.0, 1.0
xfrac = [
[], # 0-point, not dfined
[], # 1-point, not dfined
[-1.0 / sqrt(3), + 1.0 / sqrt(3)], # 2-point order
[-sqrt(3.0/5.0), 0.0, +sqrt(3.0/5.0)], # 3-point order
[-0.861136311594052, -0.339981043584856,
+0.339981043584856, +0.861136311594052], # 4-point order
[-0.906179845938663992797626878299, -0.538469310105683091036314420700, 0.0,
+0.538469310105683091036314420700, +0.906179845938663992797626878299] # 5-point order
]
weight = [
[], # 0-point, not dfined
[], # 1-point, not dfined
[1.0, 1.0], # 2-point order
[5.0/9.0, 8.0/9.0, 5.0/9.0], # 3-point order
[0.347854845137453, 0.652145154862546,
0.652145154862546, 0.347854845137453], # 4-point order
[0.236926885056189087514264040719, 0.478628670499366468041291514835, 0.568888888888888888888888888888,
+0.478628670499366468041291514835, +0.236926885056189087514264040719] # 5-point order
]
[ドキュメント]
def integ_GL(np, func, xmin, xmax):
"""
概要: ガウス・ルジャンドル法を用いて定積分を計算する。
詳細説明: この関数は、指定された積分点数 (`np`) に基づくガウス・ルジャンドルの求積公式を適用し、
`xmin` から `xmax` までの区間で `func` の定積分を数値的に計算します。
標準的なガウス・ルジャンドルの定義域 `[-1, 1]` から `[xmin, xmax]` への変数変換を行い、
重みと積分の和を計算します。ガウス点と重みはグローバル変数 `xfrac` および `weight` から参照されます。
:param np: 積分点数 (int)。2点、3点、4点、5点が有効です。
:param func: 積分対象の関数オブジェクト (callable)。xを引数にとり、数値を返す関数です。
:param xmin: 積分区間の下限 (float)
:param xmax: 積分区間の上限 (float)
:returns: 計算された定積分の近似値 (float)
"""
px = xfrac[np]
pw = weight[np]
S = 0.0
n = len(px)
xh = (xmin + xmax) / 2.0
h = xmax - xmin
for j in range(n):
x = xh + px[j] * h / 2.0;
S += pw[j] * func(x) * h;
S *= 0.5
return S
#===================
# main routine
#===================
[ドキュメント]
def main():
"""
概要: ガウス・ルジャンドル法による数値積分を実行し、結果を出力するメインルーチン。
詳細説明: この関数は、まず積分対象関数の解析値と解析積分値を出力します。
その後、異なる積分点数（3点、4点、5点）に対して `integ_GL` 関数を呼び出し、
それぞれの数値積分結果と使用されたガウス点および重みを出力します。
`xmin` と `xmax` はモジュールレベルで定義されたグローバル変数を使用します。
"""
print("Numerical integration using by Gauss-Legendre method")
fx = func(xmin)
Fx = F(xmax) - F(xmin)
print("")
print("Analytical values:")
print(" f({})={}".format(xmin, fx))
print(" integ(f)[{}, {}]={}".format(xmin, xmax, Fx))
print("")
for np in [3, 4, 5]:
print("{}-point formulat".format(np))
px = xfrac[np]
pw = weight[np]
print(" x=", px)
print(" w=", pw)
S = integ_GL(np, func, xmin, xmax)
print(" S={}".format(S))
print("")
if __name__ == '__main__':
main()