はじめに
同一直線上に存在しない3点を通る円周を描くことは可能である。つまり三角形には必ず外接円が存在する。そこで、$xy$座標上に三角形ABCをおいた場合、『3点から等距離の点を外心』とする外接円を描くことができる。あるいは『3辺の垂直2等分線は外心で交わる』ことを利用しても外接円を描くことができる。
そこで、今回は上記の2種類の方法を用いて外接円を描写することができるか挑戦する。
具体的には、前半では推定する外心の座標を変数として、その点が『3点から等距離』になる場合において値が最小になるような目的関数を与え、それを最適化することで外心の座標の近似値を推定する。
後半では、『3辺の垂直2等分線は外心で交わる』ことから、三角形の代表の2辺の垂直2等分線を考え、交点を代数的に算出する。
最終的に、前半と後半の内容で以下のような、三角形の外接円を描写することができることを確かめる。
再急降下法による近似解
導入
$xy$座標上の点$A(x_1,y_1),B(x_2,y_2),C(x_3,y_3)$を結んだ$\Delta ABC$を考える。ここで、推定される外心位置を$D(x_0,y_0)$とおく。このとき、以下の目的関数を考える。
f(x_0,y_0)=(AD-BD)^2+(BD-CD)^2+(AD-CD)^2
この目的関数$f(x_0,y_0)$は$AD=BD=CD(=R)$となるときに最小値0となる。
この、目的関数の最適化（最小値探索）を再急降下法というアルゴリズムを用いて行う。
アルゴリズムの概要は以下の記事を参照されたい。
最急降下法
ここでも簡潔に最急降下法アルゴリズムを説明しておこう。
$z=f(x,y)$の最小値を求めることを考える。
解空間の最小値を取る極値が一つ解空間の谷が一つの場合に有効な最小値の推定手法である。$(a,b)$に初期値を与え、以下のような式を用いて$(a,b)$を更新していく。つまり、$(a_k,b_k)$を算出する。
a_{k+1}=a_{k}-\alpha(\frac{\delta }{\delta a_{i}}f(a_{k},b_{k}))
b_{k+1}=b_{k}-\beta(\frac{\delta }{\delta b_{k}}f(a_k,b_k))
ただし$\alpha,\beta$は$(a,b)$の学習率である。
これは、ボールが坂に下り落ちる様に似ている。偏微分係数は各方向の坂の傾きと捉えることができる。
ただし、今回のプログラムでは簡単のため、$x_0,y_0$を推定するとき、共通の学習率$\alpha$を用いた。
プログラム
上記の考察より以下の外心座標推定プログラムを作成した。具体的なアルゴリズムは、適当な初期座標を与え目的関数を求め、最適化していくというものである。
目的関数の推移
まずは、試行回数における目的関数の値の推移を以下のプログラムでグラフ化した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.1
#微小分
delta=1.0*10**-5
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#外心の推定値の初期値
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()
##半径の推定
# r0=((x0-x1)**2+(y0-y1)**2)**0.5
# #円の描画(極座標)
# theta=np.linspace(0,math.pi*2,100)
##パラメータ表示
# x=x0+r0*np.cos(theta)
# y=y0+r0*np.sin(theta)
##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(x,y,color='blue')
##アスペクト比を等しく
# plt.axis('equal')
# plt.savefig('3点を通る外接円の推定.png')
# plt.show()
これを実行すると以下のようなグラフが出力される。
このように、数十回程度で値はほぼ最適化されて0に漸近する。
推定外接円
上記によって最適化された外心の推定値を用いて推定半径を求め、外接円を描写する。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#再急降下法の学習率
alpha=0.1
#微小分
delta=1.0*10**-5
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#外心の推定値の初期値
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()
#半径の推定
r0=((x0-x1)**2+(y0-y1)**2)**0.5
#円の描画(極座標)
theta=np.linspace(0,math.pi*2,100)
#パラメータ表示
x=x0+r0*np.cos(theta)
y=y0+r0*np.sin(theta)
#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(x,y,color='blue')
#アスペクト比を等しく
plt.axis('equal')
plt.savefig('3点を通る外接円の推定.png')
plt.show()
これを実行すると、以下のような推定された外接円を描写することができる。
このように、かなりいい精度で外接円を描写することができる。
厳密解
上記の方法の場合、あくまで推定値を用いているので精度をよくするためには十分な試行回数を要する。また、あくまで、近似値を求めるだけなので、後半では厳密解を求め、それを用いるものとする。
導入
$\Delta ABC$について$AB$の中点を$E$とし$BC$の中点を$F$とする。
DE \perp AB,DF \perp BC
が成立する。ゆえに以下のことが成立する。
\overrightarrow{DE}\cdot \overrightarrow{AB}=\overrightarrow{DF}\cdot \overrightarrow{BC}=0
したがって、
\begin{pmatrix}
x_0-\frac{x_1+x_2}{2} \\ y_0-\frac{y_1+y_2}{2}
\end{pmatrix}
\cdot
\begin{pmatrix}
x_3-x_1 \\ y_3-y_1
\end{pmatrix}
=0
同様のことを、$AB,AC$でも行う。
\begin{equation}
\left\{ \,
\begin{aligned}
& (x_3-x_1)x_0+(y_3-y_1)y_0 = \frac{1}{2}(x_3^2-x_1^2+y_3^2-y_1^2) \\
& (x_2-x_1)x_0+(y_2-y_1)y_0 = \frac{1}{2}(x_2^2-x_1^2+y_2^2-y_1^2) \\
\end{aligned}
\right.
\end{equation}
ここで、文字の置き換えを用いて上記の連立方程式を以下のようにあらわす。
(紛らわしいが、ここではA,B,Cは点ではないので注意する。このような変数名を用いたのは、後に述べるプログラミングの変数として使いやすいからである。)
\begin{equation}
\left\{ \,
\begin{aligned}
& Ax_0+By_0 = C \\
& Dx_0+Ey_0 = F \\
\end{aligned}
\right.
\end{equation}
この連立方程式をクラメルの公式を用いて解く。
\Delta=
\begin{vmatrix}
A & B \\
D & E
\end{vmatrix}
=AE-BD
\Delta_1=
\begin{vmatrix}
C & B \\
F & E
\end{vmatrix}
=CE-BF
\Delta_2=
\begin{vmatrix}
A & C \\
D & F
\end{vmatrix}
=AF-CD
となることから、
x_0=\frac{\Delta_1}{\Delta}
y_0=\frac{\Delta_2}{\Delta}
となる。
プログラム
これをプログラムに反映させると以下のようになる。これは、初期値を用いるといった手法ではなく、代数解析的に一発で値を厳密的に求めることができる。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#三角形の座標
x1=1
y1=2
x2=-3
y2=4
x3=-5
y3=-6
#クラメルの公式の準備（行列式の計算）
A=x3-x1
B=y3-y1
D=x2-x1
E=y2-y1
C=0.5*(x3**2-x1**2+y3**2-y1**2)
F=0.5*(x2**2-x1**2+y2**2-y1**2)
#行列式の計算
delta=(A*E-B*D)
delta1=(C*E-B*F)
delta2=(A*F-C*D)
# クラメルの公式による解の計算
a=delta1/delta
b=delta2/delta
r=((x1-a)**2+(y1-b)**2)**0.5
plt.plot(x1,y1,color='red',marker='o')
plt.plot(x2,y2,color='red',marker='o')
plt.plot(x3,y3,color='red',marker='o')
#極座標方式、つまり媒介変数表示で円を描写する。
theta=np.linspace(0,math.pi*2,100)
x=a+r*np.cos(theta)
y=b+r*np.sin(theta)
plt.plot(x,y,color='blue')
plt.axis('equal')
plt.savefig('厳密解による外接円の描写.png')
plt.show()
このように、正確に外接円を描写することができた。ただし、このプログラムは、2点のx座標が等しい場合は、分母が0となるので上手く計算することができない。
まとめ
今回は、3点を通る円の描写を近似解を推定する方法と厳密解による方法の2種類の方法をそれぞれ使用することで求めることを行った。近似解は、最適値（最小値）推定のアルゴリズムの一種である再急降下法を用いて求めた。一方で、厳密解はクラメルの公式を用いることで求めたものを使用した。このように、外接円を計算で用いるのは手計算では難しいのでプログラミングを用いた。ちなみに、外接円の方程式を求める問題は、大学受験では頻出である。したがって、受験生は今回の記事等を参考にして解法を身に着けてもらいたい。
参考文献
3元連立方程式を用いた外心の求め方