基于python的数学建模---差分方程

一、递推关系——酵母菌生长模型

代码：

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]#19个数据 与i相对应
plt.title('Relationship between time and number')#创建标题
plt.xlabel('time')#X轴标签
plt.ylabel('number')#Y轴标签
plt.plot(time,number)#画图
plt.show()#显示

分析：

酵母菌数量增长有一个这样的规律:当某些资源只能支撑某个最大限度的种群数量，而不能支持种群数量的无限增长，当接近这个最大值时，种群数量的增长速度就会慢下来。

两个观测点的值差△p来表征增长速度
△p与目前的种群数量有关，数量越大，增长速度越快
△p还与剩余的未分配的资源量有关，资源越多，增长速度越快
然后以极限总群数量与现有种群数量的差值表征剩余资源量

$&nbsp;$

认为delta_p为二次函数时

代码：

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]

b ：把k(665-pn)看成一个整体

代码：

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]

模型：

为什么没有后面的b？一开始酵母菌需要有一定的数量

预测曲线：

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()

预测与实际曲线

import matplotlib.pyplot as plt
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]
time = [i for i in range(0,19)]
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.scatter(time,number,s=100,alpha=1.0,marker='o')
plt.show()

二、显式差分——热传导方程

其中，k为热传导系数，第2式是方程的初值条件，第3、4式是边值条件,热传导方程如下：

绘制初值条件函数图像（第二个式子）

from matplotlib import pylab
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['font.sans-serif'] = ['Microsoft YaHei']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题


def initialCondition(x):
    return 4.0 * (1.0 - x) * x

xArray = np.linspace(0, 1.0, 50)
yArray = np.array(list(map(initialCondition, xArray)))
pylab.figure(figsize=(12, 6))
pylab.xlabel('$x$', fontsize=15)
pylab.ylabel('$f(x)$', fontsize=15)
pylab.title(u'一维热传导方程初值条件')
pylab.plot(xArray, yArray)
plt.show()

三、马尔科夫链

构建差分方程组

代码：

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)