Posted at

[Pythonによる科学・技術計算] 1階常微分方程式の数値解法,初期値問題,数値計算

More than 1 year has passed since last update.

numpy, sympy, scipyを利用して,

1階常微分方程式を初期条件の下で解く。

問題:

$x'(t)+k x(t)=0$, 初期条件 $x(0)=1$, (k=1はパラメータ)

厳密解は$x(t)=e^{-t}$である。

from sympy import *      # sympyライブラリから全ての機能をimoprt

import numpy as np # numpyをnpという名前でimport
from scipy.integrate import odeint
import matplotlib.pyplot as plt

"""
例題: 1階常微分方程式:
dx/dt = -kxを x(0)=1の初期条件の下で解く
"""

x=Symbol('x') # 文字 'x'を変数xとして定義
y=Symbol('y') # 文字 'y'を変数yとして定義
t=Symbol('t') # 文字 't'を変数yとして定義
k=Symbol('k')

def func(x, t, k): # dx/dtを関数として定義
dxdt=-k*x
return dxdt

k = 1.0 #パラメータへ値を設定する
x0 = 1.0 # 初期条件

t=np.linspace(0,10,101) # 積分の時間刻み設定: 0から10までを101等分する

sol=odeint(func, x0, t, args=(k,)) # 数値的に微分方程式を解く。結果をsolというリストへ格納する。

#比較のための厳密解(e^-t)をリストexact_solへセット
exact_sol=[]
tt=np.linspace(0,10,101)
for i in tt:
exact_sol.append(exp(-1.0*float(i)))
#

#可視化

plt.plot(t, sol, linewidth=1,label='numerical') # 数値解
plt.plot(tt, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('t', fontsize=18)
plt.ylabel('x', fontsize=18,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


結果

t.png