Pythonでルジャンドル多項式を使ってガウス積分をする方法

  • このエントリーをはてなブックマークに追加

数値解析の分野ではガウス積分という求積手法がよく用いられます。ガウス積分はルジャンドル多項式を使って積分点と重みを算出しますが、初学者はまずこれらの概念を理解するのが難しいです。この記事では簡単な関数を題材にガウス積分を計算する方法をPythonコードと共に紹介します。

こんにちは。wat(@watlablog)です。今回は数値解析の主要技術の1つ、ガウス積分をPythonで書いてみます

動作環境

この記事ではPythonコードを挟みながら説明しますが、以下のPython環境で動作確認をしました。

Python Python 3.9.6
PyCharm (IDE) PyCharm CE 2020.1
Numpy 1.23.0
matplotlib 3.7.1

ガウス積分の概要

ガウス積分の定義式

Advertisements

ガウス積分、またはガウス求積Gaussian Quadrature)とは、関数の積分値を近似計算する手法の一つです。この手法は特定の積分区間において、関数の値を積分点ガウス点)で評価します。
関数\(f(x)\)を\([a, b]\)の積分区間でガウス積分をする定義式を(1)に示します。

\[ \int_{a}^{b} f(x)dx \approx \sum_{i=1}^{n} w_{i}f(x_{i})dx \tag{1} \]

ここで\(x_{i}\)が積分点の位置、\(w_{i}\)が\(i\)番目の積分点における重みweightings)です。\(n\)が積分点の個数で次数order)とも呼ばれます。\(n\)次ガウス積分は\(2n-1\)次の多項式の積分値を正確に求めることができる特徴を持ちます。

積分の近似計算にはその他にも台形公式やシンプソン則を用いるものがありますが、高精度かつ計算回数の少なさからくる高速性が数値解析で使われる理由となっています。

積分点の求め方:ルジャンドル多項式を利用する

式(1)のガウス積分では積分点の位置を知る必要があります。積分点はルジャンドル多項式Legendre polynomials)のゼロ点を用います(ルジャンドル多項式については以下の外部サイトがわかりやすいと思います)。

高校数学の美しい物語:ルジャンドル多項式の性質と計算

百聞は一見にしかずということで、ルジャンドル多項式を可視化してみましょう。以下のPythonコードを実行することで各次数の多項式プロットとゼロ点を可視化できます。

こちらがグラフです。ルジャンドル多項式は[-1, 1]の範囲で定義される正規直交多項式であり、次数に応じて異なる数のゼロ点を持つ特徴があります。

ルジャンドル多項式の可視化結果

変数変換

実際の被積分関数は必ずしも[-1, 1]の範囲にあるわけではないので、式(2)の関係式を使って変数変換する必要があります。

\[ x = \frac{b-a}{2}\xi + \frac{b+a}{2} \tag{2} \]

式(2)は\(\xi=\)-1の時\(x=a\)、\(\xi=\)1の時\(x=b\)となります。

実際には先にルジャンドル多項式から積分点を計算しておき、式(1)でガウス積分を計算する段階で変数変換を反映させます。[-1, 1]の範囲におけるガウス積分式は式(3)です。

\[ \int_{-1}^{1} f(\xi)dx \approx \sum_{i=1}^{n} w_{i}f(\xi_{i})dx \tag{3} \]

式(2)の変換式を考慮すると、積分の演算子は式(4)の関係を持ちます。

\[ dx = \frac{b-a}{2}d\xi \tag{4} \]

式(4)を考慮すると変数変換を考慮したガウス積分は式(5)になります。

\[ \begin{eqnarray} \int_{a}^{b} f(x)dx &=& \int_{-1}^{1} f(\frac{b-a}{2}\xi + \frac{b+a}{2})\frac{b-a}{2}d\xi\\ &\approx&\frac{b-a}{2}\sum_{i=1}^{n} w_{i}f(\frac{b-a}{2}\xi_{i} + \frac{b+a}{2})d\xi \tag{5} \end{eqnarray} \]

それではこの内容をPythonでコーディングしてみましょう!

Pythonでガウス積分するコード(イメージ図付き)

理論や証明は難しいですが、使うだけなら簡単です。以下に任意の関数をガウス積分するコードを示します。f関数のreturnに関数形状を書きます。ルジャンドル多項式はnp.polynomial.legendre.leggaussを使って求めます。積分区間a, bと次数nを設定して実行してみましょう。

まずはn=2の結果です。黄色点で積分点(ガウス点)、矩形エリアに重みを考慮した面積をプロットしています。

積分点=2の結果

出力ウィンドウには「Integral value= 2.666666666666667」と出ていますが、式(6)から理論値と一致していることが確認できます。

\[ \begin{eqnarray} \int_{0}^{2} x^{2} dx &=& \frac{x^{3}}{3}\Big|_{0}^{2}\\ &=& \frac{8}{3} = 2.66666... \tag{6} \end{eqnarray} \]

以下は次数違いの積分結果です。ガウス積分の特性から、積分値はほぼ変わりませんが矩形形状の幅寸法が重みによって変わる様子を確認できると思います。

次数違いの積分結果

まとめ

本日はここまで。
この記事ではガウス積分を使って積分値を近似計算するための関連式を紹介し、Pythonで計算する方法を紹介しました。

ガウス積分をするためにはルジャンドル多項式を使い積分点と重みを算出する必要があること、そして変数変換が必要であることを学びました。

最後はmatplotlibを使い可視化を行うことで、ガウス積分のイメージを掴むことができました。
今後は多重積分における使用等に進みたいと思います。

また、以下の記事で有限要素法を学び始めたので、有限要素法の計算に取り入れるにはどうしたら良いかも考えていきたいと思います。

Pythonで基礎から学ぶ1次元梁モデルの有限要素法

難しいと思っていたガウス積分も使うだけなら簡単でした!
Xでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

SNSでもご購読できます。

コメントを残す

*