Интерполяция функций с использованием Interpolations.jl и Dierckx.jl в Julia

Question

Я пытаюсь выполнить интерполяцию функции с помощью пакетов Interpolations.jl и Dierckx.jl в Julia.Документ 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?Функция interpolate в Interpolations.jl не принимает массив xs в качестве аргумента, а только ys, поэтому я думаю, что это возможно из-за моего "неправильного" выбора xs?

Тогда япереключился на Dierckx.jl, который принимает и xs и ys в функциях Spline1D и Spline2D.Мне показалось, что 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

Итак, что же не так с кодом выше?Что я неправильно понял о функциях этих пакетов?

[Редактировать] Я попытался построить контур для сравнения интерполированной 2D-функции, созданной 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, поскольку это решение чисто Джулии.

Как вы и ожидали, вам необходимо соответствующим образом масштабировать базовую сетку.Это так же просто, как один дополнительный вызов функции (см. масштабированные линии 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

и графически

Очевидно, что теперь интерполяция способна охватить большинство характеристикосновная функция.