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

最新推荐文章于 2022-08-26 16:53:55 发布
最新推荐文章于 2022-08-26 16:53:55 发布
阅读量836 收藏 18
点赞数 2
分类专栏： 数学建模_高教杯 文章标签： python 可视化 plot
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_45931661/article/details/119493281
版权

目录

序言：

在这里插入图片描述

微分方程解析解

在这里插入图片描述

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()
一只特立独行的猫
关注
  • 2
    点赞
  • 踩
  • 18
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
fangchengqiugeng.rar_偏微分方程_偏微分方程 数值解
07-14
偏微分方程的数值解法，包括所有源程序，请多多指教！
微分方程解析解+数值解
qq_51070956的博客
07-02 2258
动手练
Matlab与微分方程解析解（dsolve）
醉后不知天在水，满船清梦压星河
07-25 3万+
微分方程是动态系统建模的基础。本文介绍了微分方程的解析解法。 常系数线性微分方程：
疫情传播SEIR模型（python）
jkchen's Haven
05-12 1万+
假设有N个人，两两之间有接触的关系为p。 事件开始： 选中一个人变为潜伏者EEE 时间推移： 健康者SSS：若与xxx个感染者III接触，则有1−(1−w)x1-(1-w)^x1−(1−w)x的概率变为潜伏者EEE 潜伏者EEE：有β\betaβ的概率病发变为感染者III，有μ\muμ的概率变为终态RRR 感染者III：有μ\muμ的概率变为终态RRR 参数设置： N=1000,p=0.006w=0.2,β=0.5,μ=0.2N = 1000,p = 0.006\\ w=0.2,\beta=0.5,\
SI，SIS，SIR，SEIRD模型
qq_37730871的博客
08-26 7334
因为个人工作需要系统地整理SI，SIR以及SEIR模型，故对三个模型进行原理介绍以及对比。文中关于SI，SIS，SIR的所有的截图都来自西工大肖华勇老师在慕课上的分享，原视频戳SEIRD模型则来自发表在SCI上的paper，想看原文戳。...
【Matlab代码】微分方程的解析解和数值解
qq_41985559的博客
09-25 2万+
Matlab求微分方程的解析解 用法：dsolve(‘equation’,‘condition;,‘v’) 说明：（1）equation是方程式，condition是条件，v是自变量（确省为t） （2）若不带条件，则解中带积分常数。 （3）如果没有显示解，则系统尝试给出隐氏解。 （4）如果无隐氏解，则返回空符号。 格式：（1）y’表示为Dy，y’‘表示为D2y，依次类推。 （2）有多个方程或多个条件时，写多个对应的参数即可。 求微分方程y’’-5y’+6y=e^ax y=dsolve('D2y-5*Dy+6
partial-differential-equations.rar_matlab数值求解_偏微分方程_数值解 matlab
07-15
可供用matlab求解偏微分方程的数值解学习
数值方法求解常微分方程_数值方法解常微分方程_
10-01
数值方法求解常微分方程的实验例题，由matlab语言实现
DYa.rar_偏微分方程_偏微分方程 数值解_隐格式
07-14
介绍了DYa算法的具体应用,是解决偏微分方程数值解的交换隐格式的方法
偏微分方程(PDEs)的MATLAB数值解法.rar_matlab 解偏微分方程 PDF_偏微分方程_数值解法
07-14
matlab 解偏微分方程 数值解法 pdf
微分方程的求解方法
热门推荐
Spuer_Tiger
12-13 14万+
文章目录前言Ⅰ.首先介绍一些关于微分方程的概念Ⅱ.在考研范围内的微分方程有哪几类Ⅲ.微分方程的求解方法1.一阶微分方程的求解①可分离变量型的解法②齐次型的解法③一阶线性型的解法（重难点）2.二阶可降阶微分方程的求解3.高阶常系数线性微分方程的求解 前言 本文主要介绍了考研范围的微分方程的求解类型及对应的求解方法，主要内容参考自张宇《闭关修炼》，希望本文对您有所帮助。 Ⅰ.首先介绍一些关于微分方程的概念 一阶是什么： 一阶微分方程就是指只有一阶导数或微分的微分方程。（注：阶数是微分方程中含有的导数或微分的
微分方程概述及matlab代码（求解析解）
小风啊吼吼吼的博客
11-05 4690
微分方程概述及matlab代码（求解析解） 文章目录微分方程概述及matlab代码（求解析解）微分方程概述引例：导弹射击问题微分方程基本概念Matlab求微分方程的解析解关于导弹射击问题的matlab求解关于Matlab求解析解的几点注意 微分方程概述 当我们描述实际对象的某些特性随时间（空间）而演变的过程，分析它的变化规律，预测它的未来形态、研究它的控制手段时。通常要建立对象的动态模型。 在许多实际问题中，当直接导出变量之间的函数关系较为困难，但导出包含未知函数的导数或微分的关系式比较容易时，可以用建立微
matlab求解微分方程解析解
qq_45341830的博客
08-09 1万+
解析解 1求解微分方程方程组用函数“dsolve” dsolve(‘方程1’,‘方程2’…‘方程n’,‘初始条件’,‘自变量’) clc,clear %% %求解du/dt=i+u^2的通解 dsolve('Du=1+u^2','t') %% %求解微分方程组 dy^2/dx^2+4*dy/dx+29*y=0,y(0)=0,y`(0)=15;的通解 y=dsolve('D2y+4*Dy+29*...
Python解微分方程
linshen1213的博客
08-19 4万+
Python解微分方程微分方程回顾微分方程：python 解析解（SymPy）微分方程：python数值解(SciPY)微分方程组：python数值解 微分方程回顾 微分方程是用来描述某一类函数与其导数之间的关系的方程，其解是一个符合方程的函数。（解为常数值的是代数方程） 微分方程的求解是研究微分方程的重点之一，例如解是否存在，存在是否唯一等等。只有少数类型的微分方程存在解析解；无法求得解析解时，...
Matlab求微分方程的解析解（dsolve）
qq_36890370的博客
08-10 1万+
一、Matlab求解函数如下 dsolve('方程1','方程2',...,'方程n','初始条件','自变量') 注意 1）在表达微分方程时，用字母D表示微分，D2、D3分别表示二阶、三阶微分，后面跟的是求解的因变量； 2）自变量可以指定，不写时默认为t； 3）用单引号； 4）在单引号中的公式算数符号不能省略：比如‘y-Dy=2*x’，这里的乘号不能省略； 5）初始值可以不给，不给求出来就是通解 6）如果微分方程形式较为复杂...
Matlab求常微分方程组的解析解
qq_18820125的博客
03-07 1万+
Matlab求常微分方程的解析解 最近同学毕设需要求解循坏摆的微分方程，我在帮忙过程中学习了一下常微分方程的解析解和数值解的求法，在此分享。 以下讲解遵循Matlab官方文档提供的方程和写法。 一阶常微分方程求解（无初值） 方程：dy/dt=ay 代码： syms y(t) a; %syms 定义方程用的变量和未知常量(可有可无)。y(t)表示y是关于t的函数，求解t。 eqn=diff(y...
一个微分方程的解析解
weixin_30838921的博客
11-08 339
近日我在给工程硕士讲《高等工程数学》，在讲授微分方程数值解这一节课的时候。 遇到一个一阶非线性微分方程初值问题，如下： \[\begin{cases} \dfrac{\text{d}y}{\text{d}x} = y-\dfrac{2x}{y}\\ y(0)=1\\ \end{cases} \] 解析解为 \(y = \sqrt{1+2x}\)，学生问怎么求出来的？ 现解答如下： 两边同时乘以 \...
常微分方程的解析解(方法归纳)以及基于Python的二阶微分方程边值问题的数值算例实现
图灵的猫的博客
10-06 3853
常微分方程的解析解(方法归纳)以及基于Python的微分方程数值解算例实现 本文归纳常见常微分方程的解析解解法以及基于Python的微分方程数值解算例实现。 如果没有出现意外，本文将不包含解法的推导过程。 常微分方程的解析解 考虑常微分方程的解析解法，我们一般可以将其归纳为如下几类： 可分离变量的微分方程（一阶） 一阶齐次（非齐次）线性微分方程（一阶） 二阶常系数微分方程（二阶） 高阶常系数微分方程（nnn阶） 可分离变量的微分方程（一阶） 这类微分方程可以变形成如下形式： f(x)dx=g(y)dy
微分方程解法总结
linkequa的博客
06-19 8万+
1，微分方程的概念 微分方程是一种描述函数与其导数关系的数学方程，它的解通常是函数，而初等代数中方程的解通常为数值。 2，分类 2.1 常微分方程与偏微分方程 常微分方程的未知数是单一变量的函数。表达通式为： f(x,y,y′,y′′,...,y(n))=0,f(x, y, y', y'', ..., y^{(n)})=0,f(x,y,y′,y...
python如何求解微分方程_常微分方程数值解：Python求解
最新发布
06-02
Python中可以使用SciPy库中的odeint函数来求解常微分方程的数值解。 首先，需要导入SciPy库和NumPy库： ```python import numpy as np from scipy.integrate import odeint ``` 然后，需要定义一个函数来描述微分...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值