离散动力系统之特征值和特征向量的应用

一路阳光Lemon

已于 2022-04-05 22:33:01 修改

阅读量2.2k

点赞数 3

分类专栏：线性代数数学建模文章标签：抽象代数学习 python 线性代数矩阵

于 2022-04-05 22:32:32 首次发布

本文链接：https://blog.csdn.net/qq_43727222/article/details/123856877

版权

数学建模同时被 2 个专栏收录

5 篇文章 0 订阅

订阅专栏

线性代数

2 篇文章 0 订阅

订阅专栏

本文探讨了利用阶段矩阵模型预测斑点猫头鹰种群灭绝的可能性，分析了森林破坏对幼年猫头鹰存活率的影响。模型显示，如果幼年猫头鹰找到新栖息地的比例低于50%，种群将走向灭绝。然而，如果这一比例提高，猫头鹰数量将逐渐增长。通过对线性动力系统的深入研究，揭示了特征值和特征向量在预测生态系统动态中的关键作用。

摘要由CSDN通过智能技术生成

【背景引入】

1990年，在利用或滥用太平洋西北部大面积森林问题上，北方的斑点猫头鹰成为一个争论的焦点。环境保护学家试图说服联邦政府，如果采伐原始森林（长有200年以上的树木）的行为得不到制止，猫头鹰将濒临灭绝，因为猫头鹰喜好在那里居住。

而木材行业却争辩说猫头鹰不应被划为“濒临灭绝动物”，并引用一些已发表的科学报告来支持其论点，对木材行业来说，如果政府出合新的伐木限制，预计将失去30000 至100000 个工作岗位。

数学生态学家处于争论双方的中间，由于争论的双方都想游说数学生态学家，数学生态学家加快了他们对斑点猫头鹰种群的动力学研究。

