概要
周長 L
長半径 a
短半径 b
の楕円で、Lとaまたはbが分かっているときに、残された半径を求める。
環境
python 3.10.6
sympy==1.11.1
numpy==1.21.5
科学技術計算ライブラリであるscipyをインストールしておく
pip install scipy
コード
import scipy as sp
import numpy as np
def get_ellipse_another_radius(L,r):
"""
楕円の周長と半径rから残されたほうの半径を求める
L:楕円の周長
r:長半径または短半径
"""
k = 0.9 # 離心率初期値
dk = 0.001 # デルタk
# rが長半径か短半径か判定する
if 2*np.pi*r > L:
# rが長半径の場合
a = r
# 楕円の周の長さの公式
def F(k):return 4*a*sp.special.ellipe(k**2)-L
# F(k)の微分
def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
# ニュートン・ラプソン法で離心率を求める
for i in range(10):
k = k - F(k)/F_diff(k)
# 離心率の式から短半径を求める
b = a*np.lib.scimath.sqrt(1-k**2)
return b
else:
# rが短半径の場合
b = r
# 楕円の周の長さの公式
def F(k):return 4*b/np.lib.scimath.sqrt(1-k**2)*sp.special.ellipe(k**2)-L
# F(k)の微分
def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
# ニュートン・ラプソン法で離心率を求める
for i in range(10):
k = k - F(k)/F_diff(k)
# 離心率の式から長半径を求める
a = b/np.lib.scimath.sqrt(1-k**2)
return a
解説
楕円の基本公式
第二種完全楕円積分\\
E(k)=\int_{0}^{\frac{\pi}{2}}\sqrt{1-k^2\sin^2t}dt\\
楕円の周長の公式\\
L=4aE(k)\\
離心率\\
k=\sqrt{1-\frac{b^2}{a^2}}\hspace{20px},\hspace{20px}(b<a)
一方の半径rが長半径aか短半径bかで計算が変わるため、まずはrが長半径か短半径かを判定する。
それには下の図を見ると分かりやすい。
2πr > Lのときはrが長半径である。
まずはrが長半径の場合について求める。
rが長半径の場合
# rが長半径の場合
a = r
離心率は変形すると
b=a\sqrt{1-k^2}
# 離心率の式から短半径を求める
b = a*np.lib.scimath.sqrt(1-k**2)
となり、aの値は分かっているので離心率kが分かればbが求まる。
離心率kを求めるにはニュートン・ラフソン法を用いる
ニュートン・ラフソン法\\
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
詳しくは解説しないが、f(x)=0のときのxの近似値を反復計算によって得ることができる手法である。
楕円の周長の公式を次のようにする
f(k)=4aE(k)-L=0
# 楕円の周の長さの公式
def F(k):return 4*a*sp.special.ellipe(k**2)-L
第二種完全楕円積分はscipy.special.ellipeで計算できる。※m=k^2 となっているので、k^2 を代入するところだけは注意
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html
離心率kをニュートン・ラフソン法に当てはめると下記になる。
k_{n+1} = k_n-\frac{f(k_n)}{f'(k_n)}\\
k_0 = 0.9とする\\
計算を繰り返していくと離心率kの近似値が求まる。
今回は１０回繰り返してみる
k = 0.9 # 離心率初期値
dk = 0.001 # デルタk
# F(k)の微分
def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
# ニュートン・ラフロソン法で離心率を求める
for i in range(10):
k = k - F(k)/F_diff(k)
kの近似値が得られたのでbを求めることができる
# 離心率の式から短半径を求める
b = a*np.lib.scimath.sqrt(1-k**2)
rが短半径の場合
# rが長半径の場合
b = r
離心率を変形して
a=\frac{b}{\sqrt{1-k^2}}
aを、楕円の周長の公式に代入して長半径の場合と同じようにする
f(k)=\frac{4b}{\sqrt{1-k^2}}E(k)-L=0
あとは同じくニュートン・ラフソン法を用いて離心率kを計算し、長半径aが求まる。
以下省略
まとめ
楕円の周長と一方の半径rから残された半径を求める方法について記述しました。
- rが長半径なのか短半径なのか判別する。
- ニュートン・ラフソン法に楕円の周長の公式を当てはめて、離心率kの近似値を求める。
- 離心率の公式からもう一方の半径を求める。