FastGaussQuadrature.jlで数値積分しましょう

私の近況

最近FastGaussQuadrature.jlにPRを送って、高速化したりドキュメント整備したりロゴ追加したりしました。
パッケージやGauss求積への理解が深まったので、この機会にFastGaussQuadrature.jlを世に広めようと筆を執りました。

記事の方針

Gauss求積の解説から始めるので、「コード例を早く読みたいよ！」って人はさっさとスクロールしてね。

他の積分パッケージとの比較等もしたかったですが…分量が多くなったので見送りました。なので、今回の記事は「Gauss求積の理論 + FastGaussQuadrature.jlの使い方」をお送りします。

そもそもGauss求積って何なん？

「数値積分の一種で、$n$点での関数の値のみから$2n-1$次の多項式まで厳密に計算できるもの」 がGauss求積と呼ばれるものです！

任意の連続関数は多項式近似できるので、高次の多項式まで厳密に数値積分できるGauss求積は実用の観点からも便利です。数学的には$2n-1$より強い精度が実現できないことが証明されてるので、やっぱりGauss求積はすごいです。

PDF版はこちら！

上記のフローチャートでは、積分の(数値)計算の方法として

• 解析解を頑張る
• Monte Carlo法で確率的に求める
• Newton-Cotes法で多項式近似して求める
• Gauss求積で求める(10種類)

を挙げていました。
FastGaussQuadrature.jlが使えるのは最後のGauss求積の所です。

Gauss求積の理論

少し長くなりますが、Gauss求積の理論を簡単に(?)紹介します。以下が本記事での理論の目次と概略です。

• 数値積分ってどんな風に考えるの？
  • 計算資源は有限なので、有限個の点$\{x_i\}$での値$\{f(x_i)\}$しか参照できません。
  • 数値積分作用素$\tilde{I}$に線形性を仮定すると$\tilde{I}(f)=\sum_{i}\alpha_i f(x_i)$の形に限定されます。
• Lagrange補間による多項式近似・数値積分公式
  • 積分点$\{x_i\}$が与えられれば、自動的にウェイト$\{\alpha_i\}$を決定できます。$n-1$次の精度が実現できます。
• Newton-Cotes公式の戦略
  • 積分点$\{x_i\}$をどうやって決めるか？とりあえず等間隔に取るのがNewton-Cotesの方法です。ここでも$n-1$次の精度が実現できます。
• Gauss求積の戦略
  • 積分点$\{x_i\}$をどうやって決めるか？Gauss求積では直交多項式の零点(根)を使います。なんと$2n-1$次の精度が実現できます。
• 一般化：重み関数
  • $\int_D f(x)dx\approx \sum_{i}\alpha_i f(x_i)$を一般化して$\int_D f(x) w(x)dx\approx \sum_{i}\alpha_i f(x_i)$が考えられます。この$w(x)$が重み関数です。重み関数(と積分区間)に依存して積分点$\{x_i\}$とウェイト$\{\alpha_i\}$が決定できます。

数値積分ってどんな風に考えるの？

数値計算なのでどうしても近似計算になってしまいます。
計算機(コンピュータ)の都合上、幾つかの性質が従います。

仮定：被積分関数は有限回しか値を参照できない

解析学で扱うような「極限」は、数値計算では通常は計算できません。
特に重要な性質として、「被積分関数$f$の値のは有限回しか参照できない」ことが挙げられます。この参照の回数を$n$としましょう。
本記事では基本的に1次元を扱うので、$x_1<\cdots<x_n$の点列が取れると考えてOKです。

これらの点列$\{x_i\}$とそこでの値$\{f(x_i)\}$だけで数値積分$\tilde{I}(f)$を計算したいのですが、これだとやっぱり無理があって、下の図のように厳密な積分値が一致しない$I(f)<I(g)$の場合にも数値積分値が一致$\tilde{I}(f)=\tilde{I}(g)$してしまいます。

なので、そもそも積分を離散的に計算することには多少の無理があります。^[1]

公理：「関数を数値積分する操作」は線形写像

与えられた関数$f$を積分する操作を$I$としましょう。

$$I(f) = \int_D f(x)dx$$

つまり、この写像$I$は積分作用素で、(関数)$\mapsto$(積分値)という対応を与えるものです。($D$は適当な積分区間です。)

さて、本記事(Gauss求積)の目標は、「なるべく少ない($n$個)の積分点から数値積分$\tilde{I}(f)$を求めてなるべく厳密な積分$I(f)$に近いようにする」ということでした。

$$\tilde{I}(f) \approx I(f) = \int_D f(x)dx$$

