浅谈牛顿迭代法的应用、拓展及模板

sjh0626 · 2025-06-18 20:49:44 · 算法·理论

简介

牛顿迭代法可用来解非线性方程 f(x)=0 的所有根（包括复数解）。

公式

我们目前已知 f(x) 在 x_n 处的一阶泰勒展开：

然后再解一下 x：

因为已知 f(x)=0，并且我们为了迭代，还是把约等于改成等于吧，所以最后公式为：

怎么样？是不是一下子就豁然开朗了？

例题

求平方根

是的，你没听错，一个 sqrt 就能够搞定的事，居然还要手搓牛迭？

但是，这是每个初学者都要经历的，所以，继续往下看吧！

题目：给定一个非线性方程 x^2=n，求 x 的值。

有的人就要说了：这也不是等于 0 啊，你不会在骗人吧？

别急，我们观察方程，注意到可以构造为 f(x)=x^2-n=0 的方程。

所以可以做了，我们先求出 f(x) 的导数，即 f^{'}(x)=2x。

可以直接上公式了！

double f(double x){return x*x-n;}
double df(double x){return 2*x;}
double solve(){
    double x=1;
    while(fabs(f(x))>1e-6){x=x-f(x)/df(x);}
}

我们随便带入一个 n 吧，比如 2，可得 \sqrt{2} \approx1.41421356。

求多项式多解

题目 1：UVA10428 The Roots

这时有人发出了尖锐的爆鸣声：“？？？求出所有根，不是，这是人啊？”

别急，我们就需要用到一个非常好的东西：因式定理。

我们该如何证明呢？

我们假设一个多项式方程 f(x)=a_nx^n+a_{n-1}x^n+\dots+a_0=0。

那么这个多项式一定可以分解为 f(x)=(x-x_n)\times(x-x_{n-1})\times\dots\times(x-x_1)。

因为这个方程是一元 n 次方程，那么也一定有 n 个根（含复根、重根），那么 n 个根 x_1,x_2,\dots,x_n 就可以对应 (x-x_1)=0,(x-x_2)=0,\dots,(x-x_n)=0。

如果还不懂，就可以理解为当 x=x_1,x_2,\dots,x_n 时 f(x) 一定等于 0。

那么事情就变得简单了，我们可以先牛顿迭代出其中的一个根 z，然后就可以用 f(x) 去除以 (x-z)，然后继续牛顿迭代，循环往复，一定能求出 n 个根。

至于多项式除法：

随便带入一个多项式进去吧。

设 ax^3+bx^2+cx+d=(x-z)(nx^2+mx+p)=nx^3+(m-zn)x^2+(p-zm)x-zp。

可得 n=a,m=b+zn,p=c+zm。

再将 n,m,p,a,b,c,d 分别设为 b_2,b_1,b_0,a_3,a_2,a_1,a_0。

我们可得 b_i=a_{i+1}+b_{i+1}z，这就是新多项式的每一项系数了。

诶，这时就会有人问了：“f(x) 的导数咋求啊？”

这简单，我们可知幂函数 x^n 的导数为 nx^{n-1}，并且导数有一个法则：F(x)=f(x)+g(x) 可求出导数为 F^{'}(x)=f^{'}(x)+g^{'}(x)。于是每一项就正常处理就对了。

代码：

#include<bits/stdc++.h>
#define sjh0626s return
#define code 0
using namespace std;
int n,cnt=0,len;
double a[100010],temp,nx,z[100010]; 
double f(double x){
    double res=0;
    for(int i=n;i>=0;i--){
        res+=a[i]*pow(x,i);
    }
    return res;
}
double d(double x){
    double res=0;
    for(int i=n;i>=1;i--){
        res+=a[i]*i*pow(x,i-1);
    }
    return res;
} 
double newton_root(double x){
    for(int i=1;i<=150;i++){
        x-=f(x)/d(x);
    }
    return x;
}
void next_equation(double x){
    double b[100010]={0};
    for(int i=n-1;i>=0;i--)b[i]=a[i+1]+b[i+1]*x;
    for(int i=n;i>=0;i--)a[i]=b[i];
}
int main(){
    cin>>n;
    while(n!=0){
        for(int i=n;i>=0;i--){
            cin>>a[i];
        }
        cout<<"Equation "<<++cnt<<": ";
        len=0;
        for(n;n>=1;n--){
            nx=newton_root(1);
            z[++len]=nx;
            next_equation(nx);
        }
        sort(z+1,z+len+1);
        for(int i=1;i<len;i++)printf("%.4lf ",z[i]);
        printf("%.4lf\n",z[len]);
        cin>>n;
    }
    sjh0626s code;
}

题目 2：P1024 [NOIP 2001 提高组] 一元三次方程求解

这其实和上一道题一样，只是变成了求三次方程，一样，上模板！

代码：

#include<bits/stdc++.h>
#define sjh0626s return
#define code 0
#define ll long long
#define PII pair<int,int>
using namespace std;
int t=1;
double a,b,c,d;
template<class T>
class Newton{ //牛顿迭代法封装
public:
    T a[100010];
    T b[100010];
    int n;
    Newton(){}
    Newton(initializer_list<T>lis){ 
        n=lis.size()-1;
        int i=n;
        for(auto it=lis.begin();it!=lis.end();it++,i--)a[i]=*it;
    }
    Newton(vector<T>vec){
        n=vec.size()-1;
        int i=n;
        for(auto it=vec.begin();it!=vec.end();it++,i--)a[i]=*it;
    }
    T f(T x){ //原函数
        T res(0.0),px(1.0);
        for(int i=0;i<=n;i++){res=res+a[i]*px,px=px*x;}
        return res;
    }
    T df(T x){ //导数
        T res(0.0),px(1.0);
        for(int i=1;i<=n;i++)res=res+a[i]*i*px,px=px*x;
        return res;
    }
    T solve(T x){ //迭代
        while(fabs(f(x))>1e-12)x=x-f(x)/df(x);
        return x;
    }
    void change(T z){ //因式分解求 f(x)=(x-z)g(x)
        b[n]=0;
        for(int i=n-1;i>=0;i--)b[i]=a[i+1]+b[i+1]*z;
        for(int i=n;i>=0;i--)a[i]=b[i];
    }
};
set<double>ans;
void solve(){
    cin>>a>>b>>c>>d;
    Newton<double> s({a,b,c,d}); //初始化系数
    for(int i=1;i<=3;i++){
        double z=s.solve(1);
        ans.insert(z);
        s.change(z);
    }
    for(auto z:ans){
        printf("%.2lf ",z);
    }
}
int main(){
//    freopen(".in","r",stdin);
//    freopen(".out","w",stdout);
//    cin>>t;
    while(t--){
        solve();
    }
    sjh0626s code;
}

求复平面上根

题目：SP481 KMSL4B - Roots of polynomial

这时就又有热于提问的观众要问了：“主播主播，一般的我们会了，但是复数运算我们不会啊！”

欸，那我们先说复数吧！

我们要知道复数的概念，要先知道虚数单位 i，i 的定义是 i^2=-1。很好，可以讲复数了，复数有一个实部 a，还有一个虚部 b，那么它的值就是 z=a+bi。

先说加法，很简单啊，就是两边实部相加，虚部相加啊。那么值就是 z=(a+bi)+(c+di)=(a+c)+(b+d)i。

减法就不用了，和加法一样。

再说乘法吧，这也不算难，就只需要展开而已，z=(a+bi)(c+di)=ac+adi+bci-bd=(ac-bd)+(ad+bc)i。

最后再说除法吧，这咋办啊，我们要多项式除法啊，不用慌，一个待定系数，解决所有问题。

我们设 (q+pi)(c+di)=a+bi。

展开，得 (qc-pd)+(qd+pc)i=a+bi。

由此得 a=qc-pd,b=qd+pc。

再推，得 p=\frac{qc-a}{d}。

带入 b，得 b=qd+\frac{qc^2-ac}{d}=\frac{qd^2+qc^2-ac}{d}。

解一下 q：