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

# はじめに
平面上の2曲線の交点を求めたいとき、[終結式(resultant)](https://ja.wikipedia.org/wiki/%E7%B5%82%E7%B5%90%E5%BC%8F)を使えば求めることはできます。また、パラメータで表された平面曲線も終結式でパラメータ消去ができます。
でも、その計算がちょっと大変なときがあります。今回はそれをSymPyで解決してみました。

# 例1 交点を求める1[^1]
```math
\displaylines{
C:\; p(x,y) = x^3 + y^3 - xy = 0 \\
D:\; q(x,y) = 5(x^2+y^2) - 6xy -(x+y) = 0
}
```
![desmos-graph.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3766036/bbefbab1-121f-46ad-10b4-a338677878f0.png)
Cが赤色のグラフでDが青色のグラフです([グラフのリンク](https://www.desmos.com/calculator/slwir6s0c3))。

これらの交点のx座標を終結式によって求めてみようと思います。
それぞれをyについて整理すると
```math
\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)は次のようになります。
```math
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に計算させてみよう」というのがこの記事の趣旨になります。

```python: 例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
```
となります。つまり
```math
\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座標は
```math
x = 0, \frac{1}{2}
```
となることが分かりました。これらの値をpとqに代入すると、結局pとqの両方をみたす(x,y)は
```math
(x,y) = (0,0), \, (\frac{1}{2},\frac{1}{2})
```
となります。

終結式の計算をSymPyにさせることで、二曲線の交点を求める際に有用であることが分かりました。


# 例2 交点を求める2[^2]
```math
\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
}
```
![sample_curve2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3766036/4471ff00-1dff-602d-965b-3a75124a9656.png)
Cが赤色のグラフでG0が青色のグラフです([グラフのリンク](https://www.desmos.com/calculator/bnjiyjbyjc))。

これも交点のx座標を終結式によって求めてるために、それぞれをyについて整理します。
```math
\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)は次のようになります。
```math
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で行列式を計算してみます。


```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
```
となります。これは
```math
\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座標は
```math
x = -1, 0, 1
```
となるので、それぞれのxの値でpとqを両方みたす(x,y)を計算すると
```math
(x,y) = (-1,0), \, (0, 1), \, (1, 0)
```
となることが分かりました。

# 例3 パラメータ消去1
次はパラメータ表記されている曲線のパラメータを消去してみます。
よく知られた円の例からです。
```math
\displaylines{
x = \frac{1 - t^2}{1 + t^2} \\
y = \frac{2t}{1 + t^2}
} \hspace{30pt} \left( \; t \in \mathbb{R} \; \right)
```
![circle_param.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3766036/3ad4a144-621a-f7b9-851c-81c1ab8168f5.png)
[グラフのリンク](https://www.desmos.com/calculator/l8anby6lcu)

これをtについて整理して改めてp,qとすると
```math
\displaylines{
  p(x,t) = -(x+1)t^2 + (1 - x) = 0 \\
 q(y,t) = -yt^2 + 2t - y = 0
}
```
とできます。
tについての係数に着目すると終結式R(p,q)は次のようになります。
```math
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で計算していきます。
```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
```
ですので、これは
```math
\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]です。
```math
\displaylines{
x = \frac{3ct}{1+t^3} \\
y = \frac{3ct^2}{1+t^3}
} \hspace{30pt} \left( \; t \neq -1 \; \right)
```
![descartes_leaf.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3766036/4c7e71d3-d94f-f47b-b224-4a69ccfbb88a.png)
[グラフ](https://www.desmos.com/calculator/g78twpaom8)


先程と同じようにtについて整理して改めてp,qとすると
```math
\displaylines{
  p(x,t) = -xt^3 + 3ct - x = 0 \\
 q(y,t) = -yt^3 + 3ct^2 - y = 0
}
```
とできます。


終結式R(p,q)は次のようになります。
```math
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で計算すると
```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
```
です。これは
```math
\begin{align}
Res(p,q) &= 81c^4xy - 27c^3x^3 -27c^3y^3  &(実行結果の2行目) \\
  &= 27c^3(cxy - x^3 - y^3)
\end{align}
```
を意味しますので、パラメータを消去できたことが分かりました。
[比較のためのグラフ](https://www.desmos.com/calculator/2p3qqbyjcx)

# おわりに
手計算だと時間を要する行列式の計算を、要素が式の場合で実行させてみました。改めてSymPyの強力さを感じます。



[^1]: 数学のかんどころ12 平面代数曲線 p127
[^2]: 数学のかんどころ12 平面代数曲線 p138の射影曲線でz=1としてxy平面上の式にしたものです
[^3]: [Wikipedia デカルトの正葉線](https://ja.wikipedia.org/wiki/%E3%83%87%E3%82%AB%E3%83%AB%E3%83%88%E3%81%AE%E6%AD%A3%E8%91%89%E7%B7%9A)