sympyで線型代数学の問題を解く

はじめに

線型代数学では掃き出し法や正規直交化法などなかなか疲れる手計算が多いかと思います。入試などでは手計算をミスなくこなす必要があるかと思いますが、実際に線型代数学を使いたい場面では具体的な値の計算は計算機にやってほしいことがほとんどだと思います。

そこで、この記事では記号計算をやってくれるsympyを用いて線型代数学の諸問題を解く方法について書きます。

スバラシク実力がつくと評判の線形代数キャンパス・ゼミ―大学の数学がこんなに分かる!単位なんて楽に取れる!

教養の線形代数

線型代数入門 (基礎数学1)

これらを見た意味は特にないです。持っている本だというだけです。

ベクトルと空間座標の基本について

$\boldsymbol{a} = \left[ 1\ 1\ -1 \right]^T, \boldsymbol{b} = \left[ 2\ 0\ 1\right]^T$のとき、$\boldsymbol{a}$と垂直な単位ベクトル$\boldsymbol{e}$を求めよ。

あるベクトルに垂直な単位ベクトルは次の式で求めることができます。

$$ \boldsymbol{b - \frac{a \cdot b}{| a |^2} a } $$

これをsympyで計算させます。

a, b = symbols('a b')
a = Matrix([1, 1, -1])
b = Matrix([2, 0, 1])
u = b - ( (transpose(a) * b ) / a.norm()**2 ).norm() * a
u.normalized()

これを実行すると、$ \left[ \begin{matrix} \frac{5 \sqrt{42}}{42} \\ - \frac{\sqrt{42}}{42}\\ \frac{2 \sqrt{42}}{21} \end{matrix} \right]$という計算結果が得られます。

求めるべき単位ベクトルはこのベクトルと反対向きのものもあるので、答えは次のようになります。

$$ \boldsymbol{e} = \pm \left[\begin{matrix} \frac{5 \sqrt{42}}{42}\\ - \frac{\sqrt{42}}{42}\\ \frac{2 \sqrt{42}}{21} \end{matrix} \right] $$

平面$\pi: x+3y+2z+1=0$と並行な平面$\alpha$が点$(1,-1,-2)$を通るとする。 (1) 平面$\alpha$の方程式を求めよ。
(2) 平面$\alpha$と直線$L: \frac{x-2}{-2} = \frac{y}{2} = \frac{z+1}{1}$との交点$B$を求めよ。

(1) について 平面$\pi$は点$(1,3,2)$を法線ベクトルとし、明らかに点$(-1,0,0)$を通過するので次のように平面を指定できます。（もっといい感じに指定する方法があるかもしれませんが）

x, y, z = symbols('x y z')
p = Plane(Point3D(-1,0,0), normal_vector=(1,3,2))

平面$\alpha$は平面$\pi$に並行なので、平面$\pi$と同じ法線ベクトルをもち、点$(1,-1,2)$を通るので次のように求めることができます。

alpha = Plane(Point3D(1,-1,2), normal_vector=(1,3,2))
alpha.is_parallel(p)

平面$\alpha$は平面$\pi$と並行なのでis_parallelを用いて確認するとTrueとなります。

(2)について 直線$L$は点$(2, 0, -1)$および点$(2-2, 0-2, -1-1)$を通るので次のように指定できます。

l = Line3D(Point3D(2,0,-1), Point3D(2-(-2),0-2,-1-1))

直線$L$と平面$\alpha$の交点はintersectionによって求めることができます。

alpha.intersection(l)

この結果より、$(-4, -2, -2)$で平面$\alpha$と直線$L$が交わることが分かります。

2個の1次方程式 $$ x + 2 y + 3 z = 1\\ 3x + 2 y + z = -1 $$ の表す直線のベクトル表示を求めよ。

plane1 = Plane(Point3D(1,0,0), normal_vector=(1, 2, 3))
plane2 = Plane(Point3D(0,0,-1), normal_vector=(3,2,1))
plane1.intersection(plane2)

