数学建模算法-day[17]

8 差分方程

差分方程是在离散的时间点上描述研究对象动态变化规律的数学表达式。有的实际问题本身就是以离散形式出现的，也有的是将现实世界中随时间连续变化的过程离散化。差分方程与微分方程都是描述状态变化问题的机理建模方法，是同一建模问题的两种思维（离散或连续）方式。

利用差分方程建模，通常是把问题看作一个系统，考察系统状态变量的变化。即首先考察任意两个相邻位置（时间或空间），通常称为一个微元，考察状态值在这个微元内的变化，即输入、输出变化情况，分析变化的原因，进而利用自然科学中的一些相应规律，如质量守恒、动量守恒、能量守恒等公理或定律，建立微元两端状态变量之间的关联方程。其遵循的一个基本准则就是：未来值=现在值+变化值。

若记$x_2$为第$k$个时刻或位置的状态变量值，则把各个状态变量的状态值次序排列，就形成了一个有序序列{xk}k=0n\{x_k\}_{k=0}^n{xk​}k=0n​。若序列中的$x$和其前一个状态或前几个状态值$x_i(0\le i<k)$存在某种关联，则把它们的关联关系用代数方程的形式表达出来，即建立状态之间的关联方程，称为差分方程。差分方程也称为递推关系。

8.1 差分方程建模

例8.1贷款问题。在现实生活中，经常会遇到贷款问题，如购房、买车、投资等，如何根据自身偿还能力，确定合适的贷款额度及偿还期限，是每个借款人应考虑的现实问题。

问题提出
假定某消费者购房需贷款30万元，期限为30年，已知贷款年利率为5.1%，采用固定额度还款方式，问每月应还款额是多少。

问题分析
当月欠款=上月欠款的当月本息-当月还款。

基本假设
假定在还款期限内，利率保持不变。

模型建立
记$y_n$为第$n$个月的欠款总额(单位：元)；$r$为月利率,$r=\frac{0.051}{12}\times 100\%=0.425\%$；$x$为当月还款额度(单位:元)；$N$为还款期限,$N=12\times 30=360$(月)；$Q$为贷款总额(单位:元)。则数学模型为
$$\{\begin{array}{rl}& y_n=(1+r)y_{n-1}-x,n=1,2,\cdots ,N,\\ & y_0=Q.\end{array}\text{\hspace{1em}(8.1)}$$

模型求解

由式(8.1),可通过递推方法求得
$$\begin{array}{rl}{y}_n& =(1+r)y_{n-1}-x\\ & =(1+r)\left[(1+r)y_{n-2}-x\right]-x\\ & ⋮\\ & =(1+r{)}^n{y}_0-x\left[1+(1+r)+(1+r{)}^2+\cdots +(1+r{)}^{n-1}\right]\\ & =(1+r{)}^n{y}_0-x\frac{(1+r{)}^n-1}{r}.\end{array}$$
每月应还款的确定：由到期应全部还清的条件，即
$$\begin{gathered} y_N=y_{360}=0, \\ 0=y_N=(1+r{)}^{N}Q-x\frac{(1+r{)}^N-1}{r}, \end{gathered}$$
解得
$$x=\frac{(1+r{)}^{N}Qr}{(1+r{)}^N-1}=\frac{(1+0.00425{)}^{360}\times 300000\times 0.00425}{(1+0.00425{)}^{360}-1}=1628.85(\text{元}).$$
到期后累计还款额度$1628.85\times 360=586386$元。

例8.2(续例6.18) 再论美国人口增长模型。以$\Delta t=10$年作为一个时间间隔步长，记$x_k$为$t=k$时的人口数量(单位：百万)，考察从$t=k$到$t=k+1$时段内人口的变化量。若假设美国人口增长服从Logistic规律，则可建立如下所示的差分方程模型：
$$x_{k+1}-x_k=r(1-s{x}_k)x_k\Delta t,\text{\hspace{1em}(8.2)}$$
式中：$r$为美国人口的固有增长率;$s$为阻滞系数。

因而建立如下的差分方程模型：
$$\{\begin{array}{l}{x}_{k+1}=(1+10r-10sr{x}_k)x_k=\alpha x_k+\beta x_k^2,k=0,1,\cdots ,\\ x_0=\mathrm{3.9.}\end{array}\text{\hspace{1em}(8.3)}$$
式中：$\alpha =1+10r;\beta =-10rs<0$。

该模型是一个非线性一阶差分方程。下面使用线性最小二乘拟合模型式(8.3)中的位置参数$\alpha$和$\beta$.

记已知的22个人口数据分别为$x_k(k=0,1,\cdots ,21),$用$=0,1,\cdots ,21$分别表示1790年、1800年、…、2000年。把已知的22个数据代入式(8.3)中的第一式，得到关于$\alpha ,\beta$的超定线性方程组
$$\left[\begin{array}{cc}{x}_0& x_0^2\\ x_1& x_1^2\\ ⋮& ⋮\\ x_{20}& x_{20}^2\end{array}\right]\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{21}\end{array}\right],$$
解之，得到$\alpha ,\beta$的最小二乘估计值
$$\hat{\alpha }=1.2095,\hat{\beta }=-0.0004,$$
2010年人口的预测值为307.6809百万。

clc;clear
a=readmatrix('美国人口数据.xlsx')
x=a([2:2:end],:)';x=x(~isnan(x))
A=[x(1:end-1),x(1:end-1).*x(1:end-1)];
b=x(2:end);
ab=pinv(A)*b
x2010=ab(1)*x(end)+ab(2)*x(end)^2

例8.3(续例6.2) 目标跟踪问题。

问题分析
把导弹与乙舰看作两个运动的质点$P(x(t),y(t))$和$Q(\tilde{x}(t),\tilde{y}(t))$,该问题就变成两个质点随时间的运动问题。

基本假设
(1)忽略潮流对两个质点运动的阻尼作用，即始终假定导弹和乙舰以恒定速度运动。

(2)导弹运动方向自始至终都指向乙舰, 即任意时刻导弹运动轨迹曲线的切线与P、Q两点之间的割线重合。

模型建立
首先把时间等间距离散化为一系列时刻:
$t_0$<$t_1$<$t_2$<$\cdots$<$t_k$<$\cdots$,式中:$△t=t_{k+1}-t_k$为等时间步长。

记$u$为导弹运行的速度,则$u=5v_0$。进一步地,记$P(x_k,y_k)$为质点$P$在$t=t_k$时刻的位置,则$Q$点位置是$Q_k=(1,v_0{t}_k)$,如图8.1所示,从点$P_k$到$Q_k$构成的割线向量$\overrightarrow{P_k{Q}_k}$为
$$\overrightarrow{P_k{Q}_k}=(1-x_k,v_0{t}_k-y_k),$$
式中: $x_k=x(t_k);y_k=y(t_k)$。

图8.1导弹追踪示意图

由基本假设(2)可知,$P$点的运动方向始终指向$Q$,故向量$\overrightarrow{P_k{Q}_k}$ 的方向就是导弹在 $t=t_k$ 时刻运动的方向,其方向向量可由如下单位向量(方向余弦)表示:
$$e^{(k)}=(e_1^{(k)},e_2^{(k)})=\frac{\overrightarrow{P_k{Q}_k}}{||\overrightarrow{P_k{Q}_k||}},$$
其中
$$\begin{gathered} ||\overrightarrow{P_k{Q}_k||}=\sqrt{(1-x_k{)}^2+(v_0{t}_k-y_k{)}^2} \\ {e}_1^{(k)}=\frac{1-x_k}{(1-x_k{)}^2+(v_0{t}_k-y_k{)}^2} \\ \frac{v_0{t}_k-y_k}{(1-x_k{)}^2+(v_0{t}_k-y_k{)}^2} \end{gathered}$$
以时间从 $t_k$ 到 $t_{k+1}$ 作为一个微元,当运动时间从 $t_k$ 变为 $t_{k+1}$ 时,在这个微小的时间单元内，假定导弹质点的运动方向不变，则在$t=t_{k+1}=t_0+(k+1)\Delta t$时刻，$P$点的位置为$P_{k+1}(x_{k+1},y_{k+1})$,满足
$$\{\begin{array}{rl}& x_{k+1}=x_k+u{e}_1^{(k)}\Delta t,\\ & y_{k+1}=y_k+u{e}_2^{(k)}\Delta t,\\ & x_0=0,y_0=0,\end{array}\text{\hspace{1em}(8.4)}$$
式中:$u=5v_0$是导弹的运行速度;$u{e}_1^{(k)}$为导弹的速度向量在$x$轴方向的投影分量;$u{e}_2^{(k)}$为导弹的速度向量在$y$轴方向的投影分量;$x_0,y_0$为导弹的初始位置。

这是一个关于参变量(时间$\Delta t$)的差分方程组，令$k=0,1,2,\cdots ,$即可求出在一系列离散时间点上的导弹位置

模型求解
取$v_0=1,\Delta =0.00005,$计算结果如图8.2所示，即乙舰大约行驶到0.2084处时被击中，经过的时间大约为0.2084。

图8.2 导弹追踪计算结果

clc;clear
v=1;u=5*v;h=0.00005;
x=0;y=0;t=0;
plot(x,y,'.'),hold on
while x<0.99999pq=[1-x,v*t-y];x=x+u*pq(1)/norm(pq)*h;y=y+u*pq(2)/norm(pq)*h;t=t+h;plot(x,y,'.')
end
x,y
xlabel('$x$','Interpreter','latex')
ylabel('$y$','Interpreter','Latex','Rotation',0)

例8.4 商品销售量预测。某商品前5年的销售量见表8.1。现希望根据前五年的统计数据预测第6年起该商品在各季度中的销售量。

表8.1 前5年销售数据表

	第1年	第2年	第3年	第4年	第5年
第1季度	11	12	13	15	16
第2季度	16	18	20	24	25
第3季度	25	26	27	30	32
第4季度	12	14	15	15	17

从表8.1可以看出，该商品在前5年相同季节里的销售量呈增长趋势，而在同一年中销售量先增后减，第一季度的销售量最小而第3季度的销售量最大。预测该商品以后的销售情况，根据本例中数据的特征，可以用回归分析方法按季度建立4个经验公式，分别用来预测以后各年同一季度的销售量。例如，如认为第1季度的销售量大体按线性增长，可设销售量$y_t^{(1)}=at+b$,如如下Matlab程序:

clc;clear
a=readmatrix('销售数据.xlsx')
y1=a(2,:);y1=y1(~isnan(y1))';
x=[[1:5]',ones(5,1)];
z=x\y1

求得$a=z(1)=1.3,b=z(2)=9.5$。

根据$y_t^{(1)}=1.3t+9.5$，预测第6年起第一季度的销售量为$y_6^{(1)}=17.3,y_7^{(1)}=18.6,\cdots$。由于数据少，用回归分析效果不一定好。

如认为销售量并非逐年等量增长而是按前一年或前几年同期销售量的一定比例增长的，则可建立相应的差分方程模型。仍以第1季度为例，简单起见，不再引入上标，以$y_t$表示第$t$年第1季度的销售量，建立形式如下的差分公式：
$$y_t=a_1{y}_{t-1}+a_2,$$
或
$$y_t=a_1{y}_{t-1}+a_2{y}_{t-2}+a_3$$
等。

上述差分方程中的系数不一定能使所有统计数据吻合，较为合理的方法是用最小二乘法求一组总体吻合较好的数据。以建立二阶差分方程$y_t=a_1{y}_{t-1}+a_2{y}_{t-2}+a_3$为例，选取$a_1,a_2,a_3$使得对于已知观测数据$y_t(t=1,2,3,4,5)$,使
$$\sum _{t=3}^5[y_t-(a_1{y}_{t-1}+a_2{y}_{t-2}+a_3){]}^2$$
达到最小。编写Matlab程序如下：

clc;clear
a=readmatrix('销售数据.xlsx')
y1=a(2,:);y1=y1(~isnan(y1))';
b=y1(3:5);
A=[y1(2:4),y1(1:3),ones(3,1)];
z=pinv(A)*b

求得$a_1=z(1)=-1,a_2=z(2)=3,a_3=z(3)=-8$。即所求二阶差分方程为
$$y_t=-y_{t-1}+3y_{t-2}-8.$$
显然这一差分方程恰好使所有统计数据吻合，但这只是一个巧合。根据这一方程，可迭代求出以后各年第1季度销售量的预测值$y_6=21,y_7=19,\cdots$。

上述为预测各年第1季度销售量而建立的二阶差分方程，虽然其系数与前5年第1季度的统计数据完全吻合，但用于预测时预测值与事实不符。凭直觉，第6年估计值明显偏高，第七年销售量预测值甚至小于第6年。稍做分析不难看出，如分析对每一季度建立一个差分方程，则根据统计数据拟合出的系数可能会相差甚大，但对同一种商品，这种差异应当是微小的，故应给根据统计数据建立一个共用于各个季度的差分方程。为此，将季度编号为$t=1,2,\cdots ,20$,令$y_t=a_1{y}_{t-4}+a_2$或$y_t=a_1{y}_{t-4}+a_2{y}_{t-8}+a_3$等，利用全体数据来拟合，求拟合的最好的系数。以二阶差分方程为例，求$a_1,a_2,a_3$使得
$$Q(a_1,a_2,a_3)=\sum _{t=9}^{20}[y_t-(a_1{y}_{t-4}+a_2{y}_{t-8}+a_3){]}^2$$
达到最小，计算得$a_1=z(1)=0.8737,a_2=z(2)=0.1941,a_3=z(3)=0.6957$,故求得二阶差分方程
$$y_t=0.8737y_{t-4}+0.1941y_{t-8}+0.6957,$$
根据此式迭代，可求得第6年和第7年第1季度销售量的预测值为
$$y_{21}=17.5869,y_{25}=\mathrm{19.1676.}$$
还是较为可信的。

计算的Matlab程序如下:

clc;clear
a=readmatrix('销售数据.xlsx')
y=a(~isnan(a));
b=y(9:20);
A=[y(5:16),y(1:12),ones(12,1)];
z=pinv(A)*b
for t=21:25y(t)=[y(t-4),y(t-8),1]*z;
end
yhat=y([21,25]) %提取t=21,25时的预测值

例8.5 养老保险。某保险公司的一份材料指出，在每月交费200元至59岁年底，60岁开始领取养老金的约定下，男子若25岁投保，届时月养老金2282元；假定人的寿命为75岁，试求出保险公司为了兑现保险责任，每月至少应有多少投资收益率？

解 设$r$表示投保金的投资收益率，缴费期间月缴费额为$p$元，领养老金期间月领取额为$q$元，缴费的月数为$N$，到75岁时领取养老金的月数为$M$,投保人在投保后第$k$个月所交保险费及收益的累计总额为$F_k$,那么容易得到数学模型为分段表示的差分方程
$$\begin{gathered} F_{k+1}=F_k(1+r)+p,k=0,1,\cdots ,N-1, \\ {F}_{k+1}=F_k(1+r)-q,k=N,N+1,\cdots ,M-1, \end{gathered}$$
这里$p=200,q=2282,N=420,M=600$。

可推出差分方程的解(这里$F_0=F_M=0$)为
$$F_k=[(1+r{)}^k-1]\frac{p}{r},k=0,1,2,\cdots ,N,\text{\hspace{1em}(8.5)}$$
$$F_k=\frac{q}{r}[1-(1+r{)}^{k-M}],k=N+1,\cdots ,M,\text{\hspace{1em}(8.6)}$$
由式(8.5)和式(8.6)得
$$F_N=[(1+r{)}^N-1]\frac{p}{r},$$
$$F_{N+1}=\frac{q}{r}[1-(1+r{)}^{N+1-M}],$$
由于$F_{N+1}=F_N(1+r)-q$,因此得
$$\frac{q}{r}[1-(1+r{)}^{N+1-M}]=[(1+r{)}^N-1]\frac{p}{r}(1+r)-q,$$
化简得
$$(1+r{)}^M-(1+\frac{q}{p})(1+r{)}^{M-N}+\frac{q}{p}=0,$$
记$x=1+r$,代入数据得
$$x^{600}-12.41x^{180}+11.41=0.$$
利用Matlab程序，求得$x=1.0049$,因而投资收益率$r=0.49$.

计算的Matlab程序如下:

clc;clear
M=600;N=420;p=200;q=2282;
eq=@(x) x^M-(1+q/p)*x^(M-N)+q/p;
sol=fzero(eq,[1.0001,1.5])

8.2 差分方程的基本概念和理论

1．基本概念

定义 8.1 称形如
$$y_{n+k}+a_1{y}_{n+k-1}+a_2{y}_{n+k-2}+\cdots +a_k{y}_n=0\text{\hspace{1em}(8.7)}$$
的差分方程为 $k$ 阶常系数线性齐次差分方程，其中 $a_i(1⩽i⩽k)$ 为常数，$a_k\ne 0,k⩾1$ 。

定义 8.2 称方程
$${\lambda }^k+a_1{\lambda }^{k-1}+\cdots +a_{k-1}\lambda +a_k=0\text{\hspace{1em}(8.8)}$$
为差分方程(8.7)的特征方程，方程的根称为差分方程(8.7)的特征根。

定义 8.3 称形如
$$y_{n+k}+a_1{y}_{n+k-1}+a_2{y}_{n+k-2}+\cdots +a_k{y}_n=f(n)\text{\hspace{1em}(8.9)}$$
的差分方程为 $k$ 阶常系数线性非齐次差分方程，其中 $a_i(1⩽i⩽k)$ 为常数，$a_k\ne 0,k⩾1$ ， $f(n)\ne 0$ 。

定理 8.1 若 $k$ 阶差分方程(8.7)的特征方程(8.8)有 $k$ 个互异的特征根 ${\lambda }_1,{\lambda }_2,\cdots$ ， ${\lambda }_k$ ，则
$$y_n=c_1{\lambda }_1^n+c_2{\lambda }_2^n+\cdots +c_k{\lambda }_k^n\text{\hspace{1em}(8.10)}$$
是差分方程(8.7)的一个通解，其中 $c_1,c_2,\cdots ,c_k$ 为任意常数。
进一步地，若给定一组初始条件：
$$y_0=u_0,y_1=u_1,\cdots ,y_{k-1}=u_{k-1},$$
则利用待定系数法，可以确定差分方程满足初始条件的特解。

定理8.2 若 $k$ 阶差分方程(8.7)的特征方程(8.8)有 $t$ 个互异的特征根 ${\lambda }_1,{\lambda }_2,\cdots$ ， ${\lambda }_t$ ，重数依次为 $m_1,m_2,\cdots ,m_t$ ，其中 $m_1+m_2+\cdots +m_t=k$ ，则差分方程的通解为
$$y_n=\sum _{j=1}^{m_1}{c}_{1j}{n}^{j-1}{\lambda }_1^n+\sum _{j=1}^{m_2}{c}_{2j}{n}^{j-1}{\lambda }_2^n+\cdots +\sum _{j=1}^{m_t}{c}_{ij}{n}^{j-1}{\lambda }_t^n\text{\hspace{1em}(8.11)}$$
定理8.3 $k$ 阶常系数非齐次差分方程(8.9)的通解 $y_n$ 等于对应齐次差分方程的通解加上非齐次差分方程的特解，即
$$y_n=\widetilde{y_n}+y_n^{\ast },\text{\hspace{1em}(8.12)}$$
式中:$\widetilde{y_n}$为对应齐次差分方程的通解;$y_n^{\ast }$为对应非齐次差分方程的特解。

2.差分方程的平衡点及稳定性

1. 一阶线性方程的平衡点及稳定性**
   考虑一阶线性常系数差分方程，一般形式为：
   $$y_{k+1}+a{y}_k=b,k=0,1,2,\cdots \text{\hspace{1em}(8.13)}$$
   式中：$a,b$ 为常数。

称 $y^{\ast }$ 为方程 (8.13) 的平衡点，若满足 $y^{\ast }+a{y}^{\ast }=b$。求解得$y^{\ast }=\frac{b}{1+a}$。

当 $k\to \infty$ 时，若 $y_k\to y^{\ast }$，则平衡点 $y^{\ast }$ 是稳定的；否则 $y^{\ast }$ 是不稳定的。

为了理解平衡点的稳定性，可以用变量代换将方程 (8.13) 的稳定性问题转换为：
$$y_{k+1}+a{y}_k=0,k=0,1,2,\cdots \text{\hspace{1em}(8.14)}$$
的平衡点 $y^{\ast }=0$ 的稳定性问题。而对于方程 (8.14) 的解可由递推公式直接给出：
$$y_k=(-a{)}^k{y}_0,k=1,2,3,\cdots$$
所以当且仅当 $|a|<1$ 时，方程(8.14)的平衡点 $y^{\ast }=0$ 稳定（从而方程(8.13)的平衡点稳定）。

2. 一阶线性常系数差分方程组的平衡点及稳定性

对于 $n$ 维向量 $\mathbf{y}(k)$ 和 $n\times n$ 常数矩阵 $\mathbf{A}$ 构成的一阶线性常系数齐次差分方程组：
$$\mathbf{y}(k+1)+\mathbf{A}\mathbf{y}(k)=0,k=0,1,2,\cdots \text{\hspace{1em}(8.15)}$$
其平衡点 ${\mathbf{y}}^{\ast }=0$ 稳定的条件是：
$\mathbf{A}$ 的所有特征根均满足 $|{\lambda }_i|<1(i=1,2,\cdots ,n)$，即特征根均在复平面的单位圆内）。

对于 $n$ 维向量 $\mathbf{y}(k)$ 和 $n\times n$ 常数矩阵 $\mathbf{A}$ 构成的一阶线性常系数非齐次差分方程组：
$$\mathbf{y}(k+1)+\mathbf{A}\mathbf{y}(k)=\mathbf{B},k=0,1,2,\cdots \text{\hspace{1em}(8.16)}$$
其平衡点为$\boldsymbol{y}^* = (\boldsymbol{E} + \boldsymbol{A})^{-1} \boldsymbol{B}
$，其中 $\mathbf{E}$ 为 $n$ 阶单位方阵。其稳定性条件与齐次方程(8.15)相同，即$\mathbf{A}$ 的所有特征根均满足 $|{\lambda }_i|<1(i=1,2,\cdots ,n)$。

3. 二阶线性常系数差分方程的平衡点及稳定性

考察二阶线性常系数齐次差分方程：
$$y_{k+2}+a_1{y}_{k+1}+a_2{y}_k=0\text{\hspace{1em}(8.17)}$$
的平衡点($y^{\ast }=0$)的稳定性。

方程(8.17)的特征方程为
$${\lambda }^2+a_1\lambda +a_2=0,$$
记它的特征根为${\lambda }_1,{\lambda }_2$,则方程(8.17)的通解可表示为
$$y_k=c_1{\lambda }_1^k+c_2{\lambda }_2^k,\text{\hspace{1em}(8.18)}$$
式中：$c_1,c_2$为待定常数，由两个初始条件的值确定。

由式(8.18)很容易就可以得到，当且仅当
$$\mid {\lambda }_1\mid <1,\mid {\lambda }_2\mid <1$$
时，方程(8.17)的平衡点才是稳定的。

与一阶线性方程一样，非齐次方程
$$y_{k+2}+a_1{y}_{k+1}+a_2{y}_k=b$$
的平衡点的稳定性条件和方程(8.17)相同。

上述结果可以推广到$n$阶线性常系数差分方程的平衡点及其稳定性问题。即平衡点稳定的充要条件是其特征方程的根$|{\lambda }_i|<1(i=1,2,\cdots ,n)$。

4. 一阶非线性差分方程
   考察一阶非线性差分方程
   $$y_k=f(y_k)\text{\hspace{1em}(8.19)}$$
   式中：$f(y_k)$ 为已知函数。其平衡点 $y^{\ast }$ 由代数方程 $y^{\ast }=f(y^{\ast })$ 解出。

现分析 $y^{\ast }$ 的稳定性。将方程的右端在 $y^{\ast }$ 点作泰勒多项式展开，只取一阶导数项，则式 (8.19) 可近似为
$$y_{k+1}\approx f(y^{\ast })+f^{\prime }(y^{\ast })(y_k-y^{\ast })\text{\hspace{1em}(8.20)}$$
故 $y^{\ast }$ 也是近似齐次线性差分方程 (8.20) 的平衡点。从而由一阶齐次线性差分方程的平衡点稳定性理论可知，$y^{\ast }$ 稳定的充要条件为 $|f^{\prime }(y^{\ast })|<1$。

8.3 莱斯利(Leslie)种群模型

莱斯利模型是研究动物种群数量增长的重要模型，这一模型研究了种群中雌性动物的年龄分布和数量增长的规律。

在某动物种群中，仅考察雌性动物的年龄和数量。设雌性动物的最大生存年龄为 $L$ (单位:年或其他时间单位)，把 $[0,L]$ 等分为 $n$ 个年龄组，每一年龄组的长度为 $L/n$，$n$ 个年龄组分别为
$$\left[0,\frac{L}{n}\right),\left[\frac{L}{n},\frac{2L}{n}\right),\cdots ,\left[\frac{(n-1)L}{n},L\right].$$

设第 $i$ 个年龄组的生育率为 $a_i$，存活率为 $b_i(i=1,2,\cdots ,n)$，$a_i,b_i$ 均为常数，且 $a_i\ge 0$ $(i=1,2,\cdots ,n)$，$0<b_i\le 1(i=1,2,\cdots ,n-1)$。同时，设至少有一个 $a_i>0(1\le i\le n)$，即至少有一个年龄组的雌性动物具有生育能力。

利用统计资料可获得基年 ($t=0$) 该种群在各年龄组的雌性动物数量，记为 $x_i^{(0)}(i=1,2,\cdots ,n)$，为 $t=0$ 时第 $i$ 年龄组雌性动物的数量，就得到初始时刻种群数量分布向量：
$$x^{(0)}=[x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)}{]}^{\text{T}}.$$

如果以年龄组的间隔 $\frac{L}{n}$ 作为时间单位，记 $t_1=\frac{L}{n},t_2=\frac{2L}{n},\cdots ,t_k=\frac{kL}{n},\cdots$，在 $t_k$ 时第 $i$ 年龄组雌性动物的数量为$x_i^{(k)}(i=1,2,\cdots ,n),t_k$时各年龄组种群数量分布向量为
$$x^{(k)}=[x_1^{(k)},x_2^{(k)},\cdots ,x_n^{(k)}{]}^T,k=0,1,2,\cdots .\text{\hspace{1em}(8.21)}$$
随着时间的变化，由于出生、死亡以及年龄的增长，该种群中每一个年龄组的雌性动物数量都将发生变化。实际上，在 $t_k$ 时刻，种群中第 1 个年龄组的雌性动物数量应等于在 $t_{k-1}$ 和 $t_k$ 之间出生的所有雌性幼体的总和，即：
$$x_1^{(k)}=a_1{x}_1^{(k-1)}+a_2{x}_2^{(k-1)}+\cdots +a_n{x}_n^{(k-1)}\text{\hspace{1em}(8.22)}$$
同时，在 $t_k$ 时刻，第 $i+1$ 个年龄组（$i=1,2,\cdots ,n-1$）中雌性动物的数量应等于在 $t_{k-1}$ 时刻第 $i$ 个年龄组中雌性动物数量 $x_i^{(k-1)}$ 乘存活率 $b_i$，即：
$$x_{i+1}^{(k)}=b_i{x}_i^{(k-1)}(i=1,2,\cdots ,n-1)\text{\hspace{1em}(8.23)}$$

综合上述分析，由式 (8.22) 和式 (8.23) 可得到在 $t_k$ 和 $t_{k-1}$ 时各年龄组中雌性动物数量间的关系：

$$\{\begin{array}{l}{x}_1^{(k)}=a_1{x}_1^{(k-1)}+a_2{x}_2^{(k-1)}+\cdots +a_n{x}_n^{(k-1)},\\ x_2^{(k)}=b_1{x}_1^{(k-1)},\\ x_3^{(k)}=b_2{x}_2^{(k-1)},\\ ⋮\\ x_n^{(k)}=b_{n-1}{x}_{n-1}^{(k-1)}.\end{array}\text{\hspace{1em}(8.24)}$$

记矩阵

$$\mathbf{L}=\left[\begin{array}{ccccc}{a}_1& a_2& \cdots & a_{n-1}& a_n\\ b_1& 0& \cdots & 0& 0\\ 0& b_2& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & b_{n-1}& 0\end{array}\right]$$

则式 (8.24) 可写成

$${\mathbf{x}}^{(k)}=\mathbf{L}{\mathbf{x}}^{(k-1)},k=1,2,3,\cdots \text{\hspace{1em}(8.25)}$$

式中：$\mathbf{L}$ 为莱斯利矩阵。

由式 (8.25) 可得 ${\mathbf{x}}^{(1)}=\mathbf{L}{\mathbf{x}}^{(0)},{\mathbf{x}}^{(2)}=\mathbf{L}{\mathbf{x}}^{(1)}={\mathbf{L}}^2{\mathbf{x}}^{(0)},\cdots$，一般有

$${\mathbf{x}}^{(k)}=\mathbf{L}{\mathbf{x}}^{(k-1)}={\mathbf{L}}^k{\mathbf{x}}^{(0)},k=1,2,3,\cdots$$

若已知初始时种群数量分布向量 ${\mathbf{x}}^{(0)}$，则可以推算任一时刻 $t_k$ 该种群数量分布向量，并以此对该种群的总量进行科学的分析。

例 8.6 某种雌性动物的最大生存年龄为 15 年，以 5 年为一间隔，把这一动物种群分为三个年龄组 $[0,5),[5,10),[10,15]$，利用统计资料，已知 $a_1=0,a_2=4,a_3=3;b_1=0.5,b_2=0.25$。在初始时刻 $t=0$ 时，三个年龄组的雌性动物个数分别为 500, 1000, 500，则初始种群数量分布向量和莱斯利矩阵分别为
$${\mathbf{x}}^{(0)}=[500,1000,500{]}^{\text{T}},\mathbf{L}=\left[\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\right]$$

于是

$${\mathbf{x}}^{(1)}=\mathbf{L}{\mathbf{x}}^{(0)}=\left[\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\right]\left[\begin{array}{c}500\\ 1000\\ 500\end{array}\right]=\left[\begin{array}{c}5500\\ 250\\ 250\end{array}\right]$$

$${\mathbf{x}}^{(2)}=\mathbf{L}{\mathbf{x}}^{(1)}=\left[\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\right]\left[\begin{array}{c}5500\\ 250\\ 250\end{array}\right]=\left[\begin{array}{c}1750\\ 2750\\ 62.5\end{array}\right]$$

$${\mathbf{x}}^{(3)}=\mathbf{L}{\mathbf{x}}^{(2)}=\left[\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\right]\left[\begin{array}{c}1750\\ 2750\\ 62.5\end{array}\right]=\left[\begin{array}{c}11187.5\\ 875\\ 687.5\end{array}\right]$$
为了分析当 $k\to \infty$ 时，该动物种群数量分布向量的特点，先求出矩阵 $\mathbf{L}$ 的特征值与特征向量，为此计算 $\mathbf{L}$ 的特征多项式：

$$|\lambda \mathbf{E}-\mathbf{L}|=\left|\begin{array}{ccc}\lambda & -4& -3\\ -0.5& \lambda & 0\\ 0& -0.25& \lambda \end{array}\right|=\left(\lambda -\frac{3}{2}\right)\left({\lambda }^2+\frac{3}{2}\lambda +\frac{1}{4}\right),$$

可得 $\mathbf{L}$ 的特征值 ${\lambda }_1=\frac{3}{2},{\lambda }_2=\frac{-3+\sqrt{5}}{4},{\lambda }_3=\frac{-3-\sqrt{5}}{4}$，不难看出 ${\lambda }_1$ 是矩阵 $\mathbf{L}$ 的唯一正特征值，且 $|{\lambda }_1|>|{\lambda }_2|,|{\lambda }_1|>|{\lambda }_3|$。$\mathbf{L}$ 有 3 个互异特征值，因此矩阵 $\mathbf{L}$ 可相似对角化。

设矩阵 $\mathbf{L}$ 属于特征值 ${\lambda }_i(i=1,2,3)$ 的特征向量为 ${\mathbf{\alpha }}_i$。不难计算 $\mathbf{L}$ 属于特征值 ${\lambda }_1=\frac{3}{2}$ 的特征向量为 ${\mathbf{\alpha }}_1=[18,6,1{]}^{\text{T}}$，记矩阵 $\mathbf{P}=[{\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2,{\mathbf{\alpha }}_3],\mathbf{\Lambda }=\text{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$，则

$${\mathbf{P}}^{-1}\mathbf{L}\mathbf{P}=\mathbf{\Lambda }\text{或}\mathbf{L}=\mathbf{P}\mathbf{\Lambda }{\mathbf{P}}^{-1}$$

${\mathbf{L}}^k=\mathbf{P}{\mathbf{\Lambda }}^k{\mathbf{P}}^{-1}$，于是

$${\mathbf{x}}^{(k)}={\mathbf{L}}^k{\mathbf{x}}^{(0)}=\mathbf{P}{\mathbf{\Lambda }}^k{\mathbf{P}}^{-1}{\mathbf{x}}^{(0)}={\lambda }_1^k\mathbf{P}\left[\begin{array}{ccc}1& 0& 0\\ 0& {\left({\lambda }_2/{\lambda }_1\right)}^k& 0\\ 0& 0& {\left({\lambda }_3/{\lambda }_1\right)}^k\end{array}\right]{\mathbf{P}}^{-1}{\mathbf{x}}^{(0)}$$

即

$$\frac{1}{{\lambda }_1^k}{\mathbf{x}}^{(k)}=\mathbf{P}\left[\begin{array}{ccc}1& 0& 0\\ 0& {\left({\lambda }_2/{\lambda }_1\right)}^k& 0\\ 0& 0& {\left({\lambda }_3/{\lambda }_1\right)}^k\end{array}\right]{\mathbf{P}}^{-1}{\mathbf{x}}^{(0)}$$

因为 $\left|\frac{{\lambda }_2}{{\lambda }_1}\right|<1,\left|\frac{{\lambda }_3}{{\lambda }_1}\right|<1$，所以

$$\underset{k\to \infty }{\lim }\frac{1}{{\lambda }_1^k}{\mathbf{x}}^{(k)}=\mathbf{P}\text{diag}(1,0,0){\mathbf{P}}^{-1}{\mathbf{x}}^{(0)}\text{\hspace{1em}(8.26)}$$

记列向量 ${\mathbf{P}}^{-1}{\mathbf{x}}^{(0)}$ 的第一个元素为 $c$（常数），则式 (8.26) 可化为
$$\underset{k\to \infty }{\lim }\frac{1}{{\lambda }_1^k}{\mathbf{x}}^{(k)}=[{\mathbf{\alpha }}_1,{\mathbf{\alpha }}_2,{\mathbf{\alpha }}_3]\left[\begin{array}{c}c\\ 0\\ 0\end{array}\right]=c{\mathbf{\alpha }}_1$$

于是，当 $k$ 充分大时，近似地有：

$${\mathbf{x}}^{(k)}=c{\lambda }_1^k{\mathbf{\alpha }}_1=c{\left(\frac{3}{2}\right)}^k\left[\begin{array}{c}18\\ 6\\ 1\end{array}\right]$$

其中 $c=\frac{2250}{19}$。

这一结果说明，当时间充分长，这种动物中雌性的年龄分布将趋于稳定，即 3 个年龄组的数量比为 $18:6:1$。并由此可近似得到在 $t_k$ 时刻种群中雌性动物的总量，从而对整个种群的总量进行估计。

莱斯利模型在分析动物种群的年龄分布和总量增长方面有广泛应用，这一模型也可应用于人口增长的年龄分布问题。

计算的Matlab程序如下：

clc;clear
syms k positive integer
X0=[500;1000;500];
L=[0,4,3;0.5,0,0;0,0.25,0];
L=sym(L);   %转为符号矩阵
X1=L*X0,X2=L*X1,X3=L*X2
p=charpoly(L)   %计算特征多项式
r=roots(p)    %计算符号特征值
[P,D]=eig(L)   %计算符号特征向量和特征值
XL=P*diag([1,0,0])*inv(P)*X0
tc=inv(P)*X0;c=tc(1)

参考文献
[1] 司守奎,孙玺青 数学建模算法与应用第3版[M]