🐈

数値積分 (1次元)

2023/06/23に公開

Julia

tech

初めに

物理の計算において、1, 2...次元空間で定義された関数を数値的に積分する操作が頻繁に現れます。
ここでは、数値積分の仕方は導入してみます。

以下のコードは適当なプロジェクト環境を作って、その中で作業してみましょう。

参考GitHubレポジトリ

プロジェクト環境の準備

GitHub上でレポジトリを先に作成して、git cloneしてから、ローカルで以下のプロジェクト環境を作成するのが楽です。プロジェクト環境を作りたいディレクトリ内で、activateすることで、Project.tomlを作ります。

$ cd ~/workspace
$ git clone git@github.com:shinaoka/numericalintegration.git
$ cd numericalintegration
$ julia
(@v1.9) pkg> activate .
  Activating new project at `~/workspace/numericalintegration`

(numericalintegration) pkg> add FastGaussQuadrature # Gauss-Legendre法
(numericalintegration) pkg> add QuadGK # 自動適応型積分
(numericalintegration) pkg> add LinearAlgebra # 線形代数用標準ライブラリ
(numericalintegration) pkg> add Plots # グラフ描画
$ git add Project.toml
$ git commmit -m "New env"
$ git push

数値積分とは

簡単のため、閉区間での積分を考えます。

\int_a^b f(x) \mathrm{d}x

$f(x)$ が $[a, b]$ の任意の点 $x$ で数値的に評価できるということを前提に、積分の値を数値的に求めことを「数値積分」と呼んでいます。

様々な種類の数値積分法が存在しますが、その多くは、関数の積分値を

\int_a^b f(x) \mathrm{d}x \approx \sum_{n=1}^N w_n f(x_n) = S

と近似します。ここで、 $x_n$ は「注意深く選ばれた」 $N$ 個の関数の評価点 (積分点)、 $w_n$ は対応する重みです (通常は正)。
数値積分の肝は、どのように $x_n$ と $w_n$ を選ぶかにあります。

長方形近似

最も簡単な選択は、

x_n = (n-1) \frac{b-a}{N-1} + a, w_n = (b-a)/N

でしょう。これは不定積分の短冊による近似 (リーマン和)に対応します。
数値積分法の性能は、積分点 $N$ に対してどのように積分値の誤差が減少するか、で測ることができます。このルールの場合、

| S - \int_a^b f(x) \mathrm{d}x | \approx 1/N

のように、非常にゆっくりとしか誤差が減少しないことが知られています (参考: Wikipedia上の記事)。
そのため、精度の悪さが許容できる場合以外、実用的ではありません。

Gauss-Legendre法

上の方法は、被積分関数 $f(x)$ を長方形に近似することに基づいて、積分点と重みを選びました。
一方、Legendre多項式 $P_N(x)$ は、 $N$ 個の零点を持つことが知られています。高精度な積分法としてよく使われるGauss-Legendre法では、積分区間は $[-1, 1]$ にスケールした上で、 $P_N(x)$ の零点を積分点として選びます (重みの選び方の解説は省略)。

詳しい解説は、他の記事にまかせて、ここでは使い方を学びましょう。
なお、Gauss-Legendre法を $2N + 1$ 次までの多項式に適用した場合、厳密な結果を与えます。
また、実用的には、十分滑らかな関数に適用した場合、誤差は指数的に収束します。

FastGaussQuadrature.jlを使って数値積分してみましょう。
具体的には、 $f(x) = x^2$ の0から1までの積分を計算します。厳密な値は $1/3$ です。

using FastGaussQuadrature
using LinearAlgebra

N = 10 # 10点ガウス求積を使用
a = 0.0 # 積分の下限
b = 1.0 # 積分の上限

# 積分点と重みを計算
x, w = gausslegendre(N)

# 積分区間の調整
function adjust_xw(x, w, a, b)
    @assert b > a
    x_adjusted = (b - a) * 0.5 * (x .+ 1) .+ a
    w_adjusted = (0.5 * (b - a)) .* w
    return x_adjusted, w_adjusted
end

# 全てのデータが引数で渡されているか、チェックしておく
@code_warntype adjust_xw(x, w, a, b)

MethodInstance for adjust_xw(::Vector{Float64}, ::Vector{Float64}, ::Float64, ::Float64)
  from adjust_xw(x, w, a, b) @ Main ~/workspace/numericalintegration/gauss_legendre.ipynb:12
Arguments
  #self#::Core.Const(adjust_xw)
  x::Vector{Float64}
  w::Vector{Float64}
  a::Float64
  b::Float64
Locals
  w_adjusted::Vector{Float64}
  x_adjusted::Vector{Float64}
Body::Tuple{Vector{Float64}, Vector{Float64}}
1 ─       Core.NewvarNode(:(w_adjusted))
│         Core.NewvarNode(:(x_adjusted))
│   %3  = (b > a)::Bool
└──       goto #3 if not %3
2 ─       goto #4
3 ─ %6  = Base.AssertionError("b > a")::Core.Const(AssertionError("b > a"))
└──       Base.throw(%6)
4 ┄ %8  = Main.:+::Core.Const(+)
│   %9  = (b - a)::Float64
│   %10 = Base.broadcasted(Main.:+, x, 1)::Core.PartialStruct(Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(+), Tuple{Vector{Float64}, Int64}}, Any[Core.Const(+), Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(1)]), Core.Const(nothing)])
│   %11 = Base.materialize(%10)::Vector{Float64}
│   %12 = (%9 * 0.5 * %11)::Vector{Float64}
│   %13 = Base.broadcasted(%8, %12, a)::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(+), Tuple{Vector{Float64}, Float64}}
│         (x_adjusted = Base.materialize(%13))
│   %15 = Main.:*::Core.Const(*)
│   %16 = (b - a)::Float64
│   %17 = (0.5 * %16)::Float64
│   %18 = Base.broadcasted(%15, %17, w)::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(*), Tuple{Float64, Vector{Float64}}}
│         (w_adjusted = Base.materialize(%18))
│   %20 = Core.tuple(x_adjusted, w_adjusted)::Tuple{Vector{Float64}, Vector{Float64}}
└──       return %20

# 簡単なテストを実行

let # x_, w_が名前空間を汚さないようにスコープを作る
   # Test case: [-1, 1]
   x_, w_ = adjust_xw(x, w, -1.0, 1.0) 
   @assert x_ ≈ x
   @assert w_ ≈ w
   @assert sum(w_) ≈ 2.0

   # Test case: [0, 1]
   x_, w_ = adjust_xw(x, w, 0.0, 1.0) 
   @assert x_ ≈ 0.5 .* x .+ 0.5
   @assert sum(w_) ≈ 1.0
end

# 被積分関数の定義
f(x) = x^2

# 数値積分の計算 (.はbroadcast, dotはベクトル同士の内積)
x_adjusted, w_adjusted = adjust_xw(x, w, a, b)
integral = dot(w_adjusted, f.(x_adjusted))

println("Approximated integral: ", integral)

Approximated integral: 0.3333333333333334

練習問題

積分点数 $N$ を増やしたとき、誤差がどのように減少するか調べて、Plots.jlでプロットしましょう。
ただし、 $f(x)=x^2$ は $N=2$ のGauss-Legendre積分で厳密に評価できてしまうので、関数を変えて遊んでみましょう。

自動適応型積分

Work in progress...

初めに