数模竞赛-线性代数模型

特征值&特征向量

差分方程-Fibonacci数列

化为一阶差分方程组：

$$\{\begin{array}{l}\begin{array}{c}{F}_{k+1}=F_{k+1},\\ F_{k+2}=F_{k+1}+F_k,\end{array}\end{array}k=0,1,2,\cdots$$

写成矩阵形式：

$${\alpha }_{k+1}=A{\alpha }_k$$

其中：

$$A=\left[\begin{array}{c}\begin{array}{cc}0& 1\\ 1& 1\end{array}\end{array}\right],{\alpha }_k=\left[\begin{array}{c}\begin{array}{c}{F}_k\\ F_{k+1}\end{array}\end{array}\right],{\alpha }_0=\left[\begin{array}{c}\begin{array}{c}1\\ 1\end{array}\end{array}\right]$$

由递推式得：

$${\alpha }_k=A^k{\alpha }_0,k=1,2,3,\cdots$$

设特征值为${\lambda }_1,{\lambda }_2$，矩阵$P$为特征向量矩阵，即：

$$A^k=P\left[\begin{array}{c}\begin{array}{cc}{\lambda }_1^k& 0\\ 0& {\lambda }_2^k\end{array}\end{array}\right]P^{-1}$$

则：

$$a_k=A^k\left[\begin{array}{c}\begin{array}{c}1\\ 1\end{array}\end{array}\right]$$

代码：

import sympy
sympy.var('k',positive=True,integer=True)
a=sympy.Matrix([[0,1],[1,1]])
val=a.eigenvals() #求特征值 无用
vec=a.eigenvects() #求特征向量 无用
P,D=a.diagonalize() #相似对角化
ak=P@(D**k)@(P.inv())
F=ak@sympy.Matrix([1,1])
s=sympy.simplify(F[0])
print(s) #输出递推式
sm=[]
for i in range(20):
    sm.append(s.subs(k,i).n())
print(sm)

求特征根：

$$|\lambda E-A|={\lambda }^2-\lambda -1$$

俩根互异，特征根可知，通解为：

$$F_k=c_1{\lambda }_1^k+c_2{\lambda }_2^k$$

求$c_1,c_2$，即：

$$\{\begin{array}{l}\begin{array}{c}{c}_1+c_2=1,\\ c_1{\lambda }_1+c_2{\lambda }_2=1\end{array}\end{array}$$

代码：

import sympy
sympy.var('t c1 c2')
t0=sympy.solve(t**2-t-1)
eq1=c1+c2-1
eq2=c1*t0[0]+c2*t0[1]-1
s=sympy.solve([eq1,eq2])
print(s[c1],s[c2])

或者直接尝试求解通项公式：

import sympy
sympy.var('k')
y=sympy.Function('y')
f=y(k+2)-y(k+1)-y(k)
s=sympy.rsolve(f,y(k),{y(0):1,y(1):1})
print(s)

Leslie种群模型

在某动物种群中，仅考查雌性动物的年龄和数量，设雌性动物的最大生存年龄为$L$，把$[0,L]$等分为$n$个年龄组：

$$[0,{\displaystyle \frac{L}{n}}),[{\displaystyle \frac{L}{n}},{\displaystyle \frac{2L}{n}}),\cdots ,[{\displaystyle \frac{n-1}{n}}L,L]$$

设第$i$个年龄组的生育率为$a_i$，存活率为$b_i$，记$x_i^0$为$t=0$时第$i$个年龄组雌性动物的数量，得到初始时刻种群数量分布向量：

$$x^{(0)}={[x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)}]}^{\mathrm{T}}$$

则第$1$个年龄组为：

$$x_1^{(k)}=\sum _{i=1}^n{a}_i{x}_i^{(k-1)}$$

其他年龄组则生存率：

$$x_{i+1}^{(k)}=b_i{x}_i^{(k-1)}$$

记Leslie矩阵：

$$L=\left[\begin{array}{c}\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\\ ⋮& ⋮& & ⋮& ⋮\\ 0& 0& \cdots & b_{n-1}& 0\end{array}\end{array}\right]$$

一般有$x^{(k)}=L{x}^{(k-1)}=L^k{x}^{(0)}$。

例题：某动物雌性最大年龄为$15$年，$5$年为一间隔，分成$3$个年龄组，已知$a_1=0,a_2=4,a_3=3,b_1=0.5,b_2=0.25$，在$t=0$时雌性动物数量为$500,1000,500$。

即：

$$x^{(0)}={[500,1000,500]}^{\mathrm{T}},L=\left[\begin{array}{c}\begin{array}{ccc}0& 4& 3\\ 0.5& 0& 0\\ 0& 0.25& 0\end{array}\end{array}\right]$$

于是有$x^{(i)}=L{x}^{(i-1)}$，求特征值$\mid \lambda E-L\mid =0$得${\lambda }_1,{\lambda }_2,{\lambda }_3$，对应特征向量为${\alpha }_1,{\alpha }_2,{\alpha }_3$，记矩阵$P=[{\alpha }_1,{\alpha }_2,{\alpha }_3],\mathrm{\Lambda }=\mathrm{diag}({\lambda }_1,{\lambda }_2,{\lambda }_3)$，则$L=P\mathrm{\Lambda }{P}^{-1}$，于是有$x^{(k)}=L^k{x}^{(0)}$，即

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

因${\displaystyle \left|{\displaystyle \frac{{\lambda }_2}{{\lambda }_1}}\right|<1,\left|{\displaystyle \frac{{\lambda }_3}{{\lambda }_1}}\right|<1}$，有：

$$\underset{k\to \mathrm{\infty }}{lim}{\displaystyle \frac{1}{{\lambda }_1^k}}=P\mathrm{diag}(1,0,0)P^{-1}{x}^{(0)}$$

于是当$k$充分大，近似成立：

$$x^{(k)}={\displaystyle \frac{2250}{19}}{\left({\displaystyle \frac{3}{2}}\right)}^k\left[\begin{array}{c}\begin{array}{c}18\\ 6\\ 1\end{array}\end{array}\right]$$

其中$c={\displaystyle \frac{2250}{19}}$为列向量$P^{-1}{x}^{(0)}$的第一个元素。

代码：

import numpy,sympy
X0=numpy.array([500,1000,500])
L=numpy.array([[0,4,3],[0.5,0,0],[0,0.25,0]])
X1=L@X0
X2=L@X1
X3=L@X2
Ls=sympy.Matrix([[0,4,3],[sympy.Rational(1,2),0,0],[0,sympy.Rational(1,4),0]])
sympy.var('lmbda')
p=Ls.charpoly(lmbda) #计算特征多项式
w1=sympy.roots(p) #计算特征值
w2=Ls.eigenvals() #直接计算特征值
v=Ls.eigenvects() #特征向量
print(w2,v)
P,D=Ls.diagonalize() #相似对角化
Pinv=P.inv()
Pinv=sympy.simplify(Pinv)
cc=Pinv@X0
print(P,cc[0])

奇异值分解

有$A=(a_{ij}{)}_{m\times n},\mathrm{rank}(A)=r$，存在正交矩阵$U,V$，使得：

$$U^{\mathrm{T}}AV=\left[\begin{array}{c}\begin{array}{cc}\mathrm{\Sigma }& 0\\ 0& 0\end{array}\end{array}\right]$$

其中$\mathrm{\Sigma }=\mathrm{diag}{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_r,{\sigma }_1⩾{\sigma }_2⩾\cdots ⩾{\sigma }_r>0$，且${\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_r^2$为矩阵$A^{\mathrm{T}}A$的正特征值。

称下式为$A$的奇异值分解，${\sigma }_i,i=1,2,\cdots ,r$为$A$的奇异值。

$$A=U\left[\begin{array}{c}\begin{array}{cc}\mathrm{\Sigma }& 0\\ 0& 0\end{array}\end{array}\right]V^{\mathrm{T}}$$

矩阵的F范数定义为：

$${\Vert \begin{array}{c}A\end{array}\Vert }_{\mathrm{F}}^2=\mathrm{tr}\left(A^{\mathrm{T}}A\right)=\sum _{i=1}^m\sum _{j=1}^n{a}_{ij}^2=\sum _{i=1}^r{\sigma }_i^2$$

例题1：

求矩阵$A$的奇异值分解。

$$A=\left[\begin{array}{c}\begin{array}{ccc}1& 0& 1\\ 0& 1& 1\\ 0& 0& 0\end{array}\end{array}\right]$$

代码：

import numpy
a=numpy.array([[1,0,1],[0,1,1],[0,0,0]])
u,s,vt=numpy.linalg.svd(a) #即a=u@numpy.diag(s)@vt
print(u,s,vt)