では、この数値積分$\tilde{I}$はどのように定めれば良いでしょうか？
すでに見たように、数値積分においては$n$点の情報$\{(x_i, f(x_i))\}$しか使えず、これらから$\tilde{I}(f)$を作る必要があります。

ここで 積分作用素$I$が線形作用素 だったことから 数値積分作用素$\tilde{I}$も線形作用素 だと公理で要請しましょう。

有限個の点$\{x_i\}$の参照に関する操作を$F$とし、その点での値$\{f(x_i)\}$から数値積分値を作る写像を$G$とします。これらの合成写像が数値積分で、$\tilde{I}=G\circ F$になります。

さて、$\tilde{I}, F$は線形写像だったので、必然的に$G$も線形写像になります。$G$は実線形写像$\mathbb{R}^n\to\mathbb{R}$なので、$n$個の実係数を用意すれば構成できます。これらの実数を$\{\alpha_i\}$とすれば

$$\tilde{I}(f) = \sum_i \alpha_i f(x_i)$$

のように数値積分作用素を構成できることが分かります。^[2] この$\alpha_i$をウェイトと呼びます。
上記を整理すると次式のような関係になります。^[3]

つまり、数値積分作用素$\tilde{I}$を定めるに当って、あとは積分点$\{x_i\}$とウェイト$\{\alpha_i\}$を探せばOKということになります。^[4]

定義：“良い数値積分”の尺度は「$m$次の多項式を厳密に計算可能か」で測る

$n$点を参照して、なるべく多くの被積分関数に対して良い近似を与えるような計算をしたい、というのは自然な要求でしょう。
では「良い近似」ってなんでしょうか？

近似の評価には色々な方法が考えられますが、思い切って 「より高次の多項式までは厳密に計算可能なほど良い近似」 と定義してしまいましょう！

任意の連続関数は多項式で近似できるので、多項式に関して精度評価することには一定の妥当性があるでしょう。より正確に述べるならば、数値積分$\tilde{I}$の性能を以下の式で評価することになります。

$$\max\{m \in \mathbb{N} \mid \forall f \in (m\text{次多項式全体の集合}), I(f)=\tilde{I}(f)\}$$

図を交えて考えましょう。線形代数の絵で描くと次のようになります。^[5] 写像$I-\tilde{I}$は「厳密な積分と数値積分の差」を返す写像で、この値は$0$になるのが理想です。しかし既に見たように、有限個の点でしか被積分関数の値を参照できないので、理想通りにはなりません。

$\operatorname{Ker}(I-\tilde{I})$は線形写像$I-\tilde{I}$の核で、この範囲が大きい程より厳密に積分できる範囲が広いということを表しています。その大きさを「どのくらい多項式が含まれているか」で測ることを考えるという訳です。^[6] 数値積分$\tilde{I}$が$m$次多項式まで厳密に積分できる性能を持つとき、 $m$次の精度を持つ と言います。

Lagrange補間による多項式近似・数値積分公式

ここでは積分点$\{x_i\}$が予め与えられた状況を想定して、ウェイト$\{\alpha_i\}$を定める問題を考えましょう。すでに見たように、このウェイトが決まれば数値積分の写像$\tilde{I}$が完全に決定されたことになります。

このウェイト決定問題に役立つのが … Lagrange補間です！

平面上に1点$(a_1, b_1)$が与えられれば、その1点を通る0次多項式(x軸に平行な直線)$p_0$が一意的に存在します。(desmosで確認)

$$p_0(x) = b_1$$

平面上に2点$(a_1, b_1), (a_2, b_2)$が与えられれば、その2点を通る1次多項式(直線)$p_1$が一意的に存在します。(desmosで確認)

$$p_1(x) = b_{1}\frac{x-a_{2}}{a_{1}-a_{2}}+b_{2}\frac{x-a_{1}}{a_{2}-a_{1}}$$

平面上に3点$(a_1, b_1), (a_2, b_2), (a_3, b_3)$が与えられれば、その3点を通る2次多項式(放物線)$p_2$が一意的に存在します。(desmosで確認)

$$p_2(x) = b_{1}\frac{(x-a_{2})(x-a_{3})}{(a_{1}-a_{2})(a_{1}-a_{3})}+b_{2}\frac{(x-a_{1})(x-a_{3})}{(a_{2}-a_{1})(a_{2}-a_{3})}+b_{3}\frac{(x-a_{1})(x-a_{2})}{(a_{3}-a_{1})(a_{3}-a_{2})}$$

より一般に、平面上に$n$点$(a_1, b_1), \dots, (a_n, b_n)$が与えられれば、その$n$点を通る$n-1$次多項式$p_{n-1}$が一意的に存在します。これがLagrange補間ってやつです。

