[Python에 의한 과학·기술 계산] 1층 상미분 방정식의 수치 해법, 초기값 문제, 수치 계산

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()

결과