再急降下法を用いた三角形の推定オイラー線の描写について

# はじめに
三角形の外心と重心と垂心は同一直線上にあることが知られている。これをオイラー線と呼ぶ。一方で、外心と垂心は、三角形の頂点の座標が与えられたとしても、厳密な値を求めることは難しく、式も複雑な形をしている、そこで、目的関数を用いた再急降下法で近似値を推定する方法により、外心と垂心の推定座標をまずは求めて、オイラー線の近似線を作成する。
結果は、以下のような図のようになり、オイラー線を近似的にではあるが描写できていることが分かる。

![垂心推定とオイラー線.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/f7bfb1f8-2225-4d10-a95b-6e9dfbf29c39.png)

ただし、厳密に三点は同一直線上に存在するが、今回は推定座標を用いたので、少し同一直線上からずれてしまった。

# オイラー線について
三角形の外心と重心と垂心は同一直線上にあることは、外心を始点ベクトルとしたベクトルの内積計算により証明することができる。証明は、高校数学の美しい物語様の以下の記事を参照されたい。

https://manabitimes.jp/math/630

# 垂心の推定に用いる目的関数について
$\Delta ABC$の垂心を$H$とおいたとき、$BC\perp AH,AC\perp BH,AB \perp CH$が成立する。
したがって、以下のことが成立する。
```math
\overrightarrow{BC}\cdot \overrightarrow{AH}=\overrightarrow{AC}\cdot \overrightarrow{BH}=\overrightarrow{AB}\cdot \overrightarrow{CH}=0
```
そこで、目的関数を以下のように定義する。

```math
f_H=(\overrightarrow{BC}\cdot \overrightarrow{AH})^2+(\overrightarrow{AC}\cdot \overrightarrow{BH})^2+(\overrightarrow{AB}\cdot \overrightarrow{CH})^2
```

この関数は$H$が垂心となるときに、最小値0を取る。この関数を最小にする座標を再急降下法により推定する。

# 垂心の近似座標を算出するプログラム
以下のようなプログラムを作成した。

```python suisin_near.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.001
#微小分
delta=1.0*10**-5
#三角形の座標

x1=1
y1=2

x2=-3
y2=4

x3=-5
y3=-6

#垂心の推定値の初期値
xH=1
yH=2
#目的関数の定義
def fH(xH,yH):
 
  AH=np.array([xH-x1,yH-y1])
  BH=np.array([xH-x2,yH-y2])
  CH=np.array([xH-x3,yH-y3])

  BC=np.array([x3-x2,y3-y2])
  AC=np.array([x3-x1,y3-y1])
  AB=np.array([x2-x1,y2-y1])
  X=(AH@BC)**2+(BH@AC)**2+(CH@AB)**2
  
  #print(X)
  return X
#垂心の推定値を求める
i=0
f_ary=[]
i_ary=[]
while fH(xH,yH)>delta:
    #微分の計算
    dfdx=(fH(xH+delta,yH)-fH(xH,yH))/delta
    dfdy=(fH(xH,yH+delta)-fH(xH,yH))/delta
    #再急降下法によるパラメータの更新
    xH=xH-alpha*dfdx
    yH=yH-alpha*dfdy
    #反復回数と目的関数の値を格納
    f_ary.append(fH(xH,yH))
    i_ary.append(i)
    i=i+1
##結果の出力
plt.xlabel('反復回数')
plt.ylabel('目的関数')
plt.plot(i_ary,f_ary)
plt.savefig('推定垂心における目的関数の推移.png')

plt.show()
    
# ##3点の座標をプロット
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
plt.plot([x1,x2,x3,x1],[y1,y2,y3,y1],color='red')
##推定した円をプロット
plt.plot(xH,yH,color='blue',marker='o')
##アスペクト比を等しく
plt.axis('equal')
plt.savefig('垂心推定.png')
plt.show()
```

 以下のような結果となった。

まず、目的関数の推移を以下に示す。


 ![推定垂心における目的関数の推移.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/88c54d76-37e3-4a8b-9285-2c5caa7e9127.png)

次に、垂心の推定値のプロット結果を示す。

![垂心推定.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/29acbfa4-3439-46ae-be67-47f22b04d12f.png)

この結果が確からしいか推定するために、以下のようなオイラー線を推定するプログラムを作成した。
ただし、外心の推定方法は、以下の記事を参照されたい。

https://qiita.com/arairuca/items/fbecad249f8d4dced499

``` eiler_delta.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.001
#微小分
delta=1.0*10**-5
#三角形の座標

x1=1
y1=2

x2=-3
y2=4

x3=-5
y3=-6

#垂心の推定値の初期値
xH=1
yH=2
#目的関数の定義
def fH(xH,yH):
 
  AH=np.array([xH-x1,yH-y1])
  BH=np.array([xH-x2,yH-y2])
  CH=np.array([xH-x3,yH-y3])

  BC=np.array([x3-x2,y3-y2])
  AC=np.array([x3-x1,y3-y1])
  AB=np.array([x2-x1,y2-y1])
  X=(AH@BC)**2+(BH@AC)**2+(CH@AB)**2
  
  #print(X)
  return X
#垂心の推定値を求める
i=0
f_ary=[]
i_ary=[]
while fH(xH,yH)>delta:
    #微分の計算
    dfdx=(fH(xH+delta,yH)-fH(xH,yH))/delta
    dfdy=(fH(xH,yH+delta)-fH(xH,yH))/delta
    #再急降下法によるパラメータの更新
    xH=xH-alpha*dfdx
    yH=yH-alpha*dfdy
    #反復回数と目的関数の値を格納
    f_ary.append(fH(xH,yH))
    i_ary.append(i)
    i=i+1
##結果の出力
# plt.xlabel('反復回数')
# plt.ylabel('目的関数')
# plt.plot(i_ary,f_ary)
# plt.savefig('推定垂心における目的関数の推移.png')

# plt.show()
    
# ##3点の座標をプロット
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
plt.plot([x1,x2,x3,x1],[y1,y2,y3,y1],color='red')
##推定した円をプロット
plt.plot(xH,yH,color='blue',marker='o')
# ##アスペクト比を等しく



#外心の推定値の初期値
x0=1
y0=2
#目的関数の定義
def f(x0,y0):
    #外心の推定値を用いて推定半径を計算
    r1=((x0-x1)**2+(y0-y1)**2)**0.5
    r2=((x0-x2)**2+(y0-y2)**2)**0.5
    r3=((x0-x3)**2+(y0-y3)**2)**0.5
    #目的関数の値を計算
    X=(r1-r2)**2+(r1-r3)**2+(r2-r3)**2
    return X

i=0
f_ary=[]
i_ary=[]
while f(x0,y0)>delta:
    #微分の計算
    dfdx=(f(x0+delta,y0)-f(x0,y0))/delta
    dfdy=(f(x0,y0+delta)-f(x0,y0))/delta
    #再急降下法によるパラメータの更新
    x0=x0-alpha*dfdx
    y0=y0-alpha*dfdy
    #反復回数と目的関数の値を格納
    f_ary.append(f(x0,y0))
    i_ary.append(i)
    i=i+1
# #結果の出力
# plt.xlabel('反復回数')
# plt.ylabel('目的関数')
# plt.plot(i_ary,f_ary)
# plt.savefig('推定外心における目的関数の推移.png')
# plt.show()

#外心の推定座標をプロット
plt.plot(x0,y0,color='blue',marker='o')

# 重心の座標の作成
xG=1/3*(x1+x2+x3)
yG=1/3*(y1+y2+y3)
# 重心をプロット
plt.plot(xG,yG,color='blue',marker='o')

plt.plot([x0,xG,xH],[y0,yG,yH],color='blue')
plt.axis('equal')
plt.savefig('垂心推定とオイラー線.png')
plt.show()
```

このプログラムを実行すると以下のようなグラフが出力される。

![垂心推定とオイラー線.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/78be8b88-1537-4319-8ecf-beeeeafe190a.png)

青点がそらぞれ、外心、重心、垂心である。また青線がオイラー線の推定直線である。

このように、厳密には少しずれてしまい同一直線上にはなさそうだが、おおよそ、ほぼ直線状にはありそうである。これは、今回のプログラムは目的関数の最適化を近似的にしているだけであり、厳密解を求めてはいないことに起因する。

# まとめ

今回は、幾何学の一大テーマの一つである、オイラー線の描写に挑戦してみた。厳密的な証明や議論は、ここでは語れないほど膨大であるが、今回は導入部だけを紹介した。また、このような幾何学には計算機による数値解析のアプローチもある。今回は、大学受験生や数学オリンピック等を目指す高校生、幾何学に興味のある学生からプログラムによる数値解析に興味のある大学生や社会人向けに書いた記事である。なので、数学的な厳密さというよりも楽しく理解できるという記事を目標として執筆した限りである。