【問題】三点の座標$p1, p2, p3$が与えられた時その3点を頂点とする3角形の外心$p$の座標を求めよ
三点が決まれば外心（3点から等距離にある点）は求められると知っていてもその計算はかなり大変です。これを楽をしてsympyにやってもらおうということです。
三点の座標から作った方程式をsympyに渡して解く
2点p,qの距離をd(p,q)と表すと \\
外心pはp1, p2, p3との距離がすべて等しいので \\
d(p,p1) = d(p,p2) = d(p,p3) \\
これを素直にsympyに渡します、関数のd2はルートを取らなくて良いように$d^2$を返します。
from sympy import *
def d2(p1, p2):
return (p1[0]-p2[0])**2 + (p1[1]-p2[1])**2
def circumcente(p1,p2,p3): # get circumcente
x, y = symbols('x y')
p = (x,y)
eq = [ Eq(d2(p,p1), d2(p,p2)),Eq(d2(p,p1), d2(p,p3))]
return tuple(solve(eq, (x, y)).values())
print(circumcente((3,5),(2,-2),(-6,2)))
# (-1, 2)
solveでは{x:1, y:2} のようにDICT形式で返るので、values()を使って答えだけをtupleで返すようにしました。
このやり方は式を与えるだけなので、方程式が解ける問題ならどんなものでもつかえます。ただし大量の問題を解かせるには効率的ではないので、一般解を出す方法を考えます。
三点の座標をsymbolで渡して一般解を求める
先ほど作った関数circumcenteに対して数値の代わりにsymbolを渡すと一般解が得られます（かなり長いですが）
x, y, x1, y1, x2, y2, x3, y3 = symbols('x y x1 y1 x2 y2 x3 y3')
p = circumcente((x1,y1), (x2,y2), (x3,y3))
print("x =", p[0])
print("y =", p[1])
# x = (x1**2*y2 - x1**2*y3 - x2**2*y1 + x2**2*y3 + x3**2*y1 - x3**2*y2 + y1**2*y2 - y1**2*y3 - y1*y2**2 + y1*y3**2 + y2**2*y3 - y2*y3**2)/(2*x1*y2 - 2*x1*y3 - 2*x2*y1 + 2*x2*y3 + 2*x3*y1 - 2*x3*y2)
# y = (-x1**2*x2 + x1**2*x3 + x1*x2**2 - x1*x3**2 + x1*y2**2 - x1*y3**2 - x2**2*x3 + x2*x3**2 - x2*y1**2 + x2*y3**2 + x3*y1**2 - x3*y2**2)/(2*x1*y2 - 2*x1*y3 - 2*x2*y1 + 2*x2*y3 + 2*x3*y1 - 2*x3*y2)
この一般解を使って関数を作る
この一般解の式を使って計算すれば高速な関数を作ることが出来ます。ただコードを見てもなぜそうなるのかは、わかりずらいですけど、、、
def circumcente1(p1,p2,p3): # get circumcente
(x1, y1), (x2, y2), (x3, y3) = p1, p2, p3
x = (x1**2*y2 - x1**2*y3 - x2**2*y1 + x2**2*y3 + x3**2*y1 - x3**2*y2 + y1**2*y2 - y1**2*y3 - y1*y2**2 + y1*y3**2 + y2**2*y3 - y2*y3**2)/(2*x1*y2 - 2*x1*y3 - 2*x2*y1 + 2*x2*y3 + 2*x3*y1 - 2*x3*y2)
y = (-x1**2*x2 + x1**2*x3 + x1*x2**2 - x1*x3**2 + x1*y2**2 - x1*y3**2 - x2**2*x3 + x2*x3**2 - x2*y1**2 + x2*y3**2 + x3*y1**2 - x3*y2**2)/(2*x1*y2 - 2*x1*y3 - 2*x2*y1 + 2*x2*y3 + 2*x3*y1 - 2*x3*y2)
return (x,y)
（開発環境：Google Colab）
この考え方はProject Euler Problem 727: Triangle of Circular Arcsを解くのに役に立ちます。