1、微分方程分類

微分方程是用來描述某一類函數(shù)與其導(dǎo)數(shù)之間關(guān)系的方程，其解是一個符合方程的函數(shù)。

微分方程按自變量個數(shù)可分為常微分方程和偏微分方程

常微分方程（ODE：ordinary differential equation）

偏微分方程（兩個以上的自變量）

2、微分方程解析解

具備解析解的ODE（常微分方程），我們可以利用SymPy庫進(jìn)行求解

以求解阻尼諧振子的二階ODE為例，其表達(dá)式為：

Demo代碼

import sympy
 
 
def apply_ics(sol, ics, x, known_params):
    free_params = sol.free_symbols - set(known_params)
    eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))]
    sol_params = sympy.solve(eqs, free_params)
    return sol.subs(sol_params)
 
 
# 初始化打印環(huán)境
sympy.init_printing()
# 標(biāo)記參數(shù)，且均為正
t, omega0, gamma = sympy.symbols("t, omega_0, gamma", positive=True)
# 標(biāo)記x是微分函數(shù)，非變量
x = sympy.Function("x")
# 用diff()和dsolve得到通解 
# ode 微分方程等號左邊的部分，等號右邊為0
ode = x(t).diff(t, 2) + 2 * gamma * omega0 * x(t).diff(t) + omega0 ** 2 * x(t)
ode_sol = sympy.dsolve(ode)
# 初始條件：字典匹配
ics = {x(0): 1, x(t).diff(t).subs(t, 0): 0}
x_t_sol = apply_ics(ode_sol, ics, t, [omega0, gamma])
sympy.pprint(x_t_sol)

運(yùn)行結(jié)果：

3、微分方程數(shù)值解

當(dāng)ODE無法求得解析解時，可以用scipy中的integrate.odeint求數(shù)值解來探索其解的部分性質(zhì),并輔以可視化，能直觀地展現(xiàn) ODE解的函數(shù)表達(dá)。

以如下一階非線性（因為函數(shù)y冪次為2）ODE為例：

現(xiàn)用odeint求其數(shù)值解

3.1 場線圖與數(shù)值解

Demo代碼

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

def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None):
    f_np = sympy.lambdify((x, y_x), f_xy, 'numpy')
    x_vec = np.linspace(x_lim[0], x_lim[1], 20)
    y_vec = np.linspace(y_lim[0], y_lim[1], 20)

    if ax is None:
        _, ax = plt.subplots(figsize=(4, 4))

    dx = x_vec[1] - x_vec[0]
    dy = y_vec[1] - y_vec[0]

    for m, xx in enumerate(x_vec):
        for n, yy in enumerate(y_vec):
            Dy = f_np(xx, yy) * dx
            Dx = 0.8 * dx**2 / np.sqrt(dx**2 + Dy**2)
            Dy = 0.8 * Dy*dy / np.sqrt(dx**2 + Dy**2)
            ax.plot([xx - Dx/2, xx + Dx/2], [yy - Dy/2, yy + Dy/2], 'b', lw=0.5)

    ax.axis('tight')
    ax.set_title(r"$%s$" %(sympy.latex(sympy.Eq(y_x.diff(x), f_xy))), fontsize=18)

    return ax

x = sympy.symbols('x')
y = sympy.Function('y')
f = x-y(x)**2

f_np = sympy.lambdify((y(x), x), f)
## put variables (y(x), x) into lambda function f.
y0 = 1
xp = np.linspace(0, 5, 100)
yp = integrate.odeint(f_np, y0, xp)
## solve f_np with initial conditons y0, and x ranges as xp.
xn = np.linspace(0, -5, 100)
yn = integrate.odeint(f_np, y0, xn)

fig, ax = plt.subplots(1, 1, figsize=(4, 4))
plot_direction_field(x, y(x), f, ax=ax)
## plot direction field of function f
ax.plot(xn, yn, 'b', lw=2)
ax.plot(xp, yp, 'r', lw=2)
plt.show()

運(yùn)行結(jié)果：

3.2 洛倫茲曲線與數(shù)值解

以求解洛倫茲曲線為例，以下方程組代表曲線在xyz三個方向上的速度，給定一個初始點，可以畫出相應(yīng)的洛倫茲曲線：

Demo代碼

import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
 
 
def dmove(Point, t, 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.001)
P1 = odeint(dmove, (0., 1., 0.), t, args=([10., 28., 3.],))
P2 = odeint(dmove, (0., 1.01, 0.), t, args=([10., 28., 3.],))
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(P1[:, 0], P1[:, 1], P1[:, 2])
ax.plot(P2[:, 0], P2[:, 1], P2[:, 2])
plt.show()

運(yùn)行結(jié)果：

4、傳染病模型

模型一：SI-Model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.25
# gamma為恢復(fù)率系數(shù)
gamma = 0
# I_0為感染者的初始人數(shù)
I_0 = 1
# S_0為易感者的初始人數(shù)
S_0 = N - I_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,I_0)


def funcSI(inivalue,_):
    Y = np.zeros(2)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0] * X[1]) / N + gamma * X[1]
    # 感染個體變化
    Y[1] = (beta * X[0] * X[1]) / N - gamma * X[1]
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSI,INI,T_range)


plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'red',label = 'Infection',marker = '.')
plt.title('SI Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

模型二：SIS model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.25
# gamma為恢復(fù)率系數(shù)
gamma = 0.05
# I_0為感染者的初始人數(shù)
I_0 = 1
# S_0為易感者的初始人數(shù)
S_0 = N - I_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,I_0)


def funcSIS(inivalue,_):
    Y = np.zeros(2)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0]) / N * X[1] + gamma * X[1]
    # 感染個體變化
    Y[1] = (beta * X[0] * X[1]) / N - gamma * X[1]
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSIS,INI,T_range)

plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'red',label = 'Infection',marker = '.')
plt.title('SIS Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

模型三：SIR model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.25
# gamma為恢復(fù)率系數(shù)
gamma = 0.05
# I_0為感染者的初始人數(shù)
I_0 = 1
# R_0為治愈者的初始人數(shù)
R_0 = 0
# S_0為易感者的初始人數(shù)
S_0 = N - I_0 - R_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,I_0,R_0)


def funcSIR(inivalue,_):
    Y = np.zeros(3)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0] * X[1]) / N
    # 感染個體變化
    Y[1] = (beta * X[0] * X[1]) / N - gamma * X[1]
    # 治愈個體變化
    Y[2] = gamma * X[1]
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSIR,INI,T_range)


plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'red',label = 'Infection',marker = '.')
plt.plot(RES[:,2],color = 'green',label = 'Recovery',marker = '.')
plt.title('SIR Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

模型四：SIRS-Model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.25
# gamma為恢復(fù)率系數(shù)
gamma = 0.05
# Ts為抗體持續(xù)時間
Ts = 7
# I_0為感染者的初始人數(shù)
I_0 = 1
# R_0為治愈者的初始人數(shù)
R_0 = 0
# S_0為易感者的初始人數(shù)
S_0 = N - I_0 - R_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,I_0,R_0)


def funcSIRS(inivalue,_):
    Y = np.zeros(3)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0] * X[1]) / N + X[2] / Ts
    # 感染個體變化
    Y[1] = (beta * X[0] * X[1]) / N - gamma * X[1]
    # 治愈個體變化
    Y[2] = gamma * X[1] - X[2] / Ts
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSIRS,INI,T_range)


plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'red',label = 'Infection',marker = '.')
plt.plot(RES[:,2],color = 'green',label = 'Recovery',marker = '.')
plt.title('SIRS Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

模型五：SEIR-Model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.6
# gamma為恢復(fù)率系數(shù)
gamma = 0.1
# Te為疾病潛伏期
Te = 14
# I_0為感染者的初始人數(shù)
I_0 = 1
# E_0為潛伏者的初始人數(shù)
E_0 = 0
# R_0為治愈者的初始人數(shù)
R_0 = 0
# S_0為易感者的初始人數(shù)
S_0 = N - I_0 - E_0 - R_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,E_0,I_0,R_0)


def funcSEIR(inivalue,_):
    Y = np.zeros(4)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0] * X[2]) / N
    # 潛伏個體變化
    Y[1] = (beta * X[0] * X[2]) / N - X[1] / Te
    # 感染個體變化
    Y[2] = X[1] / Te - gamma * X[2]
    # 治愈個體變化
    Y[3] = gamma * X[2]
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSEIR,INI,T_range)


plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'orange',label = 'Exposed',marker = '.')
plt.plot(RES[:,2],color = 'red',label = 'Infection',marker = '.')
plt.plot(RES[:,3],color = 'green',label = 'Recovery',marker = '.')

plt.title('SEIR Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

模型六：SEIRS-Model

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

# N為人群總數(shù)
N = 10000
# β為傳染率系數(shù)
beta = 0.6
# gamma為恢復(fù)率系數(shù)
gamma = 0.1
# Ts為抗體持續(xù)時間
Ts = 7
# Te為疾病潛伏期
Te = 14
# I_0為感染者的初始人數(shù)
I_0 = 1
# E_0為潛伏者的初始人數(shù)
E_0 = 0
# R_0為治愈者的初始人數(shù)
R_0 = 0
# S_0為易感者的初始人數(shù)
S_0 = N - I_0 - E_0 - R_0
# T為傳播時間
T = 150

# INI為初始狀態(tài)下的數(shù)組
INI = (S_0,E_0,I_0,R_0)


def funcSEIRS(inivalue,_):
    Y = np.zeros(4)
    X = inivalue
    # 易感個體變化
    Y[0] = - (beta * X[0] * X[2]) / N + X[3] / Ts
    # 潛伏個體變化
    Y[1] = (beta * X[0] * X[2]) / N - X[1] / Te
    # 感染個體變化
    Y[2] = X[1] / Te - gamma * X[2]
    # 治愈個體變化
    Y[3] = gamma * X[2] - X[3] / Ts
    return Y

T_range = np.arange(0,T + 1)

RES = spi.odeint(funcSEIRS,INI,T_range)


plt.plot(RES[:,0],color = 'darkblue',label = 'Susceptible',marker = '.')
plt.plot(RES[:,1],color = 'orange',label = 'Exposed',marker = '.')
plt.plot(RES[:,2],color = 'red',label = 'Infection',marker = '.')
plt.plot(RES[:,3],color = 'green',label = 'Recovery',marker = '.')

plt.title('SEIRS Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

結(jié)語

參考：

https://www.bilibili.com/video/BV12h411d7Dm
https://zhuanlan.zhihu.com/p/104091330

學(xué)習(xí)來源：

B站及其課堂PPT
對其中代碼進(jìn)行了復(fù)現(xiàn)

「文章僅作為學(xué)習(xí)筆記，記錄從0到1的一個過程」

希望對您有所幫助，如有錯誤歡迎小伙伴指正～

Python數(shù)學(xué)建模系列（五）：微分方程

1、微分方程分類

2、微分方程解析解

3、微分方程數(shù)值解

3.1 場線圖與數(shù)值解

3.2 洛倫茲曲線與數(shù)值解

4、傳染病模型

模型一：SI-Model

模型二：SIS model

模型三：SIR model

模型四：SIRS-Model

模型五：SEIR-Model

模型六：SEIRS-Model

結(jié)語

1、微分方程分類

2、微分方程解析解

4、傳染病模型