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

Trapezoidal rule

代码可能不够优雅，但我只是Haskell的学生，想先解决当前的问题，再考虑编码风格。

Answer 1

是的，确实是分数加上你有一些错误的因素（inner 分数乘以 2）——这是你的代码的固定版本：

integration :: (Double -> Double) -> Double -> Double -> Double
integration f a b = h * (f a + f b + innerSum) / 2
    where 
        h = (b - a) / 1000 
        innerPts  = map ((2*) . f . (a+)) (points (1000-1) h) 
        innerSum = sum innerPts

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

它给出了合理的近似值（到 1000 点）：

λ> integration (const 2) 1 2
2.0
λ> integration id 1 2
1.5

λ> integration (\x -> x*x) 1 2
2.3333334999999975

λ> 7/3
2.3333333333333335

Answer 2

_{注意：此答案是用文字书写的Haskell。用 .lhs 作为扩展名保存并加载到 GHCi 中以测试解决方案。}

找到罪魁祸首

首先我们来看一下integration。在其当前形式中，它仅包含函数值 f x 的总和。尽管目前因素不正确，但总体方法很好：您在网格点评估 f。但是，我们可以使用以下函数来验证是否有问题：

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

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

重新访问网格点

尽管您已经链接了这篇文章，但让我们看一下梯形规则 (public domain, original file by 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)

如果我们使用集合 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，甚至列表理解。

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

Haskell numerical integration via Trapezoidal rule results in wrong sign

haskell

numerical-methods

找到罪魁祸首

重新访问网格点

收尾

运动