前言

在阅读本文之前，请先看一道题：求解x，使得 $e^{-x}-x = 0$

如果你能用纸笔解出这道题的话，请务必留言告知我方法，谢谢。

如果你不能用纸笔解出这道题的话，那么我们可以完全借助代码来对其进行求解。

以下就是我写的求解器，它几乎能够对所有的一元函数进行自动求解。

你可以在输入框内尝试任意一元函数，按Enter键或点击求解按钮进行求解。如有bug请留言，谢谢。

它的原理有一点复杂，下面将会一一讲到。

牛顿法

牛顿法是求解的核心方法，它的维基百科的定义为：

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

简言之，牛顿法就是对x进行迭代，直至x收敛于某个很小的范围。

其迭代的公式为：$$x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}$$ 其中${f'(x_n)}$为${f(x_n)}$的导函数。

要理解牛顿法，你就要了解泰勒展开。牛顿法其实是利用了一阶的泰勒展开：$$f(x) = f(x_0)+(x -x_0)\cdot f'(x_0)$$ 对上述方程 $f(x) = 0$ 进行求解，得到： $$x = x_0 – \frac {f(x_0)}{f'(x_0)}$$ 这就是牛顿迭代公式的逼近原理。

你可以通过这幅图直观的看到x逼近的过程：

所以，对于任意的一元函数，我们都可以尝试用牛顿法来求得其近似解，一个简单伪代码如下：

声明 x0, x1, i
//当误差小于10^-9时，或者迭代步数超过10^5时，迭代结束
while (Math.abs(x0 - x1) > 10E-9) {
    //进行迭代
    x0 = x1
    x1 = x0 - (f(x0) / f'(x0))
    //记步
    ++ i
    //若步数过长，终止迭代
    if( i < 10E5) {
        break
    }
}

以下，我们将根据牛顿法，来构建自动求解器。

构建求解器

在构建求解器时，有几个关键问题需要解决：解析输入的表达式，表达函数，求导函数方程，对函数进行代入求值。

其中，最优先的一个问题是：我们怎么储存（表达）函数?

《SICP》的第一章提到了一种表示法：二分表达树。

以下是表达式((5 + z) / -8) * (4 ^ 2)的二分表达树的图例：

为什么选择这种表达方式呢？

主要是因为它是树形结构，方便递归处理节点，而我们之后求导函数其实就是用的递归思路，包括代入求值也是递归的思路。

预处理表达式

首先，我们需要预处理输入的表达式字符串。因为在数学中有一些简略或者多余的写法需要在此规范化。

例如：

• 5x表示5*x，通常都省略了乘号
• .5+x表示0.5+x，省略了小数点前的0
• -x中的-其意思是取负操作，是一个单目运算符。为了不与减法冲突，我们也需要将其预处理为(0-x)
• +5-x的加号完全没有意义，我们需要将其省去

比如，输入串 .5e^-x-2x，经过预处理后将变成 0.5e^(0-x)-2x。

自然的输入串经过预处理后，就应该是一个中缀的表达式字符串，这是人类能够自然理解的表达式形式。

但是为了将表达式储存成二叉表达树，我们还需要将中缀表达式转换成后缀表达式。

由迪杰斯特拉提出的调度场算法可以解决这个问题。

调度场算法

调度场算法基本和我们在栈 递归 汉诺塔文中提到的利用栈来计算表达式的方法类似。

它用队列表达输出的后缀表达式，利用了栈来储存操作符和函数。

其详细算法过程为（摘录自维基百科）:

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

其分支逻辑比较繁琐，需要细细阅读和小心实践。

构建二叉表达式树

假设输入表达式为：(a+b) * (c * (d+e))，经过调度场算法，我们得到a b + c d e + * *的后缀表达式。

此时我们便可以利用后缀表达式的特点，快速的构建出一颗二叉表达树来。

构建的算法很简单：

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

整个过程如果用图来说明，会十分容易理解。

1. 首先我们遍历到a，它是操作数，入栈，b也是操作数，也入栈
2. 遍历到+，它是二元操作符，弹出两个节点，然后添加至左右节点，入栈
3. 接下来三个都是操作数，将它们都入栈
4. 然后遍历到+，同2
5. 接下来是*，同2
6. 仍然是*，同2；遍历完毕，我们得到如下的二叉表达树：

我们定义二叉表达树为如下：

/**二叉表达树*/
function ExpTree(key) {
    //根节点
    this.root = null;
    //未知量的表达形式，默认为x
    this.key = key || 'x';
    //中序遍历的字符串
    this.infix = '';
    //后续遍历的字符串
    this.postfix = '';
}

定义树的节点如下：

/**二叉表达树节点*/
function ExpTreeNode(key, left, right) {
    //节点值
    this.key = key;
    //左子树，默认为空
    this.left = left || null;
    //右子树，默认为空
    this.right = right || null;
}

求值

对二叉表达树进行代入求值的算法应该很容易就能想到。

利用二叉树的递归特性，根为操作符或函数，左子树右子树是递归定义。我们只需要将左右子树的值递归求出，然后在进行操作符运算即可。

代码如下：

/**递归计算表达式的值*/
ExpTree.prototype.cal = function(node, x) {
    if(!node) {
        return;
    }
    //有子树，进行递归运算
    if(node.left || node.right) {
        //进行操作符运算
        return this.op( node.key, this.cal(node.left, x), this.cal(node.right, x) );
    }
    //否则，直接返回节点值
    switch(node.key) {
        case 'x':
            return x;
        case 'e':
            return Math.E;
        case 'pi':
            return Math.PI;
        default:
            return node.key;
    }
};

/**执行原子计算*/
ExpTree.prototype.op = function(opr, left, right) {
    switch(opr) {
        case '+':
            return left + right;
        case '-':
            return left - right;
        case '*':
            if(left == 0 || right == 0) {
                return 0;
            }
            return left * right;
        case '/':
            return left / right;
        case '^':
            return Math.pow(left, right);
        case 'sin':
            return Math.sin(right);
        case 'cos':
            return Math.cos(right);
        case 'tan':
            return Math.tan(right);
        case 'cot':
            return 1 / Math.tan(right);
        case 'arcsin':
            return Math.asin(right);
        case 'arccos':
            return Math.acos(right);
        case 'arctan':
            return Math.atan(right);
        case 'ln':
            return Math.log(right);
        case 'log':
            return Math.log(right) / Math.log(left);
        default:
            //error!
            console.log(opr, 'error!');
            return 0;
    }
};

构建导函数树

至此，我们只剩下了求解导函数的步骤。这一步也是比较复杂的操作，因为导函数的规则实在是很多。

请参见导函数的维基百科的定义。

但是重点是，我们用什么方法来求解导函数，又该如何表达导函数呢？

首先，表达导函数应该用二叉表达树来进行表示，因为可以直接对其进行代入求值，而且二叉表达树具有递归的特性；

其次，由于二叉表达树的根节点总是操作符或函数的特性，左右子树也是表达式，我们可以用递归的思路来求解导函数。

主要思路为： 对于一个二叉表达树，根节点为操作符或函数，左右子树则为操作符或函数的运算的参数。

比如根节点为+操作，左右子树为a,b：

具体的表达式则为 $a + b$，我们对加法求导法则是$(a+b)’ = a’ + b’$。我们可以发现求导也可以递归表达。即上面那颗二叉树的导函数为根为+，左子树为a的导函数的树，右子树为b的导函数的树。

所以我们可以用递归来表示求导的过程，下面是完全规则的求导代码：

ExpTree.prototype.dao = function(node) {
    if(!node) {
        return;
    }
    var t;
    switch(node.key) {
        case '+':
        case '-':
            //(l+r)' = l' + r'
            //(l-r)' = l' - r'
            t = new ExpTreeNode(node.key, this.dao(node.left), this.dao(node.right));
            break;
        case '*':
            //(l*r)' = l'*r + l*r'
            t = new ExpTreeNode('+');
            t.left = new ExpTreeNode('*', this.dao(node.left), node.right);
            t.right = new ExpTreeNode('*', node.left, this.dao(node.right));
            break;
        case '/':
            //(l/r)' = (l'*r - l*r') / (r*r)
            t = new ExpTreeNode('/');
            t.left = new ExpTreeNode('-');
            t.left.left = new ExpTreeNode('*', this.dao(node.left), node.right);
            t.left.right = new ExpTreeNode('*', node.left, this.dao(node.right));
            t.right = new ExpTreeNode('*', node.right, node.right);
            break;
        case '^':
            t = new ExpTreeNode('*');
            if(node.right._containParam(node.right)) {
                //(l^r)' = (l^r) * (ln(l)*r)'
                t.left = node;
                t.right = this.dao(new ExpTreeNode('*', new ExpTreeNode('ln', null, node.left), node.right));
            }
            else {
                //(l^c)' = (l'*c) * l^(c-1)
                t.left = new ExpTreeNode('*', this.dao(node.left), node.right);
                t.right = new ExpTreeNode('^', node.left, new ExpTreeNode('-', node.right, new ExpTreeNode(1)));
            }
            break;
        case 'sin':
            //sin(r)' = cos(r)*r'
            t = new ExpTreeNode('*');
            t.left = new ExpTreeNode('cos', null, node.right);
            t.right = this.dao(node.right);
            break;
        case 'cos':
            //cos(r)' = 0-sin(r)*r'
            t = new ExpTreeNode('-', new ExpTreeNode(0));
            t.right = new ExpTreeNode('*', new ExpTreeNode('sin', null, node.right), this.dao(node.right));
            break;
        case 'tan':
            //tan(r)' = r' / (cos(r) * cos(r))
            t = new ExpTreeNode('/', this.dao(node.right));
            t.right = new ExpTreeNode('*', new ExpTreeNode('cos', null, node.right), new ExpTreeNode('cos', null, node.right));
            break;
        case 'cot':
            //cot(r)' = 0-r'/(sin(r)*sin(r))
            t = new ExpTreeNode('-', new ExpTreeNode(0));
            t.right = new ExpTreeNode('/', this.dao(node.right));
            t.right.right = new ExpTreeNode('*', new ExpTreeNode('sin', null, node.right), new ExpTreeNode('sin', null, node.right));
            break;
        case 'arcsin':
            //arcsin(r)' = r' / (1-r*r)^0.5
            t = new ExpTreeNode('/', this.dao(node.right))
            t.right = new ExpTreeNode('^', null, new ExpTreeNode(0.5));
            t.right.left = new ExpTreeNode('-', new ExpTreeNode(1), new ExpTreeNode('*', node.right, node.right));
            break;
        case 'arccos':
            t = new ExpTreeNode('/', this.dao(node.right))
            t.right = new ExpTreeNode('^', null, new ExpTreeNode(0.5));
            t.right.left = new ExpTreeNode('-', new ExpTreeNode(1), new ExpTreeNode('*', node.right, node.right));
            t = new ExpTreeNode('-', new ExpTreeNode(0), t);
            //arccos(r)' = 0 - (r' / (1-r*r)^0.5)
            break;
        case 'arctan':
            t = new ExpTreeNode('/', this.dao(node.right))
            t.right = new ExpTreeNode('+', new ExpTreeNode(1), new ExpTreeNode('*', node.right, node.right));
            //arcsin(r)' = r' / (1+r*r)
            break;
        case 'ln':
            //ln(r)' = r'/r
            t = new ExpTreeNode('/', this.dao(node.right), node.right);
            break;
        case 'log':
            //log(l,r) = r'/(r*ln(l))
            t = new ExpTreeNode('/', this.dao(node.right));
            t.right = new ExpTreeNode('*', node.right, new ExpTreeNode('ln', null, node.left));
            break;
        case 'x':
            //x' = 1
            t = new ExpTreeNode(1);
            break;
        default:
            //常量的导数为0
            t = new ExpTreeNode(0);
            break;
    }
    return t;
};

最终计算

我们把导函数表达成了二叉表达树的形式，至此，我们已经把所有问题全部解决。

求解表达式解的过程将十分清晰：

1. 将输入的表达式构建成表达式树t，同时构建其导函数的表达式树dt
2. 利用牛顿法，进行迭代求解

开放源码

源码地址为：https://github.com/iWoz/DSinJS/blob/master/ExpressionTree/ExpTree.js。

源码集成在JavaScript的算法与数据结构系列中，托管于GitHub:

Demo托管于JsFiddle。

欢迎讨论！