ガウス求積法の導出 #

ガウス求積法は、

の形の関数を積分する場合に最適な数値積分法です。
数値積分に使用する分点の位置が制限されますが、$N$点だけ非積分関数を参照することで
(重み関数)×($2N-1$次の多項式) の被積分関数を厳密に積分することができます。

ガウス求積法は、少ない分点数で高精度の結果が得られるため科学技術計算の分野で特に効果を発揮します。本稿ではガウス求積法はどのように導かれたかについて説明します。

1. 問題設定 #

区間$x=[a,b]$で定義された重み関数$\omega(x)$と、同じ区間で定義された$M$次多項式である関数$f(x)$の積を、$N$個の点を用いて数値的に積分することを考えます。

数式で問題を表しましょう。$f(x)$は$M+1$個の既知の係数$a_m$を用いて

$$ \begin{align}\label{e1} f(x)=\sum_{m=0}^{M}a_m x^m \end{align} $$

と書けているときに、積分

$$ \begin{align}\label{e2} \int_{a}^{b}\omega(x)f(x)dx \end{align} $$

を解くことが目的となります。

数値的に積分を行う場合、コンピュータでは連続的な関数は扱えないので、被積分関数の離散的な点だけを用いて計算を試みます。つまり、適当な$N$個の点を用いて

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

となるような右辺の $x_i, \omega_i, (i=1,2,\cdots,N)$ とは何か？を求めることが本稿で考える問題となります。

ここで、$\omega_i$は “重み"と呼ばれ、それぞれの分点の位置$x_i$で評価した被積分関数の積分結果への影響度を表しています。

2. ガウス求積法 #

2.1. 離散点数と到達可能な次数の関係 #

例えば、$f(x)$の無限個の値を参照することができれば、厳密な積分を実行することができます。また、$f(x)$が一次関数であれば、始点($x=a$) と終点($x=b$) の2点における値が判明していれば厳密な積分を実施することができます。

つまり、多項式$f(x)$の次数$M$と離散点数$N$の間には関係式があることが予想できるので、まずはその関係式を導きます。

式\eqref{e2}に式\eqref{e1}を代入すると

$$ \begin{align} \int_a^b w(x)f(x)dx&=\int_a^b w(x)\sum_{m=0}^M a_m x^m dx \nonumber\\ &= \sum_{m=0}^N a_m g_m \label{e4b} \end{align} $$

と書けます。ここで、

$$ \begin{align} g_m \equiv \int_a^b w(x) x^m dx \end{align} $$

と定義しました。一方で、式\eqref{e3}の右辺に式\eqref{e1}を代入すると、

$$ \begin{align} \sum_{i=1}^N f(x_i)\omega_i &= \sum_{i=1}^N \sum_{m=0}^M a_m x_i^m\omega_i \nonumber\\ &= \sum_{m=0}^M a_m \left(\sum_{i=1}^N \omega_i x_i^m\right) \label{e5b} \end{align} $$

と書けます。 もし、任意の$a_m$について式\eqref{e3}が厳密に一致するならば、$m=0, 1, \cdots, M$のそれぞれの係数について、

$$ \begin{align} g_m = \sum_{i=1}^N\omega_i x_i^m,~~(m=0, 1,\cdots, M) \end{align} $$

が成立していなければなりません。この式は全ての$m$について成り立つため、$M+1$個の方程式が立てられます。
一方で右辺に含まれる未知の変数は、$x_i, \omega_i$のそれぞれが$N$個あるので、系の自由度としては$2N$となります。

もし両辺を一致させるような係数の組み合わせがあるならば、方程式の数$M+1$よりも系の自由度$2N$の方が多くなければならないため、

$$ \begin{align} M+1\leq 2N~~\to~~M \leq 2N-1 \end{align} $$

が満たされる必要があります。関数の評価点数を最も少なくなるように考えると、上式の等号が成立した場合で所望の関数を厳密に積分することができます。

そのため、等号が成立した時を考えて、

$$ \begin{align} M = 2N-1 \end{align} $$

が最も無駄が無いと考えます。つまり、$N$個の点を利用すれば$2N-1$次の多項式を厳密に計算できる点の位置$x_i$と$w_i$が存在することが言えます。

以上で、被積分関数を参照する点の数$N$と、到達することができる被積分関数の多項式の次数$M$の関係が求められました。ここからは、具体的に$x_i, \omega_i$はどのような値を持つか、調べていきます。

2.2. 分点の位置 #

さて、前節で点数$N$と次数$M$の関係を見ましたが、具体的な値は何一つ決まっていません。
そのため、これから具体的な分点の位置$x_i$と重み$\omega_i$を求めます。

積分

$$ \begin{align} \label{ed0} \int_a^b w(x)f(x)dx \end{align} $$

を求めるために、離散的な$N$個の点$x=x_i$における$f(x)$の値を利用して積分を実行します。つまり、

