通过梯形规则的Haskell数值积分导致错误的符号

Question

我编写了一些代码,用于使用梯形规则以数字方式集成函数.它有效,但它产生的答案有一个错误的标志.为什么会这样？

代码是:

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + partial_sum)
    where 
        h = (b - a) / 1000 
        most_parts  = map f (points (1000-1) h) 
        partial_sum = sum most_parts

points  :: Double -> Double -> [Double]
points x1 x2 
    | x1 <= 0 = []
    | otherwise = (x1*x2) : points (x1-1) x2

梯形规则

代码可能不够优雅,但我只是Haskell的学生,并希望首先处理当前的问题,然后编码风格很重要.

Answer 1

_{注意:这个答案是用文字Haskell编写的.将其保存.lhs为扩展名并将其加载到GHCi中以测试解决方案.}

找到罪魁祸首

首先,让我们来看看integration.在其当前形式中,它仅包含函数值的总和f x.即使这些因素目前不正确,整体方法也很好:您f在网格点进行评估.但是,我们可以使用以下函数来验证是否存在错误:

ghci> integration (\x -> if x >= 10 then 1 else (-1)) 10 15
-4.985

等一等.x在[10,15]中甚至不是负数.这表明您使用了错误的网格点.

网格点重新审视

即使你已经链接了这篇文章,让我们来看看梯形规则的示例性使用(公共领域,Oleg Alexandrov的原始文件):

虽然这不使用均匀网格,但我们假设6个网格点与网格距离是等距的h = (b - a) / 5.x这些点的坐标是什么？

x_0 = a + 0 * h (== a)
x_1 = a + 1 * h
x_2 = a + 2 * h
x_3 = a + 3 * h
x_4 = a + 4 * h
x_5 = a + 5 * h (== b)

如果我们使用set a = 10和b = 15(因此h = 1),我们应该最终使用[10, 11, 12, 13, 14, 15].我们来检查你的points.在这种情况下,您将使用points 5 1并最终使用[5,4,3,2,1].

而且有错误.points不尊重边界.我们可以很容易地通过解决这个问题pointsWithOffset:

> points  :: Double -> Double -> [Double]
> points x1 x2 
>     | x1 <= 0 = []
>     | otherwise = (x1*x2) : points (x1-1) x2
>
> pointsWithOffset :: Double -> Double -> Double -> [Double]
> pointsWithOffset x1 x2 offset = map (+offset) (points x1 x2)

这样一来,我们仍然可以使用当前的points定义,从生成网格点x1来0(几乎).如果我们使用integration与pointsWithOffset我们结束了

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + partial_sum)
    where 
        h = (b - a) / 1000 
        most_parts  = map f (pointsWithOffset (1000-1) h a)  
        partial_sum = sum most_parts

捆绑松散的目的

但是,这并没有考虑到您在梯形规则中使用了两次所有内点.如果我们添加因子,我们最终会

> integration :: (Double -> Double) -> Double -> Double -> Double
> integration f a b = 
>      h / 2 * (f a + f b + 2 * partial_sum)
>    --    ^^^              ^^^
>    where 
>        h = (b - a) / 1000 
>        most_parts  = map f (pointsWithOffset (1000-1) h a)  
>        partial_sum = sum most_parts

这为我们的测试函数提供了正确的值.

行使

您当前的版本仅支持1000网格点.添加一个Int参数,以便可以更改网格点的数量:

integration :: Int -> (Double -> Double) -> Double -> Double -> Double
integration n f a b = -- ...

此外,尝试写points以不同的方式,例如从去a到b,使用takeWhile和iterate,甚至一个列表理解.