微分方程_数值解_解析解_场线图_传染病模型

序言：

微分方程解析解

import numpy as np
from scipy import integrate
import sympy

#sol:通解 ics:初始条件 x:自变量 known_params:已知参数
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() 
#标记参数 postive表示参数为正
t, omega0, gamma = sympy.symbols("t, omega_0, gamma", positive=True) 
#标记微分函数x
x = sympy.Function('x') 
#给出求解方程 x(t).diff(t,2)是x关于t的二阶导数
ode = x(t).diff(t,2)+2*gamma*omega0*x(t).diff(t)+omega0**2*x(t)
#算出通解
ode_sol = sympy.dsolve(ode)
#给定初始条件 x(0)=0 x的一阶偏导在t等于0时为0
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.微分方程数值解

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import sympy
#x_lim:x的可视化范围 #ylim: y的可视化范围 ax:是否提供背景图
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)#符号表达式转隐函数
y0 = 1
xp = np.linspace(0,5,100)
yp = integrate.odeint(f_np, y0, xp)#初始y0解f_np,x范围xp
xn = np.linspace(0,-5,100)
yn = integrate.odeint(f_np, y0, xp)
fig, ax = plt.subplots(1,1,figsize=(4,4))
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.1场线图与数值解

2.2洛伦兹曲线与数值解

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

3.传染病模型

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

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

SIR 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 # S_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 = 'Recovery',marker = '.')
plt.title('SIR Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

SIRS 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为恢复率系数
Ts = 7 # Ts为抗体持续时间
I_0 = 1 # I_0为感染者的初始人数
R_0 = 0 # R_0为治愈者的初始人数
S_0 = N - I_0 - R_0 # S_0为易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, I_0, R_0) # INI为初始状态下的数组
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()

SIER Model:

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N为人群总数
beta = 0.6 # β为传染率系数
gamma = 0.1 # gamma为恢复率系数
Te = 14 # Te为疾病潜伏期
I_0 = 1 # I_0为感染者的初始人数
E_0 = 0 # E_0为潜伏者的初始人数
R_0 = 0 # R_0为治愈者的初始人数
S_0 = N - I_0 - R_0 - E_0 # S_0为易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, E_0, I_0, R_0) # INI为初始状态下的数组
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('SETR Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

SETRS Model:

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N为人群总数
beta = 0.6 # β为传染率系数
gamma = 0.1 # gamma为恢复率系数
Ts = 7 # Ts为抗体持续时间
Te = 14 # Te为疾病潜伏期
I_0 = 1 # I_0为感染者的初始人数
E_0 = 0 # E_0为潜伏者的初始人数
R_0 = 0 # R_0为治愈者的初始人数
S_0 = N - I_0 - R_0 - E_0 # S_0为易感染者的初始人数
T = 150 # T为传播时间
INI = (S_0, E_0, I_0, R_0) # INI为初始状态下的数组
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('SETRS Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()