1. 程式人生 > >用python實現解常微分方程組的簡單示例以及用odeint解常微分方程的範例

用python實現解常微分方程組的簡單示例以及用odeint解常微分方程的範例

背景:

包括兩個部分,一個是演示怎麼自己寫程式碼解常微分方程,另一部分就是示範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()

結果: