13. 差分方程模型

定义

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

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

= ( x k + 2 − x k + 1 ) − ( x k + 1 − x k ) =(x_{k+2}-x_{k+1})-(x_{k+1}-x_k) =(xk+2​−xk+1​)−(xk+1​−xk​)

= x k + 2 − 2 x k + 1 + x k . =x_{k+2}-2x_{k+1}+x_k. =xk+2​−2xk+1​+xk​.

类似有三阶差分、四阶差分……，记
Δ 3 x k = Δ ( Δ 2 x k ) , Δ 4 x k = Δ ( Δ 3 x k ) , … \Delta^3x_k=\Delta(\Delta^2x_k),\Delta^4x_k=\Delta(\Delta^3x_k),\dots Δ3xk​=Δ(Δ2xk​),Δ4xk​=Δ(Δ3xk​),…
含有未知函数 x k x_k xk​​ 的差分的方程称差分方程，差分方程中所含未知函数差分的最高阶数称差分方程的阶，如 n n n 阶差分方程形式为
F ( k , x k , Δ x k , Δ 2 x k , … , Δ n x k ) = 0 , F(k,x_k,\Delta x_k,\Delta^2x_k,\dots,\Delta^nx_k)=0, F(k,xk​,Δxk​,Δ2xk​,…,Δnxk​)=0,
或
G ( k , x k , x k + 1 , x k + 2 , … , x k + n ) = 0. G(k,x_k,x_{k+1},x_{k+2},\dots,x_{k+n})=0. G(k,xk​,xk+1​,xk+2​,…,xk+n​)=0.
若差分方程中所含未知函数和未知函数的各阶差分均为一次的，则称线性差分方程，形式为
x k + n + a 1 x k + n − 1 + a 2 x k + n − 2 + ⋯ + a n x k = f ( k ) , x_{k+n}+a_1x_{k+n-1}+a_2x_{k+n-2}+\dots+a_nx_k=f(k), xk+n​+a1​xk+n−1​+a2​xk+n−2​+⋯+an​xk​=f(k),
对应齐次方程
x k + n + a 1 x k + n − 1 + a 2 x k + n − 2 + ⋯ + a n x k = 0. x_{k+n}+a_1x_{k+n-1}+a_2x_{k+n-2}+\dots+a_nx_k=0. xk+n​+a1​xk+n−1​+a2​xk+n−2​+⋯+an​xk​=0.

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

对于 n n n 阶常系数线性齐次差分方程，特征方程为
λ n + a 1 λ n − 1 + ⋯ + a n = 0. \lambda^n+a_1\lambda^{n-1}+\dots+a_n=0. λn+a1​λn−1+⋯+an​=0.
若特征方程有 n n n 个互不相同的实根 λ 1 , λ 2 , … , λ n \lambda_1,\lambda_2,\dots,\lambda_n λ1​,λ2​,…,λn​，则通解为
c 1 λ 1 k + c 2 λ 2 k + ⋯ + c n λ n k . c_1\lambda_1^k+c_2\lambda_2^k+\dots+c_n\lambda_n^k. c1​λ1k​+c2​λ2k​+⋯+cn​λnk​.
若 λ \lambda λ 是特征方程的 m m m​ 重实根，则通解中对应的项为
c 1 ρ k cos ⁡ ( φ k ) + c 2 ρ k sin ⁡ ( φ k ) , c_1\rho^k\cos(\varphi k)+c_2\rho^k\sin(\varphi k), c1​ρkcos(φk)+c2​ρksin(φk),

ρ = α 2 + β 2 , φ = arctan ⁡ β α . \rho=\sqrt{\alpha^2+\beta^2},\varphi=\arctan\frac{\beta}{\alpha}. ρ=α2+β2 ​,φ=arctanαβ​.

若 λ = α + β i \lambda=\alpha+\beta i λ=α+βi 是特征方程的 m m m 重复根，则通解中对应的项为
( c 1 + c 2 k + ⋯ + c m k m − 1 ) ρ k cos ⁡ ( φ k ) (c_1+c_2k+\dots+c_mk^{m-1})\rho^k\cos(\varphi k) (c1​+c2​k+⋯+cm​km−1)ρkcos(φk)

+ ( c m + 1 + c m + 2 k + ⋯ + c 2 m k m − 1 ) ρ k sin ⁡ ( φ k ) . +(c_{m+1}+c_{m+2}k+\dots+c_{2m}k^{m-1})\rho^k\sin(\varphi k). +(cm+1​+cm+2​k+⋯+c2m​km−1)ρksin(φk).

例子

