前言

  最近在看神书《SICP》，刚看了第一章，虽然有些难啃，但感觉确实啃得确实“香”。说不上醍醐灌顶，但应该也是受益匪浅了。书中介绍了一些关于计算机数值求解的一些问题，这里抽取一点求平方根的算法，做个总结，希望可以便人便己。

---

一、二分法求平方根

二分法大概比较简单的一种求解的方法，它理论基础是零点存在定理。即：

通过事先确定零点所在的区间范围[a,b]（区间范围端点值异号），通过不断缩小这个区间范围，达到逐渐逼近最终解的目的。而缩小区间范围的方法就是取[a,b]的中点m，判断f(m)的正负号与f(a)、f(b)的关系，改变a的值（f(b)*f(m)<0）或者b的值（f(a）*f(m)<0)为m,每次缩减一半的解空间。循环上述操作，直到最终的范围小于所给的精度范围。

上面说的是通用的求解方程根的方法，对于平方根来说，也是适用的。代码如下：

//sqrt with binary search 
//note:l less than h
double sqrtWithBiSear(double num,double l,double h){//define function f(x) = x*x -num;auto difFunc = [num](double x)->double { return x*x - num;};//make sure that there will be one root at leastif(difFunc(l) *difFunc(h) > 0){return -1;}const double LIMIT = 0.00001; while(l<=h){double m = (l+h)/2;double difM = difFunc(m);if(fabs(difM) < LIMIT){return m;}double difL = difFunc(l);double difH = difFunc(h);if(difL* difM <= 0) {h = m;}else if( difH * difM <= 0 ) {l = m;}else {return -1;}}return -1;
}

其实对于单调递增的函数来说还可以简化一点：

//sqrt with binary search 
double sqrtWithBiSear1(double n){double h = max(1.0,n);double l = min(1.0,n);const double LIMIT = 0.00001;while(1){double guess = (l+h)/ 2;double guessSquare = guess*guess;double dif = guessSquare - n;if(fabs(dif) < LIMIT){return guess;}if(dif > 0){h = guess;}else {l = guess;}}return -1;
}

对于二分法求解来说，它的优点是对函数f(x)的要求不太高，只需要其在根所在的范围区间内连续，且在这个范围内提供一个异号的[a,b]; 但是其缺点是不能求偶数重根，也不能求复数根。由于每次砍掉一半的空间，其算法的复杂度是O(lgN)级别的，收敛的速度不算太快。
 

二、牛顿法求平方根

由于很多方程其实并不存在求根公式，而且很多时候现实中的一些应用也并不需要数学意义上的精确解，为了得到方程的近似零点，牛顿等人提出了使用迭代的方式来求这个近似的零点。

方法的思路总的来说四个字，“以直代曲”，使用切线来代替曲线。如下图所示。

下图显示了迭代的过程

通过xk，xk+1，xk+2，… 这种不断迭代的方式最终逼近至x*(当然只需要计算到指定精度的解就行)。

具体的迭代方程可以通过列切线方程，也可以通过泰勒展开式获取，最终获得的迭代方程如下图所示。

对于具体的求平方根来说，将f(x)=x^2-a=0（a>0）带入上述公式，可得

 
有了迭代方程，代码编写就不是问题了。

double sqrtWithNewdonWay(double x, double guess){//define goodEnough? functionauto goodEnough = [](double x, double guess)->bool {const double LIMIT = 0.00001;if(fabs(x - guess*guess) <= LIMIT){return true;}return false;};while(!goodEnough(x, guess)){guess = (guess + x/guess) / 2; }return guess;
}

对于牛顿法来说，它在单根附近的收敛速度较快，算法逻辑也比较简单，而且可以求代数方程的重根、复根。但是牛顿法有个问题，就是其并不是时时收敛的，和具体的函数、选取的初值都有关系。关于收敛性的讨论不在这篇博客的讨论范围之内，有兴趣的童鞋可以查阅其他资料。

 

