SymPyにより終結式で変数消去をするサンプル

はじめに

平面上の2曲線の交点を求めたいとき、終結式(resultant)を使えば求めることはできます。また、パラメータで表された平面曲線も終結式でパラメータ消去ができます。
でも、その計算がちょっと大変なときがあります。今回はそれをSymPyで解決してみました。

例1 交点を求める1¹

\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に計算させてみよう」というのがこの記事の趣旨になります。

例1のソース

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)

これを実行すると

例1の実行結果

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 交点を求める2²

\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で行列式を計算してみます。

例2のソース

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)

実行結果は

例2の実行結果

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で計算していきます。

例3のソース

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)

実行結果は

例3の実行結果

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がある場合を扱ってみたいと思います。
デカルトの葉線³です。

\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で計算すると

例4のソース

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)

実行結果は

例4の実行結果

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の強力さを感じます。

1. 数学のかんどころ12 平面代数曲線 p127 ↩

2. 数学のかんどころ12 平面代数曲線 p138の射影曲線でz=1としてxy平面上の式にしたものです ↩

3. Wikipedia デカルトの正葉線 ↩