odeintサムネ

Python SKILL

Pythonで常微分方程式を解く方法!【scipyのodeintという関数を用いる】

公開日2022/1/8 更新日2022/1/8

颯 Souのイラスト

このブログの権威性

  • 京大院卒業後、就職
  • 現在、作曲家ブロガーのフリーランス
  • 昨年、soublogが14万PV

常微分方程式をPythonで解いてみます。

つまり、

dy/dx=f(x)

の解を求めます。

解析解と、数値解を比較します。

目次は次の通りです。

例題

dy/dx=2x

初期値x=0の時y=0

解析解

dy/dx=2x

y=∫ 2x dx

y=x^2+c

初期値を入れると

c=0

よって

y=x^2

Pythonで常微分方程式の数値解を求める

・パソコンにPythonをインストールする

・pipによるパッケージのインストール

常微分方程式を解くためにscipyのodeintという関数を使います。

※Runge-Kutta法を用いる際は、odeintとは別のodeという関数を使いオプションを指定するそうです。

また、図の表示には、mtplotlibを使います。

まずはpipからダウンロードします。

python -m pip install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

ダウンロードされたら、

pip listで確認してみましょう。

下記が確認できました。

jedi 0.17.2
Jinja2 2.11.2
jsonschema 3.2.0
jupyter 1.0.0
jupyter-client 6.1.7
jupyter-console 6.2.0
jupyter-core 4.6.3
jupyterlab-pygments 0.1.2
kiwisolver 1.2.0
MarkupSafe 1.1.1
matplotlib 3.3.2
mistune 0.8.4
mpmath 1.1.0
nbclient 0.5.0
nbconvert 6.0.7
nbformat 5.0.7
nest-asyncio 1.4.1
nose 1.3.7
notebook 6.1.4
numpy 1.19.2
packaging 20.4
pandas 1.1.3
pandocfilters 1.4.2
parso 0.7.1
pickleshare 0.7.5
Pillow 7.2.0
prometheus-client 0.8.0
prompt-toolkit 3.0.7
pycparser 2.20
Pygments 2.7.1
pyparsing 2.4.7
pyrsistent 0.17.3
python-dateutil 2.8.1
pytz 2020.1
pywin32 228
pywinpty 0.5.7
pyzmq 19.0.2
qtconsole 4.7.7
QtPy 1.9.0
scipy 1.5.2
Send2Trash 1.5.0
six 1.15.0
sympy 1.6.2
terminado 0.9.1
testpath 0.4.4
tornado 6.0.4
traitlets 5.0.4
wcwidth 0.2.5
webencodings 0.5.1
widgetsnbextension 3.5.1

pipはPythonのパッケージを管理するためのツールです。

Pythonをインストールした時点で、公式で自動的にインストールされたものがあります。

それに加えて、今回pipでnumpy scipy matplotlib ipython jupyter pandas sympy noseがあることが確認できます。

ライブラリの保存場所は、

pip show []

で調べられます。

ここで[]にはライブラリ名が入ります。

つまり

pip show scipy

ですね。

・数値プログラム

さて、数値プログラムは、下記微分方程式を解くことでした。

dy/dx=2x

下記、数値プログラムを書きます。

from scipy.integrate import odeint
import matplotlib.pyplot as plt
def equation(y,x):
return 2*x
import numpy as np
x=np.arange(0,10,1) #0から10まで1刻みで解く
y0=0
y=odeint(equation,y0,x)
print(x)
print(y)
plt.plot(x,y)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()

・プログラム実行結果

まずprint(x)の結果ですが、

[0 1 2 3 4 5 6 7 8 9]

print(y)の結果ですが、


[[ 0. ]
[ 1.00000003]
[ 4.00000003]
[ 9.00000003]
[16.00000003]
[25.00000003]
[36.00000003]
[49.00000003]
[64.00000003]
[81.00000003]]

そして、

plt.plot(x,y)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()

の結果ですが、下記となります。

解析解と比較すると、

y=x^2

となっていることが確認できました。

・まとめ

Pythonで常微分方程式を数値的に解きました。

scipy.integrateを使いました。

以上、参考となれば幸いです。

Pythonに戻る

  • この記事を書いた人

颯 Sou

【経歴】京大院卒/大企業/自由を求め脱出/作曲家・ブロガーのフリーランス 【作曲】21曲リリース/AWA動画広告BGMプログラム採用実績あり/Audiostock審査通過/委嘱楽曲の提供経験あり 小学校でピアノを習い、中学校でアコギを買い、高校でエレキを始め、大学で軽音サークルに入り、社会人から作曲家でボカロPとして音楽活動をしています。 【ブログ】soublog500記事/昨年13万PV 社会人になってからブログ制作を開始。最初はボカロ活動を広めることを目的に執筆。「ボカロPになるには」というキーワードで検索順位1位に。また途中から自己啓発を発信し、「自己啓発 一覧」というキーワードで1位を獲得経験あり。SEOを考慮したブログノウハウも発信中。

-Python, SKILL
-,