Scipyで常微分方程式を解く

Scipyのodeint関数を使い、常微分方程式を解いてみる。

簡単な例として、

を、初期値

で解く。
真の解は

である。

以下がpythonスクリプト。xの0から3まで積分している。数値計算で得られた値を＋で、真の解より得られる値を実線で図示している。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt 

#solve dy/dx=y

#define derivative
def dy_dx(y,x):
    return y

#integration
x1=np.linspace(0,3, 10)
x=np.linspace(0,3, 100)
y1, infodict = integrate.odeint(dy_dx, 1., x1, full_output=True)
print infodict

# plot
plt.plot(x1,y1,'+',markersize=12)  #solusion 
plt.plot(x,np.exp(x))            #true solution 
plt.savefig("ode_ex.png")
plt.show()

参照:
LoktaVolterraTutorial
http://www.scipy.org/LoktaVolterraTutorial 
A Coupled Spring-Mass System
http://www.scipy.org/Cookbook/CoupledSpringMassSystem