$$\begin{aligned} p_{n-1}(x) = \sum_{i\in\{1,\dots,n\}}b_i\cdot \prod_{j \in \{1,\dots,n\}\setminus\{i\}}\frac{x-a_j}{a_i-a_j} \end{aligned}$$

この$p_{n-1}$は以下の性質を充たすことが簡単に確認できます。

• $x=a_k$で値$f(a_k) = b_k$を取る
• $p_{n-1}$は$n-1$次多項式

ここまで来れば、数値積分の公式まであと一歩です。点列$\{(a_i,b_i)\}$を$\{(x_i,f(x_i))\}$に置き換えて

$$\begin{aligned} p_{n-1}(x) = \sum_{i\in\{1,\dots,n\}}f(x_i)\cdot \prod_{j \in \{1,\dots,n\}\setminus\{i\}}\frac{x-x_j}{x_i-x_j} \end{aligned}$$

が得られます。これが関数$f$の近似(Lagrange補間)です。
下の図は緑色の$f$と青色の$p_3$がほぼ一致しているように見えます。

ここで数値積分$\tilde{I}(f)$を以下のように定めます。

$$\begin{aligned} \tilde{I}(f) = \int_D p_{n-1}(x) dx = \sum_{i} f(x_i) \int_D \prod_{j \in \{1,\dots,n\}\setminus\{i\}}\frac{x-x_j}{x_i-x_j} dx \end{aligned}$$

ここで$\alpha_i = \int_D \prod_{j \in \{1,\dots,n\}\setminus\{i\}}\frac{x-x_j}{x_i-x_j} dx$とおけば、当初の目標であった数値積分の公式

$$\begin{aligned} \tilde{I}(f) &= \sum_{i} \alpha_i f(x_i) \end{aligned}$$

が得られました。

重要なことは以下の4点です。

• $n-1$次の積分公式が作れた
  • つまり「$n-1$次までの多項式なら厳密に計算したるで！」ってコト
  • $n$点の積分点$\{x_i\}$が与えられている状況では、$n-1$次の積分公式($\{\alpha_i\}$の選び方)は一意的
    • $\{x_i\}$の選び方を考えなければ、最良な公式ということ
    • 一意性の証明方針：線形写像$f \mapsto I(f)-\tilde{I}(f)$が$n-1$次多項式に対して0を返すことを考えれば分かる
• $\alpha_i$の計算は被積分関数$f$に依存しない
  • つまり被積分関数ごとに計算し直す必要はないということ
• $\alpha_i$の計算は厳密に可能
  • 被積分関数は多項式なので、(展開するとかして)普通に計算できる
• $\{x_i\}$の点列の選び方は(いまのところ)どうでも良い
  • $n-1$次多項式までであれば、点列$\{x_i\}$に依存せず厳密に数値積分が可能！

Lagrange補間による数値積分をDesmosでどうぞ！

Nowton-Cotes公式の戦略

Lagrange補間で$(n-1)$次の精度を充たすウェイト$\{\alpha_i\}$が決まりましたが、積分点$\{x_i\}$はどうやって決めれば良いか分かりません。
Newton-Cotesの公式では、あまり深く考えずに、 積分点$\{x_i\}$を積分区間$[a,b]$上に等間隔に取る 戦略を採ります。

さて、Lagrange補間では有限個の点を通るように多項式近似してしまうので、十分な近似精度が得られない場合があります。
これを避けるためには、区間を適当に分割して、それぞれに対してLagrange補間による数値積分を行うのが良いでしょう。

このようにして「区間分割→Lagrange補間による積分」の公式が作れました。台形公式やSimpson公式がこれに相当します。

Gauss求積の戦略

さて、Newton-Cotes公式では積分点を等間隔に取って$(n-1)$次の精度を実現しましたが、それが最適な方法なのでしょうか？
もっと良い積分点を取って$(n-1)$次以上の多項式に対しても厳密な積分ができないでしょうか？

これを叶えるのが Gauss求積 で、 $(2n-1)$次の精度 を実現できます！(さらに、これ以上の精度を持つ数値積分が他に存在しないことが証明されています！)

大雑把な説明

数値積分をパラメータ決定問題だと思えば

$$\begin{aligned} \tilde{I}(x\mapsto x^0)&=\tilde{I}(x\mapsto x^0) \\ \tilde{I}(x\mapsto x^1)&=\tilde{I}(x\mapsto x^1) \\ &\vdots \\ \tilde{I}(x\mapsto x^{m-1})&=\tilde{I}(x\mapsto x^{m-1}) \\ \tilde{I}(x\mapsto x^m)&=\tilde{I}(x\mapsto x^m) \end{aligned}$$

を充たすような$(x_1,\dots,x_n), (\alpha_1,\dots,\alpha_n)$の$2n$変数を決める問題になります。方程式が潰れていないと仮定すれば、$2n$変数に対しては$2n$個の方程式が過不足ないことになります。つまり最大の多項式次数は$m=2n-1$になるでしょう。

では、このような$(2n-1)$次の精度を持つ積分点$\{x_i\}$を実現するにはどうすれば良いでしょうか…？
→ $n$次の直交多項式$\varphi_n$の零点を積分点とします！これがGauss求積。

直交多項式って？

関数の直交性

2つの関数$f,g$が与えられたときに

$$\langle f,g\rangle = \int_D f(x)g(x) dx$$

で内積を定めます。
ベクトルの内積

$$\langle \bm{a}, \bm{b}\rangle = \sum_i a_ib_i$$

と似た形になっていますね。
関数の内積$\langle f,g\rangle$はベクトルの内積$\langle \bm{a}, \bm{b}\rangle$と同様に双線形写像になっています。$\langle f,g\rangle=0$のときに$f$と$g$は 直交している といいます。^[7]

Fourier級数の計算で初めて関数の直交性を知った方も多いと思いますが、実は関数の直交性は数値積分にも役立つのです！
どう役に立つのか、以降の節で見ていきましょう。

直交多項式の定義と性質

任意の$n-1$次多項式と直交する$n$次多項式のことを $n$次直交多項式 と呼び、本記事では記号$\varphi_n$で表します。^[8]

$(-1,1)$区間での直交多項式の具体例を挙げましょう。

https://www.desmos.com/calculator/uedu0fwkci

ここで、それぞれの$\varphi_n$はノルムが$1$になるように正規化($\langle\varphi_n, \varphi_n\rangle=1$)されているとは限りません。^[9] なので定数倍の分だけ直交多項式$\varphi_n$には任意性があると思ってOKです。^[10]

さて、$n$次直交多項式の重要な性質として、 「ちょうど$n$個の根を積分区間に持つ」 というのがあります。代数学の基本定理から$n$個の根があるのは当たり前やろ！という話ではなくて、根の場所が積分区間に限定されているところがポイントです。

$\varphi_n$が$n$個の根を積分区間に持つことの証明
もし$\varphi_n$が積分区間に$k(<n)$個の根$(x_1,\dots,x_k)$しか持たなかったと仮定すると、$k$次多項式$f(x)=(x-x_1)\cdots(x-x_k)$と$\varphi_n$の内積$\langle f, \varphi_n \rangle$は被積分関数の符号は一定になるので積分結果(内積)は非ゼロになります。これは$\varphi_n$が直交多項式であることに矛盾。($\varphi_n$が偶数次の根(零点)を持つ場合は適当に$f$の根を間引けばOK)

なんで$2n-1$次の精度が出るの？

被積分関数$f$に対して、$n$点$\{x_i\}$でのLagrange補間を$p_{n-1}$とします。

このとき、関数$f-p_{n-1}$は$x_1,\dots,x_n$を零点に持つ関数なので

$$q(x) = \frac{f(x)-p_{n-1}(x)}{(x-x_1)\cdots(x-x_n)}$$

は極を持ちません。(零点の除去)

とくに、$f$が$(2n-1)$次の多項式であれば、割り算の結果$q(x)$は$(n-1)$次の多項式になります。
さらに、分母の$(x-x_1)\cdots(x-x_n)$は$n$次直交多項式$\varphi_n$と零点が等しい$n$次多項式なので、直交多項式の定数倍です。つまり定数$c$が存在して

$$\varphi_n(x) = c(x-x_1)\cdots(x-x_n)$$

となります。

積分しましょう！

$$\begin{aligned} \int_D f(x)-p_{n-1}(x)dx &= \frac{1}{c}\int_D\varphi_n(x) \frac{f(x)-p_{n-1}(x)}{(x-x_1)\cdots(x-x_n)} dx \\ &= \frac{1}{c}\int_D\varphi_n(x) q(x) dx = 0 \end{aligned}$$

最後の等号では、被積分関数が($n$次直交多項式)×($n-1$次多項式)になっているので直交性から$0$になることを使いました。

やったね！Gauss求積を使えば$n$点の積分点から$2n-1$次多項式まで厳密に計算できることが分かりました！

一般化：重み関数

説明の簡単のために意図的に避けていたのですが、Gauss求積には重み関数$w(x)$というものがあります。以下のような数値積分を考えます。

$$\int_D f(x) w(x) dx \approx \sum_{i}\alpha_i f(x_i)$$

これまでの積分は$w(x)=1$の定数関数と考えれば良いので、ちゃんと一般化になっていますね。

