F#实现Simpson's Rule求数值积分

       我们知道，微积分的求值是比较复杂的。一般来说，求积分有定积分和不定积分之分。不定积分需要求出具体的表达式，但被积函数非常复杂时，求解非常费劲，非常可能找不到原函数。而定积分给定了区间范围，可以利用数值方法，利用微分的本质原理来进行积分数值计算。当前，微积分的数值计算方法很多，比如常见的有如下算法：梯形规则，辛普森规则和龙贝格规则。

1 Simpson's Rule原理

      下面我们介绍一下辛普森规则如何求数值积分，积分的数值本质上就是被积函数f(x)下的面积，定积分f(x) 在[a,b]上的数值就是被积函数f(x)在x轴a和b围起来的面积。当按照x轴方向，分成非常多的份，每一个份∇x就是非常小的微分单位，比如0.00001，这样曲线的面积，就可以近似看成是梯形或者方形和三角形的和等。下面给出辛普森规则的推导公式:

    此公式也有其他的表现形式，但是基本原理相同。它也称作Simpson’s 1/3 rule 。公式来自: http://www.numericmethod.com/About-numerical-methods/numerical-integration/simpson-s-rule  。以上的公式，可以用如下示例解释（来自网址：https://assets.vividmath.com），这样更加容易理解：

2 F#Simpson's Rule实现

      结合如上的算法公式，用F#语言可实现数值积分的计算，代码如下：

module NumMethod
    let rec NInt f a b =
        let numSteps = 1000000
        let h = double(b - a) / (double(numSteps))
        let  mutable valueEven = double(0.)
        let  mutable valueOdd = double(0.)
        for i = 0 to numSteps do 
              if i % 2 <> 0 then
                    valueOdd <- valueOdd + f( a + h * double(i))
              else 
                    valueEven <-  valueEven + f(a + h * double(i))
        (f(a)+f(b)+4.0*valueOdd+2.0*valueEven)*(h/3.0)

     这里的微分步长numSteps可以根据需要进行调整，理论上，越大应该精确度越高，但是性能越差。算法的本质就是在区间[a,b]上去划分成N个小区域，然后计算每个小区域的面积，而小区域的面积的一边虽然是曲线f(x),但由于间距很小，可用直线近似替代。最后将每个小区域的面积相加即可完成数值积分的近似计算。

3 测试

      下面给定一些函数来进行定积分数值的计算测试，代码如下：

    let test = 
        let f1 x = 2. + cos(2. * sqrt(x))
        let z = NInt f1 0.1 2.0
        printfn "%f" z   //3.169783 //3.1697776595193
        /////////////////////////
        let f2 x = 2.* x
        let z = NInt f2 1. 2.0
        printfn "%f" z  //3.000004  // 3 

     由测试可知，在一定的精确度下，数值积分的结果与实际的积分结果还是比较接近的。