概要

周長 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の近似値を求める。
離心率の公式からもう一方の半径を求める。

参考文献