Return to SymPy 演習

2. SymPy で方程式を解く

SymPy で方程式（等式）を定義し，（主に代数）方程式を解く。関連して，Python のデータ型であるタプル (tuple) と辞書 (dict) について簡単に説明します。

SymPy の import

In [1]:
from sympy.abc import *
from sympy import *

init_printing()

方程式

SymPy は，式の計算だけでなく，代数方程式を解くことができます。例を以下に示します。

$x^2 + 2 x – 2 = 0$ を $x$ について解きます。

方程式の定義 Eq()

方程式とは，未知変数を含む等式のことです。

方程式は等式と同じく Eq(左辺, 右辺) のように定義します。（= は SymPy では代入を表します。）

In [2]:
eq = Eq(x**2 + 2*x - 2, 0)
eq
Out[2]:
$\displaystyle x^{2} + 2 x – 2 = 0$
方程式の左辺 .lhs，右辺 .rhs

方程式 eq の左辺は eq.lhs (left-hans-side)，右辺は eq.rhs (right-hand-side) のようにして取り出すことができます。

In [3]:
eq.lhs
Out[3]:
$\displaystyle x^{2} + 2 x – 2$
In [4]:
eq.rhs
Out[4]:
$\displaystyle 0$

方程式を解く solve()

方程式を未知変数について解く書式は solve(方程式, 未知変数) です。

たとえば，方程式 eq を x について解く場合は solve(eq, x)：

In [5]:
solve(eq, x)
Out[5]:
$\displaystyle \left[ -1 + \sqrt{3}, \ – \sqrt{3} – 1\right]$

方程式が 左辺 = 0 の形（つまり右辺がゼロ）になっている場合は，solve(方程式の左辺, 未知変数) としても同じです。

In [6]:
solve(x**2 + 2*x - 2, x)
Out[6]:
$\displaystyle \left[ -1 + \sqrt{3}, \ – \sqrt{3} – 1\right]$

3 次方程式 $x^3 – 8 = 0$ を解きます。solve(eq, x) の場合の x は省略して solve(eq) と書くこともできるようです。

In [7]:
sols = solve(x**3 - 8)
sols
Out[7]:
$\displaystyle \left[ 2, \ -1 – \sqrt{3} i, \ -1 + \sqrt{3} i\right]$

解が複数個ある場合は，solve() は解のリストを返します。

念のため，3番目の答え（Python はゼロからはじまることに注意）を3乗して展開すると…

In [8]:
sols[2]**3
Out[8]:
$\displaystyle \left(-1 + \sqrt{3} i\right)^{3}$

expand() でかっこを展開して計算させると…

In [9]:
expand(_)
Out[9]:
$\displaystyle 8$

連立方程式を解く solve([], [])

連立方程式を解く際は，方程式と変数をリストにして solve() に入れます。

solve([eq1, eq2, ..., eqn], [x1, x2, ..., xn])

例として
\begin{eqnarray}
x^2 + y^2 &=& 5 \\
x + y &=& 1
\end{eqnarray}

を $x, y$ について解きます。

In [10]:
eq1 = Eq(x**2 + y**2, 5)
eq2 = Eq(x + y, 1)

display(eq1, eq2)
solve([eq1, eq2], [x, y])
$\displaystyle x^{2} + y^{2} = 5$
$\displaystyle x + y = 1$
Out[10]:
$\displaystyle \left[ \left( -1, \ 2\right), \ \left( 2, \ -1\right)\right]$

方程式の右辺を移項して

\begin{eqnarray}
x^2 + y^2 -5 &=&0 \\
x + y -1 &=&0
\end{eqnarray}

の形にすれば，各方程式の左辺のみを使って以下のように書いてもよいです。

In [11]:
sols12 = solve([x**2 + y**2 - 5, x + y - 1], [x, y])
sols12
Out[11]:
$\displaystyle \left[ \left( -1, \ 2\right), \ \left( 2, \ -1\right)\right]$

2組の解がある場合，solve() は解 $(x, y)$ の組をリストで返します。

In [12]:
type(sols12)
Out[12]:
list

display() による表示

上のコードでは確認の意味を込めて display() で方程式を表示させています。

print() による表示と display() による表示を比較してみてください。display() で表示させたほうが，綺麗に形が整えられて表示されるので，確認しやすいと思います。

In [13]:
print(eq1)
display(eq1)

Eq(x**2 + y**2, 5)
$\displaystyle x^{2} + y^{2} = 5$

display() を使わなくても，eq1 と書けば表示されますが，以下のように1つのセル内に複数行書くと，最後の行しか表示されません。

In [14]:
eq1
eq2
Out[14]:
$\displaystyle x + y = 1$

display() を使うと…

In [15]:
display(eq1)
display(eq2)
$\displaystyle x^{2} + y^{2} = 5$
$\displaystyle x + y = 1$

