吴世涵的博客
0%

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

发表于 分类于
本文字数： 8.1k 阅读时长 ≈ 7 分钟
差分方程模型是数学建模中常见的模型之一，本文主要介绍差分方程模型的定义、常系数线性齐次差分方程的解法、常系数线性非齐次差分方程的解法、典型模型等内容，并给出了 Python 代码实现。

13. 差分方程模型

定义

设函数 ， 取非负整数，称改变量 为函数 的差分，也称函数 的一阶差分，记 。

一阶差分的差分称二阶差分 ，即

类似有三阶差分、四阶差分……，记 含有未知函数 ​ 的差分的方程称差分方程，差分方程中所含未知函数差分的最高阶数称差分方程的阶，如 阶差分方程形式为 或 若差分方程中所含未知函数和未知函数的各阶差分均为一次的，则称线性差分方程，形式为 对应齐次方程

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

对于 阶常系数线性齐次差分方程，特征方程为 若特征方程有 个互不相同的实根 ，则通解为 若 是特征方程的 ​ 重实根，则通解中对应的项为

若 是特征方程的 重复根，则通解中对应的项为

例子

求斐波那契数列 的解。

有特征方程 解得 互异，通解为 又 ，解得 ，代入得

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

若 为齐次方程的通解， 为非齐次方程的特解，则非齐次方程的通解为 。

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

典型模型

Leslie 模型

考虑对人口进行分组，为简单起见只考虑女性人口，将女性人口按年龄划分称 个年龄组。将时间离散为时段

记时段 ​ 第 ​ 年龄组的种群数量为 ​，繁殖率为 ，死亡率为 ， 称存活率。

第一年龄组种群数量是时段 各年龄组繁殖数量之和，即 时段第 年龄组的种群数量是 时段第 年龄组存活下来的数量，即 记 时段种群各年龄组的分布向量

记 则 当 ​​ 时段各年龄组人数 ​ 已知，就可以求得 ​ 时段按年龄组的分布向量 由此算出种群总量。

离散阻滞增长模型

又称离散 Logistic 模型。 其中 ​ 为固有增长率，​ 为最大容量。上式可化为 当 ​ 时，增长率为 ​，随着 增加，增长率逐渐减小，直到 时，增长率为 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

令 表示第 代植物基因型为 AA, Aa, aa 的植物占植物总数的百分率，令 表示初始分布。

选用 AA 型植物与其他植物结合。由于 AA 与 AA 结合必为 AA，AA 与 Aa 结合为 AA 的可能性为 1/2，AA 与 aa 结合不可能为 AA，得 同理有

记 则有 即 在这种情况下，。

若选用相同基因型植物相结合，则 此时

Python 代码

解差分方程

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
#! /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)

输出如下：

给出通解和部分数值解。

Leslie 模型

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

由题意，出生率为 ，成活率为 ，则 有 又 可以求解，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#! /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()

输出如下：

1
2
3
4
5
6
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

染色体遗传模型

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
#! /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()

输出如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
----- 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