1 遞推關(guān)系-酵母菌生長(zhǎng)模型

差分方程建模的關(guān)鍵在于如何得到第n組數(shù)據(jù)與第n+1組數(shù)據(jù)之間的關(guān)系

如圖所示我們用培養(yǎng)基培養(yǎng)細(xì)菌時(shí)，其數(shù)量變化通常會(huì)經(jīng)歷這四個(gè)時(shí)期。

這個(gè)模型針對(duì)前三個(gè)時(shí)期建一個(gè)大致的模型：

調(diào)整期
對(duì)數(shù)期
穩(wěn)定期

根據(jù)已有的數(shù)據(jù)進(jìn)行繪圖：

import matplotlib.pyplot as plt 
time = [i for i in range(0,19)] 
number = [9.6,18.3,29,47.2,71.1,119.1,174.6,257.3, 
          350.7,441.0,513.3,559.7,594.8,629.4,640.8,  
          651.1,655.9,659.6,661.8] 
plt.title('Relationship between time and number')#創(chuàng)建標(biāo)題 
plt.xlabel('time')#X軸標(biāo)簽 
plt.ylabel('number')#Y軸標(biāo)簽 
plt.plot(time,number)#畫圖 
plt.show()#顯示

分析：

酵母菌數(shù)量增長(zhǎng)有一個(gè)這樣的規(guī)律:當(dāng)某些資源只能支撐某個(gè)最大限度的種群數(shù)量，而不能支持種群數(shù)量的無(wú)限增長(zhǎng)，當(dāng)接近這個(gè)最大值時(shí)，種群數(shù)量的增長(zhǎng)速度就會(huì)慢下來(lái)。

兩個(gè)觀測(cè)點(diǎn)的值差△p來(lái)表征增長(zhǎng)速度
△p與目前的種群數(shù)量有關(guān)，數(shù)量越大，增長(zhǎng)速度越快
△p還與剩余的未分配的資源量有關(guān)，資源越多，增長(zhǎng)速度越快
然后以極限總?cè)簲?shù)量與現(xiàn)有種群數(shù)量的差值表征剩余資源量

模型：

「接下來(lái)計(jì)算模型表達(dá)式求系數(shù)k（已知的情況下）」當(dāng)把該式子看成「二次曲線」進(jìn)行擬合時(shí)：

import numpy as np
import matplotlib.pylab as plt
p_n = [9.6,18.3,29,47.2, 71.1,119.1, 174.6, 
      257.3, 350.7, 441.0, 513.3, 559.7, 594.8, 629.4,
      640.8, 651.1, 655.9, 659.6]
delta_p = [8.7, 10.7,18.2,23.9, 48,55.5,
          82.7, 93.4, 90.3, 72.3, 46.4,35.1,
          34.6, 11.4, 10.3,4.8,3.7,2.2]
plt.plot(p_n,delta_p)


poly = np.polyfit(p_n, delta_p, 2)
z = np.polyval(poly,p_n)
print(poly)

plt.plot(p_n, z)
plt.show()
# [-8.01975671e-04  5.16054679e-01  6.41123361e+00]
# k = -8.01975671e-04

運(yùn)行結(jié)果

當(dāng)看成「一次曲線」進(jìn)行擬合時(shí)

?
將 (665 - pn) * pn 看出一個(gè)整體那么整個(gè)式子就是一個(gè)一次函數(shù)
?

import numpy as np
import matplotlib.pylab as plt
p_n = [9.6,18.3,29,47.2, 71.1,119.1, 174.6, 
      257.3, 350.7, 441.0, 513.3, 559.7, 594.8, 629.4,
      640.8, 651.1, 655.9, 659.6]
delta_p = [8.7, 10.7,18.2,23.9, 48,55.5,
          82.7, 93.4, 90.3, 72.3, 46.4,35.1,
          34.6, 11.4, 10.3,4.8,3.7,2.2]

p_n = np.array(p_n)
x= (665 - p_n) * p_n
plt.plot(x,delta_p)

ploy = np.polyfit(x,delta_p,1)
print(ploy)
z = np.polyval(ploy,x)

plt.plot(x,z)
plt.show()

# [ 0.00081448 -0.30791574] 
# k = 0.00081448

運(yùn)行結(jié)果

預(yù)測(cè)曲線：

import matplotlib.pyplot as plt 
p0 = 9.6 
p_list = [] 
for i in range(20): 
    p_list.append(p0) 
    p0 = 0.00081448*(665-p0)*p0+p0 
plt.plot(p_list) 
plt.show()

運(yùn)行結(jié)果：

2 顯示差分-熱傳導(dǎo)方程

一維熱傳導(dǎo)方程為：

?
其中，k為熱傳導(dǎo)系數(shù)，「第2式是方程的初值條件」，第3、4式是邊值條件,熱傳導(dǎo)方程如下：
?

繪制初值條件函數(shù)圖像（第二個(gè)式子）

import numpy as np
import matplotlib.pylab as plt

x = np.linspace(0,1,100)
y = 4 * x * (1 - x)

plt.plot(x,y)
plt.show()

運(yùn)行結(jié)果

注：CAL庫(kù)沒有安裝上，以下代碼未運(yùn)行，照著PPT寫了一遍（安裝了一天裝上了但是運(yùn)行一直報(bào)錯(cuò)）

