【備忘録】PythonのSympyで(n元連立)(m次)方程式を解く【代数】

Sympyのインストール

Sympyを用いるので、下記等でSympyをインストールしておきます。

pip install sympy

以下のようにして解くことが出来ます。なお、この程度であれば numpy を用いて行列によって解いた方が良いかもしれません。

import sympy

x = sympy.Symbol('x')

y = sympy.Symbol('y')

eq1 = 2*x + 3*y - 16

eq2 = x + 7*y - 19

print(sympy.solve([eq1, eq2]))

# {x: 5, y: 2}

sympy.Symbolによる宣言は一度に行うことが出来ます。

例えば、以下のように10個の変数を一度に定義して、10元連立1次方程式を解くことが出来ます。

import random

import numpy as np

# 変数を一度に宣言

xs = [sympy.Symbol(f'x_{i}') for i in range(10)]

print(xs)

# [x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9]

xs = np.array(xs)

# 連立方程式の係数を生成する関数

def generate_coefficients():

	cs = [random.randint(-10, 10) for i in range(10)]

	cs = np.array(cs)

    return cs

# 連立方程式を生成

eqs = [generate_coefficients() @ xs - random.randint(-10, 10) for j in range(10)]

eqs

# [-8*x_0 + 3*x_1 - 4*x_2 - 8*x_3 + 3*x_4 - 5*x_5 + 6*x_6 - x_7 - 5*x_8 + x_9 + 5,

#  -6*x_0 + 7*x_1 - 9*x_3 - 7*x_4 + 5*x_5 - 5*x_6 - 6*x_7 - 6*x_8 + 7*x_9 - 9,

#  -9*x_0 + 5*x_1 - 3*x_2 + 7*x_3 + 6*x_4 + 8*x_5 + 2*x_6 - 8*x_7 - 8*x_8 + 8*x_9 + 10,

#  3*x_0 + x_1 + 6*x_2 + 8*x_3 - 10*x_4 - x_5 + x_6 - 5*x_7 + 2*x_8 + 10*x_9 + 3,

#  -10*x_0 - 3*x_1 - 10*x_2 + 7*x_3 - 10*x_4 + x_5 + x_6 - 7*x_8 + 4*x_9 - 4,

#  4*x_0 + 2*x_1 + 10*x_2 + 8*x_3 - 7*x_4 - 2*x_5 - 8*x_6 - 8*x_7 + 5*x_8 - 7*x_9 - 8,

#  -9*x_0 + x_1 + 6*x_2 + 3*x_3 + x_4 - x_5 + 7*x_6 + x_7 + 2*x_8 + x_9 + 6,

#  -2*x_0 + 4*x_1 - 3*x_2 - 10*x_3 + x_4 - 3*x_5 - 5*x_6 + 8*x_7 + 2*x_8 + 2*x_9 + 3,

#  4*x_0 - 6*x_1 - 6*x_2 + 2*x_3 + 2*x_4 - 4*x_5 - 10*x_6 - 9*x_7 + 2*x_8 + 8*x_9 + 9,

#  5*x_0 - 10*x_1 - x_2 - x_3 + x_4 - 9*x_5 + 10*x_6 + 3*x_7 + 9*x_8 + 4*x_9]

# 連立方程式を解く

print(sympy.solve(eqs))

# {x_0: 1062532487/18833392579, x_1: 20225194863/37666785158,

# x_2: -52540300375/37666785158, x_3: -13310639315/18833392579,

# x_4: -25553010347/18833392579, x_5: 58772463907/37666785158,

# x_6: 14802784526/18833392579, x_7: -18943314733/18833392579,

# x_8: 38000823649/18833392579, x_9: -21735923399/18833392579}

2次方程式となると numpy では代用しにくいく、sympy の価値が発揮されると思われます。

以下のようにして2次方程式を解くことが出来ます。

x = sympy.Symbol('x')

eq1 = x**2 + 5*x + 6

print(sympy.solve([eq1]))

# [{x: -3}, {x: -2}]

以下のように、実数解がなくとも虚数単位を用いて解いてくれます。

x = sympy.Symbol('x')

eq1 = x**2 + 4*x + 8

print(sympy.solve([eq1]))

# [{x: -2 - 2*I}, {x: -2 + 2*I}]