算法基础 - 求解非线性方程(牛顿法)

1. 前言

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         前文介绍了如何使用“高斯消元法”求解线性方程组。

本文秉承有始有终的态度，继续介绍“非线性方程”的求解算法。

本文将介绍 2 个非线性方程算法：

• 牛顿迭代法。

• 二分迭代法。

牛顿迭代法（Newton’s method）又称为牛顿-拉夫逊方法（Newton-Raphson method），是拉夫逊和牛顿同时提出来的一种在实数域和复数域上近似求解方程的方法。

为何说是近似求解方程？

因为对于多数方程式，因不存在求根公式，或者说无法或很难找到标准的可以直接套用的模板公式。因而求精准解非常困难。
即使如牛顿大神提出的方法，也只是近似求解的算法，甚至需要满足某种收敛条件的方程式才能使用牛顿迭代法求解。

下面将具体介绍这 2 种求解算法。

2. 牛顿迭代算法

下面将通过一系列演示图，直观告诉大家牛顿迭代法的算法思想。算法中，牛顿用到了微积分相关的知识。

所以，在阅读下文时，需要具备微积分的认知。

牛顿迭代算法求解方程式的过程，有点类似福尔摩斯探案。通过蛛丝马迹，先合理的预测，然后根据推理逻辑，让预测离真相近一点、再近一点……一直到找到或接近真相。

实事告诉我们，不是所有的预测都能找到真相。同理，基于预测的牛顿迭代法也不一定总是能找到方程式的解，看完下面的演示流程，你将明白。

假设现有一非线性方程式 f(x)，其在平面坐标轴上的曲线图案如下。

所谓求解，指求其与横坐标轴相交时的点的 x 值.

现在，看看牛顿是如何使用微积分思想找到这个解的。只能说，牛逼人的思想非我等凡人能比拟。

2.1. 基本思想

• 在横坐标上找一 x0 点（也称预测点），并绘制 (x0,f(x0)) 点与曲线相交的切线。切线和横坐轴相交于 x1。

• 再绘制(x1,f(x1))点与曲线的切线，此时，切线与横轴相交于x2，继续绘制出(x2,f(x2))与曲线的交点……如此迭代，直到切线与横坐标轴的交点与曲线和横坐标的交点重合，此交合点便是曲线的解。是不是很简单，为什么是牛顿发现的，而不是我？

• x0的选择并不完全是任意的，也应该有基本的推理依据。预测点是关键，如果与真实值相差太远，则迭代次数会很大。理论上，只要预测点给的好，且此方程式满足牛顿迭代算法的前提条件，无论迭代多少次，解必能找出来，无非就是时间的长短。

2.2. 如何求解 x1

现在的问题转向到如何通过已知的 x0 值计算出 x1 的值？是否存在一个标准的公式？

现在就是微积分上场的时候，请屏住呼吸！真相将昭然若揭。

• 在 x0 和 x1 之间选择任一点 x ， 从此点向上绘制垂直线，假设与切线相交的位置的纵坐标值为y 。并绘制如下箭头所指的三角形：

• 三角形为直角三角形，学过三角函数的都知道，会存在如下的关系。

现在轮到微积分知识上场，它告诉我们，其中的 tanθ 就是切线与曲线的斜率。根据微积分原理，斜率即是x0在曲线上的导数，可以根据导函数计算出来，即：tanθ=f’(x0)。太完美了，如此公式可演变如下：

k = △y / △x , tanθ

继续化丽的转身后，它便如涅槃重生一样，破茧成如下人见人爱的模样：

• 因切线与横坐标轴相交的位置y=0，从而便可以求得 x1的值：

𝑥1=𝑥0 − 𝑓(𝑥0) / 𝑓′(𝑥0)

推导方式1：

依据上式： y = f(x0) + f’(𝑥0)(𝑥 - 𝑥0)

所以： y = f(x0) + f’(𝑥0)𝑥 - f’(𝑥0)𝑥0

又因为： 切线与横坐标轴相交的位置y=0，所以： 0 = f(x0) + f’(𝑥0)𝑥 - f’(𝑥0)𝑥0

所以： 𝑥 = ( f’(𝑥0)𝑥0 - f(x0) ) / f’(𝑥0) = 𝑥0 - f(x0) / f’(𝑥0)

所以： 𝑥1 = 𝑥0 - f(x0) / f’(𝑥0)

推导方式2：

因为：

△x = X0 - X1

△y = Y0 - Y1 = Y0

k = △y / △x = f'(X0)

所以：

△x = △y  /  f'(X0)

即：

X0 - X1 =  △y  /  f'(X0) = f(X0) /  f'(X0)

所以：

X1 =  X0 -  f(X0)  /  f'(X0)

也通常写作：

 
𝑥1 = 𝑥0 − 𝑓(𝑥0) / 𝑓′(𝑥0) 
 

或者 ：

