はじめに
平面上の2曲線の交点を求めたいとき、終結式(resultant)を使えば求めることはできます。また、パラメータで表された平面曲線も終結式でパラメータ消去ができます。
でも、その計算がちょっと大変なときがあります。今回はそれをSymPyで解決してみました。
例1 交点を求める11
\displaylines{
C:\; p(x,y) = x^3 + y^3 - xy = 0 \\
D:\; q(x,y) = 5(x^2+y^2) - 6xy -(x+y) = 0
}
Cが赤色のグラフでDが青色のグラフです(グラフのリンク)。
これらの交点のx座標を終結式によって求めてみようと思います。
それぞれをyについて整理すると
\displaylines{
C:\; p(x,y) = \color{red}{1}y^3 + \color{red}{0}y^2 \color{red}{-x}y + \color{red}{x^3} = 0 \\
D:\; q(x,y) = \color{blue}{5}y^2 \color{blue}{-(6x+1)}y + \color{blue}{5x^2-x} = 0
}
となるため、pとqの終結式Res(p,q)は次のようになります。
Res(p,q) = \begin{vmatrix}
\color{red}{1} & \color{red}{0} & \color{red}{-x} & \color{red}{x^3} & 0 \\
0 & \color{red}{1} & \color{red}{0} & \color{red}{-x} & \color{red}{x^3} \\
\color{blue}{5} & \color{blue}{-6x-1} & \color{blue}{5x^2-x} & 0 & 0 \\
0 & \color{blue}{5} & \color{blue}{-6x-1} & \color{blue}{5x^2-x} & 0 \\
0 & 0 & \color{blue}{5} & \color{blue}{-6x-1} & \color{blue}{5x^2-x} \\
\end{vmatrix}
これはpとqのyに関する係数をずらしながら作成した行列式です。Res(p,q) = 0が、pとqの連立方程式でyを消去したことに相当します。
「これをSymPyに計算させてみよう」というのがこの記事の趣旨になります。
import sympy as sy
x = sy.Symbol("x")
# シルベスター行列aを作成して表示
array_a = [
1, 0, -x, x**3, 0,
0, 1, 0, -x, x**3,
5, -6*x-1, 5*x**2-x, 0, 0,
0, 5, -6*x-1, 5*x**2-x, 0,
0, 0, 5, -6*x-1, 5*x**2-x,
]
a = sy.Matrix(5, 5, array_a)
print(a)
# 終結式を求めて表示
det_a = sy.det(a)
print(det_a)
# 因数分解して交点のx座標を求める
expr = sy.factor(det_a)
print(expr)
これを実行すると
Matrix([[1, 0, -x, x**3, 0], [0, 1, 0, -x, x**3], [5, -6*x - 1, 5*x**2 - x, 0, 0], [0, 5, -6*x - 1, 5*x**2 - x, 0], [0, 0, 5, -6*x - 1, 5*x**2 - x]])
16*x**6 - 32*x**5 + 24*x**4 - 8*x**3 + x**2
x**2*(2*x - 1)**4
となります。つまり
\begin{align}
Res(p,q) &= 16x^6 - 32x^5 + 24x^4 - 8x^3 + x^2 &(実行結果の2行目)\\
&= x^2(2x-1)^4 &(実行結果の3行目)
\end{align}
という計算がなされたことになります。
結局２つの式p,qの交点のx座標は
x = 0, \frac{1}{2}
となることが分かりました。これらの値をpとqに代入すると、結局pとqの両方をみたす(x,y)は
(x,y) = (0,0), \, (\frac{1}{2},\frac{1}{2})
となります。
終結式の計算をSymPyにさせることで、二曲線の交点を求める際に有用であることが分かりました。
例2 交点を求める22
\displaylines{
C:\; p(x,y) = x^2 + y^2 - 1 = 0 \\
G0:\; q(x,y) = 5(y - 1)(x^2 - 1) + 4(x^2 + y^2 - 1) = 0
}
Cが赤色のグラフでG0が青色のグラフです(グラフのリンク)。
これも交点のx座標を終結式によって求めてるために、それぞれをyについて整理します。
\displaylines{
C:\; p(x,y) = y^2 + 0y + x^2 - 1 = 0 \\
G0:\; q(x,y) = 4y^2 + 5(x^2 - 1)y - x^2 + 1 = 0
}
これよりpとqの終結式Res(p,q)は次のようになります。
Res(p,q) = \begin{vmatrix}
1 & 0 & x^2-1 & 0 \\
0 & 1 & 0 & x^2-1 \\
4 & 5(x^2 - 1) & -x^2 + 1 & 0 \\
0 & 4 & 5(x^2 - 1) & -x^2 + 1 \\
\end{vmatrix}
同様にPythonで行列式を計算してみます。
import sympy as sy
x = sy.Symbol("x")
# シルベスター行列aを作成して表示
array_a = [
1, 0, x**2-1, 0,
0, 1, 0, x**2-1,
4, 5*(x**2-1), -x**2+1, 0,
0, 4, 5*(x**2-1), -x**2+1,
]
a = sy.Matrix(4, 4, array_a)
print(a)
# 終結式を求めて表示
det_a = sy.det(a)
print(det_a)
# 因数分解して交点のx座標を求める
expr = sy.factor(det_a)
print(expr)
実行結果は
Matrix([[1, 0, x**2 - 1, 0], [0, 1, 0, x**2 - 1], [4, 5*x**2 - 5, 1 - x**2, 0], [0, 4, 5*x**2 - 5, 1 - x**2]])
25*x**6 - 50*x**4 + 25*x**2
25*x**2*(x - 1)**2*(x + 1)**2
となります。これは
\begin{align}
Res(p,q) &= 25x^6 - 50x^4 + 25x^2 &(実行結果の2行目)\\
&= 25x^2(x-1)^2(x+1)^2 &(実行結果の3行目)
\end{align}
という計算がなされたことになります。
これより２つの式p,qの交点のx座標は
x = -1, 0, 1
となるので、それぞれのxの値でpとqを両方みたす(x,y)を計算すると
(x,y) = (-1,0), \, (0, 1), \, (1, 0)
となることが分かりました。
例3 パラメータ消去1
次はパラメータ表記されている曲線のパラメータを消去してみます。
よく知られた円の例からです。
\displaylines{
x = \frac{1 - t^2}{1 + t^2} \\
y = \frac{2t}{1 + t^2}
} \hspace{30pt} \left( \; t \in \mathbb{R} \; \right)
これをtについて整理して改めてp,qとすると
\displaylines{
p(x,t) = -(x+1)t^2 + (1 - x) = 0 \\
q(y,t) = -yt^2 + 2t - y = 0
}
とできます。
tについての係数に着目すると終結式R(p,q)は次のようになります。
Res(p,q) = \begin{vmatrix}
-(x+1) & 0 & 1 - x & 0 \\
0 & -(x+1) & 0 & 1 - x \\
-y & 2 & -y & 0 \\
0 & -y & 2 & -y \\
\end{vmatrix}
先程と同じように、これをPythonで計算していきます。
x = sy.Symbol("x")
y = sy.Symbol("y")
# シルベスター行列aを作成して表示
array_a = [
-x-1, 0, 1-x, 0,
0, -x-1, 0, 1-x,
-y, 2, -y, 0,
0, -y, 2, -y,
]
a = sy.Matrix(4, 4, array_a)
print(a)
# 終結式を求めて表示
det_a = sy.det(a)
print(det_a)
実行結果は
Matrix([[-x - 1, 0, 1 - x, 0], [0, -x - 1, 0, 1 - x], [-y, 2, -y, 0], [0, -y, 2, -y]])
4*x**2 + 4*y**2 - 4
ですので、これは
\begin{align}
Res(p,q) &= 4x^2 + 4y^2 - 4 &(実行結果の2行目) \\
&= 4(x^2 + y^2 -1)
\end{align}
となり、Res(p,q)=0が円の方程式になっていることも確認できます。
終結式の計算をSymPyにさせることで、パラメータ消去にも使えることが分かりました。
例4 パラメータ消去2
最後にパラメータの消去なのですが、定数cがある場合を扱ってみたいと思います。
デカルトの葉線3です。
\displaylines{
x = \frac{3ct}{1+t^3} \\
y = \frac{3ct^2}{1+t^3}
} \hspace{30pt} \left( \; t \neq -1 \; \right)
先程と同じようにtについて整理して改めてp,qとすると
\displaylines{
p(x,t) = -xt^3 + 3ct - x = 0 \\
q(y,t) = -yt^3 + 3ct^2 - y = 0
}
とできます。
終結式R(p,q)は次のようになります。
Res(p,q) = \begin{vmatrix}
-x & 0 & 3c & -x & 0 & 0 \\
0 & -x & 0 & 3c & -x & 0 \\
0 & 0 & -x & 0 & 3c & -x \\
-y & 3c & 0 & -y & 0 & 0 \\
0 & -y & 3c & 0 & -y & 0 \\
0 & 0 & -y & 3c & 0 & -y\\
\end{vmatrix}
これをPythonで計算すると
import sympy as sy
x = sy.Symbol("x")
y = sy.Symbol("y")
c = sy.Symbol("c")
# シルベスター行列aを作成して表示
array_a = [
-x, 0, 3*c, -x, 0, 0,
0, -x, 0, 3*c, -x, 0,
0, 0, -x, 0, 3*c, -x,
-y, 3*c, 0, -y, 0, 0,
0, -y, 3*c, 0, -y, 0,
0, 0, -y, 3*c, 0, -y,
]
a = sy.Matrix(6, 6, array_a)
print(a)
# 終結式を求めて表示
det_a = sy.det(a)
print(det_a)
実行結果は
Matrix([[-x, 0, 3*c, -x, 0, 0], [0, -x, 0, 3*c, -x, 0], [0, 0, -x, 0, 3*c, -x], [-y, 3*c, 0, -y, 0, 0], [0, -y, 3*c, 0, -y, 0], [0, 0, -y, 3*c, 0, -y]])
81*c**4*x*y - 27*c**3*x**3 - 27*c**3*y**3
です。これは
\begin{align}
Res(p,q) &= 81c^4xy - 27c^3x^3 -27c^3y^3 &(実行結果の2行目) \\
&= 27c^3(cxy - x^3 - y^3)
\end{align}
を意味しますので、パラメータを消去できたことが分かりました。
比較のためのグラフ
おわりに
手計算だと時間を要する行列式の計算を、要素が式の場合で実行させてみました。改めてSymPyの強力さを感じます。