タプル (tuple)：Python のデータ型

解のリスト sols12 の各要素は以下のように ( ) の中に複数の（今の場合は2個の）要素がカンマ (,) で区切られて並んでいます。

リストと似ていますが，[ ] の中に要素が並ぶリストと違い，( ) の中に要素が並ぶデータ型を「タプル (tuple)」と呼びます。（( ) なしでもタプル。）

In [16]:
for sol in sols12:
    print(sol)

(-1, 2)
(2, -1)

type() で確認します。

In [17]:
type(sols12[0]), type(sols12[1])
Out[17]:
(tuple, tuple)
参考：タプル (tuple) と リスト (list) の違い

リストとは異なり，タプルは変更ができません。

リストの場合は，例えば lista の2番目の要素 lista[1] の値を別の値に置き換えることができます。

In [18]:
# リストの定義

lista = [1, 2, 3]
type(lista)
Out[18]:
list
In [19]:
# lista の2番目の要素。Python はゼロ始まりだから...

lista[1]
Out[19]:
$\displaystyle 2$
In [20]:
# lista の2番目の要素の値を 20 に変更

lista[1] = 20
lista
Out[20]:
$\displaystyle \left[ 1, \ 20, \ 3\right]$

しかし，タプルの場合は，例えば tupleb の2番目の要素 tupleb[1] の値を別の値に置き換えようとすると…

In [21]:
# タプルの定義

tupleb = 1, 2, 3
type(tupleb)
Out[21]:
tuple
In [22]:
# tupleb の2番目の要素。
# リストの場合と同じように参照できる。

tupleb[1]
Out[22]:
$\displaystyle 2$
In [23]:
# tupleb の2番目の要素の値を変更しようとすると...

tupleb[1] = 20

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-6574e716f335> in <module>
      1 # tupleb の2番目の要素の値を変更しようとすると...
      2 
----> 3 tupleb[1] = 20

TypeError: 'tuple' object does not support item assignment

エラーとなります。

複数変数への一時的代入

.subs()

念のため，solve() で得られた答えを eq1 の左辺 eq1.lhs に代入して確かに答えが $5$ になっているか，確認してみます。

.subs() による一時的代入は既に説明済みですが，複数の変数（今の場合は x と y）に値 x1, y1 を一時的に代入するには以下のように

.subs({x:x1, y:y1})

または，

.subs([(x, x1), (y, y1)])

In [24]:
print('(x, y) =', sols12[0])
(x1, y1) = sols12[0]
print('eq1 の左辺 =', eq1.lhs.subs({x:x1, y:y1}))

(x, y) = (-1, 2)
eq1 の左辺 = 5
In [25]:
print('(x, y) =', sols12[0])
(x1, y1) = sols12[0]
print('eq1 の左辺 =', eq1.lhs.subs([(x, x1), (y, y1)]))

(x, y) = (-1, 2)
eq1 の左辺 = 5

もう1組の解 (2, -1) についても確認してみて。

In [ ]:
Subs()

一時的な代入を Subs(expr, variables, values).doit() でやってみます。

Subs() を使う場合は変数を tuple にできて

Subs(expr, (x, y), (x1, y1)).doit() のように書きます。

In [26]:
print('(x, y) =', sols12[0])
print('eq1 の左辺 =', Subs(eq1.lhs, (x, y), sols12[0]).doit())

(x, y) = (-1, 2)
eq1 の左辺 = 5

もう1組の解 (2, -1) についても確認してみて。

In [ ]:
例題：鶴亀算
例題
ツルとカメが合わせて8匹、足の数が合わせて26本であるとき、ツルとカメは何匹（何羽）いるか。ただしツルの足は2本、カメの足は4本である。
… とある。ツルを何匹とは数えないだろうから，
ツルとカメの個体数の合計は 8。足の数の合計が 26本。
 ツルは何羽，亀は何匹？
ツルの個体数を tsuru，カメの個体数を kame として連立方程式にして解く。
なお，最初に
from sympy.abc import *

としていると，1文字変数については宣言する必要がないが，2文字以上の変数については，以下のように宣言してやる必要がある。
In [27]:
var('tsuru kame')

Out[27]:
$\displaystyle \left( tsuru, \ kame\right)$
In [28]:
# 個体数の合計

eq1 = Eq(tsuru + kame, 8)

In [29]:
# 足の本数の合計

eq2 = Eq(2*tsuru + 4*kame, 26)

In [30]:
# 確認のため，display() で方程式を表示

display(eq1, eq2)

$\displaystyle kame + tsuru = 8$
$\displaystyle 4 kame + 2 tsuru = 26$
solve() で連立方程式を解きます。
In [31]:
sols = solve([eq1, eq2], [tsuru, kame])
print(sols)

{kame: 5, tsuru: 3}