被積分関数が何らの関数$w(x)$の積で書かれている場合には、$(fw)(x)$ではなく$f(x)$の値だけから数値積分できるという訳です。これまでと同様に、$f$が$2n-1$次以下の多項式であれば数値積分は厳密な積分に一致します。

このような重み関数$w(x)$を付け加えるには、内積の定義を変更すればOKです。

$$\langle f,g\rangle = \int_D f(x)g(x) w(x)dx$$

「内積$\langle\cdot,\cdot\rangle$の定義が変わる → 直交多項式$\varphi_n$が変わる → 直交多項式の零点$\{x_i\}$が変わる → ウェイト$\{\alpha_i\}$が変わる」という順序で積分公式が変更できます。

このようにして作った積分公式も、これまでと同様に$n$点の積分点から$2n-1$次の精度を持ちます。例えば、$(0,\infty)$上で$(4x^3-2x^2-x+4)\exp(-x)$を数値積分するには2点とれば十分です。

では、与えられた重み関数$w(x)$に対する積分点$\{x_i\}$とウェイト$\{\alpha_i\}$はどのように計算すれば良いでしょうか？

ここで役立つのがFastGaussQuadrature.jlです！
このパッケージを使えば、いくつかの重み関数に対して、$n$点の積分点$\{x_i\}$とウェイト$\{\alpha_i\}$を計算することができます。

Juliaコード例

ようやく本題です。FastGaussQuadrature.jlの使い方です！

FastGaussQuadrature.jlの目標(公式ドキュメントより)

• Juliaでの最も速いGauss求積の計算
  • あらかじめ計算した$\{x_i\}, \{\alpha_i\}$をテーブルに保持していて、それを返す訳ではないです。
  • たとえば$100000000$点の計算でも数秒程度で終わります。
• 「Gauss求積の計算は(計算資源的に)大変」という人々の先入観を打ち砕く
  • 繰り返しますが、$100000000$点の計算が数秒程度で終わります。
  • 積分点の数によってアルゴリズムを変更しています。
  • たとえば、適当に初期値を決めてNewton法で収束計算したりしてます。

計算例

実行前に

using LinearAlgebra
using FastGaussQuadrature

しておきましょう。

gausslegendre(n)

• なまえ：Gauss-Legendre求積
• 重み関数：$w(x)=1$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1} (2x^4-3x^3+x-8) dx = -\frac{76}{5}$$

julia> f(x) = 2x^4-3x^3+x-8
f (generic function with 1 method)

julia> x, α = gausslegendre(3)  # 4次多項式は3点でOK
([-0.7745966692414834, 0.0, 0.7745966692414834], [0.5555555555555556, 0.8888888888888888, 0.5555555555555556])

julia> dot(α, f.(x))
-15.2

julia> -76/5
-15.2

滑らかな関数の積分①

$\sin(2x)$を積分しましょう。

$$\int_{6}^{2} \sin(2x) dx = \frac{\cos(12)-\cos(4)}{2}$$

julia> f(x) = sin(2x)
f (generic function with 1 method)

julia> x, α = gausslegendre(100)  # とりあえず100点
([-0.7745966692414834, 0.0, 0.7745966692414834], [0.5555555555555556, 0.8888888888888888, 0.5555555555555556])

julia> x = 2*(x.+2); α = 2*α;  # 置換積分

julia> dot(α, f.(x))
-0.7487487897980528

julia> -(cos(12)-cos(4))/2
-0.748748789798052

十分に精度でてますね。積分点の数でどれくらい精度出るか見てみましょう

julia> using Plots; plotly()
Plots.PlotlyBackend()

julia> Δ(n) = ((x, α) = gausslegendre(n); x = 2*(x.+2); α = 2*α; dot(α, f.(x))+(cos(12)-cos(4))/2)
Δ (generic function with 1 method)

julia> plot([abs(Δ(n)) for n in 1:13], yscale=:log10, yticks=[0.1^n for n in -1:16], legend=false, xticks=1:13)

単調減少で誤差が減って、13個の積分点をとればeps()付近に落ち着くことが分かりました。ただ、このように精度良く積分できるかは被積分関数の“複雑さ”によって変わります。
とくに、Gauss求積では「被積分関数が多項式近似しやすいか」に精度が大きく依存します。例えば$C^1$でない被積分関数だと次のような結果になります。

滑らかでない関数の積分②

$|x-3.2|$を積分しましょう。

$$\int_{2}^{6}|x-3.2|dx = 4.64$$

julia> f(x) = abs(x-3.2)
f (generic function with 1 method)

julia> Δ(n) = ((x, α) = gausslegendre(n); x = 2*(x.+2); α = 2*α; dot(α, f.(x))-4.64)
Δ (generic function with 1 method)

julia> Δ(10)
-0.016203507648055115

julia> plot([abs(Δ(n)) for n in 1:30], yscale=:log10, yticks=[0.1^n for n in -1:16], legend=false, xticks=1:30)

全然収束しないですね。もっと$10^8$くらいまで積分点を取ってみましょうか

julia> ns = [10^i for i in 1:8];

julia> plot(ns,[abs(Δ(n)) for n in ns], yscale=:log10, yticks=[0.1^n for n in -1:16], legend=false, xscale=:log10,xticks=ns)

積分点の数を$n=10^8$まで取ってもeps()辺りには到達しません。

この例のような導関数の不連続点を持つ関数を積分する際には、適当に積分区間を分割して、それぞれの区間で数値積分をするのが良いです。

gausschebyshev(n,1)

• なまえ：第一種Gauss-Chebyshev求積
• 重み関数：$w(x)=1/\sqrt{1-x^2}$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}\frac{x^{4}+2x^{3}-x+3}{\sqrt{1-x^{2}}}dx = \frac{27}{8}\pi$$

julia> f(x) = x^4 + 2x^3 -x + 3
f (generic function with 1 method)

julia> x, α = gausschebyshev(3,1)
([-0.8660254037844387, 6.123233995736766e-17, 0.8660254037844387], [1.0471975511965976, 1.0471975511965976, 1.0471975511965976])

julia> dot(α, f.(x))
10.602875205865551

julia> 27π/8
10.602875205865551

gausschebyshev(n,2)

• なまえ：第二種Gauss-Chebyshev求積
• 重み関数：$w(x)=\sqrt{1-x^2}$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} f(x)\sqrt{1-x^2} dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3) \sqrt{1-x^{2}}dx = \frac{25}{16}\pi$$

julia> f(x) = x^4 + 2x^3 -x + 3
f (generic function with 1 method)

julia> x, α = gausschebyshev(3,2)
([-0.8660254037844387, 6.123233995736766e-17, 0.8660254037844387], [1.0471975511965976, 1.0471975511965976, 1.0471975511965976])

julia> dot(α, f.(x))
4.908738521234051

julia> 25π/16
4.908738521234052

gausschebyshev(n,3)

• なまえ：第三種Gauss-Chebyshev求積
• 重み関数：$w(x)=\sqrt{(1+x)/(1-x)}$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} f(x)\sqrt{\frac{1+x}{1-x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3) \sqrt{(1+x)/(1-x)}dx = \frac{29}{8}\pi$$

julia> f(x) = x^4 + 2x^3 -x + 3
f (generic function with 1 method)

julia> x, α = gausschebyshev(3,3)
([-0.6234898018587335, 0.22252093395631445, 0.9009688679024191], [0.3379547635663544, 1.0973322242791113, 1.7063056657443274])

julia> dot(α, f.(x))
11.388273369263

julia> 29π/8
11.388273369263

gausschebyshev(n,4)

• なまえ：第四種Gauss-Chebyshev求積
• 重み関数：$w(x)=\sqrt{(1-x)/(1+x)}$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} f(x)\sqrt{\frac{1-x}{1+x}} dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3) \sqrt{(1-x)/(1+x)}dx = \frac{25}{8}\pi$$

julia> f(x) = x^4 + 2x^3 -x + 3
f (generic function with 1 method)

julia> x, α = gausschebyshev(3,4)
([-0.900968867902419, -0.22252093395631434, 0.6234898018587336], [1.7063056657443274, 1.0973322242791113, 0.3379547635663543])

julia> dot(α, f.(x))
9.817477042468102

julia> 25π/8
9.817477042468104

gaussjacobi(n)

• なまえ：Gauss-Jacobi求積
• 重み関数：$w(x)=(1-x)^\alpha (1+x)^\beta$
• 積分区間：$(-1,1)$

$$\int_{-1}^{1} f(x)(1-x)^\alpha (1+x)^\beta dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

Gauss-Jacobi求積はGauss-Legendere求積、第一種Gauss-Chebyshev求積、第二種Gauss-Chebyshev求積、第三種Gauss-Chebyshev求積、第四種Gauss-Chebyshev求積を一般化したものになります。具体的には以下のような対応になります。

• $(\alpha,\beta)=(0,0)$：Gauss-Legendere求積に一致
• $(\alpha,\beta)=(-1/2,-1/2)$：第一種Gauss-Chebyshev求積
• $(\alpha,\beta)=(1/2,1/2)$：第二種Gauss-Chebyshev求積
• $(\alpha,\beta)=(-1/2,1/2)$：第三種Gauss-Chebyshev求積
• $(\alpha,\beta)=(1/2,-1/2)$：第四種Gauss-Chebyshev求積

多項式の積分

$$\int_{-1}^{1}x^{4} (1-x)^{1/3}(1+x)^{-1/3}dx = \frac{268}{729\sqrt{3}}\pi$$

julia> f(x) = x^4
f (generic function with 1 method)

julia> x, α = gaussjacobi(3, 1/3, -1/3)
([-0.8616781426495602, -0.1482344829762509, 0.6765792922924777], [1.0648615132183943, 0.9751785456893696, 0.37835909340452684])

julia> dot(α, f.(x))
0.6668014123659403

julia> 268π/729(√3)
0.6668014123659401

gausshermite(n)

• なまえ：Gauss-求積
• 重み関数：$w(x)=\exp(-x^2)$
• 積分区間：$(-\infty, \infty)$

$$\int_{-1}^{1} f(x)\exp(-x^2) dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3)\exp\left(-x^{2}\right)dx = \frac{15}{4}\sqrt{\pi}$$

julia> f(x) = x^4+2x^3-x+3
f (generic function with 1 method)

julia> x, α = gausshermite(3)
([-1.2247448713915892, -8.881784197001252e-16, 1.2247448713915892], [0.29540897515091974, 1.181635900603676, 0.29540897515091974])

julia> dot(α, f.(x))
6.646701940895687

julia> 15(√π)/4
6.646701940895684

gausslaguerre(n)

• なまえ：Gauss-Laguerre求積
• 重み関数：$w(x)=\exp(-x)$
• 積分区間：$(0, \infty)$

$$\int_{-1}^{1} f(x)\exp(-x) dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3)\exp(-x)dx = 38$$

julia> f(x) = x^4+2x^3-x+3
f (generic function with 1 method)

julia> x, α = gausslaguerre(3)
([0.4157745567834814, 2.2942803602790467, 6.2899450829374794], [0.7110930099291746, 0.27851773356923976, 0.010389256501586142])

julia> dot(α, f.(x))
38.00000000000007

gausslaguerre(n,α)

• なまえ：一般化Gauss-Laguerre求積
• 重み関数：$w(x)=x^\alpha\exp(-x)$
• 積分区間：$(0, \infty)$

$$\int_{-1}^{1} f(x)x^\alpha\exp(-x) dx \approx \sum_{i=1}^{n} w_i f(x_i)$$

多項式の積分

$$\int_{-1}^{1}(x^{4}+2x^{3}-x+3)x^{1/3}\exp(-x)dx = \frac{5455}{81}\Gamma(4/3)$$

julia> using SpecialFunctions

julia> f(x) = x^4+2x^3-x+3
f (generic function with 1 method)

julia> x, α = gausslaguerre(3,1/3)
([0.5804651213170846, 2.6321997616344683, 6.787335117048454], [0.5902326136158121, 0.29043933488068396, 0.012307563072753214])

julia> dot(α, f.(x))
60.138311550743964

julia> 5455gamma(4/3)/81
60.138311550743886

gaussradau(n)

• なまえ：Gauss-Radau求積
• 重み関数：$w(x)=1$
• 積分区間：$[-1, 1)$

重み関数はGauss-Legendre求積と全く同じですが、積分点の左端が常に$-1$になるように拘束されています。そのため、$2n-2$次の精度しか出ませんが、左端の点での情報が重要な場合には有効です。

julia> x, α = gaussradau(3)
([-1.0, -0.2898979485566356, 0.6898979485566357], [0.2222222222222222, 1.024971652376843, 0.7528061254009346])
julia> f(x) = x^4;

julia> I = dot(α, f.(x));

julia> I ≈ 2/5
true

gausslobatto(n)

• なまえ：Gauss-Lobatto求積
• 重み関数：$w(x)=1$
• 積分区間：$[-1, 1]$

重み関数はGauss-Legendre求積と全く同じですが、積分点の両端が常に$-1, 1$になるように拘束されています。そのため、$2n-3$次の精度しか出ませんが、両端の点での情報が重要な場合には有効です。

julia> x, α = gausslobatto(4)
([-1.0, -0.4472135954999579, 0.4472135954999579, 1.0], [0.16666666666666666, 0.8333333333333333, 0.8333333333333333, 0.16666666666666666])

julia> f(x) = x^4;

julia> I = dot(α, f.(x));

julia> I ≈ 2/5
true

精度の比較 (gausslegendre, gaussradau, gausslobatto)

gausslegendre, gaussradau, gausslobattoで精度が違うとすでに述べました。最後にこれを確認しましょう。

