【数学建模笔记 13】数学建模的差分方程模型

13. 差分方程模型

定义

设函数 $x_k=x(k)$，$k$ 取非负整数，称改变量 $x_{k+1}-x_k$ 为函数 $x_k$ 的差分，也称函数 $x_k$ 的一阶差分，记 $\Delta x_k$。

一阶差分的差分称二阶差分 ${\Delta }^2{x}_k$，即
$${\Delta }^2{x}_k=\Delta (\Delta x_k)=\Delta x_{k+1}-\Delta x_k$$

$$=(x_{k+2}-x_{k+1})-(x_{k+1}-x_k)$$

$$=x_{k+2}-2x_{k+1}+x_k.$$

类似有三阶差分、四阶差分……，记
$${\Delta }^3{x}_k=\Delta ({\Delta }^2{x}_k),{\Delta }^4{x}_k=\Delta ({\Delta }^3{x}_k),\dots$$
含有未知函数 $x_k$​ 的差分的方程称差分方程，差分方程中所含未知函数差分的最高阶数称差分方程的阶，如 $n$ 阶差分方程形式为
$$F(k,x_k,\Delta x_k,{\Delta }^2{x}_k,\dots ,{\Delta }^n{x}_k)=0,$$
或
$$G(k,x_k,x_{k+1},x_{k+2},\dots ,x_{k+n})=0.$$
若差分方程中所含未知函数和未知函数的各阶差分均为一次的，则称线性差分方程，形式为
$$x_{k+n}+a_1{x}_{k+n-1}+a_2{x}_{k+n-2}+\cdots +a_n{x}_k=f(k),$$
对应齐次方程
$$x_{k+n}+a_1{x}_{k+n-1}+a_2{x}_{k+n-2}+\cdots +a_n{x}_k=0.$$

常系数线性齐次差分方程的解法

对于 $n$ 阶常系数线性齐次差分方程，特征方程为
$${\lambda }^n+a_1{\lambda }^{n-1}+\cdots +a_n=0.$$
若特征方程有 $n$ 个互不相同的实根 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$，则通解为
$$c_1{\lambda }_1^k+c_2{\lambda }_2^k+\cdots +c_n{\lambda }_n^k.$$
若 $\lambda$ 是特征方程的 $m$​ 重实根，则通解中对应的项为
$$c_1{\rho }^k\cos (\phi k)+c_2{\rho }^k\sin (\phi k),$$

$$\rho =\sqrt{{\alpha }^2+{\beta }^2},\phi =\arctan \frac{\beta }{\alpha }.$$

若$\lambda =\alpha +\beta i$ 是特征方程的 $m$ 重复根，则通解中对应的项为
$$(c_1+c_{2}k+\cdots +c_m{k}^{m-1}){\rho }^k\cos (\phi k)$$

$$+(c_{m+1}+c_{m+2}k+\cdots +c_{2m}{k}^{m-1}){\rho }^k\sin (\phi k).$$

例子