三、不动点法求平方根

不动点法求解方程的根，其基本原理是将f(x) = 0转化为 一个等价的x = g(x), 或者说x(n+1) = g(x(n)) 。通过上述方程得到逐渐逼近一个解x0,使得x0 = g（x0）。则f(x0) = 0，即x0是方程f(x)的零点。

上述表述的不明确？下面是从【5】中“偷”的图。

还不理解？

去看看【5】中的内容吧，里面还有说不动点的几何意义。适合小阅一番。

有一点注意的是，构造出来的迭代函数g(x)不一定收敛，需要满足一定的条件，才会收敛。具体的条件这里也就不说了（因为我也不是很懂╮(╯▽╰)╭），请各位学霸同学自行查阅资料，谢谢。

 
上面是通用的求解方程根的解法，对于具体的求解平方根来说，可以构造一个和上述牛顿法一样的迭代公式，即

事实证明，它是收敛的。（如果构造一个g（x） = a/x ,它就不是收敛的）

实现的代码如下：

//fixed points of function
double sqrtWithFixPointMethod(double x,double guess){//define improveGuessauto improveGuess = [](double x,double guess)->double {return (guess+x/guess)/2;       // iterative formula};//define goodEnough? functionauto goodEnough = [](double v1, double v2)->bool {const double LIMIT = 0.00001;if(fabs(v1-v2) < LIMIT) {return true;}return false;};while(!googEnough(guess,nextGuess)){guess = nextGuess;nextGuess = improveGuess(x,guess);}return nextGuess;
}

对于不动点求解方程来说，它的收敛性也是个问题。构造不同的迭代函数g(x)可能会有不同的收敛性，而且要求g(x)在区间[a,b]上的函数值也在此区间内，对给定的初值要求较高。关于不动点法收敛性更详细的讨论可以参考【6】中所述。

 

四、更抽象的方式

上面是三种求平方根的实现代码，在SICP中还介绍了如何用更抽象的方式来定义函数。这儿也用cpp实现了一部分，供君一览。

//abstract version of binary serach method
double biSearMethod(std::function<double (double)> f,double l, double h){const double LIMIT = 0.00001;while(l<=h){double mid= (l+h)/ 2;double f_m = f(mid);if(fabs(f_m) < LIMIT){return mid;}if( f_m > 0){h = mid;}else {l = mid;}}}

//f is transform function, guess is the first guess value
double fixedPointMethod(std::function<double (double)> f, double guess){// define goodEnoughFunc with lambda expressionauto goodEnoughFunc = [](double v1, double v2) -> bool {const double LIMIT = 0.00001;if(fabs(v1-v2) < LIMIT) {return true;}return false;};// define improveGuessFunc with lambda expressionauto improveGuessFunc = [=](double guess)->double {return f(guess);};double nextGuess = improveGuessFunc(guess);while(!goodEnoughFunc(guess,nextGuess)){guess = nextGuess;nextGuess = improveGuessFunc(guess);}return nextGuess;}

//Derivative functions
std::function<double(double)> deriva(std::function<double(double)> f) {const double dx = 0.00001;return [=](double x)->double { return ((f(x+dx) - f(dx))  / dx);};
}std::function<double(double)> newtonTransfrom(std::function<double(double)> g){return [=](double x)->double { return  x - ( g(x) / (deriva(g)) (x) );     //define:xn+1 = xn - g(x) / d(g(x))};
}double newtonMethod(std::function<double(double)> h,double guess) {return fixedPointMethod(newtonTransfrom(h),guess);
}

---

参考

【1】《Structure and Interpretation of Computer Programs》
【2】方程求根：二分法–不动点迭代–牛顿法–弦截法
【3】牛顿迭代法求平方根
【4】如何通俗易懂地讲解牛顿迭代法求开方（数值分析）？
【5】不动点迭代法—单变量非线性方程近似根matlab求解
【6】不动点迭代及其收敛性