FastGaussQuadrature.jlで数値積分しましょう
私の近況
最近FastGaussQuadrature.jlにPRを送って、高速化したりドキュメント整備したりロゴ追加したりしました。
パッケージやGauss求積への理解が深まったので、この機会にFastGaussQuadrature.jlを世に広めようと筆を執りました。
記事の方針
Gauss求積の解説から始めるので、「コード例を早く読みたいよ!」って人はさっさとスクロールしてね。
他の積分パッケージとの比較等もしたかったですが…分量が多くなったので見送りました。なので、今回の記事は「Gauss求積の理論 + FastGaussQuadrature.jlの使い方」をお送りします。
そもそもGauss求積って何なん?
「数値積分の一種で、
任意の連続関数は多項式近似できるので、高次の多項式まで厳密に数値積分できる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公式の戦略
- 積分点
をどうやって決めるか?とりあえず等間隔に取るのがNewton-Cotesの方法です。ここでも\{x_i\} 次の精度が実現できます。n-1
- 積分点
- Gauss求積の戦略
- 積分点
をどうやって決めるか?Gauss求積では直交多項式の零点(根)を使います。なんと\{x_i\} 次の精度が実現できます。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\}
-
数値積分ってどんな風に考えるの?
数値計算なのでどうしても近似計算になってしまいます。
計算機(コンピュータ)の都合上、幾つかの性質が従います。
仮定:被積分関数は有限回しか値を参照できない
解析学で扱うような「極限」は、数値計算では通常は計算できません。
特に重要な性質として、「被積分関数
本記事では基本的に1次元を扱うので、
これらの点列
なので、そもそも積分を離散的に計算することには多少の無理があります。[1]
公理:「関数を数値積分する操作」は線形写像
与えられた関数
つまり、この写像
さて、本記事(Gauss求積)の目標は、「なるべく少ない(
では、この数値積分
すでに見たように、数値積分においては
ここで 積分作用素
有限個の点
さて、
のように数値積分作用素を構成できることが分かります。[2] この
上記を整理すると次式のような関係になります。[3]
つまり、数値積分作用素
m 次の多項式を厳密に計算可能か」で測る
定義:“良い数値積分”の尺度は「
では「良い近似」ってなんでしょうか?
近似の評価には色々な方法が考えられますが、思い切って 「より高次の多項式までは厳密に計算可能なほど良い近似」 と定義してしまいましょう!
任意の連続関数は多項式で近似できるので、多項式に関して精度評価することには一定の妥当性があるでしょう。より正確に述べるならば、数値積分
図を交えて考えましょう。線形代数の絵で描くと次のようになります。[5] 写像
Lagrange補間による多項式近似・数値積分公式
ここでは積分点
このウェイト決定問題に役立つのが … Lagrange補間です!
平面上に1点
平面上に2点
平面上に3点
より一般に、平面上に
この
-
で値x=a_k を取るf(a_k) = b_k -
はp_{n-1} 次多項式n-1
ここまで来れば、数値積分の公式まであと一歩です。点列
が得られます。これが関数
下の図は緑色の
ここで数値積分
ここで
が得られました。
重要なことは以下の4点です。
-
次の積分公式が作れたn-1 - つまり「
次までの多項式なら厳密に計算したるで!」ってコトn-1 -
点の積分点n が与えられている状況では、\{x_i\} 次の積分公式(n-1 の選び方)は一意的\{\alpha_i\} -
の選び方を考えなければ、最良な公式ということ\{x_i\} - 一意性の証明方針:線形写像
がf \mapsto I(f)-\tilde{I}(f) 次多項式に対して0を返すことを考えれば分かるn-1
-
- つまり「
-
の計算は被積分関数\alpha_i に依存しないf - つまり被積分関数ごとに計算し直す必要はないということ
-
の計算は厳密に可能\alpha_i - 被積分関数は多項式なので、(展開するとかして)普通に計算できる
-
の点列の選び方は(いまのところ)どうでも良い\{x_i\} -
次多項式までであれば、点列n-1 に依存せず厳密に数値積分が可能!\{x_i\}
-
Lagrange補間による数値積分をDesmosでどうぞ!
Nowton-Cotes公式の戦略
Lagrange補間で
Newton-Cotesの公式では、あまり深く考えずに、 積分点
さて、Lagrange補間では有限個の点を通るように多項式近似してしまうので、十分な近似精度が得られない場合があります。
これを避けるためには、区間を適当に分割して、それぞれに対してLagrange補間による数値積分を行うのが良いでしょう。
このようにして「区間分割→Lagrange補間による積分」の公式が作れました。台形公式やSimpson公式がこれに相当します。
Gauss求積の戦略
さて、Newton-Cotes公式では積分点を等間隔に取って
もっと良い積分点を取って
これを叶えるのが Gauss求積 で、
大雑把な説明
数値積分をパラメータ決定問題だと思えば
を充たすような
では、このような
→
直交多項式って?
関数の直交性
2つの関数
で内積を定めます。
ベクトルの内積
と似た形になっていますね。
関数の内積
Fourier級数の計算で初めて関数の直交性を知った方も多いと思いますが、実は関数の直交性は数値積分にも役立つのです!
どう役に立つのか、以降の節で見ていきましょう。
直交多項式の定義と性質
任意の
ここで、それぞれの
さて、
もし
2n-1 次の精度が出るの?
なんで被積分関数
このとき、関数
は極を持ちません。(零点の除去)
とくに、
さらに、分母の
となります。
積分しましょう!
最後の等号では、被積分関数が(
やったね!Gauss求積を使えば
一般化:重み関数
説明の簡単のために意図的に避けていたのですが、Gauss求積には重み関数
これまでの積分は
被積分関数が何らの関数
このような重み関数
「内積
このようにして作った積分公式も、これまでと同様に
では、与えられた重み関数
ここで役立つのがFastGaussQuadrature.jl
です!
このパッケージを使えば、いくつかの重み関数に対して、
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)
多項式の積分
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
滑らかな関数の積分①
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求積では「被積分関数が多項式近似しやすいか」に精度が大きく依存します。例えば
滑らかでない関数の積分②
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)
全然収束しないですね。もっと
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)
積分点の数をeps()
辺りには到達しません。
この例のような導関数の不連続点を持つ関数を積分する際には、適当に積分区間を分割して、それぞれの区間で数値積分をするのが良いです。
gausschebyshev(n,1)
- なまえ:第一種Gauss-Chebyshev求積
- 重み関数:
w(x)=1/\sqrt{1-x^2} - 積分区間:
(-1,1)
多項式の積分
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)
多項式の積分
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)
多項式の積分
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)
多項式の積分
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)
Gauss-Jacobi求積はGauss-Legendere求積、第一種Gauss-Chebyshev求積、第二種Gauss-Chebyshev求積、第三種Gauss-Chebyshev求積、第四種Gauss-Chebyshev求積を一般化したものになります。具体的には以下のような対応になります。
-
:Gauss-Legendere求積に一致(\alpha,\beta)=(0,0) -
:第一種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求積(\alpha,\beta)=(1/2,-1/2)
多項式の積分
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)
多項式の積分
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)
多項式の積分
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)
多項式の積分
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求積と全く同じですが、積分点の左端が常に
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求積と全く同じですが、積分点の両端が常に
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になる多項式を用意して調べると以下のようになります。
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
横向きが多項式次数
ちゃんと
-
gausslegendre(n)
: 次の精度2n-1 -
gaussradau(n)
: 次の精度2n-2 -
gausslobatto(n)
: 次の精度2n-3
になることが確認できましたね。
「端点の情報を積分に込めたい」というのは何となく理解できますが、精度を落としてまでやりたい場面はまだ筆者には分かっていないところです… (知ってる人いましたら教えてください! コメント頂きました!そちらも参照してください。)
まとめ
- 数値積分は楽しい!
- Lagrange補間すごい!
- 直交多項式えらい!
- Gauss求積は高精度!
- FastGaussQuadrature.jlは数値積分に便利!
- 計算が速い!
- 高精度に数値積分できる!
参考文献
- FastGaussQuadrature.jlのドキュメント
- 数値解析入門
-
点参照して数値積分 →N_1 点参照して数値積分 → …」のような反復計算(N_2 )を繰り返して収束判定する方法もありますが、今回は積分点は予め固定しているものとします。 ↩︎N_1<N_2<\cdots -
この形
はRiemann積分の有限和(短冊状のやつ)を考えれば明らかでしょう。しかし、ここで重要なのは「有限回の値の参照」と「\sum_i \alpha_i f(x_i) は線形写像」の仮定のみからこの式\tilde{I} が得られたということです。 ↩︎\sum_i \alpha_i f(x_i) -
は適当な関数空間です。可積分関数全体とか、連続関数全体とかになりますが、今回は深く立ち入りたくないので詳しく描いていません。 ↩︎X -
複素数値関数なども考えられますが、多くの場合は実数値関数の積分に帰着できるはずです。なので本記事では特に断らない限りは実数値とします。 ↩︎
-
ここでは
次多項式全体は線形空間になると仮定しています。例えばm は1次関数ですが、2x+1 と考えれば2次関数とも見做せます。 ↩︎0x^2+2x+1 -
周期関数であれば、多項式ではなく「
の線型結合までは厳密に計算可能」でも良いと思いますが、ここではとりあえず多項式を考えることにしてます。 ↩︎\cos(mx), \sin(mx) -
内積が定義されるためには、定義域となる関数空間を明示する必要がありますが、議論の簡単のために避けています。分かっている人はHilbert空間
とかを考えてください。 ↩︎L^2 -
次の多項式全体からなる線形空間の基底としてn が取れますが、これらは(正規)直交基底にはなっていないことに注意。 ↩︎\{1, x, x^2, \dots, x^{n}\} -
関数解析の文脈では、他の幾何学のように「長さ(ノルム)」それ自体が重要になることは少ないように思います。むしろノルムの構造から誘導される位相だったり、内積の構造から誘導される「直交性」だったりが重要なようです。 ↩︎
-
なので文献によって適当な定数倍だけ定義が変わってたりします。 ↩︎
Discussion
Gauss求積の記事,大変ありがたいです.
Gauss-Lobattoは両端点の積分点が-1と1に固定されていますが,これは適応積分(adaptive quadrature)のときに嬉しいと思います(ここでは,所望の精度が出るまで,区間を分割して各区間で数値積分することを繰り返すアルゴリズムの総称を適応積分と言っています).Gauss-Lobattoを用いた適応積分では,両端点の積分点が固定されていることで,区間の端点の被積分間数値を使い回せます.このため,Gauss-Legendreを適応積分で用いたときと比べて,被積分関数の計算回数が少なくなる可能性があると思います.
コメントありがとうございます!
確かにそれで計算回数を減らせますね!あらかじめ(連続な)被積分関数の(有限個の)導関数不連続点が分かっている状況でも、[Gauss-Radau] - [Gauss-Lobatto] - … - [Gauss-Lobatto] - [Gauss-Radau]のようにすれば僅かに被積分関数の計算回数を減らせそうです。