牛顿法与方程近似根

LoR_Zhong · 2022-11-13 10:56:02 · 算法·理论

引入

一般来说，求高次方程的精确根十分困难，而四次以上的方程甚至根本不存在求根公式。因此，求出高次方程的近似根便极为重要。

前置芝士

1.导数（derivative）

如何求一个二次函数 y=x^2 在某一点 P_0(x_0,y_0) 处的切线 P_0T 的斜率 k_0？

如果用初中生的思维，那你的做法一定是这样：

设切线 P_0T:y=k_0x+b，则：

由上式减下式得：

由判别式得：

从而：

将点 P_0(x_0,y_0) 代入 y=k_0x+b 得：

将 b=y_0-x_0k_0 代入 k_0^2=-4b 得：

最后用求根公式：

将 y_0=x_0^2 代入得：

这是初中生的做法，那什么是高中生的做法呢？

如下：

令 f(x)=x^2，则：

与前面的一比，是不是明显简洁多了？而这种做法，就是我所说的导数。

百度百科对导数的定义如下：

设函数 y=f(x) 在点 x_0 的某个邻域内有定义，当自变量 x 在 x_0 处有增量 \Delta x，(x_0+\Delta x) 也在该邻域内时，相应地函数取得增量 \Delta y=f(x_0+\Delta x)-f(x_0)；如果 \Delta y 与 \Delta x 之比当 \Delta x\to0 时极限存在，则称函数 y=f(x) 在点 x_0 处可导，并称这个极限为函数 y=f(x) 在点 x_0 处的导数，记作 f'(x_0) 或 y'\vert _{x=x_0} 或 \dfrac{\operatorname{d}\!y}{\operatorname{d}\!x}\vert _{x=x_0}\, .

另外，如果函数 f(x) 在点 x_0 处存在导数简称函数 f(x) 在点 x_0 处可导，否则不可导。若函数 f(x) 在区间 (a,b) 内每一点都可导，就称函数 f(x) 在区间 (a,b) 内可导。

来自小白的惶恐 简单地说，就是表示函数曲线在点 P_0(x_0,f(x_0)) 处的切线的斜率。对，就是我们刚刚求的那个。

随堂练习：求函数 f(x)=x^4-\dfrac{2\ln x}{x^3}+\ln(2x^2)-x 在 x=e 处的导数 \cdots\cdots 开个玩笑，求函数 g(x)=\sqrt{x} 在 x=3 处的导数。

算好了吗？现在公布答案：

开玩笑的那个有同学算了吗？公布答案：

或许有同学会说：诶，等等，这怎么算的啊。那，是时候把导数的运算法则教给你了。

1. \left[f(x)\pm g(x)\right]'=f'(x)\pm g'(x)\,;
2. \left[f(x)g(x)\right]'=f'(x)g(x)+f(x)g'(x)\,;
3. \left[\dfrac{f(x)}{g(x)}\right]'=\dfrac{f'(x)g(x)-f(x)g'(x)}{\left[g(x)\right]^2}\quad(g(x)\ne0)\,.

什么？难道我们问的真的不是一个函数的导数怎么算么？ 当然，还有一些基本初等函数的导数公式，它们都可以直接使用：

1. 若 f(x)=c\,\,(c 为常数 )，则 f'(x)=0\,;

2. 若 f(x)=x^\alpha\,\,(\alpha\in\mathbb{Q} ，且 \alpha\ne0\,)，则 f'(x)=\alpha x^{\alpha-1}\,;

3. 若 f(x)=\sin x，则 f'(x)=\cos x\,;

4. 若 f(x)=\cos x，则 f'(x)=-\sin x\,;

5. 若 f(x)=a^x\,\,(a>0 ，且 a\ne1\,)，则 f'(x)=a^x\ln a\,; 特别地，若 f(x)=e^x，则 f'(x)=e^x\,;

6. 若 f(x)=\log_a x\,\,(a>0，且 a\ne1\,)，则 f'(x)=\dfrac{1}{x\ln a}\,; 特别地，若 f(x)=\ln x，则 f'(x)=\dfrac{1}{x}\,.

2.高阶导数（higher derivative）

先举个简单的栗子：求函数 f(x)=x^2 的二阶导数。

我们已经知道 f(x) 的一阶导数是 f'(x)=2x，我们只要再对这个一阶导数求导即可，即 f''(x)=\left[f'(x)\right]'=2\,.

定义：

