3. 常微分方程式(scipy)

3.1. 数値積分(scipy.integrate.quad)

次の積分を計算する.

\[\int_{0}^{1} 3x^2 dx = 1\]
from scipy.integrate import quad

def f(x, a):
    return a*x*x
a=3.0
min=0.0
max=1.0
res, err = quad(f, min, max, args=(a,))
print("res", res,"err", err)

実行結果は

('res', 1.0, 'err', 1.1102230246251565e-14)

となる.

3.2. 一階の常微分方程式(scipy.integrate.odeint)

一階の微分方程式

\[\frac{dx(t)}{dt} = -a x\]

を解く.一般解は初期条件を \(x(0)=x_0\) として

\[x(t) = x_0 e^{-at}\]

である.

from scipy.integrate import odeint
import numpy as np
import pylab

def diff_op(x, time,a):
    return -a*x

x0 = 2.0
a = 2.0

time = np.linspace(0,4,40)
traject = odeint(diff_op, x0 ,time, args=(a,))

pylab.plot(time, traject, '-k', linewidth=3)
pylab.xlabel('t', fontsize=24)
pylab.ylabel('x', fontsize=24,rotation='horizontal')
pylab.show()
fig1 fig2
左図:normal plot 右図: semi-log plot