カテゴリー

Python基礎

【sympy&方程式】pythonで方程式とか微積分を解く

2021年7月11日

こんな人にオススメ

pythonで方程式とかって解ける?

もしできるならこれで課題は楽ちんだぜ!

ということで、今回はpythonのsympyを使って方程式を解いたり微分を計算してみる。執筆者が学生時代にこれを知っていたら多少は計算が楽になったのかなと思いつつ、でも結局は手計算の方が早いんじゃないかとも思いつつ。

複雑な式になるとその分、コードも複雑になるがちゃんと書き方が合っていれば自動で計算してくれるのでだいぶ楽になるだろう。

python環境は以下。

  • Python 3.9.4
  • sympy 1.8
  • numpy 1.20.3
  • plotly 4.14.3
  • plotly-orca 3.4.2
  • astropy 4.2.1
スポンサーリンク
スポンサーリンク

運営者のメガネです。TwitterInstagramも運営してます。

自己紹介はこちらから、お問い合わせはこちら。

運営者メガネ

下準備

import sympy
import sys
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
import astropy.units as u
import sympy.physics.units as units

sys.path.append('../../')
import plotly_layout_template as template

まずは下準備として今回使用するモジュールをimport。途中でplotlyを使用してグラフを作成するのでplotly用のテンプレートも同時にimport。テンプレートに関しては以下参照。

【随時更新 備忘録】plotlyのグラフ即席作成コード

続きを見る

変数や式を定義(Symbol, subs)

まずは変数や式を定義する。今まではx=10とかにしてxには値が入っていたが、変数として扱うことで値を決定せずxのまま使用することが可能。

変数を定義

# 変数を設定
x = sympy.Symbol('x')
y = sympy.Symbol('y')

print(x)
# x
print(type(x))
# <class 'sympy.core.symbol.Symbol'>
print(y)
# y
print(type(y))
# <class 'sympy.core.symbol.Symbol'>

変数を定義。変数っていうのは以下の式にあるxとかaとかのこと。

$y\ =\ a\ \times\ x\ +10$

今まではprint(x)とするとxが定義されていないのでエラーが出たが、変数として定義してあげることでxのまま使用することが可能。

今回はxyを中心に使用し、必要になった際に適宜、変数を追加する。

式を定義

eq = x ** 3 + x + y + 10
print(eq)
# x**3 + x + y + 10
print(type(eq))
# <class 'sympy.core.add.Add'>

eq2 = x * y
print(eq2)
# x*y

変数を定義できたら、この変数を使用して式を書いてみる。そうすると変数部分はそのままにした式が表示される。

なお、jupyter notebookなどを使用している場合は、sympy.init_printing()を初めに動かし、printを使用せずに出力するかdisplay(eq)のようにするとキレイなな数式で表現することが可能。

特殊な関数を使用する場合はsympy.(関数名)というようにすると使用可能。以下ではsin関数を使用。

# sin関数

# # 後にかけたい数値を置くと以下のproblemが出る
# # Operator "*" not supported for "None"Pylancere(reportOptionalOperand
# sin = sympy.sin(x) * 2

# 先にかけたい数値を置けばproblemはでない
sin = 2 * sympy.sin(x)
print(sin)
# 2*sin(x)
print(type(sin))
# <class 'sympy.core.mul.Mul'>

なお、式の後に掛け算を書くとPylanceから怒られた。これはVScodeに入っている拡張機能か何かで、エラーではないがずっと赤線を引かれるのが嫌なので先に2をかけた。

変数に値を代入

eq = x ** 3 + x + y + 10
print(eq)
# x**3 + x + y + 10

# xに1を代入
subs = eq.subs(x, 1)
print(subs)
# y + 12
print(type(subs))
# <class 'sympy.core.add.Add'>

変数を設定し、式にその変数を使用したら続いては代入。好きな変数に値を代入することが可能。上の例ではx1を代入。自動で式が整理される。

また、yxを含む式を代入することも可能。この場合は既に定義されたxと計算されて整理される。また、式を代入する際に掛け算などがある場合は自動でカッコをつけて計算してくれる。

# yにx + 3を代入
subs = eq.subs(y, x + 3)
print(subs)
# x**3 + 2*x + 13

eq = 2 * x + y + 10
# xに(y + 3)を代入
subs = eq.subs(x, y + 3)
print(subs)
# 3*y + 16

虚数を代入することも可能。虚数はsympy.Iとすることで使用することが可能。

eq = x ** 3 + x ** 2 + 2 * x + 10
# 虚数はsympy.Iで使用可能
subs = eq.subs(x, sympy.I)
print(subs)
# 9 + I

式の展開、因数分解、整理(expand, factor, simplify)

式の定義と代入が完了したら次は展開、因数分解、そして式の整理。中一とか小六とかの内容。

式を展開

# 自動では展開されない
eq = (x + 3)**2
print(eq)
# (x + 3)**2
print(type(eq))
# <class 'sympy.core.power.Pow'>

# 式を展開
expand = sympy.expand(eq)
print(expand)
# x**2 + 6*x + 9
print(type(expand))
# <class 'sympy.core.add.Add'>

式は自動では展開してくれない。なので、展開してほしかったらsympy.expandとする。

因数分解

eq = x ** 2 + 6 * x + 9
print(eq)
# x**2 + 6*x + 9

