F#实现二分法求方程的根

      在不少工程问题中，我们需要求某一个方程f(x)的根，对于一元一次方程或者二次方程，我们可以利用公式非常方便的求出根。而对于高次方程或非线性方程，求根往往就没有现成的求根公式可以套用，而需要借助数值方法来进行求解。一般来说，对于一个连续的函数来说，如果存在两个点a和b，有f(a)和f(b)异号，那么理论上就存在一个根介于[a,b]之间，设为c，那么f(c)=0。举例来说，如果求x^3-1 这个函数，由于f(0) = -1 ，而f(2)=7,即f(0)*f(2) <0 ,那么说明在[0,2]之间必定存在一个根，根据图形，结合经验，可以判断出根为1。

    那么这个根如何通过计算机程序来实现呢？那么首先需要了解一下求根算法：

1 二分法（Bisection Method）

      在数据结构和算法中，有一种算法也叫二分法，他可以用于快速查找。类似的，二分法也可以用于方程求根上。在https://x-engineer.org/官网上，有一篇解说二分法的文章，感觉非常详细，其中的核心过程如下图所示：

      在二分法当中，每一次迭代，搜索区间的长度就会缩减一半，这是一个指数级的缩减。因此，即使给定的初始搜索区间比较大，但经过二分法多次迭代后也可以迅速缩减，得到一个非常精准的结果。它和泰勒级数一样，除了能得到一个足够精确的值之外，还能得到误差的范围。详细的过程可参考:
 https://x-engineer.org/undergraduate-engineering/advanced-mathematics/numerical-methods/the-bisection-method-for-root-finding

2 二分法实现

      根据上面的二分法算法过程，这里选择F#语言进行二分法求根程序实现。通过总结多种二分法算法描述，将二分法算法总结如下:
第一步：给出一个初始搜索区间[a,b]和精确度eps，如0.0001，使得f(a)*f(b)< 0 。
第二步：计算中间值 c=(a+b)/2 。
第三步：判断闭区间的长度b-a和精确度eps的关系，当闭区间长度>e的时候，返回第二步计算c。当闭区间长度<e的时候，停止迭代。输出结果c。
第四步：求出f(c)的值，如果f(a)*f(c)< 0 ，那么在[a,c]区间进行搜索；否则在[c,b]区间进行搜索。
     这个过程F#代码如下：

    let rec roots ( f , x1 , x2 , eps ) = 
        let  m = (x1 + x2) / 2.
        if (abs(x2 - x1 )< eps ) then
            m
        else 
            if ( f(x1) * f(m) < 0. ) then
                 roots(f,x1,m,eps)
            else
                roots(f,m,x2,eps)

3 求根验证

      为了验证二分法算法是否可以正确工作，下面给出函数进行验证:

    let testRoots = 
        let f x = x ** 2. - 4.
        let r = roots(f,-3.,3.,0.001) 
        printfn "%O" r //-1.9998779296875

        let f x = x ** 2. - 4.
        let r = round(roots(f,0.,3.,0.001))
        printfn "%O" r //2

        let f x = sin(x ** 2.) - 0.7
        let r = roots(f,-3.,3.,0.001) 
        printfn "%O" r //0.8807373046875

     从上述测试可知，x ^2 - 4 = 0  的根，如果在[-3,3]上虽然有2个根，但是程序只能求出一个根 -1.9998779296875 ，理论解析解为 - 2 。如果将搜索范围给定到[0,3]，那么可以求出另外一个根，这里用round函数进行处理，四舍五入到 2 。最后一个sin(x^2) - 0.7 = 0 的根，在[-3,3]区域求出一个根0.8807373046875 。实际上，通过绘图可知，在这个区间会有多个根:

     因此，具体的根，需要结合图形来给定初始搜索空间，这往往决定着求出的是哪个根。