用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()

結果：

 

版權聲明：本文爲博主原創文章，未經博主容許不得轉載。 https://blog.csdn.net/u011702002/article/details/78118857