# 因数分解
factor = sympy.factor(eq)
print(factor)
# (x + 3)**2
print(type(factor))
# <class 'sympy.core.power.Pow'>

eq = x ** 2 * y + x ** 2 + 6 * x * y + 6 * x + 9 * y + 9
print(eq)
# x**2*y + x**2 + 6*x*y + 6*x + 9*y + 9
factor = sympy.factor(eq)
print(factor)
# (x + 3)**2*(y + 1)

展開の反対。ごちゃついた式をカッコでくくってスッキリさせるにはsympy.factorとする。

式を整理

eq1 = (x + 1) ** 2 + (x - y) * (x + y) + 2 * (5 * x + y)
eq2 = - (x + 5) ** 2 + y * (y - 2) + 25
eq = eq1 + eq2
print(eq)
# 10*x + y*(y - 2) + 2*y + (x + 1)**2 - (x + 5)**2 + (x - y)*(x + y) + 25

# 式の整理
simplify = sympy.simplify(eq)
print(simplify)
# x**2 + 2*x + 1
print(type(simplify))
# <class 'sympy.core.add.Add'>

# 因数分解してさらに整理
print(sympy.factor(simplify))
# (x + 1)**2

因数分解のようにカッコを使用せず、単に同じ次元の変数同士を計算することも可能。sympy.simplyで整理する。もちろん、因数分解できるのであればsympy.factorでさらに整理することも可能。

方程式を解く(solve)

式の展開や因数分解が完了したら、お次は方程式を解く。ここでは簡単なものを用いるが、もちろん高次元方程式も解くことが可能。注意点は(式) = 0の形式で書く必要があるということ。

$y\ =\ ax\ +\ b$

なら

$ax\ +\ b\ -\ y$

と入力する必要がある。

1次方程式

# (式)= 0として扱い、解を計算

# 1次方程式
eq = x + 5
solve = sympy.solve(eq)
print(solve)
# [-5]

# 複数変数ある場合はlistの中にdictの形で出力される
eq = x - y + 1
solve = sympy.solve(eq)
print(solve)
# [{x: y - 1}]
print(type(*solve))
# <class 'dict'>

方程式を解くにはsympy.solveとすればいい。出力はlistで複数の変数がある場合はdictで自動で選んだ変数に対する解を出力する。どの変数に対して解くかを指定するにはsympy.solveの第2引数にその変数名を入れる。

# 解を計算したい変数を指定することが可能
solve = sympy.solve(eq, x)
print(solve)
# [y - 1]

# 解を計算したい変数を指定することが可能
solve = sympy.solve(eq, y)
print(solve)
# [x + 1]

2次方程式

eq = (x + 1) ** 2
solve = sympy.solve(eq)
print(solve)
# [-1]
print(type(solve))
# <class 'list'>

eq = (x + 3) * (x - 1)
solve = sympy.solve(eq)
print(solve)
# [-3, 1]

# 解に虚数が出る時はIで表示
eq = x ** 2 + x + 1
solve = sympy.solve(eq)
print(solve)
# [-1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]

二次方程式も基本的には同じ。小数は使用されないので、小数ではなく分数やsqrt表記で表される。

2変数の連立方程式

# 式は=0の形に変形する必要がある
eq1 = 3 * x + 5 - y
eq2 = 10 * x - 2 - y
solve = sympy.solve([eq1, eq2])
print(solve)
# {x: 1, y: 8}
print(type(solve))
# <class 'dict'>

2変数の連立方程式の場合はsympy.solvelistからtupleで式を入れたらいい。出力はdict。また、2乗とかになって解が複数個ある場合はdictlistで表現される。

eq1 = 2 * x ** 2 + 3 * x - 4 - y
eq2 = 5 * x - y
solve = sympy.solve((eq1, eq2))
print(solve)
# [{x: -1, y: -5}, {x: 2, y: 10}]
print(type(solve))
# <class 'list'>

以下にこれらの式のグラフを示す。それぞれのグラフの交点が解に相当する。

関連コンテンツ

スポンサーリンク

Amazonのお買い物で損したない人へ

1回のチャージ金額通常会員プライム会員
¥90,000〜2.0%2.5%
¥40,000〜1.5%2.0%
¥20,000〜1.0%1.5%
¥5,000〜0.5%1.0%

Amazonギフト券にチャージすることでお得にお買い物できる。通常のAmazon会員なら最大2.0%、プライム会員なら2.5%還元なのでバカにならない。

ゲットしたポイントは通常のAmazonでのお買い物に使えるからお得だ。一度チャージしてしまえば、好きなタイミングでお買いものできる。

なお、有効期限は10年だから安心だ。いつでも気軽にAmazonでお買い物できる。

Amazonチャージはここから出来るで

もっとお得なAmazon Prime会員はこちらから

30日間無料登録

執筆者も便利に使わせてもらってる

スポンサーリンク

  • この記事を書いた人

メガネ

独学でpythonを学び天文学系の大学院を修了。 ガジェット好きでMac×Android使い。色んなスマホやイヤホンを購入したいけどお金がなさすぎて困窮中。 元々、人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂ったために今も人生めちゃくちゃ楽しい。 pythonとガジェットをメインにブログを書いていますので、興味を持たれましたらちょこちょこ訪問してくだされば幸いです🥰。 自己紹介→変わって楽しいの繰り返し

-Python基礎
-