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

常微分方程式をPythonで解いた数値解の図。

・あいさつ

【経歴】
京大院卒業・航空宇宙工学を学ぶ/大企業メーカーでエンジニア/
自由を求め脱出/作曲家⇄ブロガー⇄Webデザイン/
STEAM教育の大切さに気づきSTEAM+M Lab オンライン研究所を設立/作曲家で研究所所長兼任
【経過】
20曲リリース/AWA動画広告BGMプログラム採用実績あり
soublog417記事/月1万5千PV達成

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

つまり、

$\frac{dy}{dx}=f(x)$

の解を求めます。

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

例題

$\frac{dy}{dx}=2x$

初期値$x=0$の時$y=0$

解析解

$\frac{dy}{dx}=2x$

$y=\int 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

ですね。

・数値プログラム

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

$\frac{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()

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

常微分方程式をPythonで解いた結果の図。

解析解と比較すると、

$y=x^2$

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

・まとめ

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

scipy.integrateを使いました。

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

Pythonに戻る

スポンサードリンク

TechAcademyとは?
プログラミングやアプリ開発を学べる600社、30,000名を超える教育実績があるオンラインスクール。
1人では続かない方のための短期集中プログラム「オンラインブートキャンプ」を開催。
現役のプロのサポートと独自の学習システムで短期間で成長できます。
最短4週間で未経験からプロを育てるオンライン完結のスクールです、どこかに通う必要なく、自宅でもプログラミングやアプリ開発を学ぶことができます。

人工知能に関するスキルを習得したいなら下記。
機械学習のライブラリを使って実装を行います。
実務に近い形で学習し、社会でも通用するスキルが身につきます。
⓵回帰モデルの作成
②クラスタリングの実装
③Amazonレビューを評価分析
④手書き数字の画像認識

機械学習の基礎を習得したいなら下記。
期間内で4つの機械学習プログラムを開発します。
⓵画像を解析して分類
②データから花の形を分類
③住宅価格の分析と予想
④ビットコインの価格変動を予想
実務に近い形で学習することで、社会でも通用するスキルが身につきます。

統計学の基礎・データ分析の手法を習得するなら下記。

Pythonのライブラリを使って実装を行います。
実務に近い形で学習することで、社会でも通用するスキルが身につきます。

⓵区間推定、仮説検定による母集団の検証
②住宅価格の予測
③化学合成実験のデータを分散分析
④サッカーの勝敗予測