python數(shù)學(xué)建模（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()    #初始化打印環(huán)境
t, omega0, gamma = sympy.symbols("t, omega_0, gamma", positive = True)    #標(biāo)記參數(shù)掺冠，且均為正
x = sympy.Function('x')    #標(biāo)記x是微分函數(shù)，非變量
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)}$

微分方程數(shù)值解

場(chǎng)線圖與數(shù)值解

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)     #符號(hào)表達(dá)式轉(zhuǎn)隱函數(shù)
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的場(chǎng)線圖
ax.plot(xn, yn, 'b', lw = 2)
ax.plot(xp, yp, 'r', lw = 2)
plt.show()

洛倫茲曲線與數(shù)值解

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

洛倫茲曲線很好地展示了非線性的直觀形態(tài)泳姐，同時(shí)也是展示混沌系統(tǒng)的經(jīng)典例子形娇。這個(gè)例子告訴我們息堂，混沌系統(tǒng)可以是一個(gè)確定系統(tǒng)，不一定是隨機(jī)過(guò)程，同時(shí)它有初值敏感的特征

傳染病模型

SI-Model（只感染不恢復(fù)）

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N為人群總數(shù)
beta = 0.25 # β為傳染率系數(shù)
gamma = 0 # gamma為恢復(fù)率系數(shù)
I_0 = 1 # I_0為感染者的初始人數(shù)
S_0 = N- I_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, I_0) # INI為初始狀態(tài)下的數(shù)組

def funcSI(inivalue,_):
    Y = np.zeros(2)
    X = inivalue
    Y[0] = -(beta*X[0]*X[1])/N+gamma*X[1]#易感個(gè)體變化
    Y[1] = (beta*X[0]*X[1])/N-gamma*X[1]#感染個(gè)體變化
    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（會(huì)恢復(fù)）

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N為人群總數(shù)
beta = 0.25 # β為傳染率系數(shù)
gamma = 0.05 # gamma為恢復(fù)率系數(shù)
I_0 = 1 # I_0為感染者的初始人數(shù)
S_0 = N- I_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, I_0) # INI為初始狀態(tài)下的數(shù)組

def funcSIS(inivalue,_):
    Y = np.zeros(2)
    X = inivalue
    Y[0] = -(beta*X[0])/N*X[1]+gamma*X[1]#易感個(gè)體變化
    Y[1] = (beta*X[0]*X[1])/N-gamma*X[1]#感染個(gè)體變化
    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（研發(fā)出疫苗读存，可以治療）

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N為人群總數(shù)
beta = 0.25 # β為傳染率系數(shù)
gamma = 0.05 # gamma為恢復(fù)率系數(shù)
I_0 = 1 # I_0為感染者的初始人數(shù)
R_0 = 0 # R_0為治愈者的初始人數(shù)
S_0 = N - I_0 - R_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, I_0, R_0) # INI為初始狀態(tài)下的數(shù)組

def funcSIR(inivalue,_):
    Y = np.zeros(3)
    X = inivalue
    Y[0] = -(beta*X[0] *X[1])/N #易感個(gè)體變化
    Y[1] = (beta*X[0]*X[1])/N-gamma*X[1]#感染個(gè)體變化
    Y[2] = gamma*X[1] #治愈個(gè)體變化
    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（抗體有時(shí)效）

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N為人群總數(shù)
beta = 0.25 # β為傳染率系數(shù)
gamma = 0.05 # gamma為恢復(fù)率系數(shù)
Ts = 7 # Ts為抗體持續(xù)時(shí)間
I_0 = 1 # I_0為感染者的初始人數(shù)
R_0 = 0 # R_0為治愈者的初始人數(shù)
S_0 = N - I_0 - R_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, I_0, R_0) # INI為初始狀態(tài)下的數(shù)組

def funcSIRS(inivalue,_):
    Y = np.zeros(3)
    X = inivalue
    Y[0] = -(beta*X[0]*X[1])/N+X[2]/Ts #易感個(gè)體變化
    Y[1] = (beta*X[0]*X[1])/N-gamma*X[1]#感染個(gè)體變化
    Y[2] = gamma*X[1]-X[2]/Ts #治愈個(gè)體變化
    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為人群總數(shù)
