🌊

Gauss型積分公式について

2025/01/22に公開

数値積分の一種であるGauss型積分公式の説明. こちらの記事を見て, ふと調べてみたことをまとめたもの.

0. 概要

積分を数値的に計算する手法の一つであるGauss型積分公式は, 何をやっているのかいまいちピンとこないという人が多く, 精度はいいのにあまり使われてない印象. 通常使われていると思われる Newton-Cotes の積分公式との比較を交えつつ Gauss型積分公式について説明する.

参考文献

1. 数値積分の考え方

積分できない被積分関数(f(x))を積分できる, より簡単な関数で近似して積分する, というのが基本的な考え方.

1.1 Newton-Cotes の積分公式の考え方

よく使用される Newton-Cotes の積分公式では, 積分区間内に等間隔に設定した標本点から被積分関数(f(x))に対する近似関数(f_n(x), 多項式)を作成し, それを積分することによって求める積分(I)の近似値とする.

I=\int_a^b f(x) dx \sim \int_a^b f_n (x) dx

ここで f_n(x) は, 標本点を (x_i,f(x_i)) としたとき, f_n(x_i)=f(x_i) (標本点を通るという条件)を満たす, 高々 n-1 次の多項式である. f_n(x) は多項式なので, 簡単に積分できる. 一般には標本点を増やせば近似の精度は向上すると考えられるが, 実際には元の関数とまったく異なる関数になる場合がある(ルンゲの現象).

図は f(x)=1/(1+25x^2) について, 区間 [-1,1] で, 標本点を0.2間隔で11個(n=11)とって近似曲線 f_n(x) を作成したものである.

このような現象を避けるため, 通常は例えば n=3 (2次関数による近似, Simpsonの公式) 程度にとどめておき, 代わりに積分区間を細分化して近似精度を向上させる(複合公式).

Simpson の方法(複合公式)では図にあるように実線の間の区間を2次関数で近似する. 区間の両端と区間の中点(図の破線)を標本点として選び, 3点を通る曲線(2次関数)を近似曲線とする.

一方, 標本点の取り方を工夫する(等間隔から変更する)ことによってもルンゲの現象を避けることができる. Gauss型積分公式では標本点を直交多項式の零点に取ることによりルンゲの現象を回避している.

1.2 Gauss型積分公式の考え方

Gauss型積分公式は, Newton-Cotes の積分公式(複合公式)と比較すると, 積分区間を細分化しないかわりに標本点を等間隔ではなく直交多項式の零点にとるところに特徴がある. これについては後で触れることにして, まずは被積分関数を直交多項式を使って近似することを考える.

Gauss型積分公式では直交多項式 \{p_k(x)\} を使って, 以下のように被積分関数を近似する.

\begin{align*} f_n(x) &= c_0 \cdot p_0(x) + c_1 \cdot p_1(x) + \cdots + c_{n-1}\cdot p_{n-1}(x) \\ &=\sum_{k=0}^{n-1} c_k \cdot p_k(x) \end{align*}

ここで p_k(x)k 次の直交多項式であり以下の関係式を満たす(内積を (,) で表す. また, 説明のわかりやすさのため, 本質を損なわない範囲で式を簡略化(定数係数の省略)している).

\begin{align*} (p_k,p_l) &=\int_a^b p_k(x)p_l(x)dx \\ &=\delta_{kl} \\ &= \begin{cases} 1 &\text{for } k=l \\ 0 &\text{for } k\neq l \end{cases} \end{align*}

(ベクトルの直交関係と類似の関係が成立しており, 直交多項式といわれる.)
特に k=0 の場合は0次の多項式, すなわち定数関数である(p_0(x)=\mu_0).

近似関数を作るのに直交多項式を使うと, 積分を劇的に簡単にすることができる.

\begin{align*} I&=\int_a^b f(x) dx \\ &\sim \int_a^b f_n (x) dx \\ &=\int_a^b \sum_{k=0}^{n-1} c_k \cdot p_k(x) dx \\ &=\int_a^b \dfrac{1}{p_0}\sum_{k=0}^{n-1} c_k \cdot p_0 \cdot p_k(x) dx \\ &=\dfrac{1}{\mu_0}\sum_{k=0}^{n-1} c_k \int_a^b p_0 \cdot p_k(x) dx \\ &=\dfrac{1}{\mu_0}\sum_{k=0}^{n-1} c_k (p_0, p_k) \\ &=\dfrac{1}{\mu_0}\sum_{k=0}^{n-1} c_k \delta_{0k} \\ &=\dfrac{c_0}{\mu_0} \end{align*}

すなわち, 近似式における直交多項式の係数がわかれば積分の近似値が計算できる. n個の係数 \{c_k\} を定めるためには n 個の条件, すなわち n個の標本点 (x_i,f(x_i)) を通る, という条件を課す必要がある. ここで, Gauss型積分公式の場合には, 標本点を p_n(x)n個の零点にとる(これがGauss型積分公式の重要なポイント).

標本点をこのように取ると, (もともとの直交性とは別に)選点直交性といわれる以下に示す関係が成立する(内積と区別するため, <,\!>_nという記号を使用している. 添字 n により p_n(x) の零点を使用していることを表す).

\begin{align*} <p_k,p_l\!>_n&:=\sum_{i=1}^n a_i \cdot p_k(x_i)\cdot p_l(x_i) \\ &=\delta_{kl} \\ a_i^{-1} &:= \sum_{k=0}^{n-1} p_k(x_i)^2 \end{align*}(上の和は標本点に関する和($i$に関する和)

これを使って実際に \{c_k\} を求めてみよう. 近似式が (x_i,f(x_i)) を通るので, f_n(x_i)=f(x_i) であり, これを代入すると,

f(x_i)=\sum_{l=0}^{n-1} c_l \cdot p_l(x_i)

両辺にa_i p_k(x_i) をかけて i についての和を取ると, 選点直交性より,

\begin{align*} <p_k,f\!>_n&=\sum_{l=0}^{n-1} c_l \cdot <p_k,p_l\!>_n \\ &=\sum_{l=0}^{n-1} c_l\delta_{kl} \\ &=c_k \end{align*}

すなわち,

\begin{align*} c_k &= <p_k,f\!>_n \\ &= \sum_{i=1}^n a_i \cdot p_k(x_i)\cdot f(x_i) \\ \end{align*}

k=0 とすると,

c_0 = \mu_0 \sum_{i=1}^n a_i \cdot f(x_i)

これを先ほどのIの近似式に代入すると, 積分I に関する Gauss型の積分公式が得られる.

\begin{align*} I &\sim \sum_{i=1}^n a_i \cdot f(x_i) \\ a_i^{-1} &:= \sum_{k=0}^{n-1} p_k(x_i)^2 \end{align*}

(前にも記載したとおり, ここでの式変形等はわかりやすさを優先して本質を損なわない程度に式を簡略化していることに注意.)

ポイントをまとめると, 以下のとおり.

  • 直交多項式で近似式を作成(直交性により積分を劇的に簡略化できる)
  • 標本点を直交多項式の零点にとる(ルンゲの現象を回避)
  • 近似式の係数は選点直交性により容易に計算可能
  • 最終的な積分公式は a_if(x_i) の積を合計したものであり, 非常に簡単な形をしている

結局, 標本点(x_i)とウェイト(a_i)がわかれば積分を計算可能で, さらにすばらしいことに, これらのパラメータは被積分関数とは無関係に定めることができる. 一方うれしくない点としては, 標本点は数値的に求める必要があり, 異なる n に対してその都度すべて求めなおす必要がある. しかし, Gauss型積分のライブラリでは n に応じて x_ia_i 事前に計算して保持しており, 積分を高速に計算できる.

1.3 Gauss型積分公式の精度

上で見たように被積分関数の近似式は標本点がn個の場合n-1次の多項式であるので, 当然被積分関数が n-1次の多項式であれば積分の値を exact に求めることができる. しかし, 驚くべきことに, 実際に exact に求まるのは被積分関数が2n-1次の多項式までである.

近似多項式の係数 c_kn個, 標本点が n個あり, パラメータの個数だけ見れば 2n-1次の多項式の係数の数と一致しているので, おかしくはないだろうと思いつつも近似多項式はあくまで n-1次式なので, 自明ではないだろう.

被積分関数を2n-1次式であると仮定する.

f(x) = s_{2n-1}(x)

これをn次の直交多項式p_n(x)で割り算して,

s_{2n-1}(x) = q_{n-1}(x)p_n(x)+r_{n-1}(x)

これを積分して,

\begin{align*} I&=\int_a^b f(x)dx \\ &=\int_a^b s_{2n-1}(x)dx \\ &=\int_a^b q_{n-1}(x)p_n(x)dx + \int_a^b r_{n-1}(x)dx \\ &=(q_{n-1},p_n)+\int_a^b r_{n-1}(x)dx \\ &=\int_a^b r_{n-1}(x)dx \quad (\because (q_{n-1},p_n)=0) \\ &=\sum_{i=1}^n a_i \cdot r_{n-1}(x_i) \quad (※)\\ &=\sum_{i=1}^n a_i \cdot (q_{n-1}(x_i)p_n(x_i)+r_{n-1}(x_i)) \quad (\because p_n(x_i)=0)\\ &=\sum_{i-1}^n a_i \cdot s_{2n-1}(x_i) \end{align*}

(※ r_{n-1}n-1次式なので近似関数は被積分関数に一致し, Gauss型積分公式で正確な値が求まることに注意.)

式変形の途中で近似は行っていないので, 結局以下が言えたことになる.

I=\sum_{i-1}^n a_i \cdot s_{2n-1}(x_i)

これはn個の標本点によるGauss型積分公式で 2n-1次の多項式まで正確に計算できることを意味している.

証明の途中で重要なのは, 標本点が p_n(x) の零点であることを利用して元の関数 s_{2n-1}(x) に対する公式適用の形に戻していることろ(p_n(x) を使って r_{n-1}(x) の積分に帰着させているところは単純に式変形のテクニックなので, Gauss型積分公式に特有のものとまではいえない). これにより, わざわざ r_{n-1}(x) に書き換えずとも, もとの関数 s_{2n-1}(x) に Gauss型積分の公式を適用すれば n個の標本点で 2n-1 の精度が得られる, というありがたい仕組みになっている.

また, s_{2n-1}(x_i)=r_{n-1}(x_i) であるから, 前節の手続きで計算していた近似式は実際には s_{2n-1}(x) を近似しているのではなく, どちらかといえば r_{n-1}(x) を近似していることがわかる.

2. 数値例

2.1 5次関数の積分

上の精度のところで述べたように標本点がn個の Gauss型積分公式では 2n-1次の多項式まで正確に計算することができる, ということなので, n=3 として, 2n-1=5次関数を正確に積分できるのか, 実際に計算してみよう.

以下のような被積分関数を考えよう(まえにNewton補間の記事を書いたときに適当に作った関数).

\begin{align*} \tilde{f}(x)&=0.5+(x-0.02)(x-0.25)(x-0.5)(x-0.75)(x-0.98) \\ I&=\int_0^1 \tilde{f}(x) dx \end{align*}

これを Gauss-Legendre の積分公式(Legendre多項式を使った Gauss型積分公式)によって計算する. Legendre多項式は以下の内積によって定義される(前節の簡略化された例とは異なり定数係数\lambda_kが現れることに注意. また Gauss-Legendre では w(x)=1 であるが, 一般には w(x)\neq 1).

\begin{align*} (p_k,p_l)&=\int_{-1}^1 p_k(x)p_l(x)w(x)dx, \quad w(x)=1 \\ &=\lambda_k\delta_{kl} \\ \lambda_k &= \dfrac{2}{2k+1} \end{align*}

最初のいくつかの Legendre多項式を示すと以下のようになる.

\begin{align*} p_0(x)&=1 \\ p_1(x)&=x \\ p_2(x)&=\dfrac{1}{2}(3x^2-1) \\ p_3(x)&=\dfrac{1}{2}(5x^3-3x) \end{align*}

また選点直交性は以下のとおり.

\begin{align*} <p_k,p_l\!>_n &= \sum_{k=0}^{n-1}a_i p_k(x_i)p_l(x_i) \\ &=\lambda_k\delta_{kl} \\ a_i^{-1}&=\sum_{k=0}^{n-1}\dfrac{p_k(x_i)^2}{\lambda_k} \end{align*}

(ここでも \lambda_k が登場していることに注意.)

まず, Gauss-Legendre の積分公式の積分区間は [-1,1] なので, 積分区間を調整する. t=2x-1 として書き換えると,

\begin{align*} I&=\int_0^1\tilde{f}(x) dx \\ &=\int_{-1}^1 \tilde{f}\Big(\dfrac{t+1}{2}\Big) \dfrac{dt}{2} \\ &=\dfrac{1}{2}\int_{-1}^1 \Big( 0.5+\big(\dfrac{t}{2}+0.48\big)\big(\dfrac{t}{2}+0.25\big)\dfrac{t}{2}\big(\dfrac{t}{2}-0.25\big)\big(\dfrac{t}{2}-0.48\big) \Big) dt \\ &=\int_{-1}^1 \Big( \dfrac{1}{4}+\dfrac{t}{4}\big(\dfrac{t^2}{4}-0.48^2\big)\big(\dfrac{t^2}{4}-0.25^2\big) \Big) dt \end{align*}

ここで改めて積分変数をtからxに変更して, 以下の積分を得る.

I=\int_{-1}^1 \Big( \dfrac{1}{4}+\dfrac{x}{4}\big(\dfrac{x^2}{4}-0.48^2\big)\big(\dfrac{x^2}{4}-0.25^2\big) \Big) dx

後の記載のため, 被積分関数を改めて f(x) としておく.

\begin{align*} f(x) &:= \dfrac{1}{2}\tilde{f}\Big(\dfrac{x+1}{2}\Big) \\ &=\dfrac{1}{4}+\dfrac{x}{4}\big(\dfrac{x^2}{4}-0.48^2\big)\big(\dfrac{x^2}{4}-0.25^2\big) \end{align*}

念の為, f(x) もプロットしておくと以下のとおり(\tilde{f}(x) とグラフの形は同じだが, 縦軸横軸の値が異なっている点に注意).

2.1.1 解析解

まずは解析的に積分を計算してみよう. I の第2項は奇関数であるから積分には寄与しない. したがって,

\begin{align*} I&=\int_{-1}^1 \dfrac{1}{4} dx \\ &=\dfrac{1}{2} \end{align*}

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型積分公式では積分を以下の形で計算する.

I\sim \sum_{i=1}^n a_i f(x_i)

ここで, x_i は標本点, a_i はウェイトであり, 被積分関数 f とは無関係に決まる. 既存の Gauss型積分のライブラリでは, これらのパラメータを事前に用意しておくことにより, 積分計算の高速化を図っている. ここではこのようなライブラリの一つである FastGaussQuadrature.jl を使って上の積分を計算してみよう(コードは抜粋のみ, 全体は記事の最後に示す).

FastGaussQuadrature.jl を使った積分の計算
prg.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 に積分の値を求めてみる. まずは n=3 のときの零点を求める. p_3(x)=0 より,

\begin{align*} p_3(x)&=\dfrac{5x^3-3x}{2}=0\\ x&=0,\pm\sqrt{\dfrac{3}{5}} \end{align*}

これを x_1=-\sqrt{\dfrac{3}{5}}=-0.77459\cdots, x_2=0, x_3=\sqrt{\dfrac{3}{5}}=0.77459\cdots としよう. これらの零点に関する関数値は, 計算が大変なので, いったん, f(x_1), f(x_2), f(x_3) とする.

julia では以下のようなイメージになる.

標本点を求める
prg.jl(抜粋)
# 例えば n=3
x=basis.(Legendre,n) |> roots |> real

次に n-1次以下の Legendre 多項式(p_k(x))を用意しよう. ついでに内積の値\lambda_k = (p_k,p_k) も計算しておく. これらも以下のようなイメージで準備できる.

Legendre多項式と内積の係数を求める
prg.jl(抜粋)
# 例えば n=3
p=basis.(Legendre,0:(n-1))
λ=[SP.innerproduct(Legendre,p[k+1],p[k+1]) for k in 0:(n-1)]

これで, x_i, p_k(x), \lambda_k がそろったのでウェイト a_i=\{\sum_{k=0}^{n-1}p_k(x_i)^2/\lambda_k\}^{-1} を計算できる.

積分のウェイトを求める
prg.jl(抜粋)
# 手作業での計算: 必要なパラメータの算出
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

積分の計算には必要ないが近似関数も求めておこう. n=3 なので, f_n(x)=c_0 + c_1 p_1(x)+c_2 p_2(x) とし, c_k=1/\lambda_k \cdot \sum_{i=1}^{n}a_i p_k(x_i) f(x_i) である.

近似多項式の係数を求める
prg.jl(抜粋)
# 標本点によって決る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)]

(後の都合で, ここでは係数 c_k が関数 f の関数になるように定義していることに注意.)

ここまでの結果を用いて積分を計算し, 近似関数をプロットしてみよう.

手作業による積分計算
prg.jl(抜粋)
# 手作り: 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 を使用したときと(ほぼ)同じになっていることがわかる. また, 近似関数をプロットすると以下のとおり.

プロットを見るとお世辞にも近似されているとは言えない形をしている. n=3 であるので, n-1=2次関数になるはずであるが, 被積分関数が定数項を除き奇関数であるため, 偶数次の係数が0になっているだろうということは想像できる. 実際係数を確認してみると, 2次の係数は(ほぼ)0になっている.

ulia> fn_coef(3)(fifth)
3-element Vector{Float64}:
  0.25
 -0.0017587499999999895
 -6.93889390390723e-17

ちなみに, 被積分関数を p_3(x) で割り算してみると,

(手計算では厳しいので, WolframAlphaを利用した. 商と余りではなぜか勝手に全体が2倍されている.)

これからわかるように, 上で求めた近似曲線は被積分関数を p_3(x) で割った余りになっている. これは積分精度の部分で述べたとおりである.

5次関数に対する近似曲線が直線, という一見直観に反する結果になっているが, 積分の近似という意味ではこれで正しく近似できている.

2.2 コーシー・ローレンツ関数の積分

2つ目の例として, ルンゲの現象のときに使用した以下の関数(コーシー・ローレンツ関数)を考えよう.

\begin{align*} g(x)&= \dfrac{1}{1+25x^2}\\ I & = \int_{-1}^{1}g(x) dx \end{align*}

積分区間は最初から [-1,1] で考える.

2.2.1 解析解

この関数も解析的に積分することができる. y=5x とおき, さらに y=\tan\theta とおく.

\begin{align*} I&=\int_{-5}^{5}\dfrac{1}{1+y^2}\dfrac{dy}{5}\\ &=\dfrac{1}{5}\int_{\tan^{-1}(-5)}^{\tan^{-1}(5)}\dfrac{1}{1+\tan^2\theta}\dfrac{d\theta}{\cos^2\theta}\\ &=\dfrac{1}{5}\int_{\tan^{-1}(-5)}^{\tan^{-1}(5)}d\theta \\ &=\dfrac{1}{5}(\tan^{-1}(5)-\tan^{-1}(-5))\\ &=0.54936\cdots \end{align*}

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積分の方がはるかに早く真値に収束していることがわかる. 標本数が n=2^7=128 程度でマシンイプシロン程度(10^{-16})まで誤差が小さくなっている.

3. サンプルプログラム

上の記事で計算や図のプロットに使用したプログラム.

サンプルプログラム
prg.jl
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次の直交多項式の零点)とウェイト(上の記載ではa_i)は被積分関数とは関係なく事前に計算できる(FastGaussQuadrature.jl などでは事前に大量に準備されており, 1億個の標本点でも数秒で計算が終わるらしい).

一方, デメリットとしては以下のような点が挙げられる.

  • 直交多項式の零点とウェイトを外部のライブラリに依存する場合, その正確性が不安
  • Newton-Cotes公式と比較するとロジックがわかりにくい

1つ目のデメリットは実際に過去に誤植により多くのリソースが無駄になったということがあったらしい(数値解析の p.94 参照). 最近ではいわゆるリーナスの法則により間違いが放置されることは少ないと思われるが, このようなリスクがあることは心にとめておくべきだろう.

2つ目については, 本記事が少しでも役に立てば幸いである. ただ, 自分でこのような記事を書いていてなんですが, どこかのSF作家が言っていたように,

十分に発達した数学は魔法と見分けがつかない

という印象はぬぐえません(これ気に入ったので座右の銘にします(笑)).

Discussion