pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

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

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

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

pianofisica.hatenablog.com



代数方程式

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

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

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

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

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

1次方程式

 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}}

が出力されます。

2次方程式


 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

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

3次方程式


 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}  }

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

4次方程式

 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)

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

5次方程式

 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)

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

方程式の数値解


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

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

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

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

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

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

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

連立1次方程式



2変数の場合

2つの変数  x y に対する連立1次方程式

 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

を得ます。

3変数の場合

同様に、3つの変数  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)} }}

と計算されます。

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

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

pianofisica.hatenablog.com

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

連立方程式


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

# 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)

より数列が

 0, ~1, ~1, ~2, ~3, ~5, ~8, ~13, ~21, ~34, ~\dots

と求まります。

ちなみに、この漸化式によって定まる数列をFibonacci(フィボナッチ)数列といいます。

整数列が与えられたとき、その数列がどのような式で表されるのか知りたいときに便利なWebページ

The On-Line Encyclopedia of Integer Sequences® (OEIS®)

を紹介しておきます。数列の一部分を取り出して検索バーに順番通りに入力して検索すると、データベースと照合して、その数列の他の項や一般式、漸化式などの候補を示してくれます。試しに、上で求めた数列をコピー&ペーストして検索するとFibonacci数の説明が出てきます:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34 - OEIS


Python(SymPy)を使ってもう少し高等な計算を取り扱っている記事が

pianofisica.hatenablog.com

で、Bernoulli数(ベルヌーイ数)を母関数のTaylor展開から求めています。

またさらに(数学的にもPython的にも)発展した話題として、Bernoulli多項式を解説した記事が

pianofisica.hatenablog.com

です。SymPyをNumPy、Matplotlibと組み合わせて使う方法については

pianofisica.hatenablog.com

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


また、本記事で述べた事柄は、Maximaで同様のことをしている過去記事

pianofisica.hatenablog.com

の中のMaximaの入力を、PythonのSymPyライブラリの書法で置き換えたものになっています。

あわせて読んでみると、PythonMaximaの違いを見比べられるかと思います。


キーワードPython、SymPy、NumPy、JupyterNotebook