2点$\boldsymbol{a, b}$を通る直線のベクトル表示は$(1-t)\boldsymbol{a} + t \boldsymbol{b}$で求めることができます。上のコードの結果が$(-5,9,-4)$として得られるので、次のようにベクトル表示を求めます。

t = Symbol('t')
(1-t) * Matrix([-1,1,0]) + t * Matrix([-5,9,-4])

係数部分を取り出し媒介変数としてまとめてしまえば次のベクトル表示に変形できます。 $$ \begin{bmatrix} x\\ y\\ z \end{bmatrix} = \begin{bmatrix} -1\\ 1\\ -4 \end{bmatrix} + t \begin{bmatrix} 1\\ -2\\ 1 \end{bmatrix} $$

行列と行列式および連立1次方程式の解

行列の計算

行列の積$ \begin{bmatrix} -1 & 2 & 1 \ 3 & 2 & 1 \end{bmatrix} \begin{bmatrix} 4 & -1 & 1 \\ 2 & 2 & -3\\ 1 & 1 & 0 \end{bmatrix} $を求めよ。

この問題は特にsympyでやらなくても計算できますが行列を定義すれば演算することができます。

A = Matrix(([-1,2,1],[3,2,1]))
B = Matrix(([4,-1,1],[2,2,-3],[1,1,0]))
A * B

行列式

行列式$ \begin{vmatrix} 2 - \lambda & 3 & -1 \\ 2 & 1 - \lambda & 1 \\ 1 & -1 & 4 - \lambda \end{vmatrix} $を求めよ。

このような記号を含んだ問題でも行列式を求めることができます。

A, B, lmd = symbols('A B lambda')
A = Matrix(([1,2,1], [2,-1,1], [3,1,2]))
B = Matrix(([2-lmd,3,-1], [2,1-lmd,1], [1,-1,4-lmd]))

これを実行すると、 $$ \lambda^{3} + 7 \lambda^{2} - 10 \lambda - 8 $$ を得ることができます。

逆行列と連立1次方程式

$A = \begin{bmatrix} 2 & 1 & 0 \\ 1 & -1 & 2 \\ -1 & 0 & -1 \end{bmatrix} $の逆行列を余因子行列と行列式から求めよ。

逆行列$A^{-1}$は余因子行列$\tilde{A}$と行列式$|A|$を用いて次のように書くことができる。 $$A^{-1}=\frac{\tilde{A}}{|A|}$$

なお、余因子行列はadjugateで求めることができます。したがって、次のように書けば逆行列を求めることができます。

A = Symbol('A')
A = Matrix(([2,1,0], [1,-1,2], [-1,0,-1]))
A.adjugate() / A.det()

当然、逆行列が欲しい場合は直接invを用いて計算することができます。

A.inv()

連立方程式 $ \begin{cases} x_1 - x_2 + 2 x_3 = 8 & \\ 2 x_1 + 3 x_2 + x_3 = 5 & \\ - x_1 + 4 x_2 + 4 x_3 = 1 \end{cases} $を解け。

• solveを使って解く方法(1)

x1, x2, x3 = symbols('x1 x2 x3')
solve([x1-x2+2*x3-8, 2*x1+3*x2+x3-5, -x1+4*x2+4*x3-1],[x1,x2,x3])

• solveを使って解く方法(2)

solve([Eq(x1-x2+2*x3-8), Eq(2*x1+3*x2+x3-5), Eq(-x1+4*x2+4*x3-1)],[x1,x2,x3])

• linsolveを使って解く方法

system = Matrix(([1,-1,2,8], [2,3,1,5],[-1,4,4,1]))
linsolve(system, (x, y, z))

solveを使う方法では、それぞれの方程式を右辺0となるように移項し、変数を指定します。linsolveで解く場合には拡大係数行列$\boldsymbol{[A|b]}$をつくり、変数名を指定します。

正規直交基底と固有値・固有ベクトルおよび対角化

正規直交基底

