楕円の周の長さの求め方
1. 理論などどうでもよいという場合
タイトルは不穏当だが，要するに，パッケージなりを使用して，簡単に答えを得たいという場合である。
ChatGPT3.5 に聞いてみると，「Julia には Elliptic に ellipke，Python には scipy.special に elliptic_perimeter_with_ellipkがあるよ」とご丁寧にもプログラムを添えて教えてくれた。
ただ，あとになってわかったが，あいも変わらず ChatGPT3.5 は不十分な答えを返してくれたのであった。
計算結果を見てのけぞった。別の計算手段から得られた数値よりずっと大きい値が返ってくる。
私の聞き方も悪かったのだが，最初は「第二種楕円積分を使って...」で，次に「第二種完全楕円積分を使って...」と聞くと，同じ関数を引数として離心率と離心率の二乗を使ったプログラム例を提示してきた。そして，どちらも結果が異なる。
更にあとになってわかったが，ChatGPT3.5 が答えたのは第一種楕円積分と第一種完全楕円積分の方だった。
Julia では同じ関数名 ellipke が第一種と第二種の結果を返す（引数として離心率の二乗を与えるのが「完全ではないほう」，離心率を与えるのが「完全なほう」）。
Python では第一種のほうは ellipe，第二種のほうは ellipk と名前が違う。
- 第二種楕円積分による計算。
# import Pkg
# Pkg.add("Elliptic")
using Elliptic
function elliptic_circumference(a, b)
m = 1.0 - (b / a)^2 # 楕円の離心率を計算
C = 4 * a * ellipke(m^2)[2] # 楕円の周長を計算(第二種楕円積分）
return C
end
elliptic_circumference(2.5, 1.5) |> println
elliptic_circumference(0.5, 0.5) |> println
elliptic_circumference(2, 1) |> println
elliptic_circumference(1, 0.5) |> println
13.9484180877679
3.1415926535897936
10.547776863956969
5.2738884319784844
- 第二種完全楕円積分による計算
# import Pkg
# Pkg.add("Elliptic")
using Elliptic
function elliptic_circumference(a, b)
m = 1.0 - (b / a)^2 # 楕円の離心率を計算
C = 4 * a * ellipke(m)[2] # 楕円の周長を計算(第二種完全楕円積分）
return C
end
elliptic_circumference(2.5, 1.5) |> println
elliptic_circumference(0.5, 0.5) |> println
elliptic_circumference(2, 1) |> println
elliptic_circumference(1, 0.5) |> println
12.763499431699064
3.1415926535897936
9.688448220547674
4.844224110273837
2. 区分求積法
高校数学の美しい物語 で解説されている，区分求積法により求める方法である。
- 第二種楕円積分による計算
「完全ではない」ということであるが，とりあえず。
また，サイトによっては cos() ではなく sin() が使われていたりする（?）ので，注意が必要か。
using SymPy
@syms a, b, e, θ
(a, b) = (2.5, 1.5)
e = (a^2 - b^2)/a^2
4a*integrate(sqrt(1 - e^2*cos(θ)^2), (θ, 0, PI/2)).evalf()
$13.9484180877679$
- 第二種完全楕円積分による計算
k を使って定義されることもあるが，要するに cos() の係数が e^2 ではなく e である。
using SymPy
@syms a, b, e, θ
(a, b) = (2.5, 1.5)
#k = sqrt((a^2 - b^2)/a^2)
#4a*integrate(sqrt(1 - k^2*cos(θ)^2), (θ, 0, PI/2)).evalf()
e = (a^2 - b^2)/a^2
4a*integrate(sqrt(1 - e*cos(θ)^2), (θ, 0, PI/2)).evalf()
# 12.7634994316990631 後に二重階乗を使って逐次計算した値と同じになる。
$12.7634994316991$
- 離心率を使わないで書くこともできる
同じことをやっている。
using SymPy
@syms a, b, θ
(a, b) = (2.5, 1.5)
4integrate(sqrt(a^2*cos(θ)^2 + b^2*sin(θ)^2), (θ, 0, PI/2)).evalf()
$12.7634994316991$
3. 逐次計算
次に，区分求積法による計算は時間がかかることがあるので，高校数学の美しい物語 で最初に解説されているアルゴリズムに基づいてプログラムを書いた。
二重階乗も Combinatorics のものを使う。効率的に収束計算を行う（打ち切る）ために（?）小細工をしているが，最近のコンピュータなら不要であろう。
using Combinatorics
function elliptic_circumference1(a, b; maxiter=500, epsilon=1e-20)
c(t::Integer) = t == 0 ? 1 : doublefactorial(2t - 1)/doublefactorial(2t)
ϵ = sqrt(1 - b^2/a^2)
s = big"0"
for t = 0:max(10, maxiter)
delta = c(t)^2*ϵ^(2t)/(1 - 2t)
abs(delta) < epsilon && break
s += delta
end
s*a*2pi
end;
elliptic_circumference1(2.5, 1.5, maxiter=70) |> println
elliptic_circumference1(0.5, 0.5, maxiter=10) |> println
elliptic_circumference1(2, 1, maxiter=70) |> println
elliptic_circumference1(1, 0.5) |> println
12.76349943169906315639125132370200335575162431876546420693435001200081018193559
3.141592653589793115997963468544185161590576171875
9.688448220549662535946252084867017448101001516662105882180110266650046045936008
4.844224110273838165479726408797042955563087446355597419855518805638773196072963
4. 近似計算
4.1. Gauss-Kummer の公式
高校数学の美しい物語 に式が示されている。
楕円の離心率が 0 に近い場合に使える（円の離心率は 0 である）。
function elliptic_circumference2(a, b)
h = (a - b)/(a + b)
return pi*(a + b)*(1 + (h/2)^2 + h^4/64 + h^6/256)
end;
elliptic_circumference2(2.5, 1.5) |> println
elliptic_circumference2(0.5, 0.5) |> println
elliptic_circumference2(2, 1) |> println
elliptic_circumference2(1, 0.5) |> println
12.763499129827382
3.141592653589793
9.688445901297804
4.844222950648902
4.2. シュリニヴァーサ・アイヤンガー・ラマヌジャンの近似式
elliptic_circumference3(a, b) = π*(3(a + b) - sqrt((3a + b)*(a + 3b)))
elliptic_circumference3(2.5, 1.5) |> println
elliptic_circumference3(0.5, 0.5) |> println
elliptic_circumference3(2, 1) |> println
elliptic_circumference3(1, 0.5) |> println
12.763493196879272
3.141592653589793
9.688421097671288
4.844210548835644
4.3. シュリニヴァーサ・ラマヌジャンの近似式（2）
もう少し良い近似式
elliptic_circumference4(a, b) = π * (a + b)*(1 + 3*((a - b)/(a + b))^2/(10 + sqrt(4 - 3*((a - b)/(a + b))^2)))
elliptic_circumference4(2.5, 1.5) |> println
elliptic_circumference4(0.5, 0.5) |> println
elliptic_circumference4(2, 1) |> println
elliptic_circumference4(1, 0.5) |> println
12.763499431394381
3.141592653589793
9.688448216130086
4.844224108065043
4.4. 独自に（?）考案された近似式
function elliptic_circumference5(a, b)
α = (a + b)/2
β = sqrt(a*b)
# return pi*sqrt(2a^2 - b^2)(15015a^12 - 45045a^10*b^2 + 56761a8*b4 - 38447a^6*b^6 + 14724a^4*b^8 - 3008a^2*b^10 + 256b^12)/128(-2a^2 + b^2)^6
return pi*sqrt(2α^2 - β^2)*(15015α^12 - 45045α^10*β^2 + 56761α^8*β^4 - 38447α^6*β^6 + 14724α^4*β^8 - 3008α^2*β^10 + 256β^12)/(128(-2α^2 + β^2)^6)
end;
elliptic_circumference5(2.5, 1.5) |> println
elliptic_circumference5(0.5, 0.5) |> println
elliptic_circumference5(2, 1) |> println
elliptic_circumference5(1, 0.5) |> println
12.763629588373758
3.141592653589793
9.689229354206223
4.844614677103111
4.5. 関孝和の近似公式
離心率が 0 〜 0.9 の範囲では相対誤差が 1% 未満。
elliptic_circumference6(a, b) = 2sqrt(4(a - b)^2 + pi^2*a*b)
elliptic_circumference6(2.5, 1.5) |> println
elliptic_circumference6(0.5, 0.5) |> println
elliptic_circumference6(2, 1) |> println
elliptic_circumference6(1, 0.5) |> println
12.807968848195266
3.141592653589793
9.744579786153679
4.8722898930768395
4.6. 詳細不明の公式
あまり良い近似公式ではない。
elliptic_circumference7(a, b) = 2pi*sqrt((a^2 + b^2)/2)
elliptic_circumference7(2.5, 1.5) |> println
elliptic_circumference7(0.5, 0.5) |> println
elliptic_circumference7(2, 1) |> println
elliptic_circumference7(1, 0.5) |> println
12.95311834341519
3.141592653589793
9.934588265796101
4.967294132898051
5. 性能評価
離心率の異なる楕円の周長をそれぞれの方法によって求めて，違いを観察してみる。
methods = [elliptic_circumference1,
elliptic_circumference2,
elliptic_circumference3,
elliptic_circumference4,
elliptic_circumference5,
elliptic_circumference6,
]
m = length(methods)
bs = 0.01:0.01:1
n = length(bs)
e = @. sqrt(1 - bs^2)
cf = Array{Float64}(undef, (n, m));
for (i, b) in enumerate(bs)
for (j, f) in enumerate(methods)
cf[i, j] = f(1, b)
end
end
using Plots
pyplot(size=(500, 375), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(xlabel="短半径 b（長半径 a = 1）", ylabel="周長")
for j in 1:m
plot!(bs, cf[:, j]./cf[:, 1], label=false)
annotate!(-0.01, cf[1, j]./cf[1, 1], text(string(j), 10, :red, :right, :vcenter))
end
plot!(xlims=(-0.04, 1.01))
6．オンラインで計算できるサイト
いくつかの計算サイトで， a = 2.5, b = 1/5 のときの結果を計算してみた
-
楕円の円周 - 高精度計算サイト
使用公式不明
12.76348224780646841536 -
楕円の周の長さの計算
シュリニヴァーサ・アイヤンガー・ラマヌジャン
12.763 -
楕円円周電卓
2pi*sqrt((a^2 + b^2)/2) あまりよい近似公式ではない
12.95311834510651698971352818625473397176450860978523886454709379329709676760229789599996521252061537 -
楕円の面積・周長・楕円率・離心率
ラマヌジャンの近似式
12.7634931969 -
WolframAlpha
4integrate(sqrt(2.5^2cos(θ)^2 + 1.5^2sin(θ)^2), (θ, 0,pi/2))
12.763499431699064233089331002495145695979749424396093278642581743... -
楕円の面積と周長
ラマヌジャンの近似式 6次まで
12.763
おまけに
Appendix. 二重階乗
自分で定義してもよいが，Conbinatorics パッケージに doublefactorial(n::Integer) が用意されているので，それを使う。
function my_doublefactorial(number::Integer)
fact = big"1"
for m in iseven(number)+1:2:number
fact *= m
end
return fact
end;
my_doublefactorial.([0, 1, 2, 3,4, 5, 6, 10, 100])
9-element Vector{BigInt}:
1
1
2
3
8
15
48
3840
34243224702511976248246432895208185975118675053719198827915654463488000000000000
using Combinatorics
doublefactorial.([0, 1, 2, 3,4, 5, 6, 10, 100])
9-element Vector{BigInt}:
1
1
2
3
8
15
48
3840
34243224702511976248246432895208185975118675053719198827915654463488000000000000
Appendix.2.
Python による計算例
import numpy as np
import scipy.special
def elliptic_perimeter(a, b):
e = (1.0 - (b / a) ** 2)
# 第二種楕円積分を計算
return 4 * a * scipy.special.ellipe(e**2)
print(elliptic_perimeter(2.5, 1.5))
print(elliptic_perimeter(0.5, 0.5))
print(elliptic_perimeter(2, 1))
print(elliptic_perimeter(1, 0.5))
13.9484180877679
3.141592653589793
10.547776863956969
5.2738884319784844
import numpy as np
def elliptic_perimeter2(a, b):
e = (1.0 - (b / a) ** 2)
# 第二種完全楕円積分を計算
return 4 * a * scipy.special.ellipe(e)
print(elliptic_perimeter2(2.5, 1.5))
print(elliptic_perimeter2(0.5, 0.5))
print(elliptic_perimeter2(2, 1))
print(elliptic_perimeter2(1, 0.5))
12.763499431699065
3.141592653589793
9.688448220547675
4.844224110273838