$$ \begin{align} \label{ed1} (x_1, f_1), ~(x_2, f_2), \cdots , (x_N, f_N) \end{align} $$

の点が既知であると考えます。ここで、$f_i\equiv f(x_i)$と定義しました（この時点でまだ$x_i$は決まっていません）。$w(x)$については、使用する直交多項式の直交性から自動的に決まる既知の関数なので、いかなる点$x$についても既知です。

ここからかなり天下り的に考えていく項目が多くなります。出来る限り説明を試みますが、完全ではないことに注意してください。

まず、新たな関数$F(x), g(x)$を定義します。
$F(x)$は、式\eqref{ed1}の$N$個の点を全て通る$N-1$次多項式であり、一意に決まる関数として定義します。具体的には、ラグランジュ多項式¹を用いて

$$ \begin{align} \label{ed2} F(x)\equiv \sum_{k=1}^{N} f_k \left(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\right) \end{align} $$

と定義します。更に$f(x)$と$F(x)$の差の関数

$$ \begin{align}\label{ed3} g(x)\equiv w(x) [f(x)-F(x)] \end{align} $$

を定義します。ここでそれぞれの関数について次数を確認しておきます。

• $f(x)$: $2N-1$次多項式
• $F(x)$: $N-1$次多項式

ここでラグランジュ多項式$F(x)$は、点$x=x_i$において$F(x_i)=f(x_i),~~(i=1,2,\cdots,N)$となるので、$g(x)$は

$$ \begin{align} \label{ed4} g(x_i) = 0, \hspace{1em}(i=1,2,\cdots,N) \end{align} $$

の性質を持ちます。$x$の定義域で$g(x)$の積分を考えてみると

$$ \begin{align} \label{ed5} \int_{a}^b g(x)dx = \int_{a}^b w(x)f(x)dx-\int_{a}^b w(x)F(x)dx \end{align} $$

となります。$F(x)$は$M$次多項式なので、式\eqref{e3}のように離散的に表せるとすれば

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

となります。

今、右辺の$f(x)$は$2N-1$次であり、$F(x)$は$N$個の点で表現されているので、前節の議論から式\eqref{ed6}をゼロにする$(x_i, \omega_i)$の組み合わせが存在することが言えます。
もし式\eqref{ed6}をゼロにするように定めることができれば、当初の目的である$w(x)f(x)$の積分を厳密に実施することが言えます。

$g(x)$は式\eqref{ed3}で表現されており、また式\eqref{ed4}の性質を持つことが分かっているので、$N$個の$x_i$について、

$$ \begin{align} g(x)&\equiv w(x)\cdot(x-x_1)(x-x_2)\cdots(x-x_N)G(x)\nonumber\\ &=w(x)G(x)\prod_{i=1}^N (x-x_i) \end{align} $$

と書き換えることができます。ここで$G(x)$を新たに$N-1$次多項式として定義しました。
$(x-x_1)(x-x_2)\cdots(x-x_N)$は$N$次多項式であり、$w(x)$を除いた次数はあっています。

$g(x)$の積分がゼロになることを考えていたので積分は下記の通りになります。

$$ \begin{align} \label{c1} \int_{a}^b g(x)dx = \int_{a}^b w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)G(x) dx=0 \end{align} $$

$G(x)$は$N-1$次の任意の多項式なので、$N-1$次までの直交多項式$p_n(x)$の線形結合で展開することができます。つまり

$$ \begin{align}\label{c2} G(x)=w(x)\sum_{j=0}^{N-1}\alpha_j x^j=w(x) \sum_{j=0}^{N-1}\beta_j p_j(x) \end{align} $$

と書けるような$\alpha_j, \beta_j$が存在します。また、

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_N) \end{align} $$

は$N$次なので、$N$次以下の多項式の線形結合で書くことができます。

一般に$N$次の任意の多項式を展開する場合は、$N$次以下の全ての直交多項式の線形結合で書く必要がありますが、点$x_1, x_2, \cdots, x_N$を$p_N(x)=0$を満たす$N$個の点を選ぶとします。
このようにすると、$N$次だけの直交多項式だけで記述できて、

$$ \begin{align}\label{c3} (x-x_1)(x-x_2)\cdots(x-x_N)=c p_N(x) \end{align} $$

とする定数$c$を見つけることができます。定数$c$を見つけることもできますが、後々キャンセルされるので具体的な値は必要ありません。ここの条件から分点の位置が決まります。

式\eqref{c1}に式\eqref{c2}, \eqref{c3}を代入すると

