点击前往我的blog获取更好阅读体验

0 前置知识&&写在前面

• 学习牛顿迭代法需要能熟练地掌握求导QwQ

• 要想看懂牛顿迭代法的二次收敛证明需要一定的高数基础QwQ（看不懂也无所谓，会用就行）

• 学习泰勒公式需要能熟练地掌握求导并对无穷级数有一定了解QwQ

• • 都不会怎么办？没有问题，会用就行NOIP不考QwQ

• 求零点、不动点、函数近似值、极值、最值……基本可以认为是一类问题，所以本文以讨论如何 求零点的近似值 为主

• 因为笔者太弱所以没有遗传算法、模拟退火之类的玄学算法

• 笔者是\displaystyle\sum狂魔，要是有式子不理解请展开QwQ

对于低于5次的多项式方程，我们有通用的公式解法求零点的精确值

对于一些特殊高次多项式方程（例如可以因式分解的或者满足一些特定形式的方程）的和一些特殊的超越方程，我们也有方法求零点的精确值

但是其余的情况呢？

目前来说我们只能求近似值QwQ（而且在实际应用中，精确值往往也会被转换成近似值）

1 二分法 (int\long long\etc.)

这个方法可以说相当常见了（括号里是此法经常被应用于的数据类型）QwQ

对于在区间[l,r]内单调、连续且有f(l)\cdot f(r)<0成立的f(x)，做如下操作：

1. 计算mid=\frac{l+r}{2}
2. 若f(l)\cdot f(mid)<0，则令r=mid，否则令l=mid
3. 如果达到预定精度，跳转到4，否则跳转到1
4. 结束

循环次数：

\lceil\log_2 {\frac{r-l}{\epsilon}}\rceil

附程序：

double BinarySearch(double l, double r, double EPS) {
  double mid;
  while(r-l > EPS) {
    mid = l + (r - l)/2; //防溢出
    if(F(mid)*F(l) <= 0) r = mid;
    else l = mid;
  }
  return mid;
}

2 0.618法\优选法 (float\double\long double)

可能有的读者没听说过，不过这个方法其实还是很常用的QwQ

这个方法更常用于求单峰函数最值，所以这里以求单峰函数最值为例讲解

先证明一下它的优越性（过程摘自人教版高中数学选修4-7 优选法与试验设计初步）（没错真的有这本选修QwQ）

为了使每次去掉的区间有一定的规律性，我们这样来考虑：每次舍去的区间占舍去前的区间的比例数相同

下面进一步分析如何按上述两个原则确定合适的试点。如图2-1 设第1试点、第2试点分别为x_1和x_2，x_2<x_1且x_1,x_2关于[a,b]的中心对称，即x_2-a=b-x_1

（图2-1，由GeoGebra生成）

显然， 不论点x_2(或点x_1)是好点还是差点，由对称性， 舍去的区间长度都等于b-x_1，不妨设x_2是好点，x_1是差点，于是舍去(x_1,b]。再在存优范围[a,x_1]内安排第3次试验，设试点为x_3，x_3与x_2关于[a,x_1]的中心对称（如图2-2所示）。

（图2-2，由GeoGebra生成）

点x_3应在点x_2左侧，因为如果点x_3在点x_2的右侧，那么当x_3是好点，x_2是差点时，要舍去区间[a,x_2]，而它的长度与上次舍去的区间(x_1,b]的长度相同，违背成比例舍去的原则。于是，不论点x_3（或点x_2）是好点还是差点，被舍去的区间长度都等于x_1-x_2，按成比例舍去的原则，我们有等式

其中，左边是第一次舍去的比例数，右边是第二次舍去的比例数， 对式(1)变形，得 $$1-\frac{b-x_1}{b-a}=1-\frac{x_1-x_2}{x_1-a}

即

式(2)两边分别是两次舍弁后的存优范围占舍弃前全区间的比例数，设每次舍弃后的存优范围占舍弃前全区间的比例数为$t$，即 $$\frac{x_1-a}{b-a}=t$$(3) 则由$b-x_2=x_1-a$可得 $$\frac{x_2-a}{b-a}=1-t$$(4) 由式(2)得 $$\frac{x_1-a}{b-a}=\frac{\frac{x_2-a}{b-a}}{\frac{x_1-a}{b-a}}

把(3)与(4)代入(5)，得

t=\frac{1-t}{t}

即

t^2+t-1=0

解得t_1=\frac{-1+\sqrt{5}}{2}，t_2=\frac{-1-\sqrt{5}}{2}，其中t_1为对本问题有意义的根，这就是黄金分割常数，用\phi表示（注：原文用\omega表示）

一句话概括就是在缩小区间后可以只计算一个试点坐标，从而保证最优

操作流程如下（和二分法类似）

1. 计算mid1=l\cdot\phi+r\cdot(1-\phi)，mid2=l\cdot(1-\phi)+r\cdot\phi
2. 若f(l)\cdot f(mid1)>0，则令l=mid1,mid1=mid2,mid2=l\cdot(1-\phi)+r\cdot\phi，否则令r=mid2,mid2=mid1,mid1=l\cdot\phi+r\cdot(1-\phi)
3. 如果达到预定精度，跳转到4，否则跳转到2 （注意这里跳转到2）
4. 结束

附程序：

const double PHI=0.61803399; //phi记得用精确一点的值，防止卡精度
const double 1mPHI=0.38196601;
double GoldSearch(double l, double r, double EPS) {
    double mid1 = l+1mPHI*(r-l), mid2 = l+PHI*(r-l);
    while(r-l > EPS) {
        if(F(mid1) < F(mid2)) {
            l = mid1;
            mid1 = mid2;
            mid2 = l+PHI*(r-l);
        } else {
            r = mid2;
            mid2 = mid1;
            mid1 = l+1mPHI*(r-l);
        }
    }
    return (mid1+mid2)/2;
}

读者们可以在洛谷P3382中测试一下(～o￣3￣)～（虽然是三分法模板但也可以用优选法\二分法做哦）

关于这个还有一个类似方法：斐波那契法。有兴趣的读者可以查阅相关资料才不是笔者不想写_(:3」∠)_

3 泰勒公式

先讲这个是为了为下文牛顿迭代法二次收敛的证明做铺垫，不想看证明的可以略过QwQ（不过还是推荐了解一下，挺有趣的）

这里假定函数f(x)在x_0处有任意阶导数

我们可以很容易地求出多项式和类指数函数的近似值，但是像三角函数、对数函数这样的我们又该如何求近似值呢

对了，就是用泰勒公式QwQ

泰勒公式的想法很简单，就是构造一个多项式函数g(x)=\displaystyle\sum_{k=0}^n{a_kx^k}，使得它与函数f(x)在x_0处的原函数值和各阶导数均相等，即

f(x_0)=g(x_0) f^\prime(x_0)=g^\prime(x_0) f^{\prime\prime}(x_0)=g^{\prime\prime}(x_0) ... f^{(n)}(x_0)=g^{(n)}(x_0)

（f^{(n)}(x)就代表f(x)的n阶导数，f^{(0)}(x)=f(x)）

因为

g^{(m)}(x)=\displaystyle\sum_{k=m}^n{\frac{k!}{(k-m)!}a_kx^{k-m}}

于是我们就得到了（解n元一次方程组，很简单的）

g(x)=\displaystyle\sum_{k=0}^n{\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k}

当n\rightarrow\infty时，我们可以认为f(x)=g(x)

而当n有一个确定的值时，f(x)就可以写成g(x)+R_n(x)了

其中R_n(x)是余项，它有好几种不同的写法，这里选用拉格朗日余项。因为它不仅在下文中牛顿迭代法的二次收敛证明中出现了，而且形式和前面的几项非常相近

R_k(x)=\frac{f^{(k+1)}(\xi_L)}{(k+1)!}(x-x_0)^{k+1}

其中\xi_L在x和x_0之间

要是了解过拉格朗日中值定理的读者应该很快就能领会

当n\rightarrow\infty时，有（泰勒级数）

\displaystyle\sum_{k=0}^\infty{\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k}

特别地，当x_0=0时，有（麦克劳林级数）

\displaystyle\sum_{k=0}^\infty{\frac{f^{(k)}(0)}{k!}}x^k

另外注意应用麦克劳林级数并且x在某个范围之外时，得到的结果可能是发散的（就是\infty，这个不展开讲，有兴趣的读者可以去学习无穷级数相关知识）

附上Wikipedia的动图

对证明感兴趣的读者可以自行查阅相关资料

不能理解无所谓反正NOIP不考，下面给出几个常见的泰勒级数

e^x=\displaystyle\sum_{k=0}^\infty{\frac{x^k}{k!}} \sin x=\displaystyle\sum_{k=0}^\infty{(-1)^k\frac{x^{2k+1}}{(2k+1)!}} \cos x=\displaystyle\sum_{k=0}^\infty{(-1)^k\frac{x^{2k}}{(2k)!}}

（有上面三个式子你就可以证明欧拉公式之e^{i\theta}=\cos\theta+i\sin\theta，e^{i\pi}+1=0了）

\ln{(1+x)}=\displaystyle\sum_{k=1}^\infty{(-1)^{k+1}\frac{x^k}{k}} \frac{1}{1-x}=\displaystyle\sum_{k=0}^\infty{x^k} (1+x)^m=\displaystyle\sum_{k=0}^\infty{\binom{m}{k}x^k}

程序略QwQ

4 牛顿迭代法

先说说这个方法的过程

1. 随便确定一个数x_0
2. 求在f(x_0)处的切线l:[y-f(x_0)]=f^\prime(x_0)(x-x_0)
3. 求切线l的零点x_1

稍加计算便得到了

x_1=x_0-\frac{f(x_0)}{f^\prime(x_0)}

既然是迭代，那么自然就有

x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)}

其中x_n代表第n次迭代

附上Wikipedia的动图