$a_1 = [1\ 1\ 0\ 0]^T, a_2 = [0\ 1\ 1\ 0]^T, a_3 = [0\ 1\ 1\ 1]^T, a_4 = [1\ 2\ 0\ 1]^T$で与えられる$\mathbb{R}^4$の基底${a_1, a_2, a_3, a_4}$を正規直交基底に変換せよ。

基底を行列としてGramSchmidtに与えれば正規直交基底を求めることができる。このときデフォルトの状態では正規化を行わないので正規直交基底を求めたい場合はTrueを指定する必要がある。

A = [Matrix([1,1,0,0]), Matrix([0,1,1,0]), Matrix([0,1,1,1]), Matrix([1,2,0,1])]
GramSchmidt(A, True) # 正規化をしない場合はFalseのままでよい

正規直交基底は次のようになります。 $$ \left [ \left[\begin{matrix}\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\\ 0\\ 0 \end{matrix}\right], \quad \left[\begin{matrix}- \frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\0\end{matrix}\right], \quad \left[\begin{matrix}0\\0\\0\\1\end{matrix}\right], \quad \left[\begin{matrix} - \frac{\sqrt{3}}{3}\\ \frac{\sqrt{3}}{3}\\ - \frac{\sqrt{3}}{3}\\ 0\end{matrix}\right] \right ]$$

固有値と固有ベクトル

行列$ \begin{bmatrix} 3 & 2 \\ -2 & 3 \end{bmatrix} $の固有値と固有ベクトルを求めよ。

固有値はeigenvals、固有ベクトルはeigenvectsでもとめることができる。

A = Matrix(([3, 2], [-2,3]))
A.eigenvals()
A.eigenvects()

これを実行すると固有値$3-i$に対して固有ベクトル$\begin{bmatrix}i\\ 1\end{bmatrix}$、固有$3+2i$に対して固有ベクトル$\begin{bmatrix}-i\\ 1\end{bmatrix}$が求まります。

対角化

行列$A = \begin{bmatrix} 3 & 4 \\ -1 & -2 \end{bmatrix} $を対角化し、$Aⁿ$を求めよ。

A = Matrix(([3, 4], [-1,-2]))
P,D = A.diagonalize() # 対角行列D = P^{-1}AP

対角化できる行列の場合はdiagonalizeを用いて対角行列$D$と正則行列$P$を同時に求めることができます。

対角行列$D=P^{-1}AP$を$n$乗すると$D^{n}=P^{-1}A^{n}P$となるので$Aⁿ=PD^{n}P^{-1}$となります。

n=Symbol('n')
P * D**n * P.inv()

これを実行すると次の行列が得られる。 $$ \left[\begin{matrix}- \frac{\left(-1\right)^{n}}{3} + \frac{4}{3} 2^{n} & - \frac{4 \left(-1\right)^{n}}{3} + \frac{4}{3} 2^{n}\\ \frac{\left(-1\right)^{n}}{3} - \frac{2^{n}}{3} & \frac{4 \left(-1\right)^{n}}{3} - \frac{2^{n}}{3}\end{matrix}\right] $$

ジョルダン標準形

行列 $ A = \begin{bmatrix} 0 & 1 & -1 \\ 2 & 1 & -2 \\ 1 & 4 & -4 \end{bmatrix} $ をジョルダン標準形に変換せよ。

ジョルダン標準形はまあよく分かってないのですが、jordan_formで求めることができるみたいです。

A = Symbol('A')
A = Matrix([[0,1,-1],[2,1,-2],[1,4,-4]])
(P, J) = A.jordan_form()

Pが変換行列で、Jがジョルダン標準形らしいです。

おわりに

この記事では、線型代数学の演習問題をsympyで解く方法について述べました。 やはり$n$や$\lambda$といった記号をそのまま計算できるのは便利だと思いました。 また、計算結果をすぐに得ることができるので「この数値を少し変えたらどうなるだろう？」といったことをすぐに知ることができるのも便利ではないかと。 あと、sympyはjupyter notebookで使うと最高だと思います。