Python (SymPy) で方程式・連立方程式を解く、数列を求める

今回はプログラミング言語 Pythonを使って方程式・連立方程式を解いてみたいと思います。数式処理ライブラリSymPyを使って代数的に厳密に解く方法をみていきます。また、漸化式から定まる数列について、その各項を求める方法もみていきます。

SymPyのごく基本的な使い方については以下の記事を参照してください：

pianofisica.hatenablog.com

代数方程式
方程式の数値解
連立１次方程式
- ２変数の場合
- ３変数の場合
連立方程式
微分方程式
漸化式

代数方程式

多項式（変数と定数について有限回の和と積をとったもの）同士を等号で結んだ方程式を代数方程式といいます。代数方程式の解として表される数を代数的数といいます。したがって整数、有理数、そして $\sqrt{2}$ などの無理数は代数的数です。

代数的数でない数は超越数と呼ばれます。同様に、代数的でない方程式のことも超越方程式と呼ばれます。超越方程式の例としては

$1-\log x=0 \qquad \sin x=0$

などがあげられます。これらの方程式の解となるNapier数（自然対数の底） $e$ 、円周率 $\pi$ は超越数です。

以下では代数方程式をPythonの数式処理ライブラリSymPyを使って解く方法をみていきます。

１次方程式

$ax+b=0$

をPythonで解くには

import sympy
sympy.var('x, a, b')
Sol1=sympy.solve (a*x+b, x)
sympy.init_printing()
display(Sol1[0])  # sympy.solve(eq, var) で方程式 eq=0 を変数 varについて解いた根を表示

と入力します。すると解

$\displaystyle{x=-\frac{b}{a}}$

が出力されます。

２次方程式

$ax^2+bx+c=0$

をPythonで解くには

import sympy
sympy.var('x, a, b, c')
Sol2=sympy.solve (a*x**2+b*x+c, x)
sympy.init_printing()
display(Sol2)

と入力します。解として

$\displaystyle{x={-b-{\sqrt{b^2-4ac}}\over{2a}} ,\qquad x={-b+{\sqrt{b^2-4ac}}\over{2a}}}$

が出力されます。高校数学でもおなじみの公式です。

以下、その都度sympyと書くのが面倒なので

import sympy as sp

として関数の呼び出しを少し簡略化しておきます。

３次方程式

$ax^3+bx^2+cx+d=0$

をPythonで解くには

import sympy as sp
sp.init_printing()
sp.var('x, a, b, c, d')
Sol3=sp.solve (a*x**3+b*x**2+c*x+d, x)

display(Sol3[0])

display(Sol3[1])

display(Sol3[2])

と入力します。すると、解

$\displaystyle{x=r^{1/3}-\frac{s}{r^{1/3}}-\frac{b}{3a},\qquad x=\left(\frac{-1\mp i\sqrt{3}}{2}\right)r^{1/3}-\left(\frac{-1\pm i\sqrt{3}}{2}\right)\frac{s}{r^{1/3}}-\frac{b}{3a}\\ r={{\left(27a^2d^2+\left(4b^3-18abc\right)d+4ac^3-b^2c^2\right)^{1/2}}\over{2\,3^{3/2}a^2}}+{{{{bc}}-{{3ad}}}\over{6a^2}}-{{b^3}\over{27a^3}} \\ s=-\frac{b^2}{9a^2}+\frac{c}{3a} }$

が求まります。これが３次方程式の解の公式ということになります。

４次方程式

$ax^4+bx^3+cx^2+dx+e=0$

をPythonで解くには

import sympy as sp
sp.init_printing()
sp.var('x, a, b, c, d, e')
Sol4=sp.solve (a*x**4+b*x**3+c*x**2+d*x+e, x)
display(Sol4)

と入力すればよいのですが、さすがに複雑なのでここでは載せずにおきます。

５次方程式

$ax^4+bx^3+cx^2+dx+e=0$

は一般に代数的に解けないことが知られています。実際

import sympy as sp
sp.init_printing()
sp.var('x, a, b, c, d, e, f')
Sol5=sp.solve (a*x**5+b*x**4+c*x**3+d*x**2+e*x+f, x)
display(Sol5)

とすると何も出力されません。

方程式の数値解

これまで１変数の代数方程式を代数的数として解く方法をみてきましたが、一般の５次以上の方程式は特殊な場合を除いて原理的には解けないものの、応用上、数値的に解がわかれば十分という場合もあります。

裳華房の『数値計算』（柳田・中木・三村）では、数値計算の基本的な手法を紹介しています：

<br />