辞書 (dict)：Python のデータ型
同じ solve() で連立方程式を解いても，この場合の連立方程式の解 sols には，タプルではなく，{ } の中に 変数名 (key): 数値 (value) のような形でデータが格納されています。このようなデータ型を Python では「辞書 (dict)」と呼びます。
type() で確認します。
In [32]:
type(sols)

Out[32]:
dict
辞書 sols の中の，変数名 (key) kame や tsuru とペアになっている値 (value) を表示するには，以下のように。
In [33]:
print('カメは', sols[kame], '匹')
print('ツルは', sols[tsuru], '羽')

カメは 5 匹
ツルは 3 羽

念のため，得られた解を方程式 eq1 の左辺に代入して，確かに $8$ になっていることを確かめます。
In [34]:
eq1

Out[34]:
$\displaystyle kame + tsuru = 8$
今の場合，解 sols は辞書型で与えられていますから…
In [35]:
type(sols)

Out[35]:
dict
方程式 eq1 の左辺 eq1.lhs に一時的に代入する .subs() の引数には sols をそのまま入れます。
In [36]:
eq1.lhs.subs(sols)

Out[36]:
$\displaystyle 8$
参考：出力される解のデータ型の違い
SymPy の solve() で連立方程式を解くと，得られた解のデータ型は
  • 解が1組のみの場合は，辞書 (dict)
  • 解が複数組ある場合は，タプル (tuple) のリスト (list)
となるようです。辞書かタプルかによって，得られた解の参照の仕方が異なりますので少しややこしい。
1 変数の方程式の解が複数ある例。解は数値・変数のリスト (list) で出力される。
In [37]:
solve(x**2 + 2*x - 2, x)

Out[37]:
$\displaystyle \left[ -1 + \sqrt{3}, \ – \sqrt{3} – 1\right]$
In [38]:
solve(x**2 + 2*b*x + c, x)

Out[38]:
$\displaystyle \left[ – b – \sqrt{b^{2} – c}, \ – b + \sqrt{b^{2} – c}\right]$
2 変数の連立方程式の解が 1 組だけの例。解は辞書 (dict) 型で出力される。
In [39]:
var('tsuru kame')
solve([Eq(tsuru + kame, 8), 
       Eq(2*tsuru + 4*kame, 26)], [tsuru, kame])

Out[39]:
$\displaystyle \left\{ kame : 5, \ tsuru : 3\right\}$
2 変数の連立方程式の解が 2 組ある場合。解はタプル (tuple) を要素としたリストで出力される。
In [40]:
solve([Eq(x**2 + y**2, 5), Eq(x + y, 1)], [x, y])

Out[40]:
$\displaystyle \left[ \left( -1, \ 2\right), \ \left( 2, \ -1\right)\right]$
解のデータ型を辞書に統一する dict = True
solve() のオプション dict = True をつけると，解が全て（解が1個だけであっても）辞書のリストで出力されます。
解のデータ型が統一されると，解の参照の仕方も .subs() で統一できるので，ややこしさが減ります。
In [41]:
solve(x**2 + 2*x - 2, x, dict=True)

Out[41]:
$\displaystyle \left[ \left\{ x : -1 + \sqrt{3}\right\}, \ \left\{ x : – \sqrt{3} – 1\right\}\right]$
In [42]:
solve(x**2 + 2*b*x + c, x, dict=True)

Out[42]:
$\displaystyle \left[ \left\{ x : – b – \sqrt{b^{2} – c}\right\}, \ \left\{ x : – b + \sqrt{b^{2} – c}\right\}\right]$
In [43]:
var('tsuru kame')
solve([Eq(tsuru + kame, 8), 
       Eq(2*tsuru + 4*kame, 26)], 
      [tsuru, kame], dict=True)

Out[43]:
$\displaystyle \left[ \left\{ kame : 5, \ tsuru : 3\right\}\right]$
In [44]:
solve([Eq(x**2 + y**2, 5), Eq(x + y, 1)], 
      [x, y], dict=True)

Out[44]:
$\displaystyle \left[ \left\{ x : -1, \ y : 2\right\}, \ \left\{ x : 2, \ y : -1\right\}\right]$
○練習：鶴と亀と蟻
鶴と亀と蟻，個体数の合計は 10。足の数は全部で 34 本。蟻は亀より 1 匹少ない。鶴，亀，蟻，それぞれ何羽・何匹?
鶴と亀と蟻の個体数をそれぞれ tsuru, kame, ari として連立方程式をたてて solve します。蟻の足は何本かわかりますよね?
かつて線形代数の授業でこの問題を出したとき，真っ先に出た質問が「先生，蟻の足は何本ですか？」でした… orz
これに懲りて，次年度からは問題を一新！「蟻」のかわりに，「かぶと虫」にしたら（なぜか）正答率が上がりました！
In [ ]: