解一元方程的迭代法

引言:整理该博客的起因源于最近做的一道算法题:

题目描述
ax3+bx2+cx1+dx0=0 这样的一个一元三次方程。
给出该方程中各项的系数(a,b,c,d a,b,c,d均为实数)，并约定该方程存在三个不同实根(根的范围在-100至100之间)，且根与根之差的绝对值之差≥1。要求由小到大依次在同一行输出这三个实根(根与根之间留有空格)，并精确到小数点后2位。
输入格式
一行，4个实数A,B,C,D
输出格式
一行，3个实根，并精确到小数点后2位。

输入输出样例
输入
1 -5 -4 20
输出
-2.00 2.00 5.00

提示：记方程f(x)=0，若存在2个数x1和x2，且x1<x2，f(x1)f(x2)<0，则在(x1,x2)之间一定有一个根。

解析

当f(x)是不超过四次的多项式石,方程的解可用公式来求得,例如此题可采用卡丹公式,盛金公式求解.这里不做讨论.五次以上的多项式不可用公式求解,通常采用迭代法.

解法1: 枚举法

优点:简单粗暴
缺点:效率低

#include <iostream>
#include <cstdio>
using namespace std;
int main()
{
   double a,b,c,d;
   scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
   for(double i=-100;i<=100;i+=0.001)
   {
      double j=i+0.001;
      double y1=a*i*i*i+b*i*i+c*i+d;
      double y2=a*j*j*j+b*j*j+c*j+d;
      if(y1>=0&&y2<=0||y1<=0&&y2>=0)
      {
         double x=(i+j)/2;
         printf("%.2lf ",x);
      }
   }
}

解法2:二分法

在迭代过程中确定一个区间序列 { I ( k ) } \{{I^{(k)}}\} {I(k)},使每个区间都包含方程的某个解$x^{\ast }$,且区间长度趋向于零.
二分法特点:一定收敛
x(-100≤x≤100)。由于根与根之差的绝对值≥1，因此设定搜索区间[x1，x2]，其中x1=x，x2=x+1。若
(1) f(x1)=0，则确定x1为f(x)的根；
(2) f(x1)*f(x2)>0，则确定根x不在区间[x1，x2]内，设定[x2，x2+1]为下一个搜索区间
(3) f(x1)*f(x2)<0，则确定根x在区间[x1，x2]内。
如果确定根x在区间[x1，x2]内的话（f(x1)*f(x2)<0），采用二分法，将区间[x1，x2]分成左右两个子区间：左子区间[x1，x]和右子区间[x，x2],
如果f(x1)*f(x)≤0，则确定根在左区间[x1，x]内，将x设为该区间的右指针（x2=x），继续对左区间进行对分；如果f(x1)*f(x)>0，则确定根在右区间[x，x2]内，将x设为该区间的左指针（x1=x），继续对右区间进行对分；
上述对分过程一直进行到区间的间距满足精度要求为止（x2-x1<0.001）。此时确定x1为f(x)的根。

#include<cstdio>
double a,b,c,d;
double fc(double x)
{
    return a*x*x*x+b*x*x+c*x+d;
}
int main()
{
    double l,r,m,x1,x2;
    int s=0,i;
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);  //输入
    for (i=-100;i<100;i++)
    {
        l=i; 
        r=i+1;
        x1=fc(l); 
        x2=fc(r);
        if(!x1) 
        {
            printf("%.2lf ",l); 
            s++;
        }      //判断左端点，是零点直接输出。
                        
                        //不能判断右端点，会重复。
        if(x1*x2<0)                             //区间内有根。
        {
            while(r-l>=0.001)                     //二分控制精度。
            {
                m=(l+r)/2;  //middle
                if(fc(m)*fc(r)<=0) 
                   l=m; 
                else 
                   r=m;   //计算中点处函数值缩小区间。
            }
            printf("%.2lf ",r);  
            //输出右端点。
            s++;
        }
        if (s==3) 
            break;             
            //找到三个就退出大概会省一点时间
    }
    return 0;

收敛速度

二分法的收敛速度是线性的,且收敛常数q=1/2,初始长度为$L=|x^{(1)}-x^{(0)}|$,k次迭代后有$|x^{(k)}-x^{(\ast )}|\le (1/2{)}^{k}L\le \epsilon$,迭代次数k应满足$k\ge {\log }_{2}L/\epsilon$

解法3: 简单迭代法 牛顿迭代法 割线法

