在 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，它接受函数and中的xsand 。在我看来，在上面的示例中效果很好，因为我将函数插值行切换为：ysSpline1DSpline2DSpline1D

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 生成的插值二维函数和函数的实际轮廓，这产生了以下结果：

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 解决方案。

正如您所预料的，您需要适当地缩放底层网格。这就像一个额外的函数调用一样简单（请参阅包文档中的缩放 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

并以图形方式

显然，插值现在能够捕获基础函数的大多数特征。