猫头鹰的生命周期自然分为三个阶段：幼年期（1岁以前）、半成年期（1~2岁）、成年期（2岁以后) 。猫头鹰在半成年和成年期交配，开始生育繁殖，可活到20岁左右。每一对猫头鹰需要约 1000公顷(4平方英里）的土地作为自己的栖息地。生命周期的关键期是当幼年猫头鹰离开巢的时候。为生存和进入半成年期，一只幼年猫头鹰必须成功地找到一个新的栖息地安家（通常还带有一个配偶）

研究种群动力学的第 1步是建立以每年的种群量为区间的种群模型，时间为大=0.1.2..…通常可以假设在每一个生命阶段雄性和雌性的比例为 1:1，而且只计算雌性猫头鹰，第 k 年的种群量可以用向量 $x_{k} = (j_{k},s_{k},a_{k})$ 表示，其中 $j_{k},s_{k}$ 和 $a_{k}$ 分别代表雌性猫头鹰在幼年期、半成年期和成年期的数量。利用人口统计研究的实际现场数据， R.Lamberson 及其同事设计了下面的“阶段矩阵模型”

$\begin{bmatrix} j_{k+1} \\ s_{k+1} \\ a_{k+1} \end{bmatrix} = \begin{bmatrix} 0&0&0.33 \\0.18&0&0\\ 0&0.71&0.94\end{bmatrix} \begin{bmatrix} j_{k} \\ s_{k} \\ a_{k} \end{bmatrix}$ (等式1）

在这里新的幼年雌性猫头鹰在k+1年中的数量是成年雌性猫头鹰在k年里数量的0.33倍（根据每一对猫头鹰的平均生殖率市定）。此外，18%的幼年雌性猫头鹰得以生存进入半成年期，71%的半成年雌性猫头鹰和 94%的成年雌性猫头鹰生存下来被计为成年猫头鹰。

阶段矩阵模型是形式为 $x_{k+1} = Ax_{k}$ 的差分方程，这种方程被称为动力系统（或离散线性动力系统），因为它描述的是系统随时间推移的变化。

Lamberson 阶段矩阵中，18%的幼年猫头鹰生存率是受原始森林中猫头鹰数目影响最大的项目。事实上，60%的幼年猫头鹰通常生存下来后就会离开自己的巢，但是在 Lamberson 和他的同事们作研究的加州Willow Creek 地区，只有30%的幼年猫鹰在弃巢后能找到新的栖息地，其他的在寻找新家园过程中失踪了。

猫头鹰不能找到新栖息地的一个重要原因是对原始森林分散区域的砍伐加剧了原始森林的分割。当猫头鹰离开森林保护区并穿过一块滥伐地时，被捕食动物袭击的危险大增。接下来将会展示前面讨论的模型如何预测斑点猫头鹰的最终灭绝，但如果 50%的幼年猫头鹰弃巢后能找到新的栖息地，猫头鹰种群将会兴旺起来。本文的目的是剖析线性变换 $x\rightarrow Ax$ 的作用，把它分解为容易理解的元素。本文中出现的矩阵都是方阵。虽然这里讨论的主要应用是离散动力系统（包括上面的斑点猫头鹰），但有关特征值和特征向量的基本概念对纯数学和应用数学都很有用，它们出现的背景要比我们这里考虑的要广泛得多。同样，特征值还被用来研究微分方程和连续动力系统，为工程设计提供关键知识，自然地，也出现在物理和化学等领域里。

【思考一下1】

Q1：背景中讨论的模型，如何预测斑点猫头鹰的最终灭绝？

A1：在接下来的内容中，我们将讨论这个问题。

Q2: 特征值和特征方程仅被用来研究离散动力系统吗？

A2：当然不是。

还可被用来研究微分方程，连续动力系统，涉及工程设计，物理和化学等领域。

【问题描述】

捕食者-食饵系统

【例1】在加利福尼亚州的红木森林深处，作为老鼠的主要捕食者，斑点猫头鹰的食物有 80%是老
鼠。利用线性动力系统来建立猫头鹰和老鼠的自然系统模型。（其实，这个模型在某些方面与现实不符，但它能够为环境科学家们所用的更复杂的非线性模型的研究提供一个起点。）
用 $x_{k} = \begin{bmatrix} O_{k} \\ R_{k}\end{bmatrix}$ 表示在时间 $k$ （ $k$ 的单位是月），猫头鹰和老鼠的数量， $O_{k}$ 是在研究区域猫头鹰的数量， $R_{k}$ 是老鼠的数量（单位是千只）。设

$O_{k+1} = (0.5)O_{k} + (0.4)R_{k}$

$R_{k+1} = -p \cdot O_{k} + (1.1)R_{k}$ (方程组1）

【问题分析】

其中p是被指定的正参数。

第1个方程中的 $(0.5)O_{k}$ 表示，如果没有老鼠为食物，每月仅有一半的猫头鹰存活下来，而第2个方程的 $(1.1)R_{k}$ 表明如果没有猫头鹰捕食老鼠，那么老鼠的数量每月增长10%。

假如有足够多的老鼠 $(0.4)R_{k}$ 表示猫头鹰增长的数量，而负项 $(-p \cdot O_{k})$ 表示由于猫头鹰的捕食所引起的老鼠的死亡数量。（事实上，一只猫头鹰每月平均吃掉1000p只老鼠.）当p=0.104时，预测该系统的发展趋势。

【问题解答】

当p=0.104时，算出（方程组1）的系数矩阵A的特征值是 $\lambda _{1} = 1.02$ 和 $\lambda_{2} = 0.58$ 。

对应的特征向量是:

$v_{1} = \begin{bmatrix} 10\\13\end{}$ ， $v_{2} = \begin{bmatrix}5\\1 \end{}$

初试向量 $x_{0}$ 可以表示为 $x_{0} = c_{1}v_{1}+c_{2}v_{2}$ ，那么对 $k \geqslant 0$ ，

$x_{k} = c_{1}(1.02)^{k}v_{1}+c_{2}(0.58)^{k} v_{2}$

$=c_{1}(1.02)^{k}\begin{bmatrix} 10\\13\end{} + c_{2}(0.58)^{k} \begin{bmatrix} 5\\ 1 \end{}$

当 $k\rightarrow + \infty$ ， $(0.58)^{k}$ 很快趋于0。

假设 $c_{1} > 0$ ，那么对所有足够大的k， $x_{k}$ 近似等于

$c_{1}(1.02)^{k}\begin{\bmatrix}10\\13 \end{}v_{1}$

我们记为

$x_{k} \approx c_{1}(1.02)^{k}\begin{\bmatrix}10\\13 \end{}$ （近似式1）

随着k的增大，（近似式1）的近似程度会更好，故，对足够大的k，

$x_{k+1} \approx c_{1}(1.02)^{k+1}\begin{\bmatrix}10\\13 \end{} = (1.02)c_{1}(1.02)^{k}\begin{\bmatrix}10\\13 \end{} \approx 1.02x_{k}$ （近似式2）

（近似式2）表明最终 $x_{k}$ 的2个分量（猫头鹰和老鼠的数量）每月以大约 1.02 的倍数增长，即月增长率为2%。由（近似式1）， $x_{k}$ 近似于（10,13）的倍数，因此，x的2个分量之比率也近似于 10与13的比率，也就是说，对应每10只猫头鹰，大致有13000只老鼠。
例1说明了有关动力系统 $x_{k+1} =Ax_{k}$ 的两个基本事实，若A是nxn矩阵，它的特征值满足
$\left | \lambda_{1} \right | \geqslant 1$ 和 $1\geqslant \left | \lambda_{j} \right |$ ，j=2，…,n, $v_{1}$ 是 $\lambda_{1}$ 对应的特征向量，假如 $x_{0}$ 由（等式1）给出且 $c_{1} \neq 0$ ，那么对足够够大的k，

$x_{0} = c_{1}v_{1} + ...+c_{n}v_{n}$ (等式1）
$x_{k+1} \approx \lambda_{1}x_{k}$ (近似式3）

和

$x_{k} \approx c_{1}(\lambda_{1})^{k}v_{k}$ (近似式4）
式（3）和（4）的近似精度可根据需要通过取足够大的k来得到。

由式（3）， $x_{k}$ 每时段最终以近似 $\lambda_{1}$ 的倍数增长，因此， $\lambda_{1}$ 确定了系统的最终增长率。同样由式（4)，对足够大的k， $x_{k}$ 的2个分量之比近似等于 $\lambda_{1}$ 对应分量之比。

【思考一下2】

Q3：解的几何意义是什么？

A3:当A为2x2矩阵时，可以通过系统发展趋势的几何描述来补充解释代数计算。我们可以把方程 $x_{k+1} = Ax_{k}$ 看作是 $R^{2}$ 中的初始点 $x_{0}$ 被映射 $x \mapsto Ax$ 重复变换的描述。由 $x_{0},x_{1},..$ 组成的图形称为是动力系统的轨迹。

Q4：如何理解动力系统的轨迹？其本质是什么？

A4：描述线性方程如何变化。

本质：线性方程解的集合。

Q5： $x_{k+1} = Ax_{k}$ (k = 0，1，2...)（近似式3）与马尔科夫链的有什么关系？

A5：这种线性变换的等式关系可以应用于Markov Chain。

矩阵A = 移民矩阵M； $x_{0}$ = 城市和郊区的初始人口分布； $x_{k}$ ：表示k年后的人口分布。

Q6：特征值的本质是什么？特征值可能取负数吗？如果可以，会出现什么结果？

A6：本质：n次多项式的根，负数即多项式根为复数。

可以。

结果如（图1)，负数 = 复特征值。

（图1）含有复根的线性方程的解组成的轨迹

（图1）的代码展示：

#工具包导入
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import  Axes3D

#计算矩阵的特征值和特征向量
mat = np.array([[0.5,-0.6],[0.75,1.1]])
eigenvalue,featurevector = np.linalg.eig(mat)
print('特征值：',eigenvalue)
print('特征向量：',featurevector)

#计算前10次线性变换的迭代
A = [[0.5,-0.6],[0.75,1.1]]
x0 = [[2],[0]]
x1 = (np.dot(A,x0))
x2 = (np.dot(A,x1))
x3 = (np.dot(A,x2))
x4 = (np.dot(A,x3))
x5 = (np.dot(A,x4))
x6 = (np.dot(A,x5))
x7 = (np.dot(A,x6))
x8 = (np.dot(A,x7))
x9 = (np.dot(A,x8))
x10 = (np.dot(A,x9))
#打印每次的迭代值
print("前10次线性变化结果：")
print(x0)
print(x1)
print(x2)
print(x3)
print(x4)
print(x5)
print(x6)
print(x7)
print(x8)
print(x9)
print(x10)

#点的坐标，即前10次迭代
x = [2,1,-0.4,-1.64,-2.224,-1.9184,-0.84544,0.565696,1.7505536,2.23518976,1.82575002]
y = [0,1.5,2.4,2.34,1.344,-0.1896,-1.64736,-2.446176,-2.2665216,-1.18025856,0.3781079]
#给图中点做标注
txt = ['x0','x1','x2','x3','x4','x5','x6','x7','x8','x9','x10']
#显示负号
plt.rcParams['axes.unicode_minus'] = False
# 用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
plt.scatter(x, y)
for i in range(len(x)):
    plt.annotate(txt[i], xy = (x[i], y[i]), xytext = (x[i]+0.1, y[i]+0.1)) # 这里xy是需要标记的坐标，xytext是对应的标签坐标

plt.axhline(y=0.0, c='r', ls='-', lw=1)  # 垂直于y轴的参考线
plt.axvline(x=0.0, c='r', ls='-', lw=1)  # 垂直于x轴的参考线
#将离散的点依次用虚线连接
plt.plot(x, y, 'xb--')
#图片显示
plt.show()

Q7：是什么原因导致解看起来在旋转？

A7：共轭复数对导致旋转。（此处需要用复特征值知识来解答）

Q8：如果特征值不是负数，离散系统的轨迹还会出现旋转的情况吗？为什么？

A8：不会。

因为，若特征值不是复数，就不会出现共轭复数对，不会导致旋转。

eg：当 $A = \begin{bmatrix} 0.80&0 \\ 0&0.64 \end{}$ 时，动力系统 $x_{k+1} = Ax_{k}$ 的轨迹如下。

（图2）完整动力系统轨迹

【补充】

原点 = 动力系统的吸引子：所有轨迹点都趋于一点。

Q9：在什么条件下，会出现这种情况？

A9：两个特征值的绝对值都＜1

Q10：当两个特征值的绝对值都 > 1时，会出现什么情况？

A10：会出现排斥子。

eg：方程 $x_{k+1} = Ak_{k}$ 的解的若干条典型轨迹如图3，其中 $A = \begin{bmatrix}1.44&0\\ 0&1.2 \end{}$ 。

计算出矩阵A的特征值，当 $x_{0} = \begin{bmatrix}c_{1}\\c_{2} \end{}$ 时，有

$x_{k} = c_{1}(1.44)^{k} \begin{bmatrix}1\\0 \end{} + c_{2}(1.2)^{k} \begin{bmatrix} 0\\1 \end{}$ 。

结论：两项的值随着k的增大而增大，第一项的值增大的更快。

（图3）原点是排斥子（Q10 方程的解组成的轨迹图）

Q11：如何理解吸引子，排斥子？

A11：由Q8得出结论，当 $k\rightarrow +\infty$ ， $x_{k}$ 、 $(0.8)^{k}$ 和 $(0.64)^{k}$ 都趋于零，解的集合形成的轨

迹如图2，视觉上达到万物归一的效果，直观上展示了原点即吸引子。

在Q10中，除过（常数）零解， $x_{k+1} = Ak_{k}$ 的解都是无界的，离原点而去。

(图4） Q8中，当 $c_{1} = c_{2} = 1$ 时的单个轨迹图

（图4)python代码如下

#工具包导入
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import  Axes3D
import math
#当c1=c2=1
v1 = np.array([[1],[0]])
v2 = np.array([[0],[1]])
c1 = 1
c2 = 1
for k in range(7):
    xk = c1 * (math.pow(0.8, k)) * v1 + c2 * (math.pow(0.64, k)) * v2
#当c1=c2=2
c1 = 2
c2 = 2
for k in range(7):
    xk = c1 * (math.pow(0.8, k)) * v1 + c2 * (math.pow(0.64, k)) * v2

#点的坐标，即前7次迭代
x = [1.0, 0.8, 0.64, 0.512, 0.4096, 0.10737418, 0.32768, 0.262144]
y = [1.0, 0.64, 0.4096, 0.262144, 0.16777216, 0.10737418, 0.10737418, 0.06871948]
#给图中点做标注
txt = ['x0','x1','x2','x3','x4','x5','x6','x7']
#显示负号
plt.rcParams['axes.unicode_minus'] = False
# 用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
plt.scatter(x, y)
for i in range(len(x)):
    plt.annotate(txt[i], xy = (x[i], y[i]), xytext = (x[i]+0.001, y[i]+0.001)) # 这里xy是需要标记的坐标，xytext是对应的标签坐标

plt.axhline(y=0.0, c='r', ls='-', lw=1)  # 垂直于y轴的参考线
plt.axvline(x=0.0, c='r', ls='-', lw=1)  # 垂直于x轴的参考线
#将离散的点依次用虚线连接
plt.plot(x, y, 'xb--')
#图片显示
plt.show()

Q12：吸引子和排斥子可能出现在同一方程中吗？若出现，动力系统的轨迹会发生什么？

A12：可能。

当吸引子和排斥子出现在同一方程中，即方程的特征值的绝对值既有> 1的值，也有 < 1

的值。

eg：观察方程 $y_{k+1} = Dy_{k}$ 典型解的轨迹，其中 $D = \begin{bmatrix}2.0&0\\0&0.5 \end{}$ 。

计算D的特征值，若 $y_{0} = \begin{bmatrix}c_{1}\\c_{2} \end{}$ 则， $y_{k} = c_{1}2^{k}\begin{bmatrix}1\\0 \end{} + c_{2}(0.5)^k\begin{bmatrix}0\\1 \end{}$ 。

图5 Q12的动力系统的轨迹

结论：原点 = 鞍点：它在某些方向吸引解，某些方向排斥解。

图6 当 $c_{1} = c_{2} = 1$ 时，单个轨迹图

图6 python代码如下

#工具包导入
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math

#计算特征值和特征向量
mat = np.array([[2.0,0],[0,0.5]])
eigenvalue,featurevector = np.linalg.eig(mat)
print('特征值:',eigenvalue)
print('特征向量:',featurevector)
#计算前5次数的迭代
D = [[2.0,0],[0,0.5]]
x0=[[1,0],[0,1]]
x1=(np.dot(D,x0))
x2=(np.dot(D,x1))
x3=(np.dot(D,x2))
x4=(np.dot(D,x3))
x5=(np.dot(D,x4))
#打印前5次迭代结果
print('打印前6次，迭代结果：')
print(x0)
print(x1)
print(x2)
print(x3)
print(x4)
print(x5)

#画轨迹方程
#当c1=c2=1
c1=1
c2=1
#设置特征向量
v1=np.array([[1],[0]])
v2=np.array([[0],[1]])
for k in range(6):
    xk = c1 * (math.pow(2.0, k)) * v1 + c2 * (math.pow(0.5, k)) * v2
#当c1=c2=2
c1 = 2
c2 = 2
for k in range(6):
    xk = c1 * (math.pow(2.0, k)) * v1 + c2 * (math.pow(0.5, k)) * v2

#点的坐标，即前6次迭代
x = [1, 2, 4, 8, 16, 32]
y = [1,0.5,0.25,0.125,0.0625,3.125e-02]
#给图中点做标注
txt = ['x0','x1','x2','x3','x4','x5']
#显示负号
plt.rcParams['axes.unicode_minus'] = False
# 用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
plt.scatter(x, y)
for i in range(len(x)):
    plt.annotate(txt[i], xy = (x[i], y[i]), xytext = (x[i]+0.001, y[i]+0.001)) # 这里xy是需要标记的坐标，xytext是对应的标签坐标

plt.axhline(y=0.0, c='r', ls='-', lw=1)  # 垂直于y轴的参考线
plt.axvline(x=0.0, c='r', ls='-', lw=1)  # 垂直于x轴的参考线
#将离散的点依次用虚线连接
plt.plot(x, y, 'xb--')
#图片显示
plt.show()

Q13：由图1和图2，可以推测出动力系统的轨迹与什么有关吗？

A13：轨迹特征值，特征向量。

Q14：Q8中的矩阵是对角矩阵，如果是非对角矩阵，轨迹又会是什么样子？

A14：若是非对角矩阵，且特征值为复数，轨迹将发生旋转。（类似于Q6中的图)

eg:观察方程 $x_{k+1} = Ax_{k}$ 的图像，矩阵 $A = \begin{bmatrix}0.8&0.5\\-0.1&1.0 \end{}$ 。

图7 为Q14的轨迹图像，与复特征值相关的旋转

图8 当 $c_{1} = c_{2} = 1$ 时，的单个轨迹图

图8的python代码

#导入工具包
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math
#计算特征值和特征向量
#设置矩阵mat
mat = np.array([[0.8,0.5],[-0.1,1.0]])
eigenvalue,featurevector = np.linalg.eig(mat)
print('特征值：',eigenvalue)
print('特征向量：',featurevector)
#计算前10次迭代
A = [[0.8,0.5],[-0.1,1.0]]
x0= [[1-2j,1],[1+2j,1]]
x1 = (np.dot(A,x0))
x2 = (np.dot(A,x1))
x3 = (np.dot(A,x2))
x4 = (np.dot(A,x3))
x5 = (np.dot(A,x4))
x6 = (np.dot(A,x5))
x7 = (np.dot(A,x6))
x8 = (np.dot(A,x7))
x9 = (np.dot(A,x8))
x10 = (np.dot(A,x9))
#打印前10次迭代结果
print('打印前11次迭代结果：')
print(x0)
print(x1)
print(x2)
print(x3)
print(x4)
print(x5)
print(x6)
print(x7)
print(x8)
print(x9)
print(x10)
#画轨迹方程
#设置c1=c2=1
c1=1
c2=1
#设置特征向量
v1=np.array([[1-2j],[1]])
v2=np.array([[1+2j],[1]])
for k in range(10):
    xk = c1 * (math.pow(0.8, k)) * v1 + c2 * (math.pow(1, k)) * v2
#当c1=c2=2
c1 = 2
c2 = 2
for k in range(11):
    xk = c1 * (math.pow(0.8, k)) * v1 + c2 * (math.pow(1, k)) * v2

#点的坐标，即前10次迭代
x = [1-2,1.3-0.6,1.49+0.62,1.577+1.626,1.5721+2.3998,1.48933+2.93754,1.344509+3.247742,1.1541857+3.3490266,0.93470161+3.26766718,0.70140505+3.03512831, 0.46803273+2.68571386]
y = [1,0.9, 0.77, 0.621,0.4633,0.30609,0.157157,0.0227061,-0.09271247,-0.18618263,-0.25632314]
#给图中点做标注
txt = ['x0','x1','x2','x3','x4','x5','x6','x7','x8','x9','x10']
#显示负号
plt.rcParams['axes.unicode_minus'] = False
# 用来正常显示中文标签
plt.rcParams['font.sans-serif']=['SimHei']
plt.scatter(x, y)
for i in range(len(x)):
    # 这里xy是需要标记的坐标，xytext是对应的标签坐标
    plt.annotate(txt[i], xy = (x[i], y[i]), xytext = (x[i]+0.001, y[i]+0.001))

plt.axhline(y=0.0, c='r', ls='-', lw=1)  # 垂直于y轴的参考线
plt.axvline(x=0.0, c='r', ls='-', lw=1)  # 垂直于x轴的参考线
#将离散的点依次用虚线连接
plt.plot(x, y, 'xb--')
#图片显示
plt.show()

【问题结论】

现在回顾一下本文开头的问题，预测斑点猫头鹰是否会灭绝？我们使用动力系统 $x_{k+1} = Ax_{k}$ 为猫头鹰建立种群模型。在该模型中， $x_{k} = (j_{k},s_{k},a_{k})$ 的分量分别表示在时间化时幼年、半成年和成年雌性猫头型的数量，A为阶段矩阵

$A = \begin{bmatrix}0&0&0.33\\0.18&0&0\\0&0.17&0.94 \end{}$

由 pycharm算出A的特征值大约是 $\lambda _{1} = 0.98,\lambda_{2} = -0.02+0.21i,\lambda_{3} = -0.02-0.21i$ 。因为 $\left | \lambda_{2} \right |^{2} = \left | \lambda_{3} \right |^{2} =(-0.02)^{2}+(0.21)^{2} = 0.0445$ ，故三个特征值的绝对值都小于1。现在，让A作用在复向量空间 $C^{3}$ 。因为A有3个相异的特征值，故对应的了个特征向量是线性无关的，它们形成 $C^{3}$ 的一个基。记这些特征向量为 $v_{1}$ ， $v_{2}$ 和 $v_{3}$ 。那么 $x_{k+1} = Ax_{k}$ 的通解（用 $C^{3}$ 的向量表示）的形式是：

$x_{k} = c_{1}(\lambda_{1})^{k}v_{1}+c_{2}(\lambda_{2})^{k}v_{2}+c_{3}(\lambda_{3})^{k}v_{3}$ (等式2）

若初始向量 $x_{0}$ 是实向量，由于A是实矩阵，故 $x_{1} = Ax_{0}$ 也是实向量。同样，方程 $x_{k+1} = Ax_{k}$ 表明(等式2）左边的x也是实向量，尽管它表示为复向量的和。但由于所有特征值都小于1，所以(等式2）右边的每一项都趋于零向量。

因此，实序列 $\left \{ x_{k} \right \}$ 也趋于零向量。很不幸，该模型预测斑点猫头鹰最终会全部灭亡。

计算矩阵A的特征值和特征向量的python代码如下：

#导入所需工具包
import numpy as np
#设置矩阵A
mat = np.array([[0, 0, 0.33],
              [0.18, 0, 0],
              [0, 0.71, 0.94]])
#计算特征值和特征向量
eigenvalue, featurevector = np.linalg.eig(mat)
#打印特征值和特征向量
print("特征值：", eigenvalue)
print("特征向量：", featurevector)

执行结果如下：

猫头鹰还有希望吗？在矩阵A中的元素 18%源于如下事实：
尽管有 60%的幼年猫头鹰能够活下来并离巢去寻找新的栖息地，但在这60%的幼年猫头鹰中，仅有 30%的猫头鹰能活下来找到新的栖息地，森林中裸露地域的数量使得搜寻工作变得更困难和更危险，这严重影响了寻找栖息地过程中的猫头鹰的存活率。

在【背景引入】中，我们提到，若幼年猫头鹰搜寻存活率达到50%，猫头鹰种族就不会灭绝。

证明这个假设市场很容易的，首先将矩阵A（2，1）的元素换0.18换成0.3，

于是，矩阵A将变成 $A = \begin{bmatrix}0&0&0.33\\0.3&0&0\\0&0.17&0.94 \end{}$ ，计算其特征值和特征向量如下

计算python代码如下：

#导入所需工具包
import numpy as np
#设置矩阵A
mat = np.array([[0, 0, 0.33],
              [0.3, 0, 0],
              [0, 0.71, 0.94]])
#计算特征值和特征向量
eigenvalue, featurevector = np.linalg.eig(mat)
#打印特征值和特征向量
print("特征值：", eigenvalue)
print("特征向量：", featurevector)

计算结果如下：

现在矩阵A的特征值是 $\lambda_{1} = 1.01,\lambda_{2} = -0.03+0.26i,\lambda_{3} = -0.03-0.26i$ ，对于 $\lambda_{1}$ 的特征向量为 $v_{1} = (10,3,31)$ ，并设 $v_{2}$ 和 $v_{3}$ 是对应于 $\lambda_{2}$ 和 $\lambda_{3}$ 的复特征向量。此时，等式2变成

$x_{k} =c_{1}(1.01)^{k}v_{1}+c_{2}(-0.03+0.26i)^{k}v_{2}+c_{3}(-0.03-0.26i)^{k}v_{3}$ （等式3）

当 $k\rightarrow +\infty$ 时，向量 $v_{2}$ ， $v_{3}$ 趋于零，因此 $x_{k}$ 越来越接近实向量 $c_{1}(1.01)^{k}v_{1}$ ，而且，可以证明在 $x_{0}$ 的初始分解中的常量 $c_{1}$ 在 $x_{0}$ 的元素非负时是正的，因此，猫头鹰数量的长期增长率是 1.01，即猫头鹰的数量会缓慢增长。

特征向量n描述了猫头鹰在3个年龄段数量的最终分布：每31只成年猫头鹰对应大约10只幼年猫头鹰和3只半成年猫头鹰。

一路阳光Lemon

关注

3
点赞
踩
8

收藏

觉得还不错? 一键收藏
0
评论
离散动力系统之特征值和特征向量的应用

【背景引入】1990年，在利用或滥用太平洋西北部大面积森林问题上，北方的斑点猫头鹰成为一个争论的焦点。环境保护学家试图说服联邦政府，如果采伐原始森林（长有200年以上的树木）的行为得不到制止，猫头鹰将濒临灭绝，因为猫头鹰喜好在那里居住。而木材行业却争辩说猫头鹰不应被划为“濒临灭绝动物”，并引用一些已发表的科学报告来支持其论点，对木材行业来说，如果政府出合新的伐木限制，预计将失去30000 至100000 个工作岗位。数学生态学家处于争论双方的中间，...
复制链接

扫一扫

专栏目录