牛顿迭代法

牛顿迭代法快速寻找平方根

    下面这种方法可以很有效地求出根号a的近似值：首先随便猜一个近似值x，然后不断令x等于x和a/x的平均数，迭代个六七次后x的值就已经相当精确了。
    例如，我想求根号2等于多少。假如我猜测的结果为4，虽然错的离谱，但你可以看到使用牛顿迭代法后这个值很快就趋近于根号2了：

(       4  + 2/   4     ) / 2 = 2.25
(    2.25  + 2/   2.25  ) / 2 = 1.56944..
( 1.56944..+ 2/1.56944..) / 2 = 1.42189..
( 1.42189..+ 2/1.42189..) / 2 = 1.41423..
….

       
    这种算法的原理很简单，我们仅仅是不断用(x,f(x))的切线来逼近方程x^2-a=0的根。根号a实际上就是x^2-a=0的一个正实根，这个函数的导数是2x。也就是说，函数上任一点(x,f(x))处的切线斜率是2x。那么，x-f(x)/(2x)就是一个比x更接近的近似值。代入f(x)=x^2-a得到x-(x^2-a)/(2x)，也就是(x+a/x)/2。
    同样的方法可以用在其它的近似值计算中。Quake III的源码中有一段非常牛B的开方取倒函数。

牛顿迭代公式

编辑

设r是

  

的根，选取

  

作为r的初始近似值，过点

  

做曲线

  

的切线L，L的方程为

  

，求出L与x轴交点的横坐标

  

，称x1为r的一次近似值。过点

  

做曲线

  

的切线，并求该切线与x轴交点的横坐标

  

，称

  

为r的二次近似值。重复以上过程，得r的近似值序列，其中，

  

称为r的

  

次近似值，上式称为牛顿迭代公式。

用牛顿迭代法解非线性方程，是把非线性方程

  

线性化的一种近似方法。把

  

在点

  

的某邻域内展开成泰勒级数

  

，取其线性部分（即泰勒展开的前两项），并令其等于0，即

  

，以此作为非线性方程

  

的近似方程，若

  

，则其解为

  

， 这样，得到牛顿迭代法的一个迭代关系式：

  

。

已经证明，如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛。 并且，如果不为0, 那么牛顿法将具有平方收敛的性能. 粗略的说，这意味着每迭代一次，牛顿法结果的有效数字将增加一倍。[1] 

迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程，跟迭代法相对应的是直接法（或者称为一次解法），即一次性解决问题。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令（或一定步骤）重复执行，在每次执行这组指令（或这些步骤）时，都从变量的原值推出它的一个新值。

利用迭代算法解决问题，需要做好以下三个方面的工作：

一、确定迭代变量

在可以用迭代算法解决的问题中，至少存在一个可直接或间接地不断由旧值递推出新值的变量，这个变量就是迭代变量。

二、建立迭代关系式

所谓迭代关系式，指如何从变量的前一个值推出其下一个值的公式（或关系）。迭代关系式的建立是解决迭代问题的关键，通常可以使用递推或倒推的方法来完成。

三、对迭代过程进行控制

在什么时候结束迭代过程？这是编写迭代程序必须考虑的问题。不能让迭代过程无休止地执行下去。迭代过程的控制通常可分为两种情况：一种是所需的迭代次数是个确定的值，可以计算出来；另一种是所需的迭代次数无法确定。对于前一种情况，可以构建一个固定次数的循环来实现对迭代过程的控制；对于后一种情况，需要进一步分析得出可用来结束迭代过程的条件。

示例

编辑

欧几里德算法

最经典的迭代算法是欧几里德算法，用于计算两个整数a,b的最大公约数。其计算原理依赖于下面的定理：

定理：gcd(a,b) = gcd(b,a mod b)

证明：a可以表示成a = kb + r，则r = a mod b。假设d是a,b的一个公约数，则有 a%d==0,b%d==0，而r = a - kb，因此r%d==0 ，因此d是(b,a mod b)的公约数

同理，假设d 是(b,a mod b)的公约数，则 b%d==0,r%d==0 ，但是a = kb +r ，因此d也是(a,b)的公约数。

因此(a,b)和(b,a mod b)的公约数是一样的，其最大公约数也必然相等，得证。

欧几里德算法就是根据这个原理来做的，欧几里德算法又叫辗转相除法，它是一个反复迭代执行，直到余数等于0停止的步骤，这实际上是一个循环结构。其算法用C语言描述为：

1

2

3

4

5

6

7

8

9

10

11

12

13

int Gcd_2(int a,int b)/*欧几里德算法求a,b的最大公约数*/

{

if (a<=0 || b<=0)/*预防错误*/

return 0;

int temp;

while (b > 0)/*b总是表示较小的那个数，若不是则交换a，b的值*/

{

temp = a % b;/*迭代关系式*/

a = b;

b = temp;

}

return a;

}

