一.微分方程
0.微分方程分类
微分方程是用来描述某一类函数与其子数了可大个料万在其解是一个符合方程的函数。微分方程按自变量个数可分为常微分方程和偏微分方程。
sympy学习库:www.tutorialspoint.com/sympy/
1.微分方程解析解
代码如下:
import numpy as np
import sympy
# apply_ics:计算特解
# sol:通解
# ics:初始条件
# x:自变量
# known_params:参数字典
# subs:符号计算函数
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)
sympy.init_printing() #初始化打印环境
t, omega0, gamma = sympy.symbols("t, omega0, gamma", positive=True) #标记变量
x = sympy.Function('x') # 标记x是微分函数,非变量
ode = x(t).diff(t,2)+2*gamma*omega0*x(t).diff(t)+omega0**2*x(t) # 求解的微分表达式
ode_sol = sympy.dsolve(ode) # 用diff()和disolve得到通解
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)
运行结果:
2.微分方程数值解
2.1场线图与数值解
场线图代码如下:
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,y_x)两个自变量参与到f_xy规则运算(隐函数)
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)
y0 = 1
xp = np.linspace(0,5,100)
yp = integrate.odeint(f_np,y0,xp)
xn = np.linspace(0,-5,100)
yn = integrate.odeint(f_np,y0,xp)
fig,ax=plt.subplots(1,1,figsize=(4,4))
ax=plot_direction_field(x,y(x),f,ax=ax)#绘制f场线图
ax.plot(xn,yn,'b',lw=2)
ax.plot(xp,yp,'r',lw=2)
plt.show()
场线图运行结果:
2.2洛伦兹曲线与数值解
以求解洛伦兹曲线为例,以下方程组代表曲线在xyz三个方向上的速度,给定一个初始点,可以画出相应的洛伦兹曲线:
代码如下:
import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
# odeint()函数需要至少三个变量,第一个是微分方程函数,第二个是微分方程初值,第三个是微分的自变量。
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()
结果如下:
3.传染病模型
传染病模型属于传染病动力学方向
本节涉及传染病模型包括:SI、SIS、SIR、SIRS、SEIR、SEIRS共六个模型。
3.1SI-Model
代码如下:
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N为人群总数
beta = 0.25 # β为传染率系数
gamma = 0 # gamma为恢复率系数
I_0 =1 #I_0为感染者的初始人数
S_0 = N - I_0 # S_0为易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, I_0) # INI为初始状态
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()
结果如下:
3.2SIS-Model
代码如下:
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N为人群总数
beta = 0.25 # β为传染率系数
gamma = 0.05 # gamma为恢复率系数
I_0 =1 #I_0为感染者的初始人数
S_0 = N - I_0 # S_0为易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, I_0) # INI为初始状态
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()
结果如下:
3.3SIR-Model
代码如下:
import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N为人群总数
beta = 0.25 # β为传染率系数
gamma = 0.05 # gamma为恢复率系数
I_0 =1 #I_0为感染者的初始人数
R_0 = 0#R_0为易治愈者的初始人数
S_0 = N-I_0-R_0 # 易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, I_0, R_0) # INI为初始状态
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 = 'Recoveryn', marker = '.')
plt.title('SIR Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()
结果如下:
3.4SIRS-Model
代码如下:
import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt
# N为人群总数
N = 10000
# β为传染率系数
beta = 0.25
# gamma为恢复率系数
gamma = 0.05
# Ts为抗体持续时间
Ts = 7
# I_0为感染者的初始人数
I_0 = 1
# R_0为治愈者的初始人数
R_0 = 0
# S_0为易感者的初始人数
S_0 = N - I_0 - R_0
# T为传播时间
T = 150
# INI为初始状态下的数组
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()
结果如下:
3.5SEIR-Model
代码如下:
import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt
# N为人群总数
N = 10000
# β为传染率系数
beta = 0.6
# gamma为恢复率系数
gamma = 0.1
# Te为疾病潜伏期
Te = 14
# I_0为感染者的初始人数
I_0 = 1
# E_0为潜伏者的初始人数
E_0 = 0
# R_0为治愈者的初始人数
R_0 = 0
# S_0为易感者的初始人数
S_0 = N - I_0 - E_0 - R_0
# T为传播时间
T = 150
# INI为初始状态下的数组
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()
结果如下:
3.6SEIRS-Model
代码如下:
import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt
# N为人群总数
N = 10000
# β为传染率系数
beta = 0.6
# gamma为恢复率系数
gamma = 0.1
# Ts为抗体持续时间
Ts = 7
# Te为疾病潜伏期
Te = 14
# I_0为感染者的初始人数
I_0 = 1
# E_0为潜伏者的初始人数
E_0 = 0
# R_0为治愈者的初始人数
R_0 = 0
# S_0为易感者的初始人数
S_0 = N - I_0 - E_0 - R_0
# T为传播时间
T = 150
# INI为初始状态下的数组
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()
结果如下: