Julia の勉強：数値積分

QuadGK パッケージを使うと簡単に関数の定積分を求めることができます。

このパッケージには３つの関数しかないので使い方は簡単です。主に使うのは次の関数です。

quadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm)

• f は積分を行う関数です。
• a, b は積分を行う範囲です。c 以降は追加の範囲 (b～c, c～d, …) です。関数 f の変化の激しい区間は幅を狭めることでより正確な値を得ることができます。
• 残りは省略可能ですが、rtol は相対誤差の許容値、atol は絶対誤差の許容値です。
• maxevals は f の評価最大値です。
• order は積分公式の最大次数です。
• norm が指定されると E <= max(atol, rtol*norm(I)) を満たします。ただし、E は最大誤差です。

この関数の戻り値は積分値 I と最大誤差 E のペア (タプル) になります。

使用例

julia> using QuadGK
julia> f(x) = 1
f (generic function with 1 method)
julia> quadgk(f, 0, 1)
(1.0, 0.0)
julia> f(x) = x
f (generic function with 1 method)
julia> quadgk(f, 0, 1)
(0.5, 0.0)
julia> f(x) = x^2
f (generic function with 1 method)
julia> quadgk(f, 0, 1)
(0.3333333333333333, 5.551115123125783e-17)

楕円積分の例

QuadGK を使って特殊関数である楕円積分のグラフを描いてみます。

# 楕円積分
using QuadGK
using Plots

println("Start ..")
xs = Float64[]
ys = Float64[]

for m = range(-5, 0.9, step=0.01)
    f(x) = 1 / sqrt(1 - m*sin(x)^2)
    I, E = quadgk(f, 0, pi/2)
    append!(xs, m)
    append!(ys, I)
end

plot(xs, ys, show=true, framestyle=:origin, w=2)

print("> ")
read(stdin, Char)