楕円の周長の求め方

## 楕円の周の長さの求め方


## 1. 理論などどうでもよいという場合

タイトルは不穏当だが，要するに，パッケージなりを使用して，簡単に答えを得たいという場合である。

ChatGPT3.5 に聞いてみると，「Julia には Elliptic に ellipke，Python には scipy.special に elliptic_perimeter_with_ellipkがあるよ」とご丁寧にもプログラムを添えて教えてくれた。

ただ，あとになってわかったが，あいも変わらず ChatGPT3.5 は不十分な答えを返してくれたのであった。
計算結果を見てのけぞった。別の計算手段から得られた数値よりずっと大きい値が返ってくる。

私の聞き方も悪かったのだが，最初は「第二種楕円積分を使って...」で，次に「第二種完全楕円積分を使って...」と聞くと，同じ関数を引数として離心率と離心率の二乗を使ったプログラム例を提示してきた。そして，どちらも結果が異なる。

更にあとになってわかったが，ChatGPT3.5 が答えたのは第一種楕円積分と第一種完全楕円積分の方だった。

Julia では同じ関数名 ellipke が第一種と第二種の結果を返す（引数として離心率の二乗を与えるのが「完全ではないほう」，離心率を与えるのが「完全なほう」）。

Python では第一種のほうは ellipe，第二種のほうは ellipk と名前が違う。

* 第二種楕円積分による計算。


```julia
# 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


* 第二種完全楕円積分による計算


```julia
# 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. 区分求積法

[高校数学の美しい物語](https://manabitimes.jp/math/1092) で解説されている，区分求積法により求める方法である。

* 第二種楕円積分による計算

「完全ではない」ということであるが，とりあえず。
また，サイトによっては cos() ではなく sin() が使われていたりする（?）ので，注意が必要か。


```julia
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 である。


```julia
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$



* 離心率を使わないで書くこともできる

[楕円の周長の計算（区分求積法）](https://cvtech.cc/ellipse/)

同じことをやっている。


```julia
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. 逐次計算

次に，区分求積法による計算は時間がかかることがあるので，[高校数学の美しい物語](https://manabitimes.jp/math/1092) で最初に解説されているアルゴリズムに基づいてプログラムを書いた。

二重階乗も Combinatorics のものを使う。効率的に収束計算を行う（打ち切る）ために（?）小細工をしているが，最近のコンピュータなら不要であろう。


```julia
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 の公式

[高校数学の美しい物語](https://manabitimes.jp/math/1092) に式が示されている。

楕円の離心率が 0 に近い場合に使える（円の離心率は 0 である）。


```julia
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. シュリニヴァーサ・アイヤンガー・ラマヌジャンの近似式

[円周の公式（面積・円周）｜数学](http://excel-mania.com/math/t022.html)


```julia
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）

もう少し良い近似式

[楕円 - Wikipedia](https://ja.wikipedia.org/wiki/%E6%A5%95%E5%86%86)


```julia
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. 独自に（?）考案された近似式
[楕円の周長の近似式をあぶり出す](https://hyperion64.hatenadiary.org/entry/20160418/p1)


```julia
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. 関孝和の近似公式

[関孝和の楕円周を求める近似式](https://www.tcp-ip.or.jp/~n01/math/analysis/seki/seki.pdf)

離心率が 0 〜 0.9 の範囲では相対誤差が 1% 未満。


```julia
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. 詳細不明の公式

あまり良い近似公式ではない。


```julia
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. 性能評価

離心率の異なる楕円の周長をそれぞれの方法によって求めて，違いを観察してみる。



```julia
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))
```


![output_27_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/ba44ed47-2612-8c3c-0348-397d5d3691ef.png)


## 6．オンラインで計算できるサイト

いくつかの計算サイトで， a = 2.5, b = 1/5 のときの結果を計算してみた

* [楕円の円周 - 高精度計算サイト](https://keisan.casio.jp/exec/user/1495672323)
使用公式不明
12.76348224780646841536

* [楕円の周の長さの計算](http://k-ichikawa.blog.enjoy.jp/etc/HP/js/Q3/ellipse_2D.html)
シュリニヴァーサ・アイヤンガー・ラマヌジャン
12.763

* [楕円円周電卓](https://miniwebtool.com/ja/ellipse-circumference-calculator/)
2pi*sqrt((a^2 + b^2)/2)  あまりよい近似公式ではない
12.95311834510651698971352818625473397176450860978523886454709379329709676760229789599996521252061537

* [楕円の面積・周長・楕円率・離心率](https://www.math-joy-life.com/daen-area/)
ラマヌジャンの近似式
12.7634931969

* [WolframAlpha](https://www.wolframalpha.com/input?i=4integrate%28sqrt%282.5%5E2*cos%28%CE%B8%29%5E2+%2B+1.5%5E2*sin%28%CE%B8%29%5E2%29%2C+%28%CE%B8%2C+0%2Cpi%2F2%29%29&lang=ja)
4integrate(sqrt(2.5^2*cos(θ)^2 + 1.5^2*sin(θ)^2), (θ, 0,pi/2))
12.763499431699064233089331002495145695979749424396093278642581743...

* [楕円の面積と周長](https://yanban.sakura.ne.jp/cgi-bin/daende.cgi)
ラマヌジャンの近似式 6次まで
12.763

---

おまけに

## Appendix. 二重階乗
自分で定義してもよいが，Conbinatorics パッケージに doublefactorial(n::Integer) が用意されているので，それを使う。


```julia
function my_doublefactorial(number::Integer)
    fact = big"1"
    for m in iseven(number)+1:2:number
        fact *= m
    end
    return fact
end;
```





```julia
my_doublefactorial.([0, 1, 2, 3,4, 5, 6, 10, 100])
```




    9-element Vector{BigInt}:
                                                                                    1
                                                                                    1
                                                                                    2
                                                                                    3
                                                                                    8
                                                                                   15
                                                                                   48
                                                                                 3840
     34243224702511976248246432895208185975118675053719198827915654463488000000000000




```julia
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 による計算例

```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



```python
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