Gauss型積分公式について
数値積分の一種であるGauss型積分公式の説明. こちらの記事を見て, ふと調べてみたことをまとめたもの.
0. 概要
積分を数値的に計算する手法の一つであるGauss型積分公式は, 何をやっているのかいまいちピンとこないという人が多く, 精度はいいのにあまり使われてない印象. 通常使われていると思われる Newton-Cotes の積分公式との比較を交えつつ Gauss型積分公式について説明する.
参考文献
1. 数値積分の考え方
積分できない被積分関数(
1.1 Newton-Cotes の積分公式の考え方
よく使用される Newton-Cotes の積分公式では, 積分区間内に等間隔に設定した標本点から被積分関数(
ここで
図は
このような現象を避けるため, 通常は例えば
Simpson の方法(複合公式)では図にあるように実線の間の区間を2次関数で近似する. 区間の両端と区間の中点(図の破線)を標本点として選び, 3点を通る曲線(2次関数)を近似曲線とする.
一方, 標本点の取り方を工夫する(等間隔から変更する)ことによってもルンゲの現象を避けることができる. Gauss型積分公式では標本点を直交多項式の零点に取ることによりルンゲの現象を回避している.
1.2 Gauss型積分公式の考え方
Gauss型積分公式は, Newton-Cotes の積分公式(複合公式)と比較すると, 積分区間を細分化しないかわりに標本点を等間隔ではなく直交多項式の零点にとるところに特徴がある. これについては後で触れることにして, まずは被積分関数を直交多項式を使って近似することを考える.
Gauss型積分公式では直交多項式
ここで
(ベクトルの直交関係と類似の関係が成立しており, 直交多項式といわれる.)
特に
近似関数を作るのに直交多項式を使うと, 積分を劇的に簡単にすることができる.
すなわち, 近似式における直交多項式の係数がわかれば積分の近似値が計算できる.
標本点をこのように取ると, (もともとの直交性とは別に)選点直交性といわれる以下に示す関係が成立する(内積と区別するため,
これを使って実際に
両辺に
すなわち,
これを先ほどの
(前にも記載したとおり, ここでの式変形等はわかりやすさを優先して本質を損なわない程度に式を簡略化していることに注意.)
ポイントをまとめると, 以下のとおり.
- 直交多項式で近似式を作成(直交性により積分を劇的に簡略化できる)
- 標本点を直交多項式の零点にとる(ルンゲの現象を回避)
- 近似式の係数は選点直交性により容易に計算可能
- 最終的な積分公式は
とa_i の積を合計したものであり, 非常に簡単な形をしているf(x_i)
結局, 標本点(
1.3 Gauss型積分公式の精度
上で見たように被積分関数の近似式は標本点が
近似多項式の係数
被積分関数を
これを
これを積分して,
(※
式変形の途中で近似は行っていないので, 結局以下が言えたことになる.
これは
証明の途中で重要なのは, 標本点が
また,
2. 数値例
2.1 5次関数の積分
上の精度のところで述べたように標本点が
以下のような被積分関数を考えよう(まえにNewton補間の記事を書いたときに適当に作った関数).
これを Gauss-Legendre の積分公式(Legendre多項式を使った Gauss型積分公式)によって計算する. Legendre多項式は以下の内積によって定義される(前節の簡略化された例とは異なり定数係数
最初のいくつかの Legendre多項式を示すと以下のようになる.
また選点直交性は以下のとおり.
(ここでも
まず, Gauss-Legendre の積分公式の積分区間は
ここで改めて積分変数を
後の記載のため, 被積分関数を改めて
念の為,
2.1.1 解析解
まずは解析的に積分を計算してみよう.
2.1.2 Newton-Cotes の公式(Simpson の方法)
次に Newton-Cotesの公式(Simpsonの方法)を使ってみよう.
julia> simpson(fifth,-1,1,3)
0.5
被積分関数が定数項を除き奇関数なので,偶然正確な結果が得られてしまった(例が悪かったかな). 参考までに近似関数をプロットすると以下のとおり.
2.1.2 Gauss-Legendre積分(FastGaussQuadrature.jlの使用)
更に, 既存のライブラリを使って同じ結果になるかどうかを確認しよう. 前に述べたように, Gauss型積分公式では積分を以下の形で計算する.
ここで,
FastGaussQuadrature.jl を使った積分の計算
# FastGaussQuadrature.jlを使った積分の計算
function fgq(f,n,logging=false)
x,a=gausslegendre(n)
if logging
@show x; @show a; @show f.(x)
end
return dot(a,f.(x)) # dot はベクトルの内積
end # fgq
結果は以下のとおり.
julia> fgq(fifth,3,true)
x = [-0.7745966692414834, 0.0, 0.7745966692414834]
a = [0.5555555555555556, 0.8888888888888888, 0.5555555555555556]
f.(x) = [0.25136232189202845, 0.25, 0.24863767810797155]
0.5
積分値が正しく計算されていることがわかる.
2.1.3 Gauss-Legendre積分(step by step での計算)
ここではライブラリを使わず手作業で step by step に積分の値を求めてみる. まずは
これを
julia では以下のようなイメージになる.
標本点を求める
# 例えば n=3
x=basis.(Legendre,n) |> roots |> real
次に
Legendre多項式と内積の係数を求める
# 例えば n=3
p=basis.(Legendre,0:(n-1))
λ=[SP.innerproduct(Legendre,p[k+1],p[k+1]) for k in 0:(n-1)]
これで,
積分のウェイトを求める
# 手作業での計算: 必要なパラメータの算出
function gl(n) # Gauss-Legendre
p=basis.(Legendre,0:(n-1))
λ=[SP.innerproduct(Legendre,p[k+1],p[k+1]) for k in 0:(n-1)]
x=basis.(Legendre,n) |> roots |> real
a=1 ./ [dot(1 ./λ,[p[k+1](x[i])^2 for k in 0:(n-1)]) for i in 1:n]
return x,a,p,λ
end # gl
積分の計算には必要ないが近似関数も求めておこう.
近似多項式の係数を求める
# 標本点によって決るLegendre多項式由来のパラメータ
x,a,p,λ = gl(n)
# 被積分関数 f によって決まる近似多項式の係数
c(f) = [dot(a,p[k+1].(x),f.(x))/λ[k+1] for k in 0:(n-1)]
(後の都合で, ここでは係数
ここまでの結果を用いて積分を計算し, 近似関数をプロットしてみよう.
手作業による積分計算
# 手作り: FastGaussQuadrature の gausslegendre 関数相当
function gl(n) # Gauss-Legendre
p=basis.(Legendre,0:(n-1))
λ=[SP.innerproduct(Legendre,p[k+1],p[k+1]) for k in 0:(n-1)]
x=basis.(Legendre,n) |> roots |> real
a=1 ./ [dot(1 ./λ,[p[k+1](x[i])^2 for k in 0:(n-1)]) for i in 1:n]
return x,a,p,λ
end # gl
# 手作り: ガウス積分の計算(gl(n)を使った計算)
function my_gq(f,n,logging=false)
x,a,_,_ = gl(n)
if logging
@show x; @show a; @show f.(x)
end
return dot(a,f.(x))
end # my_gq
# 近似多項式の係数を与える関数を返す
function fn_coef(n)
x,a,p,λ = gl(n)
return f -> [dot(a,p[k+1].(x),f.(x))/λ[k+1] for k in 0:(n-1)]
end # fn_coef
# 近似多項式を返す
function fn(f,n,logging=false)
p=basis.(Legendre,0:(n-1))
c=fn_coef(n)
if logging
@show c(f)
end
return t -> dot(c(f),[p[k+1](t) for k in 0:(n-1)])
end # fn
# 近似多項式をプロットする
function plot_fn(f,n)
x=basis.(Legendre,n) |> roots |> real
xs=range(-1,1,400)
pal=palette(:default)
plt = plot(legend=:topright)
plot!(xs,f.(xs),lc=pal[1],label=L"f(x)")
scatter!(plt,x,f.(x),mc=pal[2],label=L"f(x_i)")
plot!(plt,xs,fn(f,n).(xs),lc=pal[3],label=L"f_n(x)")
display(plt)
# println(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
# savefig(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
end # plot_fn
julia> Prg.my_gq(Prg.fifth,3,true)
x = [-0.7745966692414833, 5.048709793414476e-29, 0.7745966692414833]
a = [0.5555555555555557, 0.8888888888888888, 0.5555555555555557]
f.(x) = [0.25136232189202845, 0.25, 0.24863767810797155]
0.5
積分値とパラメータが FastGaussQuadrature.jl を使用したときと(ほぼ)同じになっていることがわかる. また, 近似関数をプロットすると以下のとおり.
プロットを見るとお世辞にも近似されているとは言えない形をしている.
ulia> fn_coef(3)(fifth)
3-element Vector{Float64}:
0.25
-0.0017587499999999895
-6.93889390390723e-17
ちなみに, 被積分関数を
(手計算では厳しいので, WolframAlphaを利用した. 商と余りではなぜか勝手に全体が2倍されている.)
これからわかるように, 上で求めた近似曲線は被積分関数を
5次関数に対する近似曲線が直線, という一見直観に反する結果になっているが, 積分の近似という意味ではこれで正しく近似できている.
2.2 コーシー・ローレンツ関数の積分
2つ目の例として, ルンゲの現象のときに使用した以下の関数(コーシー・ローレンツ関数)を考えよう.
積分区間は最初から
2.2.1 解析解
この関数も解析的に積分することができる.
2.2.2 Newton-Cotes の公式(Simpson の方法)
次に Newton-Cotesの公式(Simpsonの方法)を使ってみよう.
julia> simpson(cauchy,-1,1,11)
0.5698340874811464
解析解と比較すると上振れしてますね. 近似関数のプロットを見てみると,
(再掲)
山の頂点付近で近似式の方が結構大きくなっているので, 感触としては妥当な結果です.
2.2.3 Gauss-Legendre積分(FastGaussQuadrature.jlの使用)
次にFastGaussQuadrature.jl を使ってみよう. なお, 標本点の個数はルンゲの現象の項に合わせて11個とする.
julia> fgq(cauchy,11)
0.5624581121773549
2.2.4 Gauss-Legendre積分(step by step での計算)
手作業で作った関数でも実行してみる. こちらも上で作ったプログラムを被積分関数を入れ替えて実行すればよい.
julia> my_gq(cauchy,11)
0.5624581121773546
FastGaussQuadrature.jlとほとんど同じ結果で解析解と比較すると上振れしてます. 近似関数のプロットを確認すると,
目分量では何とも言えないけど, 裾の方はプラスとマイナスでキャンセルしてて, やっぱり中心付近が大きくなってるのが効いてるのかな. FastGaussQuadrature.jl とほとんど同じなので, 間違ってはいないと思います.
2.2.5 Gauss-Legendre積分とNewton-Cotes の積分公式の収束性の比較
上で見たように標本点が11個だとGauss-Legendre積分でもNewton-Cotes公式でも誤差が大きい. しかし, 標本点を増やすと Gauss-Legendre積分の方が急速に真値に収束してゆく. 横軸の標本点の数, 縦軸に積分値をとって図示すると以下のようになる.
見た感じは同じように収束しているが, より詳細にみるために解析解との誤差をプロットしてみる(数値計算の教科書でよく見るプロット(笑)).
誤差でみると Gauss-Legendre積分の方がはるかに早く真値に収束していることがわかる. 標本数が
3. サンプルプログラム
上の記事で計算や図のプロットに使用したプログラム.
サンプルプログラム
module Prg
using Polynomials
using SpecialPolynomials
using FastGaussQuadrature
using LaTeXStrings
using Plots
const SP=SpecialPolynomials
# LinearAlgebra の dot 関数と同じように定義
dot(x,y) = sum([xi*yi for (xi,yi) in zip(x,y)])
dot(x,y,z) = sum([xi*yi*zi for (xi,(yi,zi)) in zip(x,zip(y,z))])
# 5次関数(積分区間[0,1])
fifth_(x)=0.5+(x-0.02)*(x-0.25)*(x-0.5)*(x-0.75)*(x-0.98)
# 積分区間調整後の関数(積分区間[-1,1])
fifth(x)=0.5*fifth_(0.5*(x+1))
# コーシー-ローレンツ(Cauchy-Lorentz)関数
cauchy(x)=1/(1+25*x*x)
# Lagrange の補間関数
lagrange(x,y)= t ->
[yi*([(t-xj)/(xi-xj) for xj in x if xj!=xi] |> prod)
for (xi,yi) in zip(x,y)] |> sum
# FastGaussQuadrature.jlを使った積分の計算
function fgq(f,n,logging=false)
x,a=gausslegendre(n)
if logging
@show x; @show a; @show f.(x)
end
return dot(a,f.(x))
end # fgq
# 手作り: FastGaussQuadrature の gausslegendre 関数相当
function gl(n) # Gauss-Legendre
p=basis.(Legendre,0:(n-1))
λ=[SP.innerproduct(Legendre,p[k+1],p[k+1]) for k in 0:(n-1)]
x=basis.(Legendre,n) |> roots |> real
a=1 ./ [dot(1 ./λ,[p[k+1](x[i])^2 for k in 0:(n-1)]) for i in 1:n]
return x,a,p,λ
end # gl
# 手作り: ガウス積分の計算(gl(n)を使った計算)
function my_gq(f,n,logging=false)
x,a,_,_ = gl(n)
if logging
@show x; @show a; @show f.(x)
end
return dot(a,f.(x))
end # my_gq
# 近似多項式の係数を与える関数を返す
function fn_coef(n)
x,a,p,λ = gl(n)
return f -> [dot(a,p[k+1].(x),f.(x))/λ[k+1] for k in 0:(n-1)]
end # fn_coef
# 近似多項式を返す
function fn(f,n,logging=false)
p=basis.(Legendre,0:(n-1))
c=fn_coef(n)
if logging
@show c(f)
end
return t -> dot(c(f),[p[k+1](t) for k in 0:(n-1)])
end # fn
# Newton-Cotes(Simpsonの方法)
function simpson(f, a, b, n)
# n は奇数で3以上と仮定
h = (b-a)/(n-1)
x = a:h:b
y = f.(x)
# x,y のインデックスの想定は 0,1,...,n だが,
# julia の配列のインデックスは1から始まる.
# そのため, 配列を参照する際に +1 を付けて調整する.
s = sum([(y[i-2+1]+4*y[i-1+1]+y[i+1]) for i in 2:2:n])
return s*h/3
end # simpton
# コーシーローレンツ関数を使った収束性の確認
# 積分値の比較
function shusoku()
ns = 11:2:31
# ns = vcat(ns,1 .+ 60:10:100)
I_nc = (n -> simpson(cauchy,-1,1,n)).(ns)
I_gl = (n -> fgq(cauchy,n,false)).(ns)
analytic = n -> (atan(5)-atan(-5))/5
pal=palette(:default)
plt = plot()
plot!(plt,ns,analytic.(ns)
,ls=:dash ,lc=:black
,xlabel=L"n",ylabel=L"I_n"
,label="Analytical solution")
plot!(plt,ns,I_nc,lc=pal[1]
,label="Newton-Cotes(Simpson)")
plot!(plt,ns,I_gl,lc=pal[3]
,label="Gauss-Legendre")
display(plt)
# println(string(nameof(var"#self#"),".png"))
# savefig(string(nameof(var"#self#"),".png"))
end # shusoku
# コーシーローレンツ関数を使った収束性の確認
# 真値との誤差の比較
function shusoku2()
ns=[2^m+1 for m in 3:13]
I_nc = (n -> simpson(cauchy,-1,1,n)).(ns)
I_gl = (n -> fgq(cauchy,n,false)).(ns)
analytic = (atan(5)-atan(-5))/5
d_nc = (t->t-analytic |> abs |> log10).(I_nc)
d_gl = (t->t-analytic |> abs |> log10).(I_gl)
log_ns = log2.(ns)
pal=palette(:default)
plt = plot(legend=:topright)
plot!(plt,log_ns,d_nc
,label="Newton-Cotes"
,lc=pal[1]
,xlabel=L"\log_2\ (n), n=2^m+1"
,ylabel=L"\log_{10}\ (I_n-I_0)")
plot!(plt,log_ns,d_gl
,label="Gauss-Legendre"
,lc=pal[3])
display(plt)
# println(string(nameof(var"#self#"),".png"))
# savefig(string(nameof(var"#self#"),".png"))
end # shusoku2
# 単純に与えられた関数をプロット
function plot_func(f,xmin,xmax,lbl)
x=range(xmin,xmax,400)
plt=plot()
plot!(x,f.(x),label=lbl)
display(plt)
# println(string(nameof(var"#self#"),"_",Symbol(f),".png"))
# savefig(string(nameof(var"#self#"),"_",Symbol(f),".png"))
end # plot_func
# lagrange 補間の近似関数をプロット
function plot_lagrange(f,n)
x = range(-1,1,400)
y = f.(x)
xn = range(-1,1,length=n)
yn = f.(xn)
fn = lagrange(xn,yn)
pal=palette(:default)
plt = plot(legend=:top)
plot!(x,y,lc=pal[1],label=L"f(x)")
plot!(plt,x,fn.(x),lc=pal[3],label=L"f_n(x)")
scatter!(plt,xn,yn,mc=pal[2],label=L"f(x_i)")
display(plt)
# println(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
# savefig(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
end # plot_lagrange
# 近似多項式をプロットする
function plot_fn(f,n)
x=basis.(Legendre,n) |> roots |> real
xs=range(-1,1,400)
pal=palette(:default)
plt = plot(legend=:topright)
plot!(xs,f.(xs),lc=pal[1],label=L"f(x)")
plot!(plt,xs,fn(f,n).(xs),lc=pal[3],label=L"f_n(x)")
scatter!(plt,x,f.(x),mc=pal[2],label=L"f(x_i)")
display(plt)
# println(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
# savefig(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
end # plot_fn
function plot_simpson(f,n)
xn = range(-1,1,length=n)
xn_o = [xn[i] for i in 1:n if isodd(i)]
xn_e = [xn[i] for i in 1:n if iseven(i)]
function fn(t)
res=0
for i in 3:2:n
if t<=xn[i]
xis = [xn[i-2],xn[i-1],xn[i]]
res = lagrange(xis,f.(xis))(t)
break
end
end
return res
end # fn
x = range(-1,1,400)
y = f.(x)
pal = palette(:default)
plt = plot()
plot!(x,y,lc=pal[1],label=L"f(x)")
plot!(ylim=ylims(),legend=(0.85,0.95))
plot!(plt ,[xn_o';xn_o']
,[zeros(Int((n-1)/2+1))';f.(xn_o)']
,lc=:black,label="")
plot!(plt ,[xn_e';xn_e']
,[zeros(Int((n-1)/2))';f.(xn_e)']
,ls=:dash,lc=:gray
,label="")
plot!(plt,x,fn.(x),lc=pal[3],label=L"f_n(x)")
scatter!(plt,xn,f.(xn),mc=pal[2],label=L"f(x_i)")
display(plt)
# println(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
# savefig(string(nameof(var"#self#"),"_",Symbol(f),"_",n,".png"))
end # plot_simpson
function main()
# Runge の現象をプロット
plot_lagrange(cauchy,11)
# Newton-Cotes(Simpson)をプロット
plot_simpson(cauchy,11)
# fifth_ 関数をプロット([0,1])
plot_func(fifth_,0,1,L"\tilde{f}(x)")
# fifth 関数をプロット([-1,1])
plot_func(fifth,-1,1,L"f(x)")
# fifth 関数の積分計算(Simpsonの方法)
println("\nNewton-Cotes (Simpson の公式)の使用:")
println("I_3=",simpson(fifth,-1,1,3))
plot_simpson(fifth,3)
# simpson(fifth,-1,1,5)
# plot_simpson(fifth,5)
# fifth 関数の積分計算(FastGaussQuadrature.jl を使用)
println("\nFastGaussQuadrature.jl の使用:")
println("I_3=",fgq(fifth,3,true))
# fifth 関数の積分計算(手作業による計算)
println("\n手作業によるGauss積分:")
println("I_3=",my_gq(fifth,3,true))
# fifth 関数を積分するための近似関数をプロット
plot_fn(fifth,3)
# Newton-Cotes (Simpson の公式)の使用:
println("\nNewton-Cotes (Simpson の公式)の使用:")
println("I_11=",simpson(cauchy,-1,1,11))
# FastGaussQuadrature.jl の使用
println("\nFastGaussQuadrature.jl の使用:")
println("I_11=",fgq(cauchy,11))
# 手作業による計算
println("\n手作業によるGauss積分:")
println("I_11=",my_gq(cauchy,11))
# cauchy 関数を積分するための近似関数をプロット
plot_fn(cauchy,11)
# 収束性の確認
shusoku()
shusoku2()
end # main
end # module Prg
using .Prg
4. まとめ
Gauss型積分公式について, 公式の導出をおこなうとともに, 具体例により計算の様子を見てきた. Gauss型積分公式には以下のような明確なメリットがある.
-
点の標本点で,n 次の精度がある(最適公式)2n-1 - 積分に必要な標本点(
次の直交多項式の零点)とウェイト(上の記載ではn )は被積分関数とは関係なく事前に計算できる(FastGaussQuadrature.jl などでは事前に大量に準備されており, 1億個の標本点でも数秒で計算が終わるらしい).a_i
一方, デメリットとしては以下のような点が挙げられる.
- 直交多項式の零点とウェイトを外部のライブラリに依存する場合, その正確性が不安
- Newton-Cotes公式と比較するとロジックがわかりにくい
1つ目のデメリットは実際に過去に誤植により多くのリソースが無駄になったということがあったらしい(数値解析の p.94 参照). 最近ではいわゆるリーナスの法則により間違いが放置されることは少ないと思われるが, このようなリスクがあることは心にとめておくべきだろう.
2つ目については, 本記事が少しでも役に立てば幸いである. ただ, 自分でこのような記事を書いていてなんですが, どこかのSF作家が言っていたように,
十分に発達した数学は魔法と見分けがつかない
という印象はぬぐえません(これ気に入ったので座右の銘にします(笑)).
Discussion