最优化方法Python计算：收敛序列极限计算

数值化地计算向量序列的极限之前，先思考如下结论
定理1 无穷序列 { x k } ⊂ R n \\{\\boldsymbol{x}_k\\}\\subset\\text{ℝ}^n {xk​}⊂Rn收敛，当且仅当对任意$\epsilon >0$存在$N\in \text{N}$，对所有$n,m>N$，
∥ x n − x m ∥ < ε . \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon. ∥xn​−xm​∥<ε.
利用本定理，可数值化地计算收敛序列 { x k } \\{\\boldsymbol{x}_k\\} {xk​}的极限：给定容错误差$\epsilon >0$，譬如取 ε = 1 0 − 6 \\varepsilon=10^{-6} ε=10−6。从一个适当的$N$值起步，譬如取 N = 1 0 3 N=10^3 N=103。从大于$N$的整数任取$n$和$m$，检测是否 ∥ x n − x m ∥ < ε \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon ∥xn​−xm​∥<ε。若是，则取$n,m$中的较大者为$k$，并以 x k \\boldsymbol{x}_k xk​作为序列极限的近似值返回并停机。否则，扩大$N$。譬如，扩大10倍，重复上述的过程。直至检测到有$n,m>N$，使得 ∥ x n − x m ∥ < ε \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon ∥xn​−xm​∥<ε。然而，这一过程不能因为没有检测到 ∥ x n − x m ∥ < ε \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon ∥xn​−xm​∥<ε而一直持续下去，所以必须设置一个最大迭代次数$i$，譬如说$i=10$（此时， N = 1 0 13 N=10^{13} N=1013）。当$i=10$时仍未检测到 ∥ x n − x m ∥ < ε \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon ∥xn​−xm​∥<ε则说明序列收敛速度非常低或根本就是一个发散序列，应停机并宣告计算无效。在Python中我们将此计算过程实现为以下的函数。

import numpy as np                                  #导入numpy
from math import nan                                #导入非数值常量nan
from scipy.stats import randint                     #导入整数均匀分布randint
def limit(x,eps=1e-6):                              #x为序列通项，eps为容错误差N=103                                         #初始步长N=1000m,n=N+randint.rvs(1,102,size=2)               #m和n大于N的随机数i=1                                             #迭代次数while i<=10 and np.linalg.norm(x(n)-x(m))>=eps: #迭代N*=10                                       #下一个步长m,n=N+randint.rvs(1,102,size=2)           #m和n大于N的随机数i+=1                                        #迭代次数自增1if i>10:                                        #N>1000*10^10return nanreturn x(max(m,n))

程序的第4~14定义函数limit。该函数有两个参数：x表示序列通项 x k \\boldsymbol{x}_k xk​，命名参数eps表示容错误差$\epsilon$，缺省值为 1 0 − 6 10^{-6} 10−6。
第5～7行执行初始化操作：第5行设置初始步长 N = 1 0 3 N=10^3 N=103为N。第6行调用randint函数（地3行导入），产生2个介于1~100之间的随机整数，加上N后赋予m和n。第7行迭代次数i初始化为1。
第8～11行的while循环进行迭代：第9行将步长N增大10倍。第10行产生2个1~100之间的整数，加上N后赋予m，n，作为大于N的序列下标。第11行将迭代次数自增1。循环往复，直至
i > 10 或 ∥ x n − x m ∥ < ε i>10\\text{或}\\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon i>10或∥xn​−xm​∥<ε
为真。
循环结束后，若在12行测得i>10，意味着 N = 1 0 13 N=10^{13} N=1013以后的元素差的范数仍大于容错误差$\epsilon$，故认为该序列不收敛，第13行返回math模块中表示非数值常量nan。否则，即
∥ x n − x m ∥ < ε \\lVert\\boldsymbol{x}_n-\\boldsymbol{x}_m\\rVert<\\varepsilon ∥xn​−xm​∥<ε
根据定理1判断序列 { x k } \\{\\boldsymbol{x}_k\\} {xk​}收敛，第14行将较大下标元素作为极限的近似值返回。
例1 考虑序列
(1) { sin ⁡ k k } \\{\\frac{\\sin{k}}{k}\\} {ksink​}；
(2) { ( k − 1 k k + 1 k ) } \\left\\{\\begin{pmatrix} \\frac{k-1}{k}\\\\\\frac{k+1}{k} \\end{pmatrix}\\right\\} {(kk−1​kk+1​​)}；
(3) { k } \\{k\\} {k}。
用上述程序中定义的函数limit判断其是否收敛，若是则计算极限的近似值。
解：(1)我们知道$\sin k$是一个有界量，而$k$是一个无穷大量，两者之比为无穷小量，即 lim ⁡ k → ∞ sin ⁡ k k = 0 \\lim\\limits_{k\\rightarrow\\infty}\\frac{\\sin{k}}{k}=0 k→∞lim​ksink​=0；
(2)根据例1.7的计算知 lim ⁡ k → ∞ ( k − 1 k k + 1 k ) = ( 1 1 ) \\lim\\limits_{k\\rightarrow\\infty}\\begin{pmatrix} \\frac{k-1}{k}\\\\\\frac{k+1}{k} \\end{pmatrix}=\\begin{pmatrix} 1\\\\1\\end{pmatrix} k→∞lim​(kk−1​kk+1​​)=(11​)；
(3)序列 { k } \\{k\\} {k}是一个发散的序列。
下列代码验证以上结果。

import numpy as np							#导入numpy
x=lambda k:np.array([np.sin(k)/k])			#序列(1)的通项
print('%.1f'%limit(x))						#计算极限
x=lambda k:np.array([(k-1)/k,(k+1)/k])		#序列(2)的通项
print(limit(x))								#计算极限
x=lambda k:k								#序列(3)的通项
print(limit(x))								#计算极限

程序的2、4、6行分别用lambda运算符定义序列(1)、(2)和(3)的通项。第3、5、7行调用limit函数分别计算各序列极限。运行程序，输出

0.0
[0.9999006 1.0000994]
nan

第1行表示序列(1)收敛于0，第2行表示序列(2)的极限向量 ( 1 1 ) \\begin{pmatrix} 1\\\\1\\end{pmatrix} (11​)的近似值，最后一行表示序列(3)无极限。