Re0微积分1：用牛顿迭代法求算术平方根

问题0：牛顿法求平方根,数学分析(zorich)3.2.7习题

假设$a>0$，定义序列:$x_{n+1}:=\frac{1}{2}\left(x_n+\frac{a}{x_n}\right)$

1. 证明在$x_1>0$的条件下$x_n\to \sqrt{a}$。
2. 估计绝对误差$e_n:=x_n-\sqrt{a}$,给出$|e_n|$的上界估计。

1. 一些基本的概念

• 序列：实际上就是一列有顺序的数。当我们定义一个序列的时候，我们定义了一个函数$a:\mathbb{N}\to \mathbb{C}$然后函数在自然数$n$处的取值为$a(n)=a_n$。 比如说常见的等差序列$1,3,5,\mathrm{7...}$,再比如$1,1/2,1/3,...$等等。
• 序列的极限：在题目中符号$x_n\to \sqrt{a}$，也就是说序列的极限为$\sqrt{a}$

在本问题的背景下，所谓序列的极限，就是与此序列有关的一个实数。当然问题要求我们去证明，这个实数就是$\sqrt{a}$。序列极限的严格意义的定义为:

极限的定义

对于实序列$x_n$,如果存在一个实数$L\in \mathbb{R}$满足，对任意的$\epsilon >0$都存在一个对应的正整数$N_{\epsilon }$使得任意$n>N_{\epsilon }$都有$|x_n-L|<\epsilon$ 那么就称$L$是序列$x_n$的极限，也可以表示为${\lim }_{n\to +\infty }{x}_n=L$

1.1 为什么需要引入这样的概念

这个分析学(微积分的推广)当中的一个基本想法有关。

分析学的一个基本想法

当我们想要研究一个复杂的对象$\text{A}$的某个性质$P$的时候，不妨想个办法找到某个由一些我们已知的并且认为相对简单的对象构成的集合$\Omega$,使得$A$可以由$\Omega$当中的元素任意逼近。那么我们可以从研究$A$的性质$P$,转而研究$\Omega$当中元素的性质$P$。

在具体的问题当中，我们或许会需要一个装置$d$去量化两个元素的性质$P$的相近的程度(假设性质$P$分别为$P_A,P_B$)，于是对象$A$可以被集合$\Omega$当中的元素任意逼近，就可以定量地表述为:$\forall \epsilon >0,\exists B\in \Omega \text{ s.t. }d(P_A,P_B)<\epsilon$ 这样的模式在分析学当中极其常见：

例子1.1.1：割圆术

假设我们现在有一个对象$S^1$为单位圆周，我们想要研究$S^1$的称之为周长的性质，即一个数量$a$.假设我们已经知道，我们可以用单位圆周的内接正多变形去任意逼近单位圆周。

1. 此处定义圆内接正多边形的集合为$\Omega =\{p_1,...,p_n,...\}$.
2. 这里的性质$P$,就是图形的周长。单位圆内接正$n$边形$p_n$的周长为$L_n=2n\sin \left(\frac{\pi }{n}\right)$
3. 此处衡量两个对象周长的相近的程度的装置自然就是实数差的绝对值，因为周长本质上还是实数。于是我们相当于得到了$\forall \epsilon >0,\exists p_n\in \Omega \text{ s.t. }|L_n-a|<\epsilon$

例子1.1.1当中，我们得到了对单位圆周的任意逼近，我们得到了什么好处？任意给定误差$\epsilon$的范围内，我可以知道$a\in (L_n-\epsilon ,L_n+\epsilon )$里面，这实际上就是古代数学家探索几百年想要搞清楚的事情，我们如何计算单位圆周的一半，即$\pi =a/2$.这里我们只要取合适的$\epsilon$以及足够大的$n$，我们就可以说:$\pi \approx 3.14$ 当然在数学里面，相应的例子不止圆的周长，还有面积，体积，方程的根等等…

不过在实际应用这个想法来达成分析学的目的还需要明确一些事情：

在分析学当中实施这个想法需要搞清楚的三件事

1. 如何找到这样一个在性质$P$的意义下，无限逼近对象$A$的集合$\Omega$?
2. 我们如何判断某个集合里面的元素会任意逼近某个对象? (收敛性的判断)
3. 如果我们的这种无限逼近的选项不止一个，例如有两个集合${\Omega }_1,{\Omega }_2$当中的元素都可以任意逼近$A$，如何对比二者逼近的品质？(误差的估计)

在本文最终要解决的问题当中，这三点皆有体现：

1. 我们想要搞清楚算术平方根$\sqrt{a}$的值(通常a是有理数，甚至正整数)，我们想到可以用有理数逼近这个值。即找到一个有理数序列$x_n$去任意逼近$\sqrt{a}$,找到这样的逼近序列的办法就是牛顿迭代法，也就是题目中定义的$x_{n+1}:=\frac{1}{2}\left(x_n+\frac{a}{x_n}\right)$对于这个迭代，理想状态下，我们只需要带入一个非零有理数$x_1$然后不断迭代，我们就可以得到$\sqrt{a}$的近似值。
2. 凭什么就说$x_n$一定要会任意逼近$\sqrt{a}$?因为$x_n$的极限就是$\sqrt{a}$,这就是为何要引入极限的概念，以及导出一些极限存在性的命题。
3. 这个逼近好不好，具体来说，给定初值$x_1$我到底需要最多几次能在给定精度下逼近$\sqrt{a}$?如何定量研究这件事？在本问题当中，我们就是要研究序列$e_n=|x_n-\sqrt{a}|$.

1.2 极限的概念的几何直观

首先在此明确，极限如果存在，那么这是一个具体的实数。比如当我说，$x_n=1/n$的极限为$0$,的时候0就是这个实数。然后这个实数对序列有特殊的意义，即满足极限定义当中的那个命题。

极限定义当中的这个命题可以有这样一个几何意义的理解：

极限的几何意义理解

任意以$L$为中心半径为$\epsilon$的实数轴上的对称开区间$I_{\epsilon }(L)=(L-\epsilon ,L+\epsilon )$的外部只有有限的，比方说不超过$N_{\epsilon }$个序列$x_n$的元素。即除开这$N_{\epsilon }$个序列的元素以外的序列剩下的无穷多个元素，全部在$I_{\epsilon }(L)$当中.

例子1.2.1

假设$y_n=\frac{(-1{)}^n}{n}$是我们要讨论极限的序列，我们遍寻所有实数轴上点(实际上在实数轴上中一目了然)，发现只有0能满足这样的几何意义。比如我们令${\epsilon }_0=1/3$,于是我们的对称区间就是$I_{{\epsilon }_0}(0)=(-1/3,1/3)$当中包含无穷多个序列的点，然而在区间外部，只有$-1,1/2,-1/3$三个序列的点。最重要的是，这个点的特殊之处在于，无论我们对称区间的半径$\epsilon$取多少，这个几何性质皆成立。

上面的例子之所以说是几何直观，我们大可以尝试把序列的点算出来，然后标注在实数轴上，这样图中会非常明显看到一个符合这样的几何性质的点。以例子1为例，mathematica的代码可以这样写:

ListPlot[Table[{(-1)^n/n, 0}, {n, 1, 100}]]

• 练习1.2.2：尝试验证$y_n=1-1/n$的极限为1.

例子1.2.3

一定要注意这里强调的是，以某个实数$L$为中心任意半径的对称区间$I_{\epsilon }(L)$的外部只有有限多个元素，而不是内部有无穷个序列的元素。这其中的差别可以从下面的例子当中看出来:

定义为$b_n=\{\begin{array}{ll}1& \text{n 是偶数}\\ -1& \text{n 是奇数}\end{array}$的序列没有极限。因为如果我们的对称区间$I_{\epsilon }(x)$的中心$x\ne 1,-1$的话，只要$\epsilon$足够小，那么其内部就不会包含任何序列$b_n$的元素，这就无法满足上面的几何性质。但是如果$x=1$,那么$I_{\epsilon }(x)$其外部始终存在无穷多个$b_n=-1$,反过来$x=-1$也是如此.因此不存在一个点满足这样的几何性质，因此该序列没有极限。

因此我们很容易想到这样一个命题：

序列极限的唯一性

如果一个序列的极限存在，那么极限是唯一的。

• 练习1.2.4：证明这个命题(使用一开始的$\epsilon -N_{\epsilon }$的定义证明)

2. 极限的基本性质

2.1 序列极限存在性的判断

定理2.1.1(单调有界原理)

单调(未必严格单调)有界的实序列极限必然存在。

回到一开始的问题，因为我们发现给具体的$a$以及初值$x_1$，只要n足够大，序列基本上是单调减少的，并且有一个明确的下界的。

n = 100;
(*这里a是2*)
(*因为收敛太快，所以为了方便观察单调性初值给得比较大，为100*)
data = RecurrenceTable[{x[k + 1] == (0.5)*(x[k] + 2/x[k]), 
    x[1] == 100}, x, {k, 1, n}];
DiscretePlot[data[[k]], {k, 1, n}]

其实从图像上就不难观察出来这个性质：

引理2.1.2

问题0当中定义的序列$x_n$当$n$足够大的时候，序列是不增的并且具有下界。

1. 下界的存在性:假设初值$x_1$为正实数,这样保证了序列$x_n>0$。因为基本不等式(算术平均值大于几何平均值，参考Power mean及其性质)会得到$x_{n+1}:=\frac{1}{2}\left(x_n+\frac{a}{x_n}\right)\ge \sqrt{a}$也就是说序列除了初值$x_1$以外，从$x_2$开始皆大于$\sqrt{a}$.
2. 单调性: 考虑$\frac{x_{n+1}}{x_n}$,这样只要此比值在某个$N$以后的任意$n>N$都小于1，那么单调减少就得到了确认。然后我们发现$\frac{x_{n+1}}{x_n}=\frac{1}{2}\left(1+\frac{a}{x_n^2}\right)\le 1$ 于是，我们意识到，如果单调有界原理成立，那么序列$x_n$的极限一定存在。

2.2 极限的四则运算

搞清楚极限的四则运算的性质，会简化我们求更为复杂的极限的过程。

命题2.2.1：极限的四则运算

假设序列两个序列的极限存在，并且$x_n\to a,y_n\to b$，那么一定有:

1. 加法：序列$x_n+y_n\to a+b$
2. 减法:序列$x_n-y_n\to a-b$
3. 乘法:序列$x_n{y}_n\to ab$
4. 除法:在$b\ne 0$的情况下，并且序列$\frac{x_n}{y_n}$有意义的情况下，$\frac{x_n}{y_n}\to \frac{a}{b}$

练习2.2.2：用极限的定义证明这四条性质。

• 其中极限的乘法性质的证明中包含了一个基本的处理技巧，即：我们可以利用三角不等式控制$|x_n{y}_n-ab|$。具体来说$\begin{array}{rl}|x_n{y}_n-ab|& \le |x_n{y}_n-x_{n}b|+|x_{n}b-ab|\\ & =|x_n|\cdot |y_n-b|+|x_n-a|\cdot |b|\end{array}$其实这样的思路也正是体现了数学中的“分而治之”的想法。此外这个技巧，在函数列的根序列的稳定性的”lemma 1.1.2”以及2. 问题2当中也重复使用了。

针对于问题0，我们的思路是极限序列的极限是存在的，那么假设$x_n\to b$.那么根据极限的四则运算$x_n+\frac{a}{x_n}\to b+\frac{a}{b}$，并且根据序列递归的定义，以及极限的唯一性，那么一定有代数方程:$2b=b+\frac{a}{b}$ 此方程的唯一正实数解就是$b=\sqrt{a}$.

于是到此，我们解决了问题0的第一部分，即$x_n\to \sqrt{a}$ 练习2.2.3：模仿差不多的思路解决问题2.

问题2.2.4

序列满足$x_{n+1}^2=x_n+1$证明对于任意初值$x_1>0,x_n>0$，此序列的极限存在并且$x_n\to \frac{\sqrt{5}+1}{2}$

解答:关于黄金比例数的迭代序列

3. 误差分析

实际上我们就是要考虑序列$e_n:=x_n-\sqrt{a}$，我们能想要知道它大概是一个怎样的增长级别，也有助于让我们判断，在给定精度下，至少要进行几次迭代才能在精度范围内，近似的得到$\sqrt{a}$。

因为我们知道$x_n$的迭代公式，这里有一个基本的想法：

基本的想法

如果我们知道两个序列$x_n,y_n$之间的一个方程$G(x_n,y_n)=0$同时知道其中一个的迭代$F(x_n,x_{n-1},...,x_1)=x_{n+1}$我们想要分析$y_n$，那么可以经由方程$G$把关于$x_n$迭代变成关于$y_n$的迭代。

• 具体的这个想法还在由素数定理得到p_n的渐近表达有所体现。

此处我们自然是可以把$x_n$的迭代通过方程$e_n=x_n-\sqrt{a}$变成$e_n$的迭代，从而得到:$\begin{array}{rl}{e}_{n+1}& =\frac{e_n^2}{2x_n}\end{array}$而前面我们分析了，如果$n>1$那么$x_n>\sqrt{a}$。于是:$|e_{n+1}|<\frac{|e_n{|}^2}{2\sqrt{a}},n>1$这说明牛顿法开根号产生的序列$x_n$是二次收敛这个速度的一个序列。

再者根据上面的不等式，我们可以迭代这个不等式最后得到: $|e_n|<2\sqrt{a}{\left(\frac{|e_1|}{2\sqrt{a}}\right)}^{2^n},n>1$ 这其实告诉我们，使用牛顿迭代法的时候，我们要先预估初始的误差$e_1$，然后把这个结果带入到绝对误差的迭代中，估计迭代n次以后的误差是多少。这表明，我们从一开始就需要对$\sqrt{a}$做出一个粗略的预估，然后把这个预估的值作为初值，同时得到绝对误差的初始值，这样才能绝对在给定精度下，需要迭代多少次。

比如说我们可以看看下面这个实际的问题：

问题3.1

对无理数$\sqrt{2}$进行有理数逼近，要求绝对误差不大于$10^{-10}$.

• 选定一个好的初值的重要性：

下面是一个Wolfram/mathematica的程序用于计算这个过程，程序会根据我们输入的精度要求，以及初值，输出最后达到这样的精度至少需要的迭代次数，以及近似值。

NewtonSqrt[a_, x0_, tol_] := Module[{x = x0, e = Abs[x0^2 - a], n = 0},
  While[e > tol,
    x = 1/2 (x + a/x);
    e = Abs[x^2 - a];
    n++;
  ];
  {n, x}
]
 
(* 使用示例 *)
(* a = 2，x0 = 1，tol = 10^-10 *)
NewtonSqrt[2, 1, 10^-10]
 

1. 如果选一个很烂的初值，比如100:

执行我们的程序:

NewtonSqrt[2, 100, 10^-10]

最后的输出是需要迭代至少$10$次，并且会得到一个非常复杂的即约分数。

2. 如果认真选一个初值，比如:$3/2$ 我们可以根据$1<\sqrt{2}<3/2$,以及我们的初值$x_1=3/2$，误差初值为$|e_1|<1/2$,根据绝对误差的迭代式子在程序运行之前预估需要迭代的次数不会超过$5$次。

NewtonSqrt[2, 3/2, 10^-10]

最后实际上我们经过了$3$次迭代，最后得到有理数的逼近结果是$\sqrt{2}\approx \frac{665857}{470832}$