通过梯形规则的 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

_{注意:这个答案是用有文化的 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 ，甚至是列表理解。