$$ \begin{align} \int_a^b g(x)dx &=\int_a^b dx w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)\cdot G(x) \nonumber\\ &=c\int_a^b dx w(x)p_N(x)\cdot \sum_{j=0}^{N-1}\beta_j p_j(x) \nonumber\\ &=c\sum_{j=0}^{N-1}\beta_j\int_a^b dx w(x) p_N(x) p_j(x) \nonumber\\ &=c h_N \sum_{j=0}^{N-1}\beta_j\delta_{N, j} \nonumber\\ &=0 \label{c4} \end{align} $$

計算できます。ここで$\delta_{i, j}$はクロネッカーのデルタを意味します。
また、$h_n$は直交多項式の直交性の右辺に現れる下記の量を意味します。

$$ \begin{align}\label{qh} \int_a^b dx w(x) p_n(x) p_m(x)=h_n \delta_{n,m} \end{align} $$

式\eqref{c4}の結果は、点$x_1, x_2, \cdots, x_N$を$p_N(x)=0$を満たす$N$個の点を選べば自動的に式\eqref{c1}が満たされる、 すなわち$\int_a^b g(x)dx = 0$になるので、この$N$個の点で式\eqref{ed0}の厳密な積分ができることを意味します。

以上より、$N$次直交多項式のゼロ点をガウス求積法における分点の位置として選ぶと

$$ \begin{gather} \label{e4} \int_a^b g(x)dx = \int_a^b w(x)f(x)dx-\sum_{i=1}^N F(x_i)\omega_i =0 \\ p_N(x)=0 \to x=x_i ,~~(i=1, 2, \cdots, N) \end{gather} $$

を満たすようにできることが分かりました。

ここまでの議論で、分点の位置$x_i$が求まりました。残るは重み$\omega_i$を具体的に決めることです。

2.3. 重み #

これまでの議論から、重み関数$w(x)$と$2N-1$次の多項式$f(x)$の積の積分は、$N$個の点を用いて

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

を厳密に積分する分点の位置と重み$x_i, \omega_i$の組があることを示しました。その後、分点の位置$x_i$は$N$次直交多項式のゼロ点、すなわち

$$ \begin{align} p_{N}(x)=0\to x=x_i ,~~(i=1, 2, \cdots, N) \end{align} $$

を満たす$N$個のゼロ点であることが分かりました。次は重み$\omega_i$を決めていきます。

式\eqref{e4}が成り立つことをみたので、$2N-1$次の$f(x)$を$N$次直交多項式のゼロ点を通る多項式の積分に置き換えることができます。 式\eqref{ed5}がゼロであることと、式\eqref{ed2}を用いて