函数 y=f(x) 的导数 y'=f'(x) 仍是 x 的函数。我们把 y'=f'(x) 的导数叫做函数 y=f(x) 的二阶导数（second derivative），记作 y'' 或 \dfrac{\operatorname{d}\!^2y}{\operatorname{d}\!x^2}，即 y''=(y')' 或 \dfrac{\operatorname{d}\!^2y}{\operatorname{d}\!x^2}=\dfrac{\operatorname{d}}{\operatorname{d}\!x}\left(\dfrac{\operatorname{d}\!y}{\operatorname{d}\!x}\right)。以此类推，n-1 阶导数的导数叫做 n 阶导数（nth derivative），记作：

或

3.泰勒级数（Taylor series）

百度百科：在数学中，泰勒级数用无限项连加式——级数来表示一个函数，这些相加的项由函数在某一点的导数求得。

定义：如果在点 x=x_0 具有任意阶导数，则幂级数

称为在点 x_0 处的泰勒级数。

我们可以用 R_n(x) 表示泰勒公式中的余项，即：

4.直线（straight line）

难的学完了再看看简单的：已知直线 y=kx+b，求它与坐标轴的交点。

嗯，有手就行：与 x 轴交点 (-\dfrac{b}{k},0)，与 y 轴交点 (0,b)。

那这个呢？已知直线 l 斜率为 k，且过点 (x_0,y_0)，求直线 l 的直线解析式。

很简单，运用一点初中知识，我们得到 l:y=kx+y_0-kx_0，那么它与 x 轴的交点的横坐标就是 \dfrac{kx_0-y_0}{k}\,.

正题

牛顿迭代法（Newton's method）

大前提：函数 f(x) 可导！

百度百科的定义：

设 r 是 f(x)=0 的根，选取 x_0 作为 r 的初始近似值，过点 (x_0,f(x_0)) 做曲线 y=f(x) 的切线 L，L:y=f(x_0)+f'(x_0)(x-x_0)，则 L 与 x 轴交点的横坐标 x_1=x_0-\dfrac{f(x_0)}{f'(x_0)}，称 x_1 为 r 的一次近似值。过点 (x_1,f(x_1)) 做曲线 y=f(x) 的切线，并求该切线与 x 轴交点的横坐标 x_2=x_1-\dfrac{f(x_1)}{f'(x_1)}，称 x_2 为 r 的二次近似值。重复以上过程，得 r 的近似值序列，其中，x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)} 称为 r 的 n+1 次近似值，上式称为牛顿迭代公式。

那为什么它可以成立呢？这要涉及到泰勒级数。用牛顿迭代法解非线性方程，是把非线性方程 f(x)=0 线性化的一种近似方法。把 f(x)=0 在点 x_0 的某邻域内展开成泰勒级数 f(x_0)=f(x_0)+f'(x_0)(x-x_0)+\dfrac{f''(x_0)}{2!}(x-x_0)^2+\cdots+\dfrac{f^{(n)}(x_0)}{n!}+R_n(x)，取其线性部分（即泰勒展开的前两项），并令其等于 0，即 f(x_0)+f'(x_0)(x-x_0)=0，以此作为非线性方程 f(x)=0 的近似方程，若 f'(x_0)\ne0，则其解为 x_1=x_0-\dfrac{f(x_0)}{f'(x_0)}，这样，得到牛顿迭代法的一个迭代关系式：x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}。

这么长一堆，真是可怕！但我们只需按照它的步骤，一步步来，我们也能求出近似根。

好，那么开始了？如果你觉得你的运算能力不错，可以心算，但有计算器的话。嘿嘿，那就 不讲武德！ 拿出来用吧。

方程：x^3+2=0，初始近似值为 x_0=-2\,.

计算器的步骤：x,-,\dfrac{\fcolorbox{black}{black}{}}{\Box},x,\text{shift},x^2\,(x^3),+,2, ↓ ,\text{shift},\int^\Box_\Box\fcolorbox{black}{black}{ }\,(\dfrac{\operatorname{d}}{\operatorname{d}\!x}\fcolorbox{black}{black}{}),(,x,\text{shift},x^2\,(x^3),+,2,), → ,x, → ,\text{sto},x\,.在这之前先设置 x 为 2，然后一直按 = 键就行了。

算好了吗？当然算不好，我还没告诉你们精确度呢，真是悲催。

所谓精确度，用希腊字母 \varepsilon 表示，当 \lvert x_{i+1}-x_i\rvert<\varepsilon 时，算法就可以停止，否则就继续算下去。

好，那现在给定精确度 \varepsilon=10^{-6}，算出来近似根是多少？

嗯~，应该是 -1.25992\cdots 吧。再开个三次方，你就会发现它确实很接近 -2，这说明牛顿迭代法是有效的不是吗？

但是你再试试这个：\sqrt[3]{x}=0，初始近似值为 x_0=-2，精确度 \varepsilon=10^{-6}\,.

啊呀呀，这下好了，不管你是用计算器还是心算，你都会发现近似根的绝对值在不断增大，可方程的根一看就是 x=0 啊。（如图）

来源：几何画板。

是时候 展现真正的技术了！ 对牛顿迭代法进行改进了！

牛顿下山法（Newton down-hill method）

要想练就绝世武功，就要忍受常人难忍受的痛

百度百科，稍微看看就行了。