積分して1になる多項式を用意して調べると以下のようになります。

$$f(x,n)=(n+1)2^{-(n+1)}(x+1)^n$$

https://www.desmos.com/calculator/mrne3tqcws

julia> f(x,p) = (p+1)*2.0^(-(p+1))*(x+1)^p
f (generic function with 2 methods)

julia> [((x, α) = gausslegendre(n); dot(α, f.(x,p)) ≈ 1) for n in 1:N, p in 0:2N]
5×11 Array{Bool,2}:
 1  1  0  0  0  0  0  0  0  0  0
 1  1  1  1  0  0  0  0  0  0  0
 1  1  1  1  1  1  0  0  0  0  0
 1  1  1  1  1  1  1  1  0  0  0
 1  1  1  1  1  1  1  1  1  1  0

julia> [((x, α) = gaussradau(n); dot(α, f.(x,p)) ≈ 1) for n in 1:N, p in 0:2N]
5×11 Array{Bool,2}:
 1  0  0  0  0  0  0  0  0  0  0
 1  1  1  0  0  0  0  0  0  0  0
 1  1  1  1  1  0  0  0  0  0  0
 1  1  1  1  1  1  1  0  0  0  0
 1  1  1  1  1  1  1  1  1  0  0

julia> [try ((x, α) = gausslobatto(n); Int(dot(α, f.(x,p)) ≈ 1)) catch; NaN end for n in 1:N, p in 0:2N]
5×11 Array{Real,2}:
 NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
   1    1    0    0    0    0    0    0    0    0    0
   1    1    1    1    0    0    0    0    0    0    0
   1    1    1    1    1    1    0    0    0    0    0
   1    1    1    1    1    1    1    1    0    0    0

横向きが多項式次数$p$で、縦向きが積分点の数$n$を表しています。
ちゃんと

• gausslegendre(n)：$2n-1$次の精度
• gaussradau(n)：$2n-2$次の精度
• gausslobatto(n)：$2n-3$次の精度

になることが確認できましたね。

「端点の情報を積分に込めたい」というのは何となく理解できますが、精度を落としてまでやりたい場面はまだ筆者には分かっていないところです… (知ってる人いましたら教えてください！ コメント頂きました！そちらも参照してください。)

まとめ

• 数値積分は楽しい！
  • Lagrange補間すごい！
  • 直交多項式えらい！
  • Gauss求積は高精度！
• FastGaussQuadrature.jlは数値積分に便利！
  • 計算が速い！
  • 高精度に数値積分できる！

参考文献

• FastGaussQuadrature.jlのドキュメント
  • https://juliaapproximation.github.io/FastGaussQuadrature.jl/dev/
• 数値解析入門
  • https://www.amazon.co.jp/dp/4781910386

脚注

1. $N_1$点参照して数値積分 → $N_2$点参照して数値積分 → …」のような反復計算($N_1<N_2<\cdots$)を繰り返して収束判定する方法もありますが、今回は積分点は予め固定しているものとします。

2. この形$\sum_i \alpha_i f(x_i)$はRiemann積分の有限和(短冊状のやつ)を考えれば明らかでしょう。しかし、ここで重要なのは「有限回の値の参照」と「$\tilde{I}$は線形写像」の仮定のみからこの式$\sum_i \alpha_i f(x_i)$が得られたということです。

3. $X$は適当な関数空間です。可積分関数全体とか、連続関数全体とかになりますが、今回は深く立ち入りたくないので詳しく描いていません。

4. 複素数値関数なども考えられますが、多くの場合は実数値関数の積分に帰着できるはずです。なので本記事では特に断らない限りは実数値とします。

5. ここでは$m$次多項式全体は線形空間になると仮定しています。例えば$2x+1$は1次関数ですが、$0x^2+2x+1$と考えれば2次関数とも見做せます。

6. 周期関数であれば、多項式ではなく「$\cos(mx), \sin(mx)$の線型結合までは厳密に計算可能」でも良いと思いますが、ここではとりあえず多項式を考えることにしてます。

7. 内積が定義されるためには、定義域となる関数空間を明示する必要がありますが、議論の簡単のために避けています。分かっている人はHilbert空間$L^2$とかを考えてください。

8. $n$次の多項式全体からなる線形空間の基底として$\{1, x, x^2, \dots, x^{n}\}$が取れますが、これらは(正規)直交基底にはなっていないことに注意。

9. 関数解析の文脈では、他の幾何学のように「長さ(ノルム)」それ自体が重要になることは少ないように思います。むしろノルムの構造から誘導される位相だったり、内積の構造から誘導される「直交性」だったりが重要なようです。

10. なので文献によって適当な定数倍だけ定義が変わってたりします。