python数学建模（3）微分方程问题

微分方程分类

微分方程解析解

import numpy as np
from scipy import integrate
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)

sympy.init_printing()    #初始化打印环境
t, omega0, gamma = sympy.symbols("t, omega_0, 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()和dsolve得到通解
ics = {x(0): 1, x(t).diff(t).subs(t, 0): 0}     #将初始条件字典匹配
x_t_sol = apply_ics(ode_sol, ics, t, [omega0, gamma])
print(sympy.latex(x_t_sol))     #以latex公式形式输出

$x{\left(t \right)} = \left(- \frac{\gamma}{2 \sqrt{\gamma - 1} \sqrt{\gamma + 1}} + \frac{1}{2}\right) e^{\omega_{0} t \left(- \gamma - \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)} + \left(\frac{\gamma}{2 \sqrt{\gamma - 1} \sqrt{\gamma + 1}} + \frac{1}{2}\right) e^{\omega_{0} t \left(- \gamma + \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)}$

微分方程数值解

场线图与数值解

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.subplot(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, xn)
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()

洛伦兹曲线与数值解

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

洛伦兹曲线很好地展示了非线性的直观形态，同时也是展示混沌系统的经典例子。这个例子告诉我们，混沌系统可以是一个确定系统，不一定是随机过程，同时它有初值敏感的特征

传染病模型

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 funcSIER(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('SIER Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

SIERS-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 funcSIERS(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('SIERS Model')
plt.legend()
plt.xlabel('Day')
plt.ylabel('Number')
plt.show()

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,456评论 5赞 477
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,370评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,337评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,583评论 1赞 273
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,596评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,572评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,936评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,595评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,850评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,601评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,685评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,371评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,951评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,934评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,167评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 43,636评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,411评论 2赞 342

python数学建模 （3）微分方程问题