N = 25 
M = 2500 
T = 1.0 
X = 1.0 
xArray = np.linspace(0,1.0,50) 
yArray = map(initialCondition, xArray) 
starValues = yArray 
U = np.zeros((N+1,M+1)) 
U[:,0]=starValues 
dx=X/N 
dt=T/N 
kappa=1.0
rho=kappa*dt/dx/dx 
for k in range(0,N): 
    for j in range(1,N): 
        U[j][k+1]=rho*U[j-1][k]+ (1.-2*rho)*U[j][k]+ rho*U[j+1][k] 
    U[0][k+1]=0. 
    U[N][k+1]=0. 
pylab.figure(figsize=(12,6)) 
pylab.plot(xArray, U[:,0]) 
pylab.plot(xArray, U[:,int(0.10/dt)]) 
pylab.plot(xArray, U[:,int(0.20/dt)]) 
pylab.plot(xArray, U[:,int(0.50/dt)]) 
pylab.xlabel(‘$x$’, fontsize=15) 
pylab.ylabel(r‘$U(\dot,\tau)$’,fontsize=15) 
pylab.title(u’一維熱傳導(dǎo)方程’，fontproperties=font) 
pylab.legend([r’$\tau=0.$’, r’$\tau=0.10$’, r’$\tau=0.20$’, r’$\tau=0.50$’],fontsize=15)

三維立體圖查看整體熱傳導(dǎo)過(guò)程：

tArray = np.linspace(0, 0.2, int(0.2/dt)+1)
xGride, tGride = np.meshgrid(xArray, tArray) 
from mpl_toolkits.mplot3d import Axes3D 
from matplotlib import cm 
fig = pylab.figure(figsize=(16,10)) 
ax = fig.add_subplot(1,1,1,projection=‘3d’) 
surface = ax.plot_surface(xGride, tGride, U[:,:int(0.2/dt)+1].T, cmap=cm.coolwarm)
ax.set_xlabel(“$x$”, fontdict={“size”:18}) 
ax.set_ylabel(r“$\tau$”, fontdict={“size”:18}) 
ax.set_zlabel(r”$U$”, fontdict={“size”:18}) 
ax.set_title(u”熱傳導(dǎo)方程 $u_\\tau = u_{xx}$”, fontproperties=font) 
fig.colorbar(surface, shrink=0.75)

3 馬爾科夫鏈-選舉投票預(yù)測(cè)

馬爾科夫鏈?zhǔn)怯删哂幸韵滦再|(zhì)的一系列事件構(gòu)成的過(guò)程：

一個(gè)事件有有限多個(gè)結(jié)果，稱為狀態(tài)，該過(guò)程總是這些狀態(tài)中的一個(gè)；
在過(guò)程的每個(gè)階段或者時(shí)段，一個(gè)特定的結(jié)果可以從它現(xiàn)在的狀態(tài)轉(zhuǎn)移到任何狀態(tài)，或者保持原狀；
每個(gè)階段從一個(gè)狀態(tài)轉(zhuǎn)移到其他狀態(tài)的概率用一個(gè)轉(zhuǎn)移矩陣表示，矩陣每行的各元素在0到1之間，每行的和為1。

實(shí)例：選舉投票預(yù)測(cè) 以美國(guó)大選為例，首先取得過(guò)去十次選舉的歷史數(shù)據(jù)，然后根據(jù)歷史數(shù)據(jù)得到選民意向的轉(zhuǎn)移矩陣。

構(gòu)建差分方程：

于是可以通過(guò)求解差分方程組，推測(cè)出選民投票意向的長(zhǎng)期趨勢(shì)

import matplotlib.pyplot as plt 
RLIST = [0.33333] 
DLIST = [0.33333] 
ILIST = [0.33333] 
for i in range(40): 
    R = RLIST[i]*0.75+DLIST[i]*0.20+ILIST[i]*0.40 
    RLIST.append(R) 
    D = RLIST[i]*0.05+DLIST[i]*0.60+ILIST[i]*0.20 
    DLIST.append(D) 
    I = RLIST[i]*0.20+DLIST[i]*0.20+ILIST[i]*0.40 
    ILIST.append(I) 
plt.plot(RLIST) 
plt.plot(DLIST) 
plt.plot(ILIST) 
plt.xlabel('Time') 
plt.ylabel('Voting percent')
plt.annotate('DemocraticParty',xy = (5,0.2)) 
plt.annotate('RepublicanParty',xy = (5,0.5)) 
plt.annotate('IndependentCandidate',xy = (5,0.25)) 
plt.show() 
print(RLIST,DLIST,ILIST)

運(yùn)行結(jié)果

最后得到的長(zhǎng)期趨勢(shì)是：

56%的人選共和黨
19%的人選民主黨
25%的人選獨(dú)立候選人

這個(gè)問(wèn)題還可以用「C-K方程」來(lái)解

import numpy as np
a = np.array([[0.75,0.05,0.20],[0.20,0.60,0.20],[0.40,0.20,0.40]])
p = np.mat(a)
for i in range(40):
    p = p*p   
print(p)

運(yùn)行結(jié)果：

結(jié)語(yǔ)

參考：

https://www.bilibili.com/video/BV12h411d7Dm
https://www.jianshu.com/p/aa477a3ebf08

學(xué)習(xí)來(lái)源：B站及其課堂PPT，對(duì)其中代碼進(jìn)行了復(fù)現(xiàn)

「文章僅作為學(xué)習(xí)筆記，記錄從0到1的一個(gè)過(guò)程」

希望對(duì)您有所幫助，如有錯(cuò)誤歡迎小伙伴指正～

Python數(shù)學(xué)建模系列（七）：差分

1 遞推關(guān)系-酵母菌生長(zhǎng)模型

2 顯示差分-熱傳導(dǎo)方程

3 馬爾科夫鏈-選舉投票預(yù)測(cè)

結(jié)語(yǔ)