$$ \begin{align} \int_a^bw(x)f(x)dx&= \int_a^b w(x)F(x)dx\nonumber\\ &= \int_a^b w(x)\sum_{k=1}^{N} f_k \Biggl(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &= \sum_{k=1}^{N} f_k \int_a^b w(x)\Biggl(\prod_{\substack{i=1,\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &\equiv \sum_{k=1}^{N}f_k \omega_k \end{align} $$

と書き変えていくことができます。つまり重みは

$$ \begin{align} \omega_k=\int_{a}^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx\label{omega0} \end{align} $$

と表されていることが分かります。位置$x_i$は既に決まっているので、式\eqref{omega0}は一意に決まります。そのためこのままでも良いのですが、積分や積和が入り組んでおり、使いにくいので式変形できないか試みていきます。

式\eqref{c3}から、

$$ \begin{align}\label{u1} (x-x_1)(x-x_2)\cdots(x-x_N)=\frac{1}{k_N} p_N(x) \end{align} $$

であることが分かっています。特に$k$番目の式についてまとめると、

$$ \begin{align}\label{uu1} \prod_{\substack{i=1\\ i\ne k}}^N (x-x_i)= \frac{cp_N(x)}{x-x_k} \end{align} $$

と書けることが分かります。また、式\eqref{u1}の両辺を微分すると

$$ \begin{align} cp’_N(x) &=\frac{\partial}{\partial x}\Bigl[(x-x_1)(x-x_2)\cdots(x-x_N)\Bigr]\nonumber\\ &=\Bigl[(x-x_2)(x-x_3)\cdots(x-x_N)\nonumber\\ &\hspace{3em}+(x-x_1)(x-x_3)\cdots(x-x_N)\nonumber\\ &\hspace{3em}+(x-x_1)(x-x_2)\cdots(x-x_{N-1})\Bigr] \nonumber\\ &=\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

となります。特に$x=x_k$ならば$j=k$の場合しか項が残らないので、和の記号は消えてしまうため

$$ \begin{align}\label{u2} cp’_N(x_k) = \prod_{\substack{i=1\\ i\ne k}}^{N} (x_k-x_i) \end{align} $$

となります。式\eqref{uu1}と\eqref{u2}を組み合わせると、

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i} &= \frac{\prod_{\substack{i=1\\ i\ne k}}^{N}{x-x_i}}{\prod_{\substack{i=1\\ i\ne k}}^{N}{x_k-x_i}}\nonumber\\ &= \frac{cp_N(x)}{x-x_k}\cdot\frac{1}{cp’_N(x_k)} \nonumber\\ &=\frac{p_N(x)}{(x-x_k)p’_N(x_k)} \end{align} $$

と計算を進めることができます。なので、重み$\omega_k$は

$$ \begin{align} \omega_k&=\int_a^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \nonumber\\ &=\int_a^b w(x)\frac{p_N(x)}{(x-x_k)p’_N(x_k)} dx \label{w1} \end{align} $$

と変形できます。

積分\eqref{w1}は、多項式の性質を用いることで計算を進めることができます。

天下り的ですが、右辺の形はクリストッフェル＝ダルブーの定理 (Christoffel–Darboux theorem) から導けそうな形であることを当たりをつけて計算します。 直交多項式（ルジャンドル多項式やラゲール多項式など）に存在するクリストッフェル＝ダルブーの定理とは、$n$次の直交多項式$p_n(x)$に存在する下記の性質です。

$$ \begin{align}\label{CDt} \sum_{j=0}^N \frac{p_j(x)p_j(y)}{h_j} = \frac{k_N}{h_N k_{N+1}}\frac{p_{N+1}(x)p_{N}(y)-p_{N}(x)p_{N+1}(y)}{x-y} \end{align} $$

ここで、$k_n$は直交多項式$p_n(x)$の$x^n$の係数を意味しています。つまり

$$ \begin{align} p_n(x)=k_n x^n+\cdots \end{align} $$

で表されるときの係数です。また$h_n$は式\eqref{qh}の通り、直交性から導かれる係数です。

式\eqref{CDt}の$y$がちょうど$x_k$の場合について考えると、関係式

$$ \begin{align} \sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} &= \frac{k_{N-1}}{h_{N-1} k_{N}}\frac{p_{N}(x)p_{N-1}(x_k)-p_{N-1}(x)p_{N}(x_k)}{x-x_k} \nonumber\\ &= \frac{k_{N-1}}{h_{N-1} k_{N}}\frac{p_{N}(x)p_{N-1}(x_k)}{x-x_k} \label{PNx} \end{align} $$

を得ます。ここで分点の位置$x_k$は$p_{N}(x_k)=0$を満たすように決めているので、これを利用しました。 続いて移項して関係式

$$ \begin{align} \frac{p_{N}(x)}{x-x_k}&=\frac{k_{N}h_{N-1} }{k_{N-1}} \frac{1}{p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} \label{w2} \end{align} $$

を得ます。

式\eqref{w2}を重みの式\eqref{w1}に代入して、

$$ \begin{align} \omega_k&=\int_{a}^b w(x)\frac{p_N(x)}{(x-x_k)p’_N(x_k)} dx\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\int_{a}^b dx w(x)\sum_{j=0}^{N-1} \frac{p_j(x)p_j(x_k)}{h_j} \nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \int_{a}^b dx w(x) p_j(x)\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \frac{1}{c_0}\int_{a}^b dx w(x) p_j(x)p_0(x)\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\sum_{j=0}^{N-1} \frac{p_j(x_k)}{h_j}\cdot \frac{h_0}{c_0}\delta_{j,0}\nonumber\\ &= \frac{k_{N}h_{N-1}}{k_{N-1} p’_N(x_k)p_{N-1}(x_k)}\label{w3} \\ &= \left[\sum_{n=0}^{N-1} \frac{p_n^2(x_k)}{h_n}\right]^{-1}\label{w3x} \end{align} $$

ここで、直交多項式のゼロ次は定数であることを使用しました($p_0(x)=c_0$)。 また、式\eqref{w3x}における最後の変形では、クリストッフェル＝ダルブーの定理における$y\to x$の極限を考えた場合の関係式

$$ \begin{align} \sum_{n=0}^{N-1} \frac{[p_n(x)]^2}{h_n} = \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} $$

と、分点の性質$p_N(x_k)=0$であることを用いて導くことができます。

実際に重みを計算する際は、式\eqref{w3}または\eqref{w3x}のどちらかを、状況によって計算しやすい方で実施します。

3. ガウス＝ロバート求積法 #

ガウス＝ロバート求積法は、端点で分点を持つように設計されたガウス求積法です。

式\eqref{c1}までは同じですが、ここから分点の位置を決める際に点$x=a,b$が選ばれるようにします。

基本的にガウス＝ロバート求積法と言った時はガウス＝ルジャンドル求積法（区間$x=[-1,1]$）を基にして両端を分点に持つようにしたものを言いますが、ルジャンドル多項式に限らず、ヤコビ多項式(Jacobi polynomials)に対して同様の考えが可能です²³。

本稿では、基本に則ってルジャンドル多項式に基づいた求積法を考えます。

3.1. 分点の制約 #

$$ \begin{align} \int_{a}^b g(x)dx = \int_{a}^b w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_N)G(x) dx=0 \end{align} $$