beta = 0.6 # β為傳染率系數(shù)
gamma = 0.1 # gamma為恢復(fù)率系數(shù)
Te = 14 # Te為疾病潛伏期
I_0 = 1 # I_0為感染者的初始人數(shù)
E_0 = 0 # E_0為潛伏者的初始人數(shù)
R_0 = 0 # R_0為治愈者的初始人數(shù)
S_0 = N - I_0 - R_0 - E_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, E_0, I_0, R_0) # INI為初始狀態(tài)下的數(shù)組
def funcSIER(inivalue,_):
    Y = np.zeros(4)
    X = inivalue
    Y[0] = -(beta*X[0]*X[2])/N #易感個(gè)體變化
    Y[1] = (beta*X[0]*X[2]/N-X[1]/Te) # 潛伏個(gè)體變化
    Y[2] = X[1]/Te-gamma*X[2]#感染個(gè)體變化
    Y[3] = gamma*X[2] #治愈個(gè)體變化
    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（考慮抗體時(shí)效）

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
N = 10000 # N為人群總數(shù)
beta = 0.6 # β為傳染率系數(shù)
gamma = 0.1 # gamma為恢復(fù)率系數(shù)
Ts = 7 # Ts為抗體持續(xù)時(shí)間
Te = 14 # Te為疾病潛伏期
I_0 = 1 # I_0為感染者的初始人數(shù)
E_0 = 0 # E_0為潛伏者的初始人數(shù)
R_0 = 0 # R_0為治愈者的初始人數(shù)
S_0 = N - I_0 - R_0 - E_0 # S_0為易感染者的初始人數(shù)
T = 150 # T為傳播時(shí)間
INI = (S_0, E_0, I_0, R_0) # INI為初始狀態(tài)下的數(shù)組
def funcSIERS(inivalue,_):
    Y = np.zeros(4)
    X = inivalue
    Y[0] = -(beta*X[0]*X[2])/N+X[3]/Ts #易感個(gè)體變化
    Y[1] = (beta*X[0]*X[2]/N-X[1]/Te) # 潛伏個(gè)體變化
    Y[2] = X[1]/Te-gamma*X[2]#感染個(gè)體變化
    Y[3] = gamma*X[2]-X[3]/Ts #治愈個(gè)體變化
    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()

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者

人面猴
序言：七十年代末谨娜，一起剝皮案震驚了整個(gè)濱河市航攒，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌趴梢，老刑警劉巖漠畜，帶你破解...
沈念sama閱讀 206,839評(píng)論 6贊 482
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件币他，死亡現(xiàn)場(chǎng)離奇詭異，居然都是意外死亡盆驹，警方通過(guò)查閱死者的電腦和手機(jī)圆丹，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,543評(píng)論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)躯喇，“玉大人辫封，你說(shuō)我怎么就攤上這事×觯” “怎么了倦微？”我有些...
開封第一講書人閱讀 153,116評(píng)論 0贊 344
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長(zhǎng)正压。經(jīng)常有香客問我欣福，道長(zhǎng)，這世上最難降的妖魔是什么焦履？我笑而不...
開封第一講書人閱讀 55,371評(píng)論 1贊 279
?港島之戀（遺憾婚禮）
正文為了忘掉前任拓劝，我火速辦了婚禮，結(jié)果婚禮上嘉裤，老公的妹妹穿的比我還像新娘郑临。我一直安慰自己，他們只是感情好屑宠，可當(dāng)我...
茶點(diǎn)故事閱讀 64,384評(píng)論 5贊 374
惡毒庶女頂嫁案：這布局不是一般人想出來(lái)的
文/花漫我一把揭開白布厢洞。她就那樣靜靜地躺著，像睡著了一般典奉。火紅的嫁衣襯著肌膚如雪躺翻。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 49,111評(píng)論 1贊 285
城市分裂傳說(shuō)
那天卫玖，我揣著相機(jī)與錄音公你，去河邊找鬼。笑死假瞬，一個(gè)胖子當(dāng)著我的面吹牛省店，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播笨触，決...
沈念sama閱讀 38,416評(píng)論 3贊 400
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼懦傍，長(zhǎng)吁一口氣：“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼！你這毒婦竟也來(lái)了芦劣？” 一聲冷哼從身側(cè)響起粗俱，我...
開封第一講書人閱讀 37,053評(píng)論 0贊 259
萬(wàn)榮殺人案實(shí)錄
序言：老撾萬(wàn)榮一對(duì)情侶失蹤，失蹤者是張志新（化名）和其女友劉穎虚吟，沒想到半個(gè)月后寸认，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體签财，經(jīng)...
沈念sama閱讀 43,558評(píng)論 1贊 300
?護(hù)林員之死
正文獨(dú)居荒郊野嶺守林人離奇死亡，尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點(diǎn)故事閱讀 36,007評(píng)論 2贊 325
?白月光啟示錄
正文我和宋清朗相戀三年偏塞，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了唱蒸。大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點(diǎn)故事閱讀 38,117評(píng)論 1贊 334
活死人
序言：一個(gè)原本活蹦亂跳的男人離奇死亡灸叼，死狀恐怖神汹，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情古今，我是刑警寧澤屁魏，帶...
沈念sama閱讀 33,756評(píng)論 4贊 324
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站捉腥，受9級(jí)特大地震影響氓拼，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜抵碟，卻給世界環(huán)境...
茶點(diǎn)故事閱讀 39,324評(píng)論 3贊 307
男人毒藥：我在死后第九天來(lái)索命
文/蒙蒙一桃漾、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧拟逮，春花似錦撬统、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,315評(píng)論 0贊 19
一樁弒父案宪摧，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽(yáng)粒竖。三九已至颅崩，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間蕊苗，已是汗流浹背沿后。一陣腳步聲響...
開封第一講書人閱讀 31,539評(píng)論 1贊 262
情欲美人皮
我被黑心中介騙來(lái)泰國(guó)打工，沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留朽砰，地道東北人尖滚。一個(gè)月前我還...
沈念sama閱讀 45,578評(píng)論 2贊 355
代替公主和親
正文我出身青樓，卻偏偏與公主長(zhǎng)得像瞧柔，于是被迫代替她去往敵國(guó)和親漆弄。傳聞我的和親對(duì)象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點(diǎn)故事閱讀 42,877評(píng)論 2贊 345

python數(shù)學(xué)建模 （3）微分方程問題