カテゴリー

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

続きを見る

変数や式を定義(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'>

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

ガジェット

2023/9/18

【デスクツアー2022下半期】モノは少なく、でも効率的に Desk Updating #0

今回はガジェットブロガーなのにデスク環境を構築していない執筆者の ...

ライフハック

2023/9/16

【Audible vs YouTube Premium】耳で聴く音声学習コンテンツを比較

ワイヤレスイヤホンが普及し耳で学習することへのハードルが格段に下 ...

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

2023/9/18

【SENNHEISER MOMENTUM True Wireless 3レビュー】全てが整ったイヤホン

今回は高音質・高機能なSENNHEISERのフラグシップ完全ワイヤレスイヤホン「SENNH ...

ライフハック

2023/3/11

【YouTube Premiumとは】メリットしかないから全員入れ

今回はYouTube Premiumを実際に使ってみてどうなのか、どんなメリット/デメリット ...

マウス

2023/8/17

【Logicool MX ERGOレビュー】疲れない作業効率重視トラックボールマウス

こんな人におすすめ トラックボールマウスの王道Logicool MX ERGOが気になるけどऩ ...

ベストバイ

2023/9/18

【ベストバイ2022】今年買って良かったモノのトップ10

2022年ベストバイ この1年を振り返って執筆者は何を買ったのか。ガジェッ& ...

スマホ

2023/1/15

【楽天モバイル×povo2.0の併用】月1,000円の保険付きデュアルSIM運用

こんな人におすすめ 楽天モバイルとpovo2.0のデュアルSIM運用って実際のとこ ...

マウス

2023/9/16

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

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

macOSアプリケーション

2022/9/30

【Chrome拡張機能】便利で効率的に作業できるおすすめの拡張機能を18個紹介する

こんな人におすすめ Chromeの拡張機能を入れたいけど、調べても同じような ...

macOSアプリケーション

2023/5/3

【Automator活用術】Macで生産性を上げる作業の自動化術

今回はMacに標準でインストールされているアプリ「Automator」を使ってできる ...

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

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

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

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

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

  • この記事を書いた人

メガネ

Webエンジニア駆け出し。独学のPythonで天文学系の大学院を修了。常時金欠のガジェット好きでM2 Pro MacBook Pro(30万円) x Galaxy S22 Ultra(17万円)使いの狂人。自己紹介と半生→変わって楽しいの繰り返しレビュー依頼など→お問い合わせ運営者情報、TwitterX@m_ten_pa、 YouTube@megatenpa、 Threads@megatenpa

-Python基礎
-