ガウス求積法の数値計算 #

1. ガウス求積法の概要 #

区間$x=[a,b]$で定義された$2N-1$次多項式$f(x)$と適当な重み関数$w(x)$の積の積分を考えます。この時、積分は特定の$N$個の離散点を評価することで厳密に計算できます。

$$ \begin{align} \int_{a}^{b}w(x)f(x)dx=\sum_{i=1}^N \omega_i f(x_i) \end{align} $$

ここで$x=[a,b]$は直交多項式の定義域,
$x_i, \omega_i$は, $N$次直交多項式 $p_N(x)$ の $p_N(x)=0$ を満たす$N$個の点$x=x_i,~~(i=1,2,\cdots, N)$と重みで、下記の数式で表現されます。

$$ \begin{gather} p_N(x)=0 \to x=x_i,\\ \omega_i= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_i)p_{N-1}(x_i)} \label{wdef} \end{gather} $$

ここで、$p’_n(x)$は$p_n(x)$の導関数$p’_n(x)=\frac{dp_n(x)}{dx}$を表し, また$h_n, k_n$は$p_n(x)$の直交性と係数から決まる下記の量です。

$$ \begin{gather} \int_a^b dx w(x) p_n(x) p_m(x)=h_n \delta_{n,m} \\ p_n(x)=k_n x^n+\cdots \end{gather} $$

2. 数値計算方法 #

2.1. 根の探索+定義 #

良く使用される方法です。

まず多項式のゼロ点を求めるために、適切な初期値を用意してゼロ点を数値的に求めます。その後、式\eqref{wdef}を用いて重みを計算する方法です。

2.1.1. ゼロ点の計算 #

多項式のゼロ点を求めます。

完全ではありませんが、多くの場合では次数$N$の漸近展開式を用いて、その第1, 2項までを考慮した式のゼロ点を解析的に解いて初期値とします。

ニュートン＝ラフソン法の詳細は省きますが、$f(x)=0$を満たす$x$を求めるために、下記の式に従って値を更新していく方法です。

$$ \begin{align} x_{n+1}=x_n-\frac{f(x_n)}{f’(x_n)} \end{align} $$

例えば、Legendre多項式$P_n(x), x=[-1, 1]$の$n$個のゼロ点を求めたいとします。$P_n(x)$のゼロ点を$x_m$と書き、更に下記の$\theta_m^{(n)}$を定義します。
$x_{n-m+1}^{(n)}=\cos\theta_m^{(n)},~~(0<\theta_1^{(n)}<\theta_2^{(n)}<\cdots<\theta_n^{(n)}<\pi)$
この時、$\theta_m^{(n)}$を求めるのであればニュートン＝ラフソン法の初期値として下記の式を用います¹。

$$ \begin{align} \theta_m^{(n)}=\frac{4m-1}{4n+1}\pi+\frac{1}{8n^2}\cot\frac{4m-1}{4n+2}\pi+O(n^{-3}) \end{align} $$

または解の下限と上限が指定できる¹

$$ \begin{align} \frac{2m-1}{2n+1}\pi \leq \theta_m^{(n)}\leq \frac{2m}{2n+1}\pi \end{align} $$

を利用して、二分法やBrent法などを用いて解きます。

その他の直交多項式については参考文献¹を参照してください。

やってみると分かりますが、初期値とニュートン法の組み合わせの方法はなかなか難しいです。十分近くないことが多く、必ず成功するとは限りません。また、上限下限が分かっていればまだやりようがありますが、そうではなく値の近似値だけが分かっている多項式の場合、精度が足らなければ取りこぼしが無いか確認するプログラムを追加しなければなりません。

2.1.2. 重みの計算 #

もし分点の位置が求められたら、次は重み\eqref{wdef}を計算します。そのために必要な係数は下記の通りです², ³。

	Sym.	$a$	$b$	$w(x)$	$h_n$	$k_n$
