2
votes

So I tried to plot a contour in Julia by interpolating a 2D function, using the following code:

using Interpolations
using Plots
gr()

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

# Interpolate the function
g_int = interpolate(g, BSpline(Quadratic(Line(OnCell()))))

# Scale the interpolated function to the correct grid 
gs_int = scale(g_int, xs, ys)

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

# Compare the real value and the interpolated value of the function at an arbitrary point
println("gs_int(3.2, 3.2) = ", gs_int(3.2, 3.2))
println("f(3.2, 3.2) = ", f(3.2, 3.2))

# Contour of the interpolated plot
p1 = contour(xs, ys, gs_int(xs, ys), fill=true)
# Real contour of the function
p2 = contour(xc, yc, f, fill=true)

plot(p1, p2)

And this obviously didn't give the correct contour, although the interpolation was seemingly correct:

a

The problem was fixed by transposing gs_int(xs, ys):

p1 = contour(xs, ys, gs_int(xs, ys)', fill=true)

Then I randomly generated some points in 2D space, and repeated the same procedures:

using DelimitedFiles
using Interpolations
using Plots
gr()

data = readdlm("./random_points.txt", Float64)

# Create a dictionary to test different orders of interpolations. 
inter = Dict("constant" => BSpline(Constant()), 
    "linear" => BSpline(Linear()), 
    "quadratic" => BSpline(Quadratic(Line(OnCell()))),
    "cubic" => BSpline(Cubic(Line(OnCell())))
)

x = range(-10, length=64, stop=10)
y = range(-10, length=64, stop=10)

v_unscaled = interpolate(data, inter["cubic"])
v = scale(v_unscaled, x, y)

# The contour of the data points
p0 = contour(x, y, data, fill=true)
display(p0)

# The contour of the interpolated function
p_int = contour(x, y, v(x,y)', fill=true)
display(p_int)

However the two contour plots don't look the same.

b

As I removed the apostrophe after v(x,y), this worked:

p_int = contour(x, y, v(x,y), fill=true)

c

Now I don't get it. When should I apply transposition, and when shouldn't I do so?

1

1 Answers

4
votes

That's because in your first example you plot a function, in the second example you plot two arrays. The two arrays don't need to be transposed as they are oriented the same way. But in the first example, the way you generate the array is transposed relative to the way Plots generates an array from the 2-d function you're passing.

When you plot a function, Plots will calculate the outcome as g = Float64[f(x,y) for y in ys, x in xs] not the other way around, like you did in your code. For a good discussion of transposes in plotting, again refer to https://github.com/JuliaPlots/Makie.jl/issues/205