これまでの導出では$(x-x_1)\cdots(x-x_N)$を多項式のゼロ点にとることによって式\eqref{c1}を満たすようにしていましたが、ここを変えます。

つまり、

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_N) \end{align} $$

の内、2つの点が端点であるとして

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_{N-3})(x-a)(x-b) \end{align} $$

として考えます。これを$N-2$次の多項式$q_{N-2}(x)$と$(x-a)(x-b)$の積を用いて

$$ \begin{align} (x-x_1)(x-x_2)\cdots(x-x_{N-3})&(x-a)(x-b)\nonumber \\ &=c q_{N-2}(x)(x-a)(x-b)\label{k3a} \end{align} $$

と書き変えられるとします。ここで$c$は適当な定数です。

$q_{N-2}(x)$は$N-2$次多項式です。これまでの議論で端点で$q_{N}(x)$がゼロにならなければならないという制約を課したため、両端の分点位置の自由度 “2” を失っています。 そのため、厳密な積分が可能な積分の次数$M$と分点数$N$の関係は、

$$ M=(2N-1) - 2 = 2N-3 $$

と変更されます。両端の位置を固定しましたが、両端における重みの自由度は減らさないので2引くだけで良いです。

この次数の減少に伴い、これまで$G(x)$は$N-1$次の任意の多項式と扱ってきましたが、これを諦めて$N-3$次の任意の多項式に次数を落とします。
つまり

$$ \begin{align} G(x)=w(x)\sum_{j=0}^{N-3}\alpha_j x^j=w(x) \sum_{j=0}^{N-3}\beta_j q_j(x) \end{align} $$

とします。すると、

$$ \begin{align} \int_a^b g(x)dx \label{k2r} &=\int_a^b dx w(x)\cdot (x-x_1)(x-x_2)\cdots(x-x_{N-3})(x-a)(x-b)\cdot G(x) \\ &=-c\int_a^b dx [w(x)(x-a)(b-x)] q_{N-2}(x)\cdot \sum_{j=0}^{N-3}\beta_j q_j(x) \nonumber\\ &=-c\sum_{j=0}^{N-3}\beta_j\int_a^b dx w’(x) q_{N-2}(x) q_j(x) \\ \end{align} $$

と書き変えられます。ここで新たな重み関数

$$ \begin{align} w’(x) = w(x)(x-a)(b-x) \end{align} $$

と置きました。もし$q_n(x)$が直交性

$$ \begin{align} \label{k4} \int_a^b dx w’(x) q_n(x) q_m(x)=h_n\delta_{n,m}, \end{align} $$

を満たすように決められていれば、

$$ \begin{align} \int_a^b g(x)dx&=c\sum_{j=0}^{N-3} h_n \beta_j \delta_{N-2,j}\\ &=0 \end{align} $$

となり、ガウス求積法の性質を満たすようにできます。

さて、以上のように式\eqref{k4}を決めましたが、直交性\eqref{k4}を満たす関数$q_{N-2}(x)$を見つけることはできるのでしょうか？ $q_{N-2}(x)$の条件として下記を性質を満たしていなければなりません。

• $q_{N-2}(x)$: $N-2$次多項式
• $q_{N-2}(a)<\infty$, $q_{N-2}(b)<\infty$
• 直交性\eqref{k4}を満たす

以降、複雑さを回避するためにルジャンドル多項式に基づくガウス＝ロバート求積法の場合のみを考えていきます⁴。

3.2. 分点の位置 #

ルジャンドル多項式の場合、$a=-1, b=1$です。 $N$次のルジャンドル多項式$P_N(x)$が次数2だけ減って直交多項式$P_{N-2}(x)$となるので、減った次数分情報を用いて両端でゼロになる条件を満たすように多項式を設計します。

つまり通常のガウス求積法では

$$ \begin{align} P_N(x)=0 \hspace{1em}\to\hspace{1em} x=x_j,\hspace{1em}(j=1,\cdots,N) \end{align} $$

を満たす$N$個のゼロ点を探してそれを分点としていましたが、これを下記のようにします。

$$ \begin{align} &P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)=0, \nonumber \\ &\hspace{3em}\to\hspace{1em} x=-1, 1, x_j,\hspace{1em}(j=1,\cdots,N-2) \label{k6} \end{align} $$

ここで、定数$u, v$は端点$x=\pm 1$で与式がゼロになるという条件から決定します。

このようにした後、式\eqref{k3a}より、$N-2$次の多項式を

$$ \begin{align} \label{k7} q_{N-2}(x)=\frac{P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)}{(1+x)(1-x)},~~q_{N-2}(\pm 1)<\infty \end{align} $$