代数方程式・連立方程式を数値的にいかにして解くかを数学的背景を基礎に解説していて、原理の面から理解したいという人にはうってつけの本だと思います。

さて、具体的な問題を考えて、例えば

$x^5+6x^4+x^3+x+1=0$

は５次方程式で、代数的に解くことはできません：

import sympy as sp
sp.init_printing()
sp.var('x')
sp.solve (x**5+6*x**4+x**3+x+1, x)

ですが数値解であれば、Pythonの数値計算ライブラリであるNumPyを用いて

import numpy
numpy.roots([1, 6, 1, 0, 1, 1])

によって

$x_1=-5.8241071 \\ x_2 =0.42084104+0.52540961\, i \\ x_3= 0.42084104-0.52540961\,i \\ x_4=-0.50878749+0.34645119\,i \\x_5=-0.50878749-0.34645119\,i$

と求めることができます。詳しくは次の記事を読んでください：

pianofisica.hatenablog.com

連立１次方程式

２変数の場合

２つの変数 $x$ と $y$ に対する連立１次方程式

$ax+by=e \\ cx+dy=f$

をPythonで解くには

import sympy as sp
sp.init_printing()
sp.var('x, y, a, b, c, d, e, f')
eq1=sp.Eq(a*x+b*y, e)
eq2=sp.Eq(c*x+d*y, f)
sp.solve ([eq1, eq2], [x, y])

と入力します。すると解

$\displaystyle{x={{de-bf}\over{ad-bc}} ,\qquad y={{af-ce}\over{ad-bc}}}$

が出力されます。例えば

$2x+y=3 \\ x+3y=4$

に対しては

sp.var('x, y')
eq3=sp.Eq(2*x+1*y, 3)
eq4=sp.Eq(1*x+3*y, 4)
sp.solve ([eq3, eq4], [x, y])

より解

$x=1, \quad y=1$

を得ます。

３変数の場合

同様に、３つの変数 $x$ 、 $y$ 、 $z$ の場合

$ax+by+cz=j \\ dx+ey+fz=k \\ gx+hy+iz=l$

は

sp.var('x, y, z, a, b, c, d, e, f, g, h, i, j, k, l')
eq1=sp.Eq(a*x+b*y+c*z, j)
eq2=sp.Eq(d*x+e*y+f*z, k)
eq3=sp.Eq(g*x+h*y+i*z, l)
sp.solve ([eq1, eq2, eq3], [x, y, z])

より

$\displaystyle{x={{j(ei-fh)-k(bi-ch)+l(bf-ce)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ y={{−𝑗(𝑑𝑖−𝑓𝑔)+𝑘(𝑎𝑖−𝑐𝑔)−𝑙(𝑎𝑓−𝑐𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ z={{𝑗(𝑑ℎ−𝑒𝑔)−𝑘(𝑎ℎ−𝑏𝑔)+𝑙(𝑎𝑒−𝑏𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)} }}$

と計算されます。

多変数の連立１次方程式は行列によっても表すことができ、

その解は係数行列の逆行列を求める操作によっても求めることができます。それは次の記事

pianofisica.hatenablog.com

で解説していますので、あわせて読んでみてください。

連立方程式

ここまでは１変数の代数方程式、多変数の連立１次方程式についてみてきましたが、多変数の連立方程式についても

# sympy.solve([eq1, eq2, ...], [var1, var2, ...]) で連立方程式 eq1, eq2, ... の解を表示

として解くことができます。

例えば、放物線 $y=ax^2+bx+c$ と直線 $y=gx+h$ の交点は連立方程式

$y=ax^2+bx+c \\ y=gx+h$

の解として求めることができますが、これをPythonで解くには

sp.var('x, y, a, b, c, g, h')
eq5=sp.Eq(y, a*x**2+b*x+c)
eq6=sp.Eq(y, g*x+h)
sp.solve ([eq5, eq6], [x, y])

と入力します。

微分方程式

未知関数とその導関数がみたす関係式として書かれる方程式を微分方程式といいます。

SymPyでの微分方程式の解析的な取り扱いについては次の記事

pianofisica.hatenablog.com

で解説しています。数値的な取り扱いについては

pianofisica.hatenablog.com

で解説しています。あわせて読んでいただけたらさいわいです。

漸化式

$F_{n+2}=F_{n+1}+F_n \qquad F_0=0 \qquad F_1 =1$

のような漸化式から定まる数列について、その各項をPythonで求めることができます：

def n(k):
    if k == 0:
        return 0
    if k == 1:
        return 1
    else:
        return n(k-1) + n(k-2) 
N_list = []
for i in range(10):
    N = n(i)
    N_list.append(N)
print(N_list)