牛顿法求解一元函数

牛顿法求解一元函数

对于一个简单的一元方程我们可以通过代数运算来求解，但是对于一个非线性的复杂一元函数例如−2x5−8x2+sin(x∗x)−2x=0$-2x^5-8x^2+sin(x\ast x)-2^x=0$这样的方程，想要通过人力计算就很难办到。
下面介绍利用牛顿法来构建的一个一元函数方程求解的程序。

牛顿法

当方程没有求根公式或者求根公式特别复杂时，利用牛顿法都可以进行迭代求解。
维基对牛顿法的定义为：

它是一种在实数域和复数域上近似求解方程的方法。方法使用函数f(x)的泰勒级数的前面几项来寻找方程 f(y)=0的根。

先理解下泰勒级数
如果一个函数是以实数或者复数为变量，并且该函数是可导的，那么这个函数就可以被展开为泰勒多项式。
有如下定理：

设 n 是一个正整数。如果定义在一个包含 a 的区间上的函数 f 在 a 点处 n+1 次可导，那么对于这个区间上的任意 x，都有：
f(x)=f(a)+f′(a)1!(x−a)+f(2)(a)2!(x−a)2+...+f(n)(a)n!(x−a)n+Rn(x)$f(x)=f(a)+\frac{f^{\prime }(a)}{1!}(x-a)+\frac{f^{(2)}(a)}{2!}(x-a{)}^2+...+\frac{f^{(n)}(a)}{n!}(x-a{)}^n+R_n(x)$
其中的多项式称为函数在a 处的泰勒展开式，剩余的Rn(x)$R_n(x)$ 是泰勒公式的余项，是(x−a)n$(x-a{)}^n$的高阶无穷小。

牛顿法就是利用一阶泰勒展开f(x)=f(a)+(x−a)⋅f′(a)$f(x)=f(a)+(x-a)\cdot f^{\prime }(a)$，对f(x)=0$f(x)=0$求解得到x=a−f(a)f′(a)$x=a-\frac{f(a)}{f^{\prime }(a)}$，利用其对x$x$进行不断的迭代，直到x$x$收敛到正确的解或者逼近正确的解。
迭代公式为：xn+1=xn−fxnf′(xn)$x_{n+1}=x_n-\frac{f{x}_n}{f^{\prime }(x_n)}$
用一段伪代码表示为：

while x<0.000000001 //精度
    x0 = x1
    x1 = x0 - F(x0)/Fdao(x1)

以上为求解方程的核心算法。

下面是对函数表达式进行处理，包括表达式的预处理、前缀表达式转化后缀表达式、表达式计算和表达式的求导及计算。

表达式的预处理

对于用户输入的表达式，我们要转化成一定的格式方便我们后续处理。
例如：

• 去除表达式中夹杂的空格
• +2x-1省略其前面的+
• 2x转化为2*x
• -x转化为(0-x)
  经过简单的预处理得到一个中缀表达式。

中缀表达式转化为后缀表达式

转化过程我们可以利用栈来作为临时容器，用队列输出后缀表达式。具体算法使用迪杰斯特拉提出的调度场算法，具体算法步骤如下（摘自维基）：

• 当还有记号可以读取时：
  
  • 读取一个记号。
  • 如果这个记号表示一个数字，那么将其添加到输出队列中。
  • 如果这个记号表示一个函数，那么将其压入栈当中。
  • 如果这个记号表示一个函数参数的分隔符（例如，一个半角逗号 , ）：
    
    • 从栈当中不断地弹出操作符并且放入输出队列中去，直到栈顶部的元素为一个左括号为止。如果一直没有遇到左括号，那么要么是分隔符放错了位置，要么是括号不匹配。
  • 如果这个记号表示一个操作符，记做o1，那么：
    
    • 只要存在另一个记为o2的操作符位于栈的顶端，并且
      
      • 如果o1是左结合性的并且它的运算符优先级要小于或者等于o2的优先级，或者
      • 如果o1是右结合性的并且它的运算符优先级比o2的要低，那么将o2从栈的顶端弹出并且放入输出队列中（循环直至以上条件不满足为止）；
    • 然后，将o1压入栈的顶端。
  • 如果这个记号是一个左括号，那么就将其压入栈当中。
  • 如果这个记号是一个右括号，那么：
    
    • 从栈当中不断地弹出操作符并且放入输出队列中，直到栈顶部的元素为左括号为止。
    • 将左括号从栈的顶端弹出，但并不放入输出队列中去。
    • 如果此时位于栈顶端的记号表示一个函数，那么将其弹出并放入输出队列中去。
    • 如果在找到一个左括号之前栈就已经弹出了所有元素，那么就表示在表达式中存在不匹配的括号。
• 当再没有记号可以读取时：
  
  • 如果此时在栈当中还有操作符：
    
    • 如果此时位于栈顶端的操作符是一个括号，那么就表示在表达式中存在不匹配的括号。
    • 将操作符逐个弹出并放入输出队列中。
• 退出算法。

假如输入的表达式为(a+b) * (c * (d+e))，经过转化得到a b + c d e + * *这样的后缀表达式。

表达式的计算

首先要解决一个问题，通过哪种数据结构来存储表达式，表达式的计算是一个不递归环的过程，我们可以用“栈+循环”来处理，但为了方便递归处理和后序的递归求导，这里用二叉树来存储后缀表达式。
下面是(a+b) * (c * (d+e))的表达式树：

该树的特点是：叶子节点为操作数，父节点为操作符，可以从叶子节点开始一直向上递归最终求解出整个表达式的结果。

构建表达式树

利用后缀表达式很容易的得出表达式树，其构建的算法如下：

• 新建一个栈，该栈保存的是一个表达树
• 从头到尾遍历后缀表达式，
  
  • 当发现操作数时，我们将其存为操作数节点，入栈；
  • 当发现操作符时，我们将其存为操作符节点，
    
    • 若是二元操作符，我们将栈顶的两个节点弹出，作为该操作符的节点左右节点，然后将该子树入栈
    • 若是一元操作符，我们将栈顶节点弹出，作为该操作符的右节点，然后将该子树入栈

最终栈底只有一个表达式树，其他情况均视为表达式错误。

求值

以根结点为操作符，左右子树为操作数，可以对二叉树进行递归求解，最终得到一个含有未知数x$x$的简化表达式。

表达式树的求导和计算

对于上面得到的二叉表达树，由于二叉表达树根结点总是操作符，左右子树总是表达式，因此可以对二叉表达树进行递归求导。求导的过程也是一个递归的构建导函数树的过程。
例如：一个子树是这样的

设计一个导函数dao(+,a,b),利用导函数规则(a+b)’=a’+b’,对该子树求导的伪代码如下：

dao(+,a,b)
   return new root(+,dao(a),dao(b))

上面是只有“+”一种情况，在具体代码中我们需要根据所有符号的导函数规则分情况处理。
最终，将导函数树进行递归计算。

至此我们得到一个表达式t$t$，一个表达式的导dt$dt$。
利用牛顿法迭代处理，最终求解出未知数。

---

完整代码的Github地址。
最终的运行效果地址点此链接
如有错误或者改进的地方，欢迎指出。