として決めます。これがガウス＝ロバート求積法で考える多項式となります。 念のためですが、式\eqref{k7}において、分子の根($x=\pm1$)を分母で割っているのでこの結果は多項式となり、$q_{N-2}(x)$も多項式となります（割った分だけ次数が減少するのでこの操作は “減次” と呼ばれます⁵）。

最後に、$q_{N-2}(x)$は直交性\eqref{k4}が満たされているのか確認します。
まず、$q_{N-2}(x)$の明らかな形を求めていきましょう。$x=\pm 1$で式\eqref{k6}がゼロになるという条件から、

$$ \begin{align} x=1 &:& 1 + u + v &= 0 \\ x=-1 &:& -1 + u - v &= 0 \end{align} $$

より、$u=0, v=-1$が導かれます。
ここで、ルジャンドル多項式の性質$P_n(1)=1, P_n(-1)=(-1)^n$を用いました。

以上から、

$$ \begin{align} \label{k8} q_{N-2}(x)=\frac{P_{N}(x)- P_{N-2}(x)}{(1-x^2)}=-\frac{2N-1}{N(N-1)}P’_{N-1}(x) \end{align} $$

となります。ここではより一般的に考えるために添え字をずらして

$$ \begin{align} \label{k8a} q_m(x)&=\frac{P_{m+2}(x)- P_{m}(x)}{(1-x^2)} \\ &=-\frac{2m+3}{(m+2)(m+1)}P’_{m+1}(x) \end{align} $$

として考えていきます。すると直交性は、

$$ \begin{align} \label{k9} &\int_{-1}^{1} (1-x^2)q_n(x)q_m(x)dx \nonumber \\ &=\int_{-1}^{1} \left[P_{n+2}(x)- P_{n}(x)\right]\cdot \left[-\frac{2m+3}{(m+2)(m+1)}P^{\prime}_{m+1}(x)\right]dx \nonumber \\ &=-\frac{2m+3}{(m+2)(m+1)}\Biggl\{ \Bigl[\bigl(P_{n+2}(x)- P_{n}(x)\bigr)P_{m+1}(x)\Bigr]_{x=-1}^{x=1} \nonumber\\ &\hspace{9em} -\int_{-1}^{1}\frac{d}{dx}\Bigl[P_{n+2}(x)- P_{n}(x)\Bigr] \cdot P_{m+1}(x) dx \Biggr\} \nonumber\\ &=\frac{(2m+3)(2n+3)}{(m+2)(m+1)}\int_{-1}^{1}P_{n+1}(x) \cdot P_{m+1}(x) dx \nonumber\\ &=\frac{(2m+3)(2n+3)}{(m+2)(m+1)}\frac{2}{2n+3}\delta_{n,m} \nonumber\\ &=\frac{2(2n+3)}{(n+2)(n+1)}\delta_{n,m} \end{align} $$

となります。

ここでルジャンドル多項式の関係式（Appendix \eqref{ap2}）を利用しました。

3.3. ガウス＝ロバート求積法で使用する多項式の性質 #

以上まで出てきたガウス＝ロバート求積法で使用する多項式$q_n(x)$の性質をまとめておきます。

• ルジャンドル多項式$P_n(x)$との関係性

$$ \begin{gather} q_n(x)=\frac{P_{n+2}(x)- P_{n}(x)}{(1-x^2)} =-\frac{2n+3}{(n+2)(n+1)}P’_{n+1}(x) \\ \int_{-1}^{1} (1-x^2)q_n(x)q_m(x)dx = \frac{2(2n+3)}{(n+2)(n+1)}\delta_{n,m} \end{gather} $$

• 対称性

$$ \begin{gather} q_{n}(-x)=(-1)^n q_n(x) \end{gather} $$

• 特別な値

$$ \begin{gather} q_{n}(1)=-\frac{2n+3}{2} \end{gather} $$

$$ \begin{gather} \frac{d}{dx}q_{n}(x)\Bigr|_{x=1}=-\frac{n(n+3)(2n+3)}{8} \end{gather} $$

※Appendix参照。

• 三項漸化式

$$ \begin{gather} \frac{n+3}{2n+5}q_{n+1}(x)=xq_n(x)-\frac{n}{2n+1}q_{n-1}(x) \end{gather} $$

• 微分

$$ \begin{gather} (1-x^2)q’_{n}(x)=-nxq_n(x)+\frac{n(2n+3)}{2n+1}q_{n-1}(x) \end{gather} $$

• 最高次の係数

$$ \begin{gather} q_n(x)=-\frac{2n+3}{n+2} \left(\frac{(2n+2)!}{2^{n+1}((n+1)!)^2}\right)x^n+\cdots \end{gather} $$

• 具体的な$q_n(x)$

3.4. 重み #

重みを求める際は、関数$c q_{N-2}(x)(x-a)(x-b)$について計算します。よって重みは