Jacobi	$P_n^{(\alpha,\beta)}(x)$	$-1$	$1$	$(1-x)^\alpha(1+x)^\beta$	$\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}$	$\frac{1}{2^n}\binom{2n+\alpha+\beta}{n}$
Legendre	$P_n(x)$	$-1$	$1$	$1$	$\frac{2}{2n+1}$	$\frac{(2n)!}{2^n(n!)^2}$
Associated Laguerre	$L_n^{(\alpha)}(x)$	$0$	$\infty$	$x^\alpha e^{-x}$	$\frac{\Gamma(\alpha+n+1)}{n!}$	$\frac{(-1)^n}{n!}$
Laguerre	$L_n(x)$	$0$	$\infty$	$e^{-x}$	$1$	$\frac{(-1)^n}{n!}$
Hermite	$H_n(x)$	$-\infty$	$\infty$	$e^{-x^2}$	$\sqrt{\pi}2^nn!$	$2^n$

また、Legendre求積法の系列ですが、端点で値が評価されるようにしたGauss-Radau, Gauss-Lobatto求積法も紹介します。

	Sym.	$a$	$b$	$w(x)$	$x_i$	$w_i$	Error
Legendre	$P_n(x)$	$-1$	$1$	$1$	$P_N(x)=0$	$\frac{2}{N P_{N-1}(x_k)P’_N(x_k)}$	$O(x^{2N})$
Radau⁴	-	$-1$	$1$	$1$	$i=1: x=-1$ $i=2,\cdots,N: \frac{P_{N-1}(x)+P_N(x)}{1+x}=0$	$\frac{1-x_i}{N^2[P’_{N-1}(x_i)]^2}$	$O(x^{2N-1})$
Lobatto⁵	-	$-1$	$1$	$1$	$i=1: x=-1$ $i=2,\cdots,N-1: P’_{N-1}(x)=0$ $i=N: x=+ 1$	$\frac{2}{N(N-1) [P_{N-1}(x_i)]^2}$	$O(x^{2N-2})$

ここで、$\alpha>-1,\beta>-1$です。

2.2. Golub-Welsch algorithm #

多項式の三項漸化式⁶と対角化を用いて直交多項式のゼロ点と重みを計算することができます。これは、Golub-Welsch algorithm と呼ばれるアルゴリズム(原論文⁷、説明⁸)で、三項漸化式を固有値問題に置き換えてから計算する方法で、精度の高い計算結果を得ることができます（日本語による説明⁹¹⁰）。

具体的にアルゴリズムの説明をしましょう。

2.2.1. ゼロ点の計算 #

直交多項式${p_n(x)}$には、三項漸化式が存在します。それは次のような形をしています⁶。

$$ \begin{gather} a_n p_{n+1}(x)=(b_n+c_n x)p_{n}(x)-d_n p_{n-1}(x),\\ (n=0,1,2,\cdots,~~p_{-1}(x)=0) \end{gather} $$

具体的な係数$a_n, b_n, c_n, d_n$は、例えば文献⁶から得られます。少し式変形して、

$$ \begin{gather} xp_{n}(x)=\frac{d_n}{c_n}p_{n-1}(x)-\frac{b_n}{c_n}p_{n}(x)+\frac{a_n}{c_n}p_{n+1}(x) \end{gather} $$

を得ます。これを行列形式として書くと、下記のように書けます。

$$ \begin{align} \label{eqxJ} x\left( \begin{array}{c} p_0(x)\\ p_1(x)\\ \vdots\\ p_{N-1}(x) \end{array} \right)= \left( \begin{array}{ccccc} -\frac{b_0}{c_0} & \frac{a_0}{c_0} & 0 & \cdots & 0\\ \frac{d_1}{c_1} & -\frac{b_1}{c_1} & \frac{a_1}{c_1} & \cdots & 0\\ \vdots&\vdots&\ddots&\vdots&\\ 0 & 0 & \cdots & -\frac{b_{N-2}}{c_{N-2}} & \frac{a_{N-2}}{c_{N-2}}\\ 0 & 0 & \cdots & \frac{d_{N-1}}{c_{N-1}} & -\frac{b_{N-1}}{a_{N-1}}\\ \end{array} \right) \left( \begin{array}{c} p_0(x)\\ p_1(x)\\ \vdots\\ p_{N-1}(x) \end{array} \right)+ \frac{a_{N-1}}{c_{N-1}} \left( \begin{array}{c} 0\\ 0\\ \vdots\\ 0\\ p_{N}(x) \end{array} \right) \end{align} $$