从上面的程序我们可以看到a，b是迭代变量，迭代关系是temp = a % b；根据迭代关系我们可以由旧值推出新值，然后循环执a = b; b = temp；直到迭代过程结束（余数为0）。在这里a好比那个胆小鬼，总是从b手中接过位置，而b则是那个努力向前冲的先锋。

斐波那契数列

还有一个很典型的例子是斐波那契（Fibonacci）数列。斐波那契数列为：0、1、1、2、3、5、8、13、21、…，即 fib⑴=0; fib⑵=1;fib(n)=fib(n-1)+fib(n-2) （当n>2时）。

在n>2时，fib(n)总可以由fib(n-1)和fib(n-2)得到，由旧值递推出新值，这是一个典型的迭代关系，所以我们可以考虑迭代算法。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

int Fib(int n) //斐波那契（Fibonacci）数列

{

if (n < 1）/*预防错误*/

return 0;

if (n == 1 || n == 2)/*特殊值，无需迭代*/

return 1;

int f1 = 1,f2 = 1,fn;/*迭代变量*/

int i;

for(i=3; i<=n; ++i)/*用i的值来限制迭代的次数*/

{

fn = f1 + f2; /*迭代关系式*/

f1 = f2;//f1和f2迭代前进，其中f2在f1的前面

f2 = fn;

}

return fn;

}

C语言代码

编辑

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

double func(double x) //函数

{

    return x*x*x*x-3*x*x*x+1.5*x*x-4.0;

}

double func1(double x) //导函数

{

    return 4*x*x*x-9*x*x+3*x;

}

int Newton(double *x,double precision,int maxcyc) //迭代次数

{

    double x1,x0;

    int k;

    x0=*x;

    for(k=0;k<maxcyc;k++)

    {

        if(func1(x0)==0.0)//若通过初值，函数返回值为0

        {

            printf("迭代过程中导数为0!\n");

            return 0;

        }

        x1=x0-func(x0)/func1(x0);//进行牛顿迭代计算

        if(fabs(x1-x0)<precision || fabs(func(x1））<precision) //达到结束条件

        {

            *x=x1; //返回结果

            return 1;

        }

        else //未达到结束条件

        {

            x0=x1; //准备下一次迭代

        }

    }

    printf("迭代次数超过预期！\n"); //迭代次数达到，仍没有达到精度

    return 0;

}

 

int main()

{

    double x,precision;

    int maxcyc;

    printf("输入初始迭代值x0:");

    scanf("%lf",&x);

    printf("输入最大迭代次数：");

    scanf("%d",&maxcyc);

    printf("迭代要求的精度：");

    scanf("%lf",&precision);

    if(Newton(&x,precision,maxcyc)==1） //若函数返回值为1

    {

        printf("该值附近的根为：%lf\n",x);

    }

    else //若函数返回值为0

    {

        printf("迭代失败！\n");

    }

    getch();

    return 0;

}

C++代码

编辑

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

//此函数是用来求一元3次方程ax^3+bx^2+cx+d=0的解

//比如 x^3-27=0，我们就可以输入1 0 0 -27，这样我们就可以得到一个解

#include<iostream>

#include<cmath>

using namespace std;

int main()

{

double diedai(double a,double b,double c,double d,double x);

double a,b,c,d;

double x=10000.0;

cout<<"请依次输入方程四个系数：";

cin>>a>>b>>c>>d;

x=diedai(a,b,c,d,x);

cout<<x<<endl;

return 0;

}

double diedai(double a,double b,double c,double d,double x)

{

while(abs(a*x*x*x+b*x*x+c*x+d)>0.000001)

{

x=x-(a*x*x*x+b*x*x+c*x+d)/(3*a*x*x+2*b*x+c);

}

return x;

}

可以得到一元3次方程3个解的程序(原创，超优化)：

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

#include<iostream>

#include<vector>

using namespace std;

 

vector<double >v;//stl vector链型数组

vector<double >::iterator it;//vector迭代器

 

int x0=5;

 

double a,b,c,d;

 

double abs(double y){    while(y<0) y=-y;    return y;}

 

double f(double x){    return a*x*x*x+b*x*x+c*x+d;}

 

double fd(double x){    return 3*a*x*x+2*b*x+c;}

 

bool u;//用来判断是否重复

 

void newton(int a1,int b1,int c1,int d1)

{   

     for(x0=-5000;x0<=5000;x0++)//在一个大区域中逐个点用牛顿法，可找出大多数3次方程所有根    

     {          

          double x1=x0;          

          while(abs(f(x1))>0.001)         

          {                  

            double x=x1;                 

             x1=x-f(x)/fd(x);          

           }          

             for( it=v.begin();it!=v.end();it++)          

             {         

                 if(abs((*it-x1))<0.01)  {u=1; break;}        

             }          

              if(u!=1&&x1<1000000000)          

              {         

                      cout<<x1<<" ";            

                      v.push_back(x1);//把已得到的解添加到vector，用于防止重复         

              }          

              u=0;   

        }

 }

 

int main(){    

     cin>>a>>b>>c>>d;     

     newton(a,b,c,d);

}