常微分方程式

SciPy常微分方程式

SciPy モジュールを使った常微分方程式の数値計算の方法を示します．

はじめに

odeint とは

Python の scipy パッケージの odeint モジュールを使うと，とても簡単に微分方程式の数値解を得ることができます．このパッケージが計算することができる方程式は，1階の常微分方程式です．1階であれば，連立微分方程式の計算も可能です．したがって，高階の常微分方程式の計算も可能と言うことです．どんな高階の微分方程式でも，1階の連立微分方程式に変換できるからです．これについては，「高階の常微分方程式」を参照願います．

odeint のアルゴリズムは，

使い方

1階の常微分方程式の数値解を得るためには，(1)常微分方程式と(2)初期条件が必要です．

とりあえず使ってみよう

まず簡単な微分方程式を odeint を使い，数値計算します．例として，

\begin{align} \ddiffA{y}{x}=\cos(x) \end{align}

の解を数値計算で求めます．初期条件は，\(y(0)=0\) とします．

常微分方程式を数値計算するプログラム (odeint)

001   # -*- coding: utf-8 -*-
002   
003   import numpy as np
004   import matplotlib.pyplot as plt
005   from scipy.integrate import odeint
006   
007   # ===========================================================
008   # 常微分方程式を解くクラス
009   # ===========================================================
010   class ODE(object):
011   
012       # -------------------------------------------------------
013       # コンストラクター
014       # -------------------------------------------------------
015       def __init__(self, diff_eq, init_con):
016   
017           self.diff_eq  = diff_eq
018           self.init_con = init_con
019   
020           
021       # -------------------------------------------------------
022       # 常微分方程式の計算
023       # -------------------------------------------------------
024       def cal_euation(self, x_min, x_max, N):
025   
026           x = np.linspace(x_min, x_max, N)             # x の配列の生成
027           y = odeint(self.diff_eq, self.init_con, x)   # 方程式の計算
028   
029           return x, y
030   
031   
032   
033   # -------------------------------------------------------
034   # 解くべき常微分方程式
035   # -------------------------------------------------------
036   def diff_eq(y, x):
037       dydx = np.cos(x)
038       return dydx
039   
040   
041   # -------------------------------------------------------
042   # プロット
043   # -------------------------------------------------------
044   def plot(x, y, x_range, y_range):
045       fig = plt.figure()
046   
047       # ----- プロットの準備 -----
048       sol = fig.add_subplot(1,1,1)
049       sol.set_xlabel("x", fontsize=20, fontname='serif')
050       sol.set_ylabel("y", fontsize=20, fontname='serif')
051       sol.tick_params(axis='both', length=10, which='major')
052       sol.tick_params(axis='both', length=5,  which='minor')
053       sol.set_xlim(x_range)
054       sol.set_ylim(y_range)
055       sol.minorticks_on()
056       sol.plot(x, y, 'b-', markersize=5)
057   
058       # ----- スクリーン表示 -----
059       fig.tight_layout()
060       plt.show()
061           
062       # ----- pdf 作成 -----
063       fig.savefig('ode_solve.pdf', orientation='portrait', \
064                   transparent=False, bbox_inches=None, frameon=None)
065       fig.clf()
066   
067   
068   
069   # -------------------------------------------------------
070   # メイン関数
071   # -------------------------------------------------------
072   if __name__ == "__main__":
073   
074       N = 1000                              # 分割数
075       min_x = 0                             # x の最小
076       max_x = 4*np.pi                       # x の最大
077       initial_condition = np.array([0])     # 初期条件
078   
079       ode = ODE(diff_eq, initial_condition)
080       x, y = ode.cal_euation(min_x, max_x, N)
081   
082       plot(x, y, (min_x, max_x), (-1.2, 1.2))
083

図1: 常微分方程式を数値計算で解いた結果

このプログラムの動作は，以下のとおりです．Matplotlib でプロットしているコマンドは，009 — 010 に書かれています．

リスト行	動作の説明
003	numpy モジュールをインポートします．これがなくても，微分方程式の数値解を計算する odeint は動作します．しかし，計算に便利なので，odeint と一緒に読み込まれることが多いです．
004	matplotlib.pyplot をインポートします．これがなくても，odeint は動作します．しかし，解の確認に便利なので一緒に使います．
015 — 018	コンストラクターです．解くべき微分方程式(コールバック関数で与えられた導関数)と初期条件を読み込んでいます．

プロット作成については，Matplotlib のページを参照ください．

いろいろな常微分方程式

連立微分方程式

カオスの研究で有名な「ローレンツ方程式」を数値計算します．解くべき微分方程式は，

