Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
2
Help us understand the problem. What are the problem?

More than 3 years have passed since last update.

posted at

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

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

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
2
Help us understand the problem. What are the problem?