こんな人にオススメ
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
下準備
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のまま使用することが可能。
今回はx
とy
を中心に使用し、必要になった際に適宜、変数を追加する。
式を定義
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'>
変数を設定し、式にその変数を使用したら続いては代入。好きな変数に値を代入することが可能。上の例ではx
に1
を代入。自動で式が整理される。
また、y
にx
を含む式を代入することも可能。この場合は既に定義された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.solve
にlist
からtuple
で式を入れたらいい。出力はdict
。また、2乗とかになって解が複数個ある場合はdict
のlist
で表現される。
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'>
以下にこれらの式のグラフを示す。それぞれのグラフの交点が解に相当する。
+ クリックでオープン
# 確認用のグラフ
def graph(equation1, equation2, save_name: str):
def scatter(x, y, name, color: str):
d = go.Scatter(
mode='lines+markers',
x=x, y=y, name=name,
marker=dict(symbol='circle'),
line=dict(dash='solid', color=color),
hovertemplate=f"{name}<br>"
+ 'x: %{x} <br>'
+ 'y: %{y} <br>'
+ "<extra></extra>",
)
return d
# 横軸の値
yoko = np.arange(-10, 10)
plot = []
zipper = zip((equation1, equation2), ('red', 'blue'))
for num, (equation, color) in enumerate(zipper, 1):
# sympy.subsだと配列を代入できないので新規で定義
func = sympy.lambdify([x, y], equation)
# xには横軸の値yokoを、yは「y=」として計算したいので0を代入
tate = func(yoko, 0)
d = scatter(x=yoko, y=tate, name=f"eq{num}", color=color)
plot.append(d)
# eq1でのtateの値
# [166 131 100 73 50 31 16 5 -2 -5 -4 1 10 23 40 61 86 115
# 148 185]
# eq2でのtateの値
# [-50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35
# 40 45]
layout = go.Layout(template=template.plotly_layout(),)
fig = go.Figure(data=plot, layout=layout)
fig.show(config=template.plotly_config())
pio.write_html(
fig, f"{save_name}.html",
auto_open=False, config=template.plotly_config()
)
pio.write_image(fig, f"{save_name}.png")
graph(equation1=eq1, equation2=eq2, save_name='simultaneous')
微分系(diff, dsolve)
ここから急に年代が上がって高校数学。微分方程式については特殊解は求めずに一般解のみとする。
微分
# 微分
eq = 2 * x ** 2 + 3 * x - 4
diff = sympy.diff(eq)
print(diff)
# 4*x + 3
print(type(diff))
# <class 'sympy.core.add.Add'>
微分はsympy.diff
で行う。ここではx
で微分した。もちろん三角関数も微分可能。
# 三角関数の微分も対応
sin = 2 * sympy.sin(x)
diff = sympy.diff(sin)
print(diff)
# 2*cos(x)
2変数の場合は変数を指定しないとエラーになる。
# 2変数の場合は変数を選択
eq = x ** 3 + x + 3 * y ** 2 + 10
# どの変数かを指定しないとエラー
diff = sympy.diff(eq)
# Since there is more than one variable in the expression, the
# variable(s) of differentiation must be supplied to differentiate x**3
# + x + 3*y**2 + 10
diff = sympy.diff(eq, x)
print(diff)
# 3*x**2 + 1
diff = sympy.diff(eq, y)
print(diff)
# 6*y
微分方程式
# fを関数として定義
f = sympy.Function('f')
print(f)
# f
print(type(f))
# <class 'sympy.core.function.UndefinedFunction'>
# 定数としてaを定義
a = sympy.Symbol('a')
# 微分方程式を解く
# Object of type "None" cannot be calledPylancere(portOptionalCall)のエラーあり
eq = sympy.diff(f(x)) - a * f(x)
print(eq)
# -a*f(x) + Derivative(f(x), x)
print(type(eq))
# <class 'sympy.core.add.Add'>
dsolve = sympy.dsolve(eq)
print(dsolve)
# Eq(f(x), C1*exp(a*x))
print(type(dsolve))
# <class 'sympy.core.relational.Equality'>
微分方程式を解くにはdiff solveという意味だろう、sympy.dsolve
を使用する。微分方程式と計算結果のEq(f(x), C1*exp(a*x))
はJupyter notebookなどでprint
を使用せずに出力すると以下のようにキレイな数式となる。
$\begin{aligned} \cfrac{d}{dx}f(x)\ &=\ af(x)\\ f(x)\ &=\ C_1e^{ax} \end{aligned}$
解が2定数を含む式の演算なら自動的に$C_2$のように定数が決められる。
eq = sympy.diff(f(x), x, 2) - a ** 2 * f(x)
print(eq)
# -a**2*f(x) + Derivative(f(x), (x, 2))
dsolve = sympy.dsolve(eq)
print(dsolve)
# Eq(f(x), C1*exp(-a*x) + C2*exp(a*x))
この時の微分方程式とその解は以下。
$\begin{aligned} \cfrac{d^2}{dx^2}f(x)\ &=\ a^2f(x)\\ f(x)\ &=\ C_1e^{-ax}\ +\ C_2e^{ax} \end{aligned}$
また、解けない(と自分は思っている)微分方程式についてはエラーが吐かれる。
# 解けるかどうかは知らんけど、エラーが出る微分方程式もある
eq = sympy.diff(f(x), x, 2) - x**2 + f(x) ** 4
print(eq)
# -x**2 + f(x)**4 + Derivative(f(x), (x, 2))
dsolve = sympy.dsolve(eq)
# raise NotImplementedError(dummy + "solve" + ": Cannot solve " + str(eq))
# NotImplementedError: solve: Cannot solve -x**2 + f(x)**4 + Derivative(f(x), (x, 2))
この時の微分方程式は以下。
$\cfrac{d^2}{dx^2}f(x)\ =\ -x^2\ +\ f^4(x)$
積分系(integrate)
微分ときたら積分と思うのは私だけだろうか。簡単にざっくりイメージは微分の逆。
不定積分
eq = x ** 2 + 6 * x + 9
integrate = sympy.integrate(eq)
# 積分定数は表記されない
print(integrate)
# x**3/3 + 3*x**2 + 9*x
print(type(integrate))
# <class 'sympy.core.add.Add'>
積分についてはsympy.integrate
で可能。不定積分であるが、積分定数(通常は$C$が使われる)は書かれないのでその点には注意。今回の積分は以下。
$\displaystyle\int dx\ x^2\ +\ 6x\ +\ 9$
もちろん三角関数も可能だが、積分結果が$\tan(x)$ではないのが気になる。
eq = 1 / (sympy.cos(x)**2)
integrate = sympy.integrate(eq)
# tan(x)にはなってくれない
print(integrate)
# sin(x)/cos(x)
定積分
eq = x ** 2 + 6 * x + 9
integrate = sympy.integrate(eq, (x, 0, 3))
print(integrate)
# 63
print(type(integrate))
# <class 'sympy.core.numbers.Integer'>
定積分を計算したい場合はsympy.integrate
の引数に初めの値と終わりの値を入れればいい。今回は0
から3
まで積分してみた。式は以下。
$\displaystyle\int_{0}^{3} dx\ x^2\ +\ 6x\ +\ 9$
計算結果は普通に足し算などが可能。
# 普通に計算可能
print(integrate + 1)
# 64
単位を付与(sympy.physics.units )
最後は単位について。というのも執筆者自身、普段からastropy
で数値に単位をつけている。単位をつけることでメートルからミリメートルのような単位換算を自動計算してくれて楽。
astropy
と単位については以下参照。
【astropy&単位】astropyで変数に単位を付与
続きを見る
【astropy&単位】astropyの単位出力書式
続きを見る
astropy
の単位だとエラー
# 変数に単位を付与することはできない
m = sympy.Symbol('m') * u.s
# m = sympy.Symbol('m') * u.s
# TypeError: unsupported operand type(s) for *: 'Symbol' and 'IrreducibleUnit'
今まで通りastropy
を使用して単位を付与してみる。まずは変数にそのまま単位を付与してみるがエラー。代入する値に単位をつけて代入しようとしてもエラー。
# 代入するときに単位を付与してもエラー
eq = 2 * x + 5
val = 3 * u.s
subs = eq.subs(x, val)
# RecursionError: maximum recursion depth exceeded in comparison
しかも変数ではない項に単位を付与してもエラー。まじですか。
# 変数ではない項に単位を付与してもエラー
eq = 2 * x + 5 * u.s
# RecursionError: maximum recursion depth exceeded in comparison
と思ったらそもそも単位ありと単位なしを同時に式に書くこと自体無理だった。
eq = 2 + 5 * u.s
# astropy.units.core.UnitConversionError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)
代入完了後に単位を付与することは可能。
# 代入後の単位の付与は可能
eq = 2 * x + 5
val = 3
subs = eq.subs(x, val) * u.s
print(subs)
# 11.0 s
しかし、これだと単位換算の恩恵を享受することができない。一応、sympy
での単位換算系はsymdim
というモジュールで可能らしいがanacondaになかったので執筆者は入れていない。以下に参考URLを書いておく。
sympy
の単位だと大丈夫
# sympyの単位を使用することが可能
val = 12 * units.m
print(val)
# 12*meter
print(type(val))
# <class 'sympy.core.mul.Mul'>
じゃあ単位をつけられないのかというとそういうわけではない。astropy
と使わなければいい。調べているとちょうどsympy
で単位を付与できる機能があるっぽいのでこれを使ってみる。
モジュールのimport
はfrom
を使う例が結構あったが、個人的にはimport sympy.physics.units as units
の方がしっくりきたのでこれでimport
。units.m
でメートル単位を使用することが可能。
sympyの単位付与の状態で代入を行うとエラーは出ない。
# sympyの単位だとエラーが起きることはない
eq = 2 * x + 5
subs = eq.subs(x, val)
print(subs)
# 24*meter + 5
print(type(subs))
# <class 'sympy.core.add.Add'>
単位換算はunits.convert_to(変換したい単位つき変数, 変換先単位)
で可能。ここではメートルをミリメートル、センチメートル、メートルのまま、そしてキロメートルに変換した。小数が使えないっぽいので数字の表記は分数となっている。
# 単位つきの値の単位変換。小数は使用されず、分数表記となる
print(units.convert_to(val, units.mm))
# 12000*millimeter
print(units.convert_to(val, units.cm))
# 1200*centimeter
print(units.convert_to(val, units.m))
# 12*meter
print(units.convert_to(val, units.km))
# 3*kilometer/250
複数の単位を組み合わせて、単位換算をすることも可能。ここではジュールJと秒sからワットWに変換した。
# ジュールJと秒sからワットWを計算
val1 = 10 * units.s
val2 = 100 * units.J
print(val1)
# 10*second
print(val2)
# 100*joule
val = val2 / val1
print(val)
# 10*joule/second
print(type(val))
# <class 'sympy.core.mul.Mul'>
print(units.convert_to(val, units.W))
# 10*watt
pythonに解いてもらう
今回はsympy
を使って変数を定義し、方程式などを解くという方法について解説した。執筆者自身、pythonを学んだのは研究室配属された大学4年生からでこの頃にはもうほとんどの単位を取得してあまり方程式を使用する機会がなかった。
また、大学院の講義でも数学系というよりは宇宙の概論などを主に学んだので今回のsympy
を使用する機会はあまりなかった。
しかし、早い時期にこの記事を見ていただいた方は、sympy
を使って方程式などを解いてもらうことで検算が可能だ(くれぐれもpythonの結果をそのまま解答にして力をつけないままなのは避けるように)。
今回のような便利な機能をまた見つけたら記事にしていく。
関連記事
【python&初級】のlistとかforとかifとかまとめ
続きを見る
【python&関数化】defとかargsとかを使って関数を作成する
続きを見る
【python3&四捨五入】四捨五入すると1.5も2.5も2.0になる問題
続きを見る
【python&csv読み込み】pythonを使ってcsvを読み込み
続きを見る