これをさらに書き変えて

$$ \begin{align}\label{eqTN} x\mathbf{p}(x)=T\mathbf{p}(x)+\frac{p_N(x)}{a_{N-1}}\mathbf{e}_N \end{align} $$

と書きます。ここで、

$\mathbf{p}(x)=(p_0(x), p_1(x), \cdots, p_{N-1}(x))^{T}$は, $N\times 1$の行列,
$\mathbf{e}(x)=(0, 0, \cdots, 0, 1)^{T}$は, $(N+1)\times 1$の行列,
$T$は, 式\eqref{eqxJ}に現れる$N\times N$の三重対角行列

を意味します。

もし、$x$が$p_N(x)=0$を満たす$x=x_j$であれば、式\eqref{eqTN}の右辺最後の項はゼロになるので、下記の方程式に書き変えられます。

$$ \begin{align}\label{eqTjN} T\mathbf{p}(x_j)=x_j\mathbf{p}(x_j) \end{align} $$

つまりこれは、三重対角行列$T$に対する固有値問題となり、$T$の固有値として$N$個のゼロ点$x_j$が導かれ、それぞれの固有値$x_j$に属する固有ベクトルとして$\mathbf{p}(x_j)=(p_0(x_j), p_1(x_j), \cdots, p_{N-1}(x))^T$が得られることを意味します。

式\eqref{eqTN}をそのまま解いても悪くないですが、対称化してから解くと更に良いです。

三項漸化式に$\kappa_n$を掛けて、更に$p’_n(x)=\kappa_n p_n(x)$となる量を導入します。すると、${p’_n}$に対する三項漸化式

$$ \begin{gather} xp’_{n}(x)=\beta_{n-1} p’_{n-1}(x)+\alpha_n p’_{n}(x)+\beta_n p’_{n+1}(x) \\ \alpha_n \equiv -\frac{b_n}{c_n},~~~\beta_n \equiv \left(\frac{a_n d_{n+1}}{c_n c_{n+1}}\right)^{1/2} \end{gather} $$

を得ます。すると、固有値問題\eqref{eqTjN}は、

$$ \begin{align}\label{eqJjNe} J\mathbf{p}’(x_j)=x_j\mathbf{p}’(x_j) \end{align} $$

となります。具体的に$J$を書けば、

$$ \begin{align} \label{eqJjN} DTD^{-1}=J= \left( \begin{array}{ccccc} \alpha_0 & \beta_0 & 0 & \cdots & 0\\ \beta_0 & \alpha_1 & \beta_1 & \cdots & 0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0 & 0 & \cdots & \alpha_{N-2} & \beta_{N-2}\\ 0 & 0 & \cdots & \beta_{N-2} & \alpha_{N-1}\\ \end{array} \right) \end{align} $$

となります。$\kappa$は漸化式によって与えられ、それは

$$ \begin{align} \kappa_{n+1}=\left(\frac{a_n c_{n+1}}{c_n d_{n+1}}\right)^{1/2} \kappa_{n},~~\kappa_0=1 \end{align} $$

と書かれます。

$J$を対角化して得られた固有ベクトルを$\mathbf{q}’(x_j)$と書きます。この$\mathbf{q}’(x_j)$は、本来欲しい多項式の値$\mathbf{p}(x_j)$とは異なります。異なる点は、２つで、それは

$\mathbf{q}’(x_j)$のそれぞれの要素が$\kappa_n$倍されている
$\mathbf{q}’(x_j)$の全体の定数倍が決まらない（式\eqref{eqJjNe}より）

という点です。

第一の点は、それぞれの要素に対して$\kappa_n$倍すればこの問題は解決します。つまり、

$$ q_n(x_j)\equiv \frac{q’_n(x_j)}{\kappa_n} $$

