Python数学建模系列（五）：微分方程
2024-05-29 17:14:51
文章目录
- 前言
- 往期文章
- 1、微分方程分类
- 2、微分方程解析解
- 3、微分方程数值解
- 3.1 场线图与数值解
- 3.2 洛伦兹曲线与数值解
- 4、传染病模型
- 模型一：SI-Model
- 模型二：SIS model
- 模型三：SIR model
- 模型四：SIRS-Model
- 模型五：SEIR-Model
- 模型六：SEIRS-Model
- 结语
前言
Hello！小伙伴！
非常感谢您阅读海轰的文章，倘若文中有错误的地方，欢迎您指出～
自我介绍 ଘ(੭ˊᵕˋ)੭
昵称：海轰
标签：程序猿｜C++选手｜学生
简介：因C语言结识编程，随后转入计算机专业，有幸拿过一些国奖、省奖…已保研。目前正在学习C++/Linux/Python
学习经验：扎实基础 + 多做笔记 + 多敲代码 + 多思考 + 学好英语！
初学Python 小白阶段
文章仅作为自己的学习笔记 用于知识体系建立以及复习
题不在多 学一题 懂一题
知其然 知其所以然！
往期文章
Python数学建模系列（一）：规划问题之线性规划
Python数学建模系列（二）：规划问题之整数规划
Python数学建模系列（三）：规划问题之非线性规划
Python数学建模系列（四）：数值逼近
1、微分方程分类
微分方程是用来描述某一类函数与其导数之间关系的方程，其解是一个符合方程的函数。
微分方程按自变量个数可分为常微分方程和偏微分方程
常微分方程（ODE：ordinary differential equation）
偏微分方程（两个以上的自变量）
2、微分方程解析解
具备解析解的ODE（常微分方程），我们可以利用SymPy库进行求解
以求解阻尼谐振子的二阶ODE为例，其表达式为：
Demo代码
import sympydef 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, omega_0, gamma", positive=True)
# 标记x是微分函数，非变量
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)
运行结果：
3、微分方程数值解
当ODE无法求得解析解时，可以用scipy中的integrate.odeint求 数值解来探索其解的部分性质,并辅以可视化，能直观地展现 ODE解的函数表达。
以如下一阶非线性（因为函数y幂次为2）ODE为例：
现用odeint求其数值解
3.1 场线图与数值解
Demo代码
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
import sympydef 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) * dxDx = 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 axx = sympy.symbols('x')
y = sympy.Function('y')
f = x-y(x)**2f_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()
运行结果：
3.2 洛伦兹曲线与数值解
以求解洛伦兹曲线为例，以下方程组代表曲线在xyz三个方向 上的速度，给定一个初始点，可以画出相应的洛伦兹曲线：
Demo代码
import numpy as np
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as pltdef dmove(Point, t, sets):p, r, b = setsx, y, z = Pointreturn 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()
运行结果：
4、传染病模型
模型一：SI-Model
import scipy.integrate as spi
import numpy as np
import matplotlib.pyplot as plt# N为人群总数
N = 10000
# β为传染率系数
beta = 0.25
# gamma为恢复率系数
gamma = 0
# I_0为感染者的初始人数
I_0 = 1
# S_0为易感者的初始人数
S_0 = N - I_0
# T为传播时间
T = 150# INI为初始状态下的数组
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 YT_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为人群总数
N = 10000
# β为传染率系数
beta = 0.25
# gamma为恢复率系数
gamma = 0.05
# I_0为感染者的初始人数
I_0 = 1
# S_0为易感者的初始人数
S_0 = N - I_0
# T为传播时间
T = 150# INI为初始状态下的数组
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 YT_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为人群总数
N = 10000
# β为传染率系数
beta = 0.25
# gamma为恢复率系数
gamma = 0.05
# 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 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 YT_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为人群总数
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] / Tsreturn YT_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为人群总数
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 YT_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为人群总数
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] / Tsreturn YT_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()
结语
参考：
- https://www.bilibili.com/video/BV12h411d7Dm
- https://zhuanlan.zhihu.com/p/104091330
学习来源：B站及其课堂PPT，对其中代码进行了复现
文章仅作为学习笔记，记录从0到1的一个过程
希望对您有所帮助，如有错误欢迎小伙伴指正～
我是 海轰ଘ(੭ˊᵕˋ)੭
如果您觉得写得可以的话，请点个赞吧
谢谢支持 ❤️
Python数学建模系列（五）：微分方程相关推荐
- Python数学建模系列（六）：蒙特卡洛算法
文章目录 前言 往期文章 1.蒙特卡洛算法 样例1 样例2 样例3 2.三门问题 3.M*M豆问题 结语 前言 Hello!小伙伴! 非常感谢您阅读海轰的文章,倘若文中有错误的地方,欢迎您指出- ...
- Python数学建模系列（八）：图论
文章目录 前言 往期文章 1 图论模型 - Dijkstra 样例1 2 图论模型-Floyd 样例2 3 机场航线设计 0.Airlines.csv数据 1.数据导入.观察变量 2.数据清洗 3.时 ...
- Python数学建模系列（四）：数值逼近
文章目录 前言 往期文章 1. 一维插值 1.1 线性插值与样条插值(B-spline) 1.2 高阶样条插值 2. 二维插值 2.1 图像模糊处理--样条插值 2.2 二维插值的三维图 3. 最小 ...
- Python数学建模系列（二）：规划问题之整数规划
文章目录 前言 整数规划 例题 方法一:分支定界法(使用scipy库) 方法二:使用pulp库进行求解 结语 前言 Hello!小伙伴! 非常感谢您阅读海轰的文章,倘若文中有错误的地方,欢迎您指出- ...
- Python小白的数学建模课-10.微分方程边值问题
小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型边值问题的建模与求解,不涉及算法推导和编程,只探讨如何使用 Pytho ...
- Python小白的数学建模课-09.微分方程模型
小白往往听到微分方程就觉得害怕,其实数学建模中的微分方程模型不仅没那么复杂,而且很容易写出高水平的数模论文. 本文介绍微分方程模型的建模与求解,通过常微分方程.常微分方程组.高阶常微分方程 3个案例手 ...
- Python小白的数学建模课-09 微分方程模型
1. 微分方程 1.1 基本概念 微分方程是描述系统的状态随时间和空间演化的数学工具.物理中许多涉及变力的运动学.动力学问题,如空气的阻力为速度函数的落体运动等问题,很多可以用微分方程求解.微分方程在 ...
- 【Python数学建模】SEIR传染病模型模型延伸-SEIDR模型（一），加入疫苗接种、政府管控、病毒变异等因素的影响
目录 一. SEIR传染病模型 二. SEIR模型的延伸--SEIDR模型 三. 模型延伸--影响因素1:疫苗接种 四. 模型延伸--影响因素2:政府管控 五. 模型延伸--影响因素3:病毒变异 写在 ...
- python可以用来数学建模吗_怎么用Python数学建模：python数据建模工具
怎么用Python数学建模 djcjfhfhhjdvjfhvfghhfgbdthhgdchfjfuivvh DSI方法在几何建模上的应用 本节叙述如何应用DSI方法来与曲面S相联系的二维图形图3.1) ...
最新文章
- 挖坑挖到cnblogs.com来...
- 查看MYSQL数据库中所有用户及拥有权限
- https HttpsURLConnection请求的单向认证
- linux内核加载卡主,请教mx6，linux3.0.35，tf卡能启动uboot但是无法加载内核问题
- LeetCode 117. Populating Next Right Pointers in Each Node II
- 【CV】使用直方图处理进行颜色校正
- Insertion Sort List（单链表插入排序）
- sql2000导出mysql_如何将sql2000的数据库导入到mysql中？
- python构建huffman树_python：哈夫曼树,PythonHuffuman
- 老挑毛 win7 linux,图解老挑毛u盘启动工具怎么重装系统
- 高等代数第3版下 [丘维声 著] 2015年版_一文搞懂代数几何发展史（一）
- 通达OA2017版工作流触发器应用实例
- 大规模额外涨薪后，Intel又准备了24亿美元，明年再涨工资！
- ABAP ALV上的红绿灯状态
- 社交之战，结局：被封？，这只是一个开始
- 基于RSA+AES 软件授权License
- P3332 [ZJOI2013]K大数查询 - 整体二分-区间修改
- 数据挖掘的过程有哪些
- android源代码下载
- excel工具栏隐藏了怎么办_Office2016基础教程第一章:认识Excel #excel #职场 #办公技巧...