大体意思即：为了保证 x_k 是在向根 r 逼近，对牛顿迭代法提出 \lvert f(x_{k+1})\rvert<\lvert f(x_k)\rvert 的条件（即函数值必须是在逼近 0），并有了新的迭代关系式：x_{k+1}=x_k-\lambda_k\dfrac{f(x_k)}{f'(x_k)}\,.

其中 \lambda_k 称为“下山因子”，满足 0<\lambda_k\le 1。当 \lambda_k 取某个值时，如果仍有 \lvert f(x_{k+1})\rvert\ge\lvert f(x_k)\rvert，就对 \lambda_k 进行减半操作，直到满足条件为止。

可参考：https://wenku.baidu.com/view/8d98b85ff78a6529647d53bd.html?fr=sogou&_wkts_=1668844957938

可以写出整数次幂非线性方程求解的代码，如下：

#include<cmath>
#include<iomanip>
#include<iostream>
using namespace std;
int _f[110],_f_[110],n;
double r,r_pre,epsilon;
double f(double x){
    double sum=_f[n];
    for(int i=n-1;i>=0;i--) sum*=x,sum+=_f[i]; //秦九韶算法 
    return sum;
}
double f_(double x){
    double sum=_f_[n-1];
    for(int i=n-2;i>=0;i--) sum*=x,sum+=_f_[i];
    return sum;
}
int main(){
    while(n<2) cin>>n; //保证是非线性方程 
    for(int i=n;i>0;i--){
        cout<<"请输入"<<i<<"次项：";
        cin>>_f[i];
        _f_[i-1]=_f[i]*i;
    }
    cout<<"请输入常数项：";
    cin>>_f[0];
    cout<<"请输入精确度：";
    cin>>epsilon;
    cout<<"请输入初始近似根：";
    cin>>r;
    do{
        r_pre=r;
        for(double lambda=1;;lambda/=2){ //牛顿下山法 
            r=r_pre-lambda*f(r_pre)/f_(r_pre);
            if(fabs(f(r))<fabs(f(r_pre))) break;
        }
    }while(fabs(r_pre-r)>=epsilon);
    cout<<fixed<<setprecision((int)log10(1/epsilon))<<r<<endl; //有效位数 
    cout<<fixed<<setprecision((int)log10(1/epsilon))<<f(r);
    return 0;
}

它甚至能解一百次方程！ 就我试验来看，它是没问题的。

后置芝士

1.复数（complex number）

上过初中的同学应该知道，方程 x^2+1=0 的 \Delta 小于 0，是没有实数根的。诶，好奇怪，为什么不说无解呢？那是因为它是有解的（这不是废话吗），只不过涉及到 \sqrt{-1} 的问题罢了。

我们令 i^2=-1，并将形如 z=a+bi（a,b 均为实数）的数称为复数，其中 a 称为实部 （real part），b 称为虚部（imaginary part），i 称为虚数单位（imaginary unit）。当 b=0 时，z 为实数；当 a=0 且 b\ne 0 时，z 为纯虚数（pure imaginary）。

这样就好了，方程 x^2+1=0 虽然在实数域没有根，但在复数域就有两个根 x_1=i,x_2=-i\,.

复数的运算法则：

1. (a+bi)\pm(c+di)=(a\pm c)+(b\pm d)i\,;
2. (a+bi)(c+di)=(ac-bd)+(bc+ad)i\,;
3. \dfrac{a+bi}{c+di}=\dfrac{(a+bi)(c-di)}{(c+di)(c-di)}=\dfrac{(ac+bd)+(bc-ad)i}{c^2+d^2}\,.

零碎知识点：

• i^{4n+1}=i,i^{4n+2}=-1,i^{4n+3}=-i,i^{4n}=1\,\,(n\in\mathbb{Z})\,;

• 复数 z_1=a+bi,z_2=c+di,z_1=z_2 当且仅当 a=c 且 b=d\,;

• 复数 z=a+bi 的模记作 \lvert z\rvert,\lvert z\rvert=\sqrt{a^2+b^2}\,;

• 复数 z=a+bi 的共轭复数记作 \overline{z},\overline{z}=a-bi\,.z+\overline{z}=2a,z\cdot\overline{z}=a^2+b^2\,;

• 复数集记作 \mathbb{C}\,;

• 百度百科

2.复数平面（complex plane）

百度百科：复数平面即是 z=a+bi，它对应的坐标为 (a,b)。其中，a 表示的是复平面内的横坐标，b 表示的是复平面内的纵坐标，表示实数 a 的点都在 x 轴上，所以 x 轴又称为“实轴”；表示纯虚数 bi 的点都在 y 轴上，所以 y 轴又称为“虚轴”。y 轴上有且仅有一个实点即为原点"0"。

复平面内的每一个点，有唯一的一个复数和它对应；反过来，每一个复数，有复平面内唯一的一个点和它对应，所以复数集 \mathbb{C} 和复平面内所有的点所成的集合是一一对应的。

3.代数基本定理（fundamental theorem of algebra）

百度百科：任何复系数一元 n\,\,(n\ge1) 次多项式方程在复数域上至少有一根。

牛顿法与复数根

同实数根一样，只需将平面直角坐标系换成复平面，就有递推式 z_{k+1}=z_k-\dfrac{f(z_k)}{f'(z_k)}，同样也可以得到复数近似根。

闭幕

文中引用资料，侵权删。