$$ \begin{align} \omega_k=\int_{a}^b w(x)\Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx \end{align} $$

で計算できます。便宜上、$x_{N-1}=a, x_{N}=b$と置きました。この段階では$w(x)=1$ですが、間違えないように下記のように消して置きます。

$$ \begin{align} \omega_k=\int_{a}^b \Biggl(\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i}\Biggr) dx\label{omega1} \end{align} $$

式\eqref{c3}から、

$$ \begin{align}\label{u1x} (x-x_1)(x-x_2)\cdots(x-x_N)=c q_{N-2}(x) (x-x_{N-1})(x-x_{N}) \end{align} $$

と書き変えることができます。これまでと同様に進めます。
$k$番目の式についてまとめると、

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N} (x-x_i)= \left\{ \begin{aligned} &\frac{cq_{N-2}(x)}{x-x_k}(x-x_{N-1})(x-x_{N}),&&(k=1,\cdots,N-2) \\ &cq_{N-2}(x)(x-x_{N}),&&(k=N-1) \\ &cq_{N-2}(x)(x-x_{N-1}),&&(k=N) \\ \end{aligned} \right.\label{uu1x} \end{align} $$

と書けることが分かります。
一方で式\eqref{u1x}の両辺を微分しますと

$$ \begin{align} (d/dx \text{ of left-hand side} ) =\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

より、

$$ \begin{align} &cq’_{N-2}(x)(x-x_{N-1})(x-x_{N}) + cq_{N-2}(x)(x-x_{N}) + cq_{N-2}(x)(x-x_{N-1}) \nonumber\\ &\hspace{18em}=\sum_{j=1}^N\prod_{\substack{i=1\\ i\ne j}}^{N}(x-x_i) \end{align} $$

となります。特に$x=x_k$ならば$j=k$の場合しか項が残らないので、和の記号は消えてしまうため

$$ \begin{align} \prod_{\substack{i=1\\ i\ne k}}^{N} (x_k-x_i) =\left\{ \begin{aligned} &cq’_{N-2}(x_k)(x_k-x_{N-1})(x_k-x_{N}), &&(k=1,\cdots,N-2) \\ &cq_{N-2}(x_k)(x_k-x_{N}),&&(k=N-1) \\ &cq_{N-2}(x_k)(x_k-x_{N-1}),&&(k=N) \end{aligned} \right.\label{u2x} \end{align} $$

となります。式\eqref{uu1x}と\eqref{u2x}を組み合わせると、

$$ \begin{align} &\prod_{\substack{i=1\\ i\ne k}}^{N}\frac{x-x_i}{x_k-x_i} \nonumber\\ &=\left\{ \begin{aligned} &\frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})}, &&(k=1,\cdots,N-2) \\ &\frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})},&&(k=N-1) \\ &\frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})}{(x_k-x_{N-1})},&&(k=N) \\ \end{aligned} \right. \end{align} $$

と計算をできます。なので、重み$\omega_k$は

• $k=1,2,\cdots,N-2$の時

$$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})} dx \end{align} $$

• $k=N-1$の時 ($x_{N-1}=-1$)

$$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})} dx \end{align} $$

• $k=N$の時 ($x_{N}=1$)

$$ \begin{align} \omega_k=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})}{(x_k-x_{N-1})} dx \end{align} $$

で求めることができます。また$q_n(x)$は直交多項式であるので、クリストッフェル＝ダルブーの定理より下記の関係式が導かれます。

$$ \begin{align} \frac{q_{N-2}(x)}{x-x_k}=\frac{k_{N-2}h_{N-3}}{k_{N-3}} \frac{1}{q_{N-3}(x_k)}\sum_{j=0}^{N-3} \frac{q_j(x)q_j(x_k)}{h_j} \label{w2x} \end{align} $$

これを利用して、それぞれの重みについて計算していくと下記のようになります。

• $k=1,2,\cdots,N-2$の時

$$ \begin{align} \omega_k&=\int_a^b \frac{q_{N-2}(x)}{(x-x_k)q’_{N-2}(x_k)}\cdot\frac{(x-x_{N-1})(x-x_{N})}{(x_k-x_{N-1})(x_k-x_{N})} dx \nonumber\\ &=\frac{1}{(1-x_k^2)q’_{N-2}(x_k)}\int_{-1}^1 \frac{q_{N-2}(x)}{(x-x_k)}\cdot(1-x^2) dx \nonumber\\ &=\frac{k_{N-2}h_{N-3}}{k_{N-3}}\frac{1}{(1-x_k^2)q’_{N-2}(x_k)q_{N-3}(x_k)} \sum_{j=0}^{N-3} \frac{q_j(x_k)}{h_j}\int_{-1}^1 (1-x^2)q_j(x)dx \nonumber\\ &=\frac{k_{N-2}h_{N-3}}{k_{N-3}}\frac{1}{(1-x^2_k)q’_{N-2}(x_k)q_{N-3}(x_k)} \\ &=\frac{2(2N-1)(2N-3)}{N(N-1)(N-2)}\frac{1}{(1-x^2_i)q’_{N-2}(x_i)q_{N-3}(x_i)}\nonumber\\ &= \frac{2}{N(N-1)P^2_{N-1}(x_i)} \end{align} $$

