ガウス求積法は、関数の数値積分を行う方法です。
N個の点で関数を参照すると、2N-1次以下の多項式で表される関数を厳密に数値積分することができる方法です。
ガウス求積法の凄さを少し述べます。
通常、3個の任意の点における値があれば、その3点を結ぶ2次関数を作ることができるので、うまい選び方をすればN個の点でN-1次の関数を厳密に表現できそうであり、それが厳密に積分できてもそこまで不思議ではありません。
ですが、ガウス求積法では分点の位置にも制限を加えることで積分の次数を倍にすることができます。
分点の位置を制限するにもかかわらず、高々二倍になるだけか…とも思われますが、これは同じ精度に達するまでに関数の評価点数を半分に減らせることと同じです。
評価点数を少なくできるというこの性質は、多次元の計算において顕著に現れます。なぜなら、1次元で半分になるので2次元では4分の1になります。
多少無理してもガウス求積法を使う価値があります。
科学技術計算ではガウス求積法を使うことをお勧めします。ガウス求積法を推奨します。
ガウス求積法を使いましょう（圧）。
本稿のPDF版はこちらのリンクからどうぞ。
https://slpr.sakura.ne.jp/qp/wp-content/uploads/2021/03/Gauss_legendre_quadrature_derivation_20210904.pdf
導出手順
導出は以下の順番で行います。
- 点数と次数の関係
- 具体的な位置と重みの計算
『点数と次数の関係』では、M個の離散的な点を用いてN次多項式を厳密に積分できるとき、NとMの関係を導出します。
『具体的な位置と重みの計算』では、具体的に、M点の位置とその重みを導出します。
ガウス求積法を利用する方は、1. が分かればよいと思います。
ガウス求積法を研究する方は、2. まで分かればよいと思います。
また、本稿ではガウス求積法の中で最も基本的な
ガウス＝ルジャンドル求積法
について記述します。他の方法も、ガウス＝ルジャンドル求積法の導出を理解できれば、おのずと分かるでしょう。
先に結論をまとめます。
\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、
を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、N個の\(x_i\)とN個の\(\omega_i\)は
から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。
点数と次数の関係
N次多項式を有限区間\([-1,1]\)で、M個の点を用いて数値的に厳密に積分することを考えます。つまり、関数\(f(x)\)がN個の既知の係数\(a_n\)を用いて
と書けている場合に、積分
を考えます。この積分を離散的な点だけを用いて計算します。つまり、適当なM個の点を用いて
となるような右辺の\(x_i, \omega_i\)の組み合わせを求めることが目標です。
式(3)をできるだけ計算していきましょう。
まず、式(3)の左辺を積分すると、
となります。
続いて、式(3)の右辺を計算すると、
と書けます。ここで、\(\omega_i\)は重みであり、適当な係数です。もしも任意の\(a_n\)について、式(3) の右辺と左辺が厳密に一致するならば、\(n=0, 1,\cdots, N\)について、
が成立します。この式は合計で\(N+1\)個の方程式の数があります。今、右辺の未知の係数は\(\omega_i, x_i, (i=1,\cdots, M)\)なので、右辺の自由度は\(2M\)個あります。
つまり、両辺を一致させるような係数の組み合わせがあるならば、
である必要があります。よって、もっとも小さくなるような\(M\)であるならば、上式の等号が成立します。できるだけ少ない点数で計算したいので、等しい場合を採用して
となります。つまり、M個の点を利用すれば\(2M-1\)次の多項式を厳密に計算できる点の位置\(x_i\)と\(w_i\)が存在します。
具体的な位置と重みの計算
具体的な分点の位置
積分
を求めるために、離散的なN個の点\(x=x_i\)における\(f(x)\)の値だけを利用します。つまり、
の点は既知とします。ここで、\(f_i\equiv f(x_i)\)と定義しました（この時点でまだ\(x_i\)は決まっていません）。
新たな関数\(f(x)\)を定義します。\(f(x)\)は、式(12)のN個の点を全て通る\(N-1\)次多項式であり、一意に決まります。具体的には、ラグランジュ多項式を用いて
と書けます。ここで、更に\(f(x)\)と\(f(x)\)の差の関数\(g(x)\)
を定義します。\(x\)の次数を確認しておくと、\(g(x), f(x)\)は\(2N-1\)次、\(f(x)\)は\(N-1\)次多項式です。
ここで、\(f(x_i)=F(x_i),~(\mbox{for all } i)\)なので、\(g(x)\)は
の性質を持ちます。積分を考えると
であり、\(f(x)\)を離散的に表すと、
です。今、\(f(x)\)は\(2N-1\)次であり、\(f(x)\)はN個の点で表現されているので、前節の議論から式(17)をゼロにするような\((x_i, \omega_i)\)の組み合わせが存在することが分かります。
今, \(g(x)\)の満たす性質として式(15)が分かっていますので、N個の\(x_i\)について、
と書き換えることができます。
ここで、\((x-x_1)(x-x_2)\cdots(x-x_N)\)はN次多項式で、\(G(x)\)は\(N-1\)次多項式です。つまり、
を満たすN個の\(x=x_i\)を探すことができれば嬉しいわけです。ここで、\(G(x)\)は\(N-1\)次多項式なので、\(N-1\)次のlegendre多項式を用いた基底関数で展開可能です。なので、
と書くことができます。また、
はN次なので、N次多項式の関数で書くことができます。
つまりN次多項式としてN次legendre多項式\(P_n(x)\)を選ぶことにします。
単にN次であるだけだと、N次以下の全てのlegendre多項式の和で書けなければなりませんが、点\(x_1, x_2, \cdots, x_N\)を\(P_N(x)=0\)を満たすN個の点を作為的に選ぶことにしましょう。
このように選ぶと、N次以下の和ではなく、N次のlegendre多項式だけで書くことができて、
とする定数\(c\)を見つけることができます。
これからの議論において定数\(c\)を具体的に決める必要はありませんが明確に決めたい場合
の関係から決まります。
式(19)に式(20), 式(22)を代入すると
となります。これはつまり、N次legendre多項式のゼロ点を求積法におけるゼロ点として選ぶと
を満たすように\(g(x)=0\)にできるということです。
ここまでの議論で、分点の位置\(x_i\)が求まりました。
残るは、まだ決定していない重み\(\omega_i\)を決めることです。
具体的な重み
上記までの議論から、\(2N-1\)次の多項式で記述される\(f(x)\)の積分は、N個の点を用いて
を厳密に積分する分点の位置と重み\(x_i, \omega_i\)の組があることを示し、分点の位置\(x_i\)は、N次legendre多項式のゼロ点、すなわち
を満たすN個のゼロ点であることが分かりました。次は重みを決めます。
式(16)より\(2N-1\)次の\(f(x)\)をN次のlegendre多項式のゼロ点を通る多項式の積分に置き換えできます。
つまり、式(13)を用いて
のように置き換えできます。つまり、重みは
ですので、これを求めることが問題となります。
積和は書き換えがやりにくいですので、意図的に消していきましょう。
式(22)から、
分かっています。意図的にk番目の式についてまとめると、
となります。また、式(35)の両辺を微分すると
となります。特に、\(x=x_k\)ならば、
式(36)と式(40)を組み合わせると、
なので、重み\(\omega_k\)は
となります。この積分を計算するために多項式の性質を用いて計算を進めましょう。
右辺の形は、クリストッフェル＝ダルブーの定理から導けそうな形であることを当たりをつけて計算します（私のような凡人には思い付きません。天下り的に使います。導けるからそれで良いです）。
直交多項式（Legendre多項式やLaguerre多項式など）にはクリストッフェル＝ダルブーの定理 (Christoffel–Darboux theorem) があります。それは、n次の直交多項式を\(P_n(x)\)と書くとき、\(P_n(x)\)には
という性質であります。ここで、\(k_n\)は、直交多項式\(P_n(x)\)の\(x^n\)の係数で、
で表されるときの係数です。この係数を具体的に求めてみましょう。そのために、
legendre多項式の性質の一つである三項漸化式
を用います。xの次数で両辺を比較すると
なので、次数を比較すれば最高次の係数について
の関係があります。よって、
を得ます。また、\(h_N\)は直交性から導かれる係数で、
を満たします。legendre多項式の場合、\(w(x)=1\)であり、\(h_n=\frac{2}{2n+1}\)です。
重みの計算に使いたいため、クリストッフェル＝ダルブーの定理に対し、\(N\to N-1\)にして、legendre多項式の場合に適用、さらに\(y=x_k\)の場合を考えて、
より、
を得ます。ここで\(P_{N}(x_k)=0\)を利用しました。
式(43)に式(56)を代入して、
となります。\(P_n(x)\)の微分は次数が一つ下がるので、\(P_{n-1}(x)\)と同じ次数になります。つまり、両者に何か関係がありそうですので頑張って見つけてみます。
微分に関する三項漸化式
を用います。特に\(n=N,~x=x_k\)の場合、
の関係があります。代入すると重みは
と求めることができます。
まとめ
\(2N-1\)次以下の多項式\(f(x)\)の積分を、N個の点だけを用いて厳密に計算できます。
つまり、
を厳密に計算する点の位置\(x_i\)と\(w_i\)が存在し、\(x_i, \omega_i\)は
から求めることができます。ここで、\(P_n(x)\)はn次ルジャンドル多項式です。
参考文献など
恐らく、ガウスが発表した論文がこれではないかと思います。ドイツ語のwikipedia（https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur）にリンクがありました。
Methodus nova integralium valores per approximationem inveniendi. In: Comm. Soc. Sci. Göttingen Math. Band 3, 1815, S. 29–76
千葉豪著, Gauss-Legendre求積公式の導出, (2016)
http://roko.eng.hokudai.ac.jp/studentadm/chiba_data/PnSn/gl5.pdf
Abramowitz and Stegun, “Handbook of Mathematical Functions.”
https://www.math.ubc.ca/~cbm/aands/page_887.htm
L. Bos et al, “An Orthogonality Property of the Legendre Polynomials”, Constr Approx (2016)
https://math.indiana.edu/promotion/files/legendrepaperinprintCA.pdf
LECTURE 16GAUSS QUADRATURE
https://coast.nd.edu/jjwteach/www/www/30125/pdfnotes/lecture16_19v17.pdf
Identities and properties for associated Legendre functions
https://www.mat.univie.ac.at/~westra/associatedlegendrefunctions.pdf
Christoffel–Darboux formula -wikipedia
https://en.wikipedia.org/wiki/Christoffel%E2%80%93Darboux_formula