用python實現解常微分方程組的簡單示例以及用odeint解常微分方程的範例
阿新 • • 發佈:2019-01-30
背景:
包括兩個部分,一個是演示怎麼自己寫程式碼解常微分方程,另一部分就是示範python怎麼呼叫解常微分方程的函式。 下面的方程組給出洛侖茲引子在三個方向上的速度,求解運動軌跡軟體:
python3.5.2部分1:(div)
程式碼:
# -*- coding: utf8 -*- import numpy as np """ 移動方程: t時刻的位置P(x,y,z) steps:dt的大小 sets:相關引數 """ def move(P,steps,sets): x,y,z = P sgima,rho,beta = sets #各方向的速度近似 dx = sgima*(y-x) dy = x*(rho-z)-y dz = x*y - beta*z return [x+dx*steps,y+dy*steps,z+dz*steps] #設定sets引數 sets = [10.,28.,3.] t = np.arange(0,30,0.01) #位置1: P0 = [0.,1.,0.] P = P0 d = [] for v in t: P = move(P,0.01,sets) d.append(P) dnp = np.array(d) #位置2: P02 = [0.,1.01,0.] P = P02 d = [] for v in t: P = move(P,0.01,sets) d.append(P) dnp2 = np.array(d) """ 畫圖 """ from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt fig = plt.figure() ax = Axes3D(fig) ax.plot(dnp[:,0],dnp[:,1],dnp[:,2]) ax.plot(dnp2[:,0],dnp2[:,1],dnp2[:,2]) plt.show()
結果:
部分2:
程式碼:
# -*- coding: utf-8 -*- import numpy as np from scipy.integrate import odeint """ 定義常微分方程,給出各方向導數,即速度 """ def dmove(Point,t,sets): """ p:位置向量 sets:其他引數 """ p,r,b = sets x,y,z = Point return np.array([p*(y-x),x*(r-z),x*y-b*z]) t = np.arange(0,30,0.01) #呼叫odeint對dmove進行求解,用兩個不同的初始值 P1 = odeint(dmove,(0.,1.,0.),t,args = ([10.,28.,3.],)) #(0.,1.,0.)是point的初值 #([10.,28.,3.],)以元祖的形式給出 point,t之後的引數 P2 = odeint(dmove,(0.,1.01,0.),t,args = ([10.,28.,3.],)) """ 畫3維空間的曲線 """ from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt fig = plt.figure() ax = Axes3D(fig) ax.plot(P1[:,0],P1[:,1],P1[:,2]) ax.plot(P2[:,0],P2[:,1],P2[:,2]) plt.show()