在 Julia 中使用 Interpolations.jl 和 Dierckx.jl 插值函数

Question

我正在尝试使用 Julia 中的 Interpolations.jl 和 Dierckx.jl 包执行函数插值。 Interpolations.jl 的文档如下：

https://github.com/JuliaMath/Interpolations.jl/blob/master/doc/Interpolations.jl.ipynb

对于 Dierckx.jl：

https://github.com/kbarbary/Dierckx.jl

所以我尝试使用不同的函数来试验插值，例如： 一个简单的代码：

using Interpolations
xs = 0:5
f(x) = x^2 * abs(sin(3*x) + cos(x))
ys = f.(xs)
f_int = interpolate(ys, BSpline(Quadratic(Line(OnCell()))))
println("f(3.2) = ", f(3.2))
println("f_int(3.2) = ", f_int(3.2))

二次插值应该是相当准确的，但结果如下：

f(3.2) = 12.007644743861604
f_int(3.2) = 2.973832923722435

那么我对 Interpolations.jl 的功能有什么误解？ Interpolations.jl 中的interpolate 函数不接受数组xs 作为参数，而只接受ys，所以我认为这可能是由于我对xs 的“不正确”选择？

然后我切换到 Dierckx.jl，它在函数 Spline1D 和 Spline2D 中同时接受 xs 和 ys。在我看来，Spline1D 在上面的示例中运行良好，因为我将函数插值的行切换为：

f_int = Spline1D(xs, ys)

然而，当我尝试 2D 时，问题又出现了：

using Dierckx
xs = 1:5
ys = 1:8
g = Float64[(3x + y ^ 2) * abs(sin(x) + cos(y)) for x in xs, y in ys]
f(x) = (3x + y ^ 2) * abs(sin(x) + cos(y))
f_int = Spline2D(xs, ys, g)
println("f(3.2, 3.2) = ", f(3.2, 3.2))
println("f_int(3.2, 3.2) = ", f_int(3.2, 3.2))

结果：

f(3.2, 3.2) = -0.6316251447925815
f_int(3.2, 3.2) = 20.578758429637535

再说一遍，上面的代码有什么问题？我对这些包的功能有什么误解？

[编辑] 我试图绘制一个轮廓来比较由 Interpolations.jl 生成的插值 2D 函数和函数的实际轮廓，这会产生以下结果：

using Interpolations
using Plots
gr()

xs = 1:0.5:5
ys = 1:0.5:8
g = Float64[(3x + y ^ 2) for x in xs, y in ys]
f(x, y) = (3x + y ^ 2)

g_int = interpolate(g, BSpline(Quadratic(Line(OnCell()))))

gs_int = scale(g_int, xs, ys)

xc = 1:0.1:5
yc = 1:0.1:5

println("gs_int(3.2, 3.2) = ", gs_int(3.2, 3.2))
println("f(3.2, 3.2) = ", f(3.2, 3.2))

p1 = contour(xs, ys, gs_int(xs, ys), fill=true)
p2 = contour(xc, yc, f, fill=true)

plot(p1, p2)

通过随机选择（x，y）的一些值，目前插值函数似乎没有问题，但是为什么插值函数的等高线图看起来如此扭曲？

Answer 1

让我专注于您使用 Interpolations.jl 的尝试，因为它是一个纯 Julia 解决方案。

如您所料，您需要适当地缩放底层网格。这就像一个额外的函数调用一样简单（参见包文档中的scaled BSplines）：

using Interpolations
xs = 0:5
f(x) = x^2 * abs(sin(3*x) + cos(x))
ys = f.(xs)
f_int = interpolate(ys, BSpline(Quadratic(Line(OnGrid()))))
sf_int = scale(f_int, xs) # new: scale the interpolation to the correct x-grid
println("f(3.2) = ", f(3.2))
println("f_int(3.2) = ", f_int(3.2))
println("sf_int(3.2) = ", sf_int(3.2)) # new: printing of the result

通过此更改，您将获得

f(3.2) = 12.007644743861604
f_int(3.2) = 2.973832923722435
sf_int(3.2) = 7.353598413214446

这更接近，但仍然很糟糕。然而，原因很简单：输入数据不足以进行良好的插值。让我们想象一下。使用当前的输入数据，我们有以下情况：

现在让我们使用更精细的输入数据网格，xs = range(0,5,length=20)。有了这个改变，我们有

f(3.2) = 12.007644743861604
f_int(3.2) = 0.6113243320846269
sf_int(3.2) = 12.002579991274903

图形化

显然，插值现在能够捕获底层函数的大部分特征。