求斐波那契数列
$$\{\begin{array}{rl}& F_n=F_{n-1}+F_{n-2},\\ & F_1=F_2=1.\end{array}$$
的解。

有特征方程
$${\lambda }^2-\lambda -1=0.$$
解得
$${\lambda }_1=\frac{1-\sqrt{5}}{2},{\lambda }_2=\frac{1+\sqrt{5}}{2}$$
互异，通解为
$$F_n=c_1(\frac{1-\sqrt{5}}{2}{)}^n+c_2(\frac{1+\sqrt{5}}{2}{)}^n.$$
又 $F_1=F_2=1$，解得 $c_1=-\frac{\sqrt{5}}{5},c_2=\frac{\sqrt{5}}{5}$，代入得
$$F_n=\frac{\sqrt{5}}{5}[(\frac{1+\sqrt{5}}{2}{)}^n-c_2(\frac{1-\sqrt{5}}{2}{)}^n].$$

常系数线性非齐次差分方程的解法

若 $x_k$ 为齐次方程的通解，$x_k^{\ast }$ 为非齐次方程的特解，则非齐次方程的通解为 $x_k+x_k^{\ast }$。

求特解一般用常数变易法，特殊情况也可采用待定系数法。

典型模型

Leslie 模型

考虑对人口进行分组，为简单起见只考虑女性人口，将女性人口按年龄划分称 $m$ 个年龄组。将时间离散为时段 $t_k,k=1,2,\dots$

记时段 $t_k$​ 第 $i$​ 年龄组的种群数量为 $x_i(k),i=1,2,\dots ,m$​，繁殖率为 ${\alpha }_i$，死亡率为 $d_i$，${\beta }_i=1-d_i$ 称存活率。

第一年龄组种群数量是时段 $t_k$ 各年龄组繁殖数量之和，即
$$x_1(k+1)=\sum _{i=1}^m{\alpha }_i{x}_i(k),$$
$t_{k+1}$ 时段第 $i+1$ 年龄组的种群数量是 $t_k$ 时段第 $i$ 年龄组存活下来的数量，即
$$x_{i+1}(k+1)={\beta }_i{x}_i(k),i=1,2,\dots ,m-1.$$
记 $t_k$ 时段种群各年龄组的分布向量

$$X(k)=\left(\begin{array}{c}{x}_1(k)\\ x_2(k)\\ ⋮\\ x_m(k)\end{array}\right),$$
记
$$L=\left(\begin{array}{ccccc}{\alpha }_1& {\alpha }_2& \dots & {\alpha }_{m-1}& {\alpha }_m\\ {\beta }_1& 0& \dots & 0& 0\\ 0& {\beta }_2& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & {\beta }_{m-1}& 0\end{array}\right),$$
则
$$X(k+1)=LX(k),k=0,1,\dots$$
当 $t_0$​​ 时段各年龄组人数 $X(0)$​ 已知，就可以求得 $t_k$​ 时段按年龄组的分布向量
$$X(k)=L^{k}X(0),k=0,1,\dots$$
由此算出种群总量。

离散阻滞增长模型

又称离散 Logistic 模型。
$$\Delta x_k=r{x}_k(1-\frac{x_k}{N}),k=0,1,2,\dots$$
其中 $r$​ 为固有增长率，$N$​ 为最大容量。上式可化为
$$\frac{x_{k+1}-x_k}{x_k}=r(1-\frac{x_k}{N}),k=0,1,2,\dots$$
当 $x=0$​ 时，增长率为 $r$​，随着 $x$ 增加，增长率逐渐减小，直到 $x=N$ 时，增长率为 0。

染色体遗传模型

如果某遗传特征由两个基因 A 和 a 控制，就有三种基因对 AA, Aa, aa。不同基因对结合，后代从父体、母体内分别取得一个基因构成新的基因对。

如父体为 AA，母体为 Aa，则后代从父体必获得 A，从母体获得 A 或 a 的概率分别为 1/2，构成新的基因对 AA 或 Aa 的概率分别为 1/2。

父体-母体/后代	AA*AA	AA*Aa	AA*aa	Aa*Aa	Aa*aa	aa*aa
AA	1	1/2	0	1/4	0	0
Aa	0	1/2	1	1/2	1/2	0
aa	0	0	0	1/4	1/2	1

令 $a_n,b_n,c_n$ 表示第 $n$ 代植物基因型为 AA, Aa, aa 的植物占植物总数的百分率，令
$$x^{(0)}=(a_0,b_0,c_0{)}^T$$
表示初始分布。

选用 AA 型植物与其他植物结合。由于 AA 与 AA 结合必为 AA，AA 与 Aa 结合为 AA 的可能性为 1/2，AA 与 aa 结合不可能为 AA，得
$$a_n=1a_{n-1}+\frac{1}{2}{b}_{n-1}+0c_{n-1},$$
同理有
$$b_n=\frac{1}{2}{b}_{n-1}+c_{n-1},$$

$$c_n=0.$$

记
$$M=\left(\begin{array}{ccc}1& 1/2& 0\\ 0& 1/2& 1\\ 0& 0& 0\end{array}\right)$$
则有
$$x^{(n)}=M{x}^{(n-1)},n=1,2,\dots$$
即
$$x^{(n)}=M^n{x}^{(0)}.$$
在这种情况下，$a_n\to 1,b_n\to 0,c_n\to 0$。

若选用相同基因型植物相结合，则
$$M=\left(\begin{array}{ccc}1& 1/4& 0\\ 0& 1/2& 0\\ 0& 1/4& 1\end{array}\right)$$
此时
$$a_n\to a_0+\frac{1}{2}{b}_0,b_n\to 0,c_n\to c_0+\frac{1}{2}{b}_0.$$

Python 代码

解差分方程

对常系数线性齐次差分方程的解法的例子求解，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-28
# @ function: 使用 sympy 求解差分方程

# %%

import numpy as np
import sympy as sp

# %%

# 定义方程
x = sp.symbols('x')
y = sp.Function('y')
f = y(x+2) - y(x+1) - y(x)

con = {
	y(1):1,
	y(2):1,
}

# 解递归方程
solve = sp.rsolve(f, y(x), con)
solve

# %%

# 画图
x1 = np.linspace(1,10,10)
y1 = []
for each in x1:
	y1.append((solve.subs(x,each).evalf()))
y1

import matplotlib.pyplot as plt

plt.plot(x1,y1)
plt.scatter(x1,y1)

输出如下：
$$-\frac{\sqrt{5}(\frac{1}{2}-\frac{\sqrt{5}}{2}{)}^x}{5}+-\frac{\sqrt{5}(\frac{1}{2}+\frac{\sqrt{5}}{2}{)}^x}{5}$$

给出通解和部分数值解。

Leslie 模型

养殖场养殖一类动物最多三年，按一岁、两岁、三岁划分为 3 个年龄组。一岁不繁殖，两岁平均一年繁殖 4 个后代，三岁平均一年繁殖 3 个后代。一岁、两岁的成活率分别为 0.5、0.25。假设开始时各年龄组动物为 1000 头，求未来动物的数量。

由题意，出生率为 $\alpha =(0,4,3{)}^T$，成活率为 $\beta =(0.5,0.25{)}^T$，则
$$L=\left(\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\right)$$
有
$$X(k+1)=LX(k),k=0,1,2,\dots$$
又
$$X(0)=(1000,1000,1000{)}^T$$
可以求解，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-28
# @ function: Lesile 模型

# %%

import numpy as np

# %%

# 源数据
alpha = np.array([0, 4, 3])
beta = np.array([.5, .25])

L = np.zeros((len(alpha), len(alpha)))
L[0, :] = alpha
for index, each in enumerate(beta):
    L[index+1, index] = each

x0 = np.array([1000, 1000, 1000])

# 迭代求解
x_list = [x0]
for i in range(5):
    x_list.append(L.dot(x_list[-1].T))
for index, x in enumerate(x_list):
    print('the {}th year: 1y: {}, 2y: {}, 3y: {}'.format(
        index, x[0], x[1], x[2]))

# %%

import matplotlib.pyplot as plt

year_list = np.array(range(len(x_list)))
y1,y2,y3 = [],[],[]
for year in year_list:
	y1.append(x_list[year][0])
	y2.append(x_list[year][1])
	y3.append(x_list[year][2])

w = .3
plt.bar(year_list-w,y1, width=w, label='1y')
plt.bar(year_list,y2, width=w, label='2y')
plt.bar(year_list+w,y3, width=w, label='3y')
plt.legend()

输出如下：

the 0th year: 1y: 1000, 2y: 1000, 3y: 1000
the 1th year: 1y: 7000.0, 2y: 500.0, 3y: 250.0
the 2th year: 1y: 2750.0, 2y: 3500.0, 3y: 125.0
the 3th year: 1y: 14375.0, 2y: 1375.0, 3y: 875.0
the 4th year: 1y: 8125.0, 2y: 7187.5, 3y: 343.75
the 5th year: 1y: 29781.25, 2y: 4062.5, 3y: 1796.875

染色体遗传模型

对例子中的两中情况求解，代码如下：

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-28
# @ function: 染色体遗传模型

# %%

import numpy as np

# %%

# 都选用 AA 结合    
M1 = np.array([[1, .5, 0],
               [0, .5, 1],
               [0, 0, 0]])

# 选用相同基因组结合
M2 = np.array([[1, .25, 0],
               [0, .5, 0],
               [0, .25, 1]])

x0 = np.array([1/3, 1/3, 1/3])

# 迭代求解
x1_list, x2_list = [x0], [x0]
for i in range(10):
	x1_list.append(M1.dot(x1_list[-1].T))
	x2_list.append(M2.dot(x2_list[-1].T))

print('----- M1 -----')
for index, each in enumerate(x1_list):
	print('iter {}: AA: {:.4f}, Aa: {:.4f}, aa: {:.4f}'.format(index, each[0], each[1], each[2]))

print('\n----- M2 -----')
for index, each in enumerate(x2_list):
	print('iter {}: AA: {:.4f}, Aa: {:.4f}, aa: {:.4f}'.format(index, each[0], each[1], each[2]))


# %%

# 情况一画图
a,b,c = [], [], []
for each in x1_list:
	a.append(each[0])
	b.append(each[1])
	c.append(each[2])

import matplotlib.pyplot as plt

it = [i for i in range(11)]
plt.plot(it,a)
plt.plot(it,b)
plt.plot(it,c)
plt.scatter(it,a, label='AA')
plt.scatter(it,b, label='Aa')
plt.scatter(it,c, label='aa')
plt.legend()

# %%

# 情况二画图
a,b,c = [], [], []
for each in x2_list:
	a.append(each[0])
	b.append(each[1])
	c.append(each[2])

import matplotlib.pyplot as plt

it = [i for i in range(11)]
plt.plot(it,a, alpha=.5)
plt.plot(it,b)
plt.plot(it,c, alpha=.5)
plt.scatter(it,a, label='AA', alpha=.5)
plt.scatter(it,b, label='Aa')
plt.scatter(it,c, label='aa', alpha=.5)
plt.legend()

输出如下：

----- M1 -----
iter 0: AA: 0.3333, Aa: 0.3333, aa: 0.3333
iter 1: AA: 0.5000, Aa: 0.5000, aa: 0.0000
iter 2: AA: 0.7500, Aa: 0.2500, aa: 0.0000
iter 3: AA: 0.8750, Aa: 0.1250, aa: 0.0000
iter 4: AA: 0.9375, Aa: 0.0625, aa: 0.0000
iter 5: AA: 0.9688, Aa: 0.0312, aa: 0.0000
iter 6: AA: 0.9844, Aa: 0.0156, aa: 0.0000
iter 7: AA: 0.9922, Aa: 0.0078, aa: 0.0000
iter 8: AA: 0.9961, Aa: 0.0039, aa: 0.0000
iter 9: AA: 0.9980, Aa: 0.0020, aa: 0.0000
iter 10: AA: 0.9990, Aa: 0.0010, aa: 0.0000

----- M2 -----
iter 0: AA: 0.3333, Aa: 0.3333, aa: 0.3333
iter 1: AA: 0.4167, Aa: 0.1667, aa: 0.4167
iter 2: AA: 0.4583, Aa: 0.0833, aa: 0.4583
iter 3: AA: 0.4792, Aa: 0.0417, aa: 0.4792
iter 4: AA: 0.4896, Aa: 0.0208, aa: 0.4896
iter 5: AA: 0.4948, Aa: 0.0104, aa: 0.4948
iter 6: AA: 0.4974, Aa: 0.0052, aa: 0.4974
iter 7: AA: 0.4987, Aa: 0.0026, aa: 0.4987
iter 8: AA: 0.4993, Aa: 0.0013, aa: 0.4993
iter 9: AA: 0.4997, Aa: 0.0007, aa: 0.4997
iter 10: AA: 0.4998, Aa: 0.0003, aa: 0.4998