を実行して$\mathbf{q}(x_j)$を得ます。このようにして$\mathbf{q}(x_j)$を固有ベクトルとして考えれば、$\mathbf{q}(x_j)$と$\mathbf{p}(x_j)$の違いは全体の定数倍だけとなります。

第二の点は、$\mathbf{q}(x_j)$と$\mathbf{p}(x_j)$の違いは定数倍だけなので、多項式の0次を合わせるようにすれば十分です。つまり、$\mathbf{p}(x_j)$は

$$ \begin{align}\label{pqrel} \mathbf{p}(x_j)=\frac{k_0}{q_0(x_j)}\mathbf{q}(x_j) \end{align} $$

または$\mathbf{q}’(x_j)$を用いて

$$ \begin{align}\label{pqrela} p_n(x_j)=\frac{k_0}{q’_0(x_j)} \kappa_0\frac{q’_n(x_j)}{\kappa_n},~~(j=0,1,\cdots,N-1) \end{align} $$

を計算すれば得られることが分かります。ここで$k_0$は下記$k_n$の$n=0$です³。

$$ p_n(x)=k_nx^n+\cdots $$

さて、ゼロ点の位置は固有値から求められますが、重みは式\eqref{wdef}に従って計算するのは少し勿体ないです。固有ベクトルの情報があるので、これを使って重みを求める式に変形しましょう。

2.2.2. 重みの計算 #

次の式\eqref{CDt}のクリストッフェル＝ダルブーの定理 (Christoffel–Darboux theorem) を用いて重みの式を変形していきます。

$$ \begin{align}\label{CDt} \sum_{\nu=0}^{N-1} \frac{p_\nu(x)p_\nu(y)}{h_\nu} = \frac{k_{N-1}}{h_{N-1} k_{N}}\frac{p_{N}(x)p_{N-1}(y)-p_{N-1}(x)p_{N}(y)}{x-y} \end{align} $$

式\eqref{CDt}で、特に$y\to x$の極限を考えます。右辺は、

$$ \begin{gather}\label{CDtlim} \frac{k_{N-1}}{h_{N-1} k_{N}}\lim_{y\to x}\left[ p_{N}(x)\frac{p_{N-1}(y)-p_{N-1}(x)}{x-y}+p_{N-1}(x)\frac{p_{N}(x)-p_{N}(y)}{x-y}\right] \\ =\frac{k_{N-1}}{h_{N-1} k_{N}}\left[-p_N(x)p’_{N-1}(x)+p_{N-1}(x)p’_N(x)\right] \end{gather} $$

となるので、関係式

$$ \begin{align} \sum_{\nu=0}^{N-1} \frac{[p_\nu(x)]^2}{h_\nu} = \frac{k_{N-1}}{h_{N-1} k_{N}}\left[-p_N(x)p’_{N-1}(x)+p_{N-1}(x)p’_N(x)\right] \end{align} $$

が導かれます。特に、$x$を$p_N(x)=0$を満たすゼロ点$x=x_j$とすると

$$ \begin{align}\label{CDnu} \sum_{\nu=0}^{N-1} \frac{[p_\nu(x_j)]^2}{h_\nu} = \frac{k_{N-1}}{h_{N-1} k_{N}}p_{N-1}(x_j)p’_N(x_j) \end{align} $$

となります。式\eqref{CDnu}の右辺は、重み\eqref{wdef}の右辺に等しいので、結果として

$$ \begin{align} \omega_j&= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_j)p_{N-1}(x_j)}=\left[\sum_{\nu=0}^{N-1} \frac{p_\nu^2(x_j)}{h_\nu}\right]^{-1}\label{wdef2} \end{align} $$

を得ます。

固有ベクトルで計算する場合、式\eqref{pqrel}から求めた固有ベクトルを用いて、式\eqref{wdef2}を計算する方が良いでしょう。つまり、

$$ \begin{align} \omega_j &=\left[\sum_{\nu=0}^{N-1} \frac{p_\nu^2(x_j)}{h_\nu}\right]^{-1}\\ &=\left(\frac{q_0(x_j)}{k_0}\right)^2 \left[\sum_{\nu=0}^{N-1} \frac{q_\nu^2(x_j)}{ h_\nu}\right]^{-1},\\ &(j=0, 1, \cdots, N-1)\nonumber \end{align} $$