求斐波那契数列
{ F n = F n − 1 + F n − 2 , F 1 = F 2 = 1. \left\{\begin{aligned} &F_n=F_{n-1}+F_{n-2},\\ &F_1=F_2=1. \end{aligned}\right. {​Fn​=Fn−1​+Fn−2​,F1​=F2​=1.​
的解。

有特征方程
λ 2 − λ − 1 = 0. \lambda^2-\lambda-1=0. λ2−λ−1=0.
解得
λ 1 = 1 − 5 2 , λ 2 = 1 + 5 2 \lambda_1=\frac{1-\sqrt5}{2},\lambda_2=\frac{1+\sqrt5}{2} λ1​=21−5 ​​,λ2​=21+5 ​​
互异，通解为
F n = c 1 ( 1 − 5 2 ) n + c 2 ( 1 + 5 2 ) n . F_n=c_1(\frac{1-\sqrt5}{2})^n+c_2(\frac{1+\sqrt5}{2})^n. Fn​=c1​(21−5 ​​)n+c2​(21+5 ​​)n.
又 F 1 = F 2 = 1 F_1=F_2=1 F1​=F2​=1，解得 c 1 = − 5 5 , c 2 = 5 5 c_1=-\frac{\sqrt5}{5},c_2=\frac{\sqrt5}{5} c1​=−55 ​​,c2​=55 ​​，代入得
F n = 5 5 [ ( 1 + 5 2 ) n − c 2 ( 1 − 5 2 ) n ] . F_n=\frac{\sqrt5}{5}[(\frac{1+\sqrt5}{2})^n-c_2(\frac{1-\sqrt5}{2})^n]. Fn​=55 ​​[(21+5 ​​)n−c2​(21−5 ​​)n].

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

若 x k x_k xk​ 为齐次方程的通解， x k ∗ x_k^* xk∗​ 为非齐次方程的特解，则非齐次方程的通解为 x k + x k ∗ x_k+x_k^* xk​+xk∗​。

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

典型模型

Leslie 模型

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

记时段 t k t_k tk​​ 第 i i i​ 年龄组的种群数量为 x i ( k ) , i = 1 , 2 , … , m x_i(k),i=1,2,\dots,m xi​(k),i=1,2,…,m​，繁殖率为 α i \alpha_i αi​，死亡率为 d i d_i di​， β i = 1 − d i \beta_i=1-d_i βi​=1−di​ 称存活率。

第一年龄组种群数量是时段 t k t_k tk​ 各年龄组繁殖数量之和，即
x 1 ( k + 1 ) = ∑ i = 1 m α i x i ( k ) , x_1(k+1)=\sum_{i=1}^m\alpha_ix_i(k), x1​(k+1)=i=1∑m​αi​xi​(k),
t k + 1 t_{k+1} tk+1​ 时段第 i + 1 i+1 i+1 年龄组的种群数量是 t k t_k tk​ 时段第 i i i 年龄组存活下来的数量，即
x i + 1 ( k + 1 ) = β i x i ( k ) , i = 1 , 2 , … , m − 1. x_{i+1}(k+1)=\beta_ix_i(k),i=1,2,\dots,m-1. xi+1​(k+1)=βi​xi​(k),i=1,2,…,m−1.
记 t k t_k tk​ 时段种群各年龄组的分布向量

X ( k ) = ( x 1 ( k ) x 2 ( k ) ⋮ x m ( k ) ) , X(k)=\begin{pmatrix} x_1(k)\\ x_2(k)\\ \vdots\\ x_m(k) \end{pmatrix}, X(k)=⎝⎜⎜⎜⎛​x1​(k)x2​(k)⋮xm​(k)​⎠⎟⎟⎟⎞​,
记
L = ( α 1 α 2 … α m − 1 α m β 1 0 … 0 0 0 β 2 … 0 0 ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 … β m − 1 0 ) , L=\begin{pmatrix} \alpha_1&\alpha_2&\dots&\alpha_{m-1}&\alpha_m\\ \beta_1&0&\dots&0&0\\ 0&\beta_2&\dots&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\dots&\beta_{m-1}&0 \end{pmatrix}, L=⎝⎜⎜⎜⎜⎜⎛​α1​β1​0⋮0​α2​0β2​⋮0​………⋱…​αm−1​00⋮βm−1​​αm​00⋮0​⎠⎟⎟⎟⎟⎟⎞​,
则
X ( k + 1 ) = L X ( k ) , k = 0 , 1 , … X(k+1)=LX(k),k=0,1,\dots X(k+1)=LX(k),k=0,1,…
当 t 0 t_0 t0​​​ 时段各年龄组人数 X ( 0 ) X(0) X(0)​ 已知，就可以求得 t k t_k tk​​ 时段按年龄组的分布向量
X ( k ) = L k X ( 0 ) , k = 0 , 1 , … X(k)=L^kX(0),k=0,1,\dots X(k)=LkX(0),k=0,1,…
由此算出种群总量。

离散阻滞增长模型

又称离散 Logistic 模型。
Δ x k = r x k ( 1 − x k N ) , k = 0 , 1 , 2 , … \Delta x_k=rx_k(1-\frac{x_k}{N}),k=0,1,2,\dots Δxk​=rxk​(1−Nxk​​),k=0,1,2,…
其中 r r r​ 为固有增长率， N N N​ 为最大容量。上式可化为
x k + 1 − x k x k = r ( 1 − x k N ) , k = 0 , 1 , 2 , … \frac{x_{k+1}-x_k}{x_k}=r(1-\frac{x_k}{N}),k=0,1,2,\dots xk​xk+1​−xk​​=r(1−Nxk​​),k=0,1,2,…
当 x = 0 x=0 x=0​ 时，增长率为 r r r​，随着 x x x 增加，增长率逐渐减小，直到 x = N x=N 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 a_n,b_n,c_n an​,bn​,cn​ 表示第 n n n 代植物基因型为 AA, Aa, aa 的植物占植物总数的百分率，令
x ( 0 ) = ( a 0 , b 0 , c 0 ) T x^{(0)}=(a_0,b_0,c_0)^T x(0)=(a0​,b0​,c0​)T
表示初始分布。

选用 AA 型植物与其他植物结合。由于 AA 与 AA 结合必为 AA，AA 与 Aa 结合为 AA 的可能性为 1/2，AA 与 aa 结合不可能为 AA，得
a n = 1 a n − 1 + 1 2 b n − 1 + 0 c n − 1 , a_n=1a_{n-1}+\frac12b_{n-1}+0c_{n-1}, an​=1an−1​+21​bn−1​+0cn−1​,
同理有
b n = 1 2 b n − 1 + c n − 1 , b_n=\frac12b_{n-1}+c_{n-1}, bn​=21​bn−1​+cn−1​,

c n = 0. c_n=0. cn​=0.

记
M = ( 1 1 / 2 0 0 1 / 2 1 0 0 0 ) M=\begin{pmatrix} 1&1/2&0\\ 0&1/2&1\\ 0&0&0 \end{pmatrix} M=⎝⎛​100​1/21/20​010​⎠⎞​
则有
x ( n ) = M x ( n − 1 ) , n = 1 , 2 , … x^{(n)}=Mx^{(n-1)},n=1,2,\dots x(n)=Mx(n−1),n=1,2,…
即
x ( n ) = M n x ( 0 ) . x^{(n)}=M^nx^{(0)}. x(n)=Mnx(0).
在这种情况下， a n → 1 , b n → 0 , c n → 0 a_n\to1,b_n\to0,c_n\to0 an​→1,bn​→0,cn​→0。

若选用相同基因型植物相结合，则
M = ( 1 1 / 4 0 0 1 / 2 0 0 1 / 4 1 ) M=\begin{pmatrix} 1&1/4&0\\ 0&1/2&0\\ 0&1/4&1 \end{pmatrix} M=⎝⎛​100​1/41/21/4​001​⎠⎞​
此时
a n → a 0 + 1 2 b 0 , b n → 0 , c n → c 0 + 1 2 b 0 . a_n\to a_0+\frac12b_0,b_n\to0,c_n\to c_0+\frac12b_0. an​→a0​+21​b0​,bn​→0,cn​→c0​+21​b0​.

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)

输出如下：
− 5 ( 1 2 − 5 2 ) x 5 + − 5 ( 1 2 + 5 2 ) x 5 -\frac{\sqrt5(\frac12-\frac{\sqrt5}{2})^x}{5}+-\frac{\sqrt5(\frac12+\frac{\sqrt5}{2})^x}{5} −55 ​(21​−25 ​​)x​+−55 ​(21​+25 ​​)x​

给出通解和部分数值解。

Leslie 模型

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

由题意，出生率为 α = ( 0 , 4 , 3 ) T \alpha=(0,4,3)^T α=(0,4,3)T，成活率为 β = ( 0.5 , 0.25 ) T \beta=(0.5,0.25)^T β=(0.5,0.25)T，则
L = ( 0 4 3 0.5 0 0 0 0.25 0 ) L=\begin{pmatrix} 0&4&3\\ 0.5&0&0\\ 0&0.25&0 \end{pmatrix} L=⎝⎛​00.50​400.25​300​⎠⎞​
有
X ( k + 1 ) = L X ( k ) , k = 0 , 1 , 2 , … X(k+1)=LX(k),k=0,1,2,\dots X(k+1)=LX(k),k=0,1,2,…
又
X ( 0 ) = ( 1000 , 1000 , 1000 ) T X(0)=(1000,1000,1000)^T 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