在 Julia 中使用 Interpolations.jl 和 Dierckx.jl 插值函数
Interpolating functions using Interpolations.jl and Dierckx.jl in Julia
我正在尝试使用 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，所以我认为这可能是由于我 "incorrect" 选择 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 生成的插值二维函数和函数的实际等高线，这产生了以下结果：
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)的值，但是为什么插值函数的等值线图看起来那么扭曲？
让我专注于您使用 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
和图形
显然，插值现在能够捕获底层函数的大部分特征。
我正在尝试使用 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，所以我认为这可能是由于我 "incorrect" 选择 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 生成的插值二维函数和函数的实际等高线，这产生了以下结果：
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)的值，但是为什么插值函数的等值线图看起来那么扭曲？
让我专注于您使用 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
和图形
显然，插值现在能够捕获底层函数的大部分特征。