对于一个未知量的非线性方程$f(x)=0$,实数$x^{\ast }$为方程的解,用迭代的办法产生一个序列 { x ( k ) } \{x^{(k)}\} {x(k)},并希望这个序列收敛于$x^{\ast }$.
主要讨论
(1)迭代公式的构造,
(2)初始近似值$x^{(0)}$的选择,
(3)迭代的收敛性
收敛的准则:
一般是第k次迭代$x^{(k)}$充分接近$x^{\ast }$或者是$|f(x^{(k)})|$充分小.

(1)简单迭代法

将原方程$f(x)=0$改写为等价方程为$x=\phi (x)$,可获得迭代形式
$$x^{(k+1)}=\phi (x^{(k)}),k=0,1,2,3,...$$

$f(x)=a{x}^3+b{x}^2+cx+d$.
$x^{(k+1)}=1/c\ast (a{x}^3+b{x}^2+d)$

(2)牛顿迭代法

牛顿迭代法的几何意义容易理解

用切线近似代替曲线y=f(x),以切线与ｘ轴的交点，即切线方程的跟作为下次近似$x^{(k＋\mathrm{１})}$

收敛性

满足以下四个条件，牛顿迭代收敛：

先对函数求导,$f^{\prime }(x)=3a{x}^2+2\ast bx+c$.
然后直接求根公式求f’(x)=0的点,也就是函数极点.
这题保证有三个不定根,所以有两个单峰.
我们分别设这两个点为p,q.
然后显然的必有三个根在[-100,p),[p,q],(q,100]三个区间内 (两极点间必定存在零点,勘根定理).
$f(x)=a{x}^3+b{x}^2+cx+d$的导函数是
$f^{\prime }(x)=3a{x}^2+2bx+c$,由$f^{\prime }(x)=0$,由此得
$x1=((-2)b-sqrt(4b^2-12ac))/(6\ast a)$;
$x2=((-2)b+sqrt(4b^2-12ac))/(6\ast a)$;
然后用牛顿迭代法

#include<iostream>
#include<cstdio>
#include<cmath>
#define eps 1e-３
using namespace std;
double x1,x2,x3,a,b,c,d;
double f(double x){return a*x*x*x+b*x*x+c*x+d;}
double df(double x){return 3*a*x*x+2*b*x+c;}
double slove(double l,double r)
{
	double x,x0=(l+r)/2;
	while(abs(x0-x)>eps)
	  x=x0-f(x0)/df(x0),swap(x0,x);
	return x;
}
int main()
{
	cin>>a>>b>>c>>d;
	double p=(-b-sqrt(b*b-3*a*c))/(3*a);
	double q=(-b+sqrt(b*b-3*a*c))/(3*a);
	x1=slove(-100,p),x2=slove(p,q),x3=slove(q,100);
	printf("%.2lf %.2lf %.2lf",x1,x2,x3);
	return 0;
}

(3)割线法

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
double f(double x,double a,double b,double c,double d)//用于计算三次函数值
{
    return a*x*x*x+b*x*x+c*x+d;
}
double solve(double a,double b,double c,double d,double x0,double x1)//割线法求解
{
    double xx[2]={x0,x1};//用xx数组滚动存储结果
    int flag=0;
    while(abs(xx[0]-xx[1])>1e-3)//控制精度
    {
        xx[flag]=xx[flag]-f(xx[flag],a,b,c,d)/((f(xx[flag],a,b,c,d)-f(xx[flag^1],a,b,c,d)))*(xx[flag]-xx[flag^1]);
        flag^=1;
    }
    return xx[0];
}

double a,b,c,d,x[3]={-1000,-1000,-1000};
int tot;

bool check(double xx)//判断是否为新的解
{
    bool flag=0;
    for(int i=0;i<3;++i)flag|=abs(xx-x[i])<0.5;
    return !flag;
}

int main()
{
    scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
    for(double i=-100;i<=100;i+=0.5)//注意这里以0.5为步长枚举。我一开始步长为一，结果有一个点答案有1.0,2.0，直接除以零GG
    {
        double xx=solve(a,b,c,d,i,i+0.5);
        if(check(xx))x[tot++]=xx;
    }
    sort(x,x+3);//不要忘记从小到大排序 
    for(int i=0;i<3;++i)printf("%.2lf ",x[i]);
}

总结

收敛速度牛顿法最快，为超线性收敛，区间法和简单迭代法为线性收敛．