利用Julia语言的集成能力

Question

我选择在我的项目中使用Julia的主要原因之一是它的速度，尤其是计算积分方面。

我想要对一个定义在区间[a,b]上的一元函数f(x)进行积分。通常情况下，Julia的quadgk函数会是一个快速而准确的解决方案。然而，我没有f(x)的解析式，只有在[a,b]区间内一些离散点xi处的函数值f(xi)，这些函数值存储在一个数组中。xi是均匀地分布在[a,b]区间上，我可以让它们之间的距离尽可能小。

一种朴素的方法是简单地定义一个函数f，该函数使用f(xi)的值进行插值，并将其提供给quadgk（并使间距尽可能小）。然而，这样做我将不知道误差大小，这很遗憾，因为QuadGK会告诉您其估计量的误差。

另一种解决方案是编写自己的函数来积分数组（例如使用梯形法则），但那样做将打败使用Julia的目的...

如何使用Julia最简单地准确地积分仅给出离散值的函数呢？

- math_lover

---

只有离散值且没有其他关于函数的信息，你不能比梯形法更好。 - Oscar Smith

@OscarSmith：我是Julia的新手。目前，我更喜欢使用Julia中已经存在的集成函数而不是编写自己的函数。显然有一个名为QuadGK的函数，但同样的问题，我猜这只接受函数，而不是值数组？ - math_lover

当你说你可以把间距缩小到任意大小时，这表明你可以访问底层函数，是这样吗？还是从某些物理过程中测量得出的？如果你对底层信号（或者说函数）有任何了解，你可能能够改进trapz。它是否平滑？是否振荡？如果你实际上正在积分f(x)^2，Trapz会在振荡方面存在一些系统性问题。 - DNF

Answer 1

由于你只有数值，而不是函数本身，因此梯形法可能是你最好的选择。包Trapz提供了这个功能(https://github.com/francescoalemanno/Trapz.jl)。然而，我认为值得尝试一下自己编写一个相当不错的实现。

function trap(A)
    return sum(A) - (A[begin] + A[end])/2
end

这对于一组1000万个浮点数只需要2.9毫秒，如果它们是Int类型，那么也是2.9毫秒。即使它们是复数，它仍然可以正常工作（并且需要8.9毫秒）。

像这样的方法是一个很好的示例，展示了在Julia中编写相当快速但仍然完全通用的代码可以多么简单。

- Oscar Smith