として計算します。

3. プログラム #

PythonでGolub-Welschアルゴリズムを実装しました。

下記URLからダウンロード・解凍してください。
https://slpr.sakura.ne.jp/sikinote/supplement/python/numeric/Integration/gauss_quadrature_abscissa_weight.zip (約4kB)

基本的には Abramowitz and Stegun., “Handbook of Mathematical Functions.” (以降A&S) を参照して、若干数値計算的には良くありませんが、極力 A&Sに現れているそのままの形でプログラムできるように実装しています。
ガウス＝ロバート求積法についてはガウス求積法の導出で導出していますので、こちらの情報を用いて実装します。

現状実装しているのは、６種であり、下記の関数を呼ぶとg_kindで指定した求積法の分点位置と重みを返します。

x, w = gauss_quadrature_abscissa_weight(n_points, g_kind="Chebyshev1")

g_kind	$a$	$b$	$w(x)$	After Scaling
Chebyshev1	$-1$	$1$	$(1-x^2)^{-1/2}$	$\int_{a}^{b} [(x-a)(b-x)]^{-1/2} f(x) dx$
Chebyshev2	$-1$	$1$	$(1-x^2)^{1/2}$	$\int_{a}^{b} [(x-a)(b-x)]^{1/2} f(x) dx$
Legendre	$-1$	$1$	$1$	$\int_{a}^{b} f(x) dx$
Lobatto	$-1$	$1$	$1$	$\int_{a}^{b} f(x) dx$
Associated Laguerre	$0$	$\infty$	$x^\alpha e^{-x}$	$\int_{0}^{\infty} x^{\alpha} e^{-\kappa x}f(x) dx$
Hermite	$-\infty$	$\infty$	$e^{-x^2}$	$\int_{-\infty}^{\infty} e^{-\kappa x^2}f(x) dx$

適当なスケーリングによって、上記表の After Scaling の積分に対応できます。

zipを解凍したままであれば、$f(x)=x^4$で指定しており、下記の実行結果が得られます。

> python main.py
Chebyshev1 1.1780972450961693
Chebyshev2 0.1963495408493613
Legendre 0.39999999999999875
Lobatto 0.39999999999999997
Laguerre 0.6163885883547535
Hermite 0.2349964007466552
>

これらは下記の積分を計算していることになります。

$$ \begin{align} \text{Chebyshev1} &: \int_{-1}^{1} [(1-x^2)]^{-1/2} x^4 dx = 3\pi/8 \approx 1.17809724509617… \\ \text{Chebyshev2} &: \int_{-1}^{1} [(1-x^2)]^{1/2} x^4 dx =\pi/16 \approx 0.196349540849362… \\ \text{Legendre} &: \int_{-1}^{1} x^4 dx = 2/5 = 0.4 \\ \text{Lobatto} &: \int_{-1}^{1} x^4 dx = 2/5 = 0.4 \\ \text{Laguerre} &: \int_{0}^{\infty} x^{-0.25} e^{-2 x} x^4 dx \approx 0.61638858835475457… \\ \text{Hermite} &: \int_{-\infty}^{\infty} e^{-2 x^2}x^4 dx = (3 \sqrt(\pi/2))/16 \approx 0.2349964007466563… \end{align} $$

4. Appendix #

4.1. ガンマ関数に関係する数値計算 #

$$ \begin{align} \Gamma(z+1)=z\Gamma(z) \end{align} $$

$$ \begin{align} \frac{\Gamma(n+b)}{\Gamma(n+a)}=\frac{\Gamma(b)}{\Gamma(a)}\cdot\prod_{k=1}^{n}\frac{n-k+b}{n-k+a} \end{align} $$

$h_n=\frac{\Gamma(\alpha+n+1)}{n!}$を計算する場合、下記の漸化式を用いると良い。

$$ \begin{gather} h_n=\frac{\alpha+n}{n}h_{n-1},~~h_0=\Gamma(\alpha+1) \end{gather} $$