• $k=N-1$の時 ($x_{N-1}=-1$)

$$ \begin{align} \omega_k&=\int_a^b \frac{q_{N-2}(x)}{q_{N-2}(x_k)}\cdot\frac{(x-x_{N})}{(x_k-x_{N})} dx \nonumber\\ &=\frac{1}{-2q_{N-2}(-1)}\int_{-1}^1 q_{N-2}(x)\cdot (x-1) dx \nonumber\\ &=\frac{1}{2 \cdot(-c_{N-2}P’_{N-1}(-1))}(-c_{N-2})\int_{-1}^1 (1-x)P’_{N-1}(x) dx \nonumber\\ &=\frac{1}{2 \cdot (-1)^N \frac{N(N-1)}{2}}\left(\Bigl[(1-x)P_{N-1}(x)\Bigr]^1_{-1} +\int_{-1}^1 P_{N-1}(x) dx \right) \nonumber\\ &=\frac{2}{N(N-1)}\frac{1}{2 \cdot (-1)^N }\bigl(-2P_{N-1}(-1) \bigr) \\ &=\frac{2}{N(N-1)} \end{align} $$

• $k=N$の時 ($x_{N}=1$)

$$ \begin{align} \omega_k=\frac{2}{N(N-1)} \end{align} $$

$k=N$の時、素直に計算しても良いですが、対称性があるので計算は不要です。 途中、$c_n$を$q_n(x)=-c_nP’_{n+1}(x)$を満たすように置いています。

3.5. まとめ #

ガウス=ロバート求積法は、$N (\geq 2)$点の離散的な関数値を用いて、$2N-3$次の多項式$f(x)$の積分を厳密に評価することが可能です。

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

• 分点の位置

区間 $x=[-1,1]$ で定義された $N-2$次の多項式

$$ \begin{align} q_{N-2}(x)=\frac{P_{N}(x)+ u P_{N-1}(x)+v P_{N-2}(x)}{(1+x)(1-x)}=0,~~q_{N-2}(\pm 1)<\infty \end{align} $$

が、$q_{N-2}(x)=0$ になるような$N-2$個の点と、$x=\pm 1$の2個の点の、合計$N$点が分点の位置となります。 $x$の昇順に番号付けをすると、

$$ \begin{gather} x_1=-1, x_N=1, \\ q_{N-2}(x)=0\to x=x_{i},~~(i=2,3,\cdots, N-1) \end{gather} $$

となります。

• 重み

$$ \begin{align} \omega_1 &= \omega_N = \frac{2}{N(N-1)}, \\ \omega_i &= \frac{2}{N(N-1)P^2_{N-1}(x_i)},~~(i=2,3,\cdots, N-1) \end{align} $$

※$N=2$の場合は、$x_1=-1, x_2=1$であり、$\omega_1=1, \omega_2=1$となるため、台形則（1次式ならば厳密）となります。

4. Appendix #

4.1. ルジャンドル多項式の関係式 #

Legendre polynomials -wikipedia

$$ \begin{align}\label{ap2} (2n+1)P_n(x)=\frac{d}{dx}\left(P_{n+1}(x)-P_{n-1}(x)\right) \end{align} $$

$$ \begin{align}\label{ap3} nP’_{n+1}(x)-(2n+1)xP’_{n}+(n+1)P’_{n-1}(x)=0 \end{align} $$

ルジャンドル多項式の微分は下記の通り書かれます⁶。

$$ \begin{align}\label{ap5} \frac{d^m}{dx^m}P_n(x) = 1 \cdot 3 \cdots (2m-1)C^{(m+\frac{1}{2})}_{n-m}(x),~~(m\leq n) \end{align} $$

ここで、$C^{(\alpha)}_{n}(x)$は Gegenbauer多項式 (またはUltraspherical多項式)で、$x=1$と対称性は下記の通りになります⁷。

$$ \begin{align}\label{ap6} C^{(\alpha)}_{n}(1)&=\binom{n+2\alpha-1}{n}\\ C^{(\alpha)}_{n}(-x)&=(-1)^n C^{(\alpha)}_{n}(x) \end{align} $$

以上より、

$$ \begin{align}\label{ap7} \frac{d^m}{dx^m}P_n(x)\Bigr|_{x=1} = 1 \cdot 3 \cdots (2m-1)\cdot \binom{n+m}{n-m},~~(m\leq n) \end{align} $$

となります。

5. 参考文献 #