\begin{align} &\ddiffA{x}{t}=-px+py& &\ddiffA{y}{t}=-xz+rx-y& &\ddiffA{z}{t}=xy-bz \end{align}

です．これを変数は \((p,\,q,\,b)=(10,\,28,\,8/3)\) で，初期値 \((x,\,y,\,z)=(0.1,\,0.1\,0.1)\) で数値計算します．odeint を使った Python のプログラム例を以下に示します．

ローレンツ方程式を数値計算するプログラム (Lorenz)

001   # -*- coding: utf-8 -*-
002   
003   import numpy as np
004   import matplotlib.pyplot as plt
005   from mpl_toolkits.mplot3d import Axes3D
006   from scipy.integrate import odeint
007   
008   # ===========================================================
009   # 常微分方程式を解くクラス
010   # ===========================================================
011   class ODE(object):
012   
013       # -------------------------------------------------------
014       # コンストラクター
015       # -------------------------------------------------------
016       def __init__(self, diff_eq, init_con):
017   
018           self.diff_eq  = diff_eq
019           self.init_con = init_con
020   
021           
022       # -------------------------------------------------------
023       # 常微分方程式の計算
024       # -------------------------------------------------------
025       def cal_euation(self, t_min, t_max, N):
026   
027           t = np.linspace(t_min, t_max, N)             # x の配列の生成
028           v = odeint(self.diff_eq, self.init_con, t)   # 方程式の計算
029   
030           return v
031   
032   
033   
034   # -------------------------------------------------------
035   # 解くべき常微分方程式
036   # -------------------------------------------------------
037   def diff_eq(v, t):
038       p = 10
039       r = 28
040       b = 8/3
041       dxdt = -p*v[0] + p*v[1]
042       dydt = -v[0]*v[2] + r*v[0] -v[1]
043       dzdt = v[0]*v[1] -b*v[2]
044       return [dxdt, dydt, dzdt]
045   
046   
047   # -------------------------------------------------------
048   # プロット
049   # -------------------------------------------------------
050   def plot(x, y, z):
051       fig = plt.figure()
052   
053       # ----- プロットの準備 -----
054       sol = fig.add_subplot(1,1,1, projection='3d')
055       sol.set_xlabel("x", fontsize=20, fontname='serif')
056       sol.set_ylabel("y", fontsize=20, fontname='serif')
057       sol.set_zlabel("z", fontsize=20, fontname='serif')
058       sol.tick_params(axis='both', length=10,  which='major')
059       sol.tick_params(axis='both', length=5,  which='minor')
060       sol.minorticks_on()
061   
062       sol.plot(x, y, z, 'b-', markersize=0)
063   
064       # ----- スクリーン表示 -----
065       fig.tight_layout()
066       plt.show()
067           
068       # ----- pdf 作成 -----
069       fig.savefig('Lorenz.pdf', orientation='portrait', \
070                   transparent=False, bbox_inches=None, frameon=None)
071       fig.clf()
072   
073   
074   
075   # -------------------------------------------------------
076   # メイン関数
077   # -------------------------------------------------------
078   if __name__ == "__main__":
079   
080       N = 10001                                # 分割数
081       min_t = 0                                # t の最小
082       max_t = 100                              # t の最大
083       initial_condition = np.array([0.1, 0.1, 0.1])  # 初期条件
084   
085       ode = ODE(diff_eq, initial_condition)
086       v = ode.cal_euation(min_t, max_t, N)
087   
088       plot(v[:, 0], v[:,1], v[:,2])
089

図2: 常微分方程式(ローレンツ方程式)を数値計算で解いた結果

高階の微分方程式

これまで計算した常備分方程式はいずれも1階でしたが，ここでは2階の微分方程式を計算します．例として，調和振動子を表す微分方程式

\begin{align} \ddiffB{2}{x}{t}=-x \end{align}

を計算します．この2階の微分方程式は，

\begin{align} &\ddiffA{x}{t}=v& &\ddiffA{v}{t}=-x& \end{align}

と書きなおすことができます．1階の連立微分方程式なので，odeint で計算することができます．初期条件

odeint の詳細

クラス scipy.integrate.odeint の引数は，次のとおりです．

scipy.integrate.odeint(
      func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None,
      mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0,
      hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5,
      printmessg=0)

func: callable(y, t0, …): ここで，計算する1階の常微分方程式を指定します．導関数を返すコールバック関数です．
y0: array: 初期条件を指定します．解となる関数の初期値です．連立微分方程式の場合，複数の初期値(ベクター)を記述します．
t: array: 解の関数の独立変数のシーケンスです．先に示した初期値は，このシーケンスの先頭と対応します．

ページ作成情報

参考資料

モジュール「pyplot.py」のソースコードが良い資料です．

更新履歴

2017年08月2日

新規作成