𝑥1 = 𝑥0 − F(𝑥0) / F′(𝑥0) 

2.3. 如何求解 x2

同理，求得 x2的值：

𝑥2 = 𝑥0 − 𝑓(𝑥1) / 𝑓′(𝑥1)

最后，可以抽象出牛顿迭代公式，即迭代法中的核心子逻辑。

$Xn+1=Xn-f(Xn)/f^{\prime }(Xn)$

2.4. 收敛性

什么是牛顿迭代算法的收敛性？

通俗理解，选定预测点后，也许中间会有偏离，或许会忽远忽近，但无论如何最终都能靠近真实解，这便是收敛性。

换一句话而言，如果通过预测点，无法收敛到真实值，则无法求出解。

如果预测点为曲线的驻点，很不幸，由此点绘制的切线不会和横坐标轴相交，是无法求方程式的解。

另，如果收敛越来越远，也不能使用牛顿迭代法。如下图所示：

怎样的方程式能使用牛顿迭代法，牛顿迭代法已经给出了答案，可自行查阅一下。

2.5. 编码实现

现在来一个具体的案例：求解如下方程式。

f(x)=x⁴-3x³+1.5x²-4

牛顿迭代法中的子逻辑需要求解函数的导函数。受限于篇幅，导函数的推导在此不负赘。

这里仅给出常见的基本函数的导函数公式，再根据导函数生成法则，直接找到求解函数的导函数：

f‘(x)=4x³-9x²+3x

牛顿迭代法可以使用 递归 和 非递归 方案实现。

因初始值为预测值，从而可能导致递归或迭代的次数会很大。前文说过，牛顿迭代法并不是一个解非线性方程式的通用算法，也就是说使用牛顿迭代法可能得不到解。

故最好在编写算法时添加如下的辅助手段：

• 保证函数在整个定义域内最好是二阶可导的。（注意：这里的二阶可导）

• 预测点会影响计算量，可限制迭代的次数，当在此限制下不能得到结果时，则增加其它的判断手段试错。

2.5.1. 递归实现

#include <iostream>
#include <cmath>
using namespace std;

/* 原函数  f(x) = x^4 - 3x^3 + 1.5x^2 - 4 */
double yfun(double x)
{
  return pow(x,4) - 3 * pow(x,3) + 1.5 * pow(x,2) - 4.0;
}

/* 导函数 f'(x)=4*x^3-9*x^2+3x */
double dfun(double x)
{
  return 4 * pow(x,3) - 9 * pow(x,2) + 3 * x;
}

/* 递归实现牛顿迭代法
* val : 预测值
* precision : 精度(误差)
* deep : 递归深度
* */
double newtonIter(double val_esti, double precision, int deep)
{
  if(deep == 0) return -1;

  //求解
  double yRes = yfun(val_esti);
  if( yRes == 0.0 || fabs(yRes) < precision )
  {
    std::cout << "deep = " << deep << ", val_esti = " << val_esti << ", precision = " << precision << ", yRes = " << yRes << std::endl;
    return val_esti; //如果找到
  }

  //根据牛顿迭代公式修正 val 值
  val_esti = val_esti - yRes / dfun(val_esti);
  std::cout << "deep = " << deep << ", val_esti = " << val_esti << ", precision = " << precision << ", yRes = " << yRes << std::endl;

  //递归
  return newtonIter(val_esti, precision, deep -1);

}

/* 非递归实现 */
double newtonIter_(double val, double jd)
{
  double res = yfun(val);

  while( !(res == 0.0 || fabs(res) < jd))
  {//根据牛顿迭代公式修正 val 值
    val = val - yfun(val) / dfun(val);
    res = yfun(val);
  }

  return val;   //如果找到
}


int main()
{
  double val_esti = 0.0;
  int deep = 0;
  double precision = 0.0;

  std::cout << "请输入预测值：" << std::endl;
  std::cin >> val_esti;

  std::cout << "递归或迭代的最大次数：" << std::endl;
  std::cin >> deep;

  std::cout << "输入精度："<< std::endl;
  std::cin >> precision;

  double val_real = newtonIter(val_esti, precision, deep);
  std::cout << "val_real = " << val_real << std::endl;

  return 0;

}

测试： 当 x=0其倒数为 0，说明为驻点，不能做为预测值。

再试着把把 2代入导函数，其导数为 2，可以作为预测值试试：

正向测试一下，把 2.64894 代入原函数，可知是符合精度要求的近似值。

2.5.2. 非递归实现

//省略……
double newtonIter_(double val,double precision,int deep) 
{
	int i=0;
	double res=0;
	while( 1 ) {
		res=yfun(val);
		if( res==0.0 || fabs(res) <precision)
			return val;
		//根据牛顿迭代公式修正 val 值
		val=val-yfun(val)/dfun(val);
		i++;
		if(i==deep)
			return -1;
	}
	
	//没有
	return -1;
}