カテゴリー

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

運営者のメガネとです。YouTubeTwitterInstagramも運営中。

自己紹介はこちらから、お問い合わせはこちらからお願いいたします。

運営者メガネ

下準備

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のグラフ即席作成コード

こんな人にオススメ 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'>

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

スイッチボット

2022/11/28

【SwitchBotロックレビュー】これからのスタンダードになりうるスマートロック

こんな人にオススメ SwitchBotからスマートロック「SwitchBotロック」が発売された ...

生活に役立つ

2022/11/28

【メガネ厳選】クソ便利に使っているサービスやアイテム達

このページでは執筆者「メガネ」が実際に使って便利だと感じているサ ...

マウス

2022/9/11

【Logicool MX ERGO vs MX Master 3】ERGOをメインにした決定的な理由

こんな疑問・お悩みを持っている人におすすめ 執筆者はLogicoolのハイエンӠ ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【ながら聴きイヤホン比較】SONY LinkBuds、ambie、BoCoはどれがおすすめ?

こんな人におすすめ 耳を塞がない開放型のイヤホンに完全ワイヤレスӟ ...

macOSアプリケーション

2022/10/15

【M1 Mac】MacBook Proに入れている便利でニッチなアプリを21個紹介する

こんな人におすすめ MacBookを購入してLINEとか必要最低限のアプリは入れた。 ...

完全ワイヤレスイヤホン(TWS)

2022/10/23

【SENNHEISER MOMENTUM True Wireless 3レビュー】高レベルでバランス型の高音質イヤホン

こんな人におすすめ SENNHEISER MOMENTUM True Wireless 3って実際のところどうなの? 評判は良い ...

完全ワイヤレスイヤホン(TWS)

2022/11/21

【SONY WF-1000XM4レビュー】神とゴミのハーフ&ハーフ

こんな人におすすめ SONYのフラグシップモデル「SONY WF-1000XM4」ってどれくらい性 ...

完全ワイヤレスイヤホン(TWS)

2022/8/19

【Nothing ear (1)レビュー】ライトな完成度、アップデートに期待

こんな人にオススメ 完全ワイヤレスイヤホン(TWS)でスケルトンボディ ...

Pythonを学びたいけど独学できる時間なんてない人へのすゝめ

執筆者は大学の研究室・大学院にて独学でPythonを習得した。

でも社会人になったら独学で行うには時間も体力もなくて大変だ。

時間がない社会人だからこそプロの教えを乞うのが効率的。

ここでは色んなタイプに合ったプログラミングスクールの紹介をする。

  • この記事を書いた人

メガネ

ベンチャー企業のWebエンジニア駆け出し。独学のPythonで天文学系の大学院を修了→新卒を1.5年で辞める→転職→今に至る。
常時金欠のガジェット好きでM1 MacBook Pro x Galaxy S22 Ultraの狂人。
人見知りで根暗だったけど、人生楽しもうと思って良い方向に狂う→人生が楽しい

ガジェットのレビューとPythonコードを記事にしています。ぜひ楽しんでください🦊
自己紹介と半生→変わって楽しいの繰り返し

-Python基礎
-