二次收敛证明（Wikipedia上的，笔者翻译QwQ）（还是那句话：看不懂无所谓，会用就行QwQ）：

根据Taylor's theorem，任何二阶导数连续的函数f(x)（设\alpha是根）都可以写成

f(\alpha)=f(x_n)+f^\prime(x_n)(\alpha-x_n)+R_1\qquad(1)

由Lagrange form of the Taylor series expansion remainder得

其中，$\xi_n$在$x_n$和$\alpha$之间。 由于$\alpha$是根，所以(1)式变为 $$0=f(\alpha)=f(x_n)+f^\prime(x_n)(\alpha-x_n)+\frac{1}{2}f^{\prime\prime}(\xi_n)(\alpha-x_n)^2\qquad(2)

(2)式两边同时除以f^\prime(x_n)，整理得

\frac{f(x_n)}{f^\prime(x_n)}+(\alpha-x_n)=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}(\alpha-x_n)^2\qquad(3)

由于

x_{n+1}=x_n-\frac{f(x_n)}{f^\prime(x_n)}\qquad(4)

代入(3)式，有

\underbrace{\alpha-x_{n+1}}_{\epsilon_{n+1}}=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}(\underbrace{\alpha-x_n}_{\epsilon_n})^2

即

\epsilon_{n+1}=\frac{-f^{\prime\prime}(\xi_n)}{2f^\prime(x_n)}\epsilon_n^2\qquad(5)

两边取绝对值，有

|\epsilon_{n+1}|=\frac{|f^{\prime\prime}(\xi_n)|}{2|f^\prime(x_n)|}\epsilon_n^2\qquad(6)

(6)式表明，如果函数满足以下条件，其为二次收敛

1. 对于所有的x\in I（I为区间[\alpha-r,\alpha+r],r\geq |\alpha-x_0|（即x_0\in I）），有f^\prime(x)\neq 0

2. 对于所有的x\in I，f^{\prime\prime}(x)连续

“足够接近”意为

$\qquad$b. 对于$C<\infty$（笔者注：这里写法好像有问题，存疑），$\frac{1}{2}|\frac{f^{\prime\prime}(x_n)}{f^\prime(x_n)}|<C|\frac{f^{\prime\prime}(\alpha)}{f^\prime(\alpha)}| \qquad$c. 对于$n\in\Bbb{N}$，$C|\frac{f^{\prime\prime}(\alpha)}{f^\prime(\alpha)}|\epsilon_n<1

当满足上述条件时，(6)式可以写为：

|\epsilon_{n+1}|\leq M\epsilon_n^2

其中

M=\displaystyle\sup_{x\in I}{\frac{1}{2}\bigg|\frac{f^{\prime\prime}(x)}{f^\prime(x)}\bigg|}

由条件3得M|\epsilon_0|<1

程序：

double Newton(double x0;int n) {
    for(int i=0; i<n; ++i) x0 = x0 - f(x0) / df(x0);
    return x0;
}

当然，上述各方法的应用范围远不止于此，有兴趣的读者可以自行查阅相关资料QwQ

5 赠品

5.1 \cos x=x的解析解

对于PhOer来说，\cos x=x这个方程应该是相当熟悉了QwQ

笔者在这里放上解析解（近似值x=0.739），详情见参考文献[2]（文献里讲的是t\sin x=x-m的解法，不过笔者太弱了，实在是看不懂QwQ）

\displaystyle\frac{\pi}{2}\exp\Bigg\{\frac{1}{\pi}\int_0^1{\frac{\arctan\bigg[\frac{(\pi x+2)\ln{\big(\frac{\sqrt{1-x^2}+1}{x}\big)}x}{x^2\ln^2{\big(\frac{\sqrt{1-x^2}+1}{x}\big)}-\pi x-1}\bigg]}{x}}\mathrm{d}x\Bigg\}

5.2 快速求\frac{1}{\sqrt{x}}

关于这个有一个相当有名的故事

为了节约篇幅，笔者把文章放在这里了 (当然直接看[5]也可以)

2018.09.18 upd:因为文章比较久远，所以在今天该算法的优势可能会没有文中说的那样明显

另外如果精度不够就多迭代几次，牛顿法是二次收敛的，一般迭代三四次就够用了

2018.10.03 upd:把测试程序补上点此 (应该在半个月前补上的QwQ)

6 后记

• 这篇文章偏数学一些，如果不能理解的话请多读几遍QwQ

• 其实可写的还有很多，限于篇幅就到此为止了（现在在后台打个字都要卡两秒QwQ）

• 因为笔者是个蒟蒻，所以如果有错误，烦请各位dalao不吝赐教

7 主要参考资料

• [1] 人教版高中数学选修4-7 优选法与试验设计初步

• [2] Siewert C E, Burniston E E. An exact analytical solution of Kepler's equation[J]. Celestial Mechanics, 1972, 6(3):294-304.

• [3] Newton's method - Wikipedia

• [4] Taylor's theorem - Wikipedia

• [5] 一个Sqrt函数引发的血案（这个是笔者能找到的最早的一篇了QwQ）