浅葱色の計算用紙

数学(広義)を扱っています。

離心率の一般化

突然ですが、円の曲率を求めるにはどうすればよいでしょうか?

半径の逆数を取ればよいですね。簡単です。

 

では、円とは限らない曲線の曲率を求めるにはどうすればよいでしょうか?

曲線上の微小距離離れた3点を取って、その3点を通る円の曲率を求めればよいですね。いわゆる微分の考え方です。

 

 

突然ですが、二次曲線の離心率を求めるにはどうすればよいでしょうか?

曲線を回転させて、標準形に変換して、数Ⅲで習った離心率の公式に代入すればよいですね。簡単です。

 

では、二次曲線とは限らない曲線の離心率を求めるにはどうすればよいでしょうか?

5点あれば二次曲線を決定できるので、曲線上の微小距離離れた5点を取って、その5点を通る二次曲線の離心率を求めればよいですね。いわゆる微分の考え方です。

しかし、私が調べた限り誰もこのことに気づいていませんでした。

 

 

そこで、私はこのことをpdfに書いてTwitterで公開しました。

もちろん、できるだけ多くの人に知ってほしいので英語で書きました。

こちらがそのpdfです。

 

drive.google.com

と言われても、英語だとわかりにくいですね。

そこで、この記事では、上のpdfに基づいてその内容を解説しようと思います。

 

上のpdfでは、任意の2次元曲線に対し、その離心率を定義し、\( y=f(x) \)に対する離心率を近似するプログラムを示し、実際にそのプログラムで\( y=x^3-3x \)の離心率を求めた結果を示しています。

 

定義

\( C \)を2次元曲線、\( P \)をその上の点とする。このとき、\( C \)の\( P \)における離心率\( e(C,P) \)を次で定義する:

\( C \)上に\( N, O, P, Q, R \)がこの順に並ぶように4点\( N, O, P, Q \)をとる。このとき、\( N, O, P, Q, R \)を通る二次曲線の離心率を\( E \)とする。\( N, O, P, Q \)を独立に\( P \)に近づけたとき、\( E \)がある値に収束するならば、その値を\( e(C,P) \)とする。

プログラム

説明はコメントでしているので、改めてここで説明する必要はないと思います。

# -*- coding: utf-8 -*-

import numpy
import sympy
import math
import matplotlib.pyplot as plt

# finds the eccentricity of the conic curve that goes through all of the following:
# (x-2e, f(x-2e)), (x-e, f(x-e)), (x, f(x)), (x+e, f(x+e)), (x+2e, f(x+2e))
def eccentricity(f, x, e):
    # Solve for the coefficients
    left = []
    right = []
    for i in [-2,-1,0,2,4]:
        local_x = x + i * e
        local_y = f(x + i * e)
        left_eq = [local_x ** 2 , local_x * local_y, local_y ** 2,
                   local_x, local_y, 1]
        left.append(left_eq)
        right.append(0)
    left.append([1,0,0,0,0,0])
    right.append(1)
    
    #print(left)
    #print(right)
    coef = numpy.linalg.solve(left, right)
    #print(coef)
   
    # Diagonalize the quadratic form
    q_form = numpy.array([[coef[0], coef[1]/2],[coef[1]/2, coef[2]]])
    # p is already orthogonal
    l, p = numpy.linalg.eig(q_form)
    pinv = numpy.linalg.inv(p)

    #print(q_form)
    #print(p)
    #print(pinv)
    #print(numpy.dot(numpy.dot(pinv, q_form), p))
    
    x1 = sympy.Symbol('x1')
    y1 = sympy.Symbol('y1')
    x2 = sympy.Symbol('x2')
    y2 = sympy.Symbol('y2')
    
    curve = numpy.dot(coef, [x1**2, x1*y1, y1**2, x1, y1, 1])
    curve_straight = curve.subs([(x1, pinv[0, 0]*x2+pinv[1, 0]*y2),
                                 (y1, pinv[0, 1]*x2+pinv[1, 1]*y2)])
    curve_straight = sympy.expand(curve_straight)
    #print(curve)
    #print(curve_straight)
    
    # Change it to "standard" form
    # s(x-a)^2 + t(y-b)^2 = 1
    x3 = sympy.Symbol('x3')
    y3 = sympy.Symbol('y3')
    x2_sq = sympy.diff(curve_straight, x2, 2).subs([(x2, 0), (y2, 0)]) / 2
    x2_li = sympy.diff(curve_straight, x2, 1).subs([(x2, 0), (y2, 0)])
    y2_sq = sympy.diff(curve_straight, y2, 2).subs([(x2, 0), (y2, 0)]) / 2
    y2_li = sympy.diff(curve_straight, y2, 1).subs([(x2, 0), (y2, 0)])
    
    #print(x2_sq, x2_li, y2_sq, y2_li)
    
    if x2_sq != 0 and y2_sq != 0: # ellipse, two lines, or hyperbola
        curve_standard = curve_straight.subs([(x2, x3 - x2_li / (x2_sq * 2)),
                                              (y2, y3 - y2_li / (y2_sq * 2))])
        curve_standard = sympy.expand(curve_standard)
        #print(curve_standard)
        
        c = curve_standard.subs([(x3, 0), (y3, 0)])
        
        if c == 0: # two lines
            return float('inf')
        else:
            curve_standard /= -c
            #print(curve_standard)
            x3_sq = sympy.diff(curve_standard, x3, 2).subs([(x3, 0), (y3, 0)]) / 2
            y3_sq = sympy.diff(curve_standard, y3, 2).subs([(x3, 0), (y3, 0)]) / 2
            #print(y3_sq, x3_sq)
            if x3_sq < 0 and y3_sq < 0: # impossible
                return float('nan')
            elif x3_sq < 0 and y3_sq > 0: # upways hyperbola
                return math.sqrt(1 - y3_sq / x3_sq)
            elif x3_sq > 0 and y3_sq < 0: # sideways hyperbola
                return math.sqrt(1 - x3_sq / y3_sq)
            else: # ellipse
                if x3_sq < y3_sq:
                    return math.sqrt(1 - x3_sq / y3_sq)
                else:
                    return math.sqrt(1 - y3_sq / x3_sq)
            
    elif x2_sq == 0 and y2_sq == 0: # line
        return float('inf')
    elif y2_sq == 0: # laid parabola
        return 1
    else: # stand parabola
        return 1

x = numpy.arange(-3, 3, 0.05)
y = []

for a in x:
    e = eccentricity(lambda x: x**3-3*x, a, 0.001)
    y.append(e)
    print(a, e)

y = numpy.array(y)

plt.plot(x, y)
plt.ylim([0,5])
plt.show()

実行結果

\( y = x^3-3x \)の、\( -3 \leq x \leq 3\)での曲率のグラフを以下に示します。

f:id:asangi_a4ac:20200219150040p:plain

y=x^3-3xの離心率のグラフ

\( x=0 \)の近くで離心率が大きくなっていますが、理由は分かりません。

\( y=x^3-3x \)のグラフは原点に関して対称で、原点の周囲の5点を通る二次曲線は重なった直線になるので、それが原因かもしれません。

これに関する調査は今後の課題とします。

 

名称

普通の(二次曲線の)離心率以外にも、何種類か離心率があるようなので、この離心率を「Ασαγγη離心率」と呼んで区別することにしました。

ギリシャ文字が入力しづらい環境のために、次のような代替表記を推奨します(上が最も推奨):

  • Asangi離心率
  • 浅葱離心率
  • アサキ゜離心率
  • アサギ離心率

しかし、次の表記は推奨されません:

  • Asagi離心率
  • アサンギ離心率

ラテン文字表記と片仮名表記で、/n/を表す文字が入るか入らないかが逆転していることに注意してください。