# My code for cubic spline interpolation plots different image everytime and not smooth

I wrote following code trying to implement a cubic spline interpolation (clamped spline). Every time I run it on jupyter notebook may plot different image. And the plotting is likely not very smooth, which is weird for cubic spline interpolation. ``````using Plots

function myspline( x::Array,y::Array,conval::Array, x_i::Float64 )
np = length(x) # the number of data points
ns = np - 1

h=Array{Float64}(undef,ns)
for i in 1:ns
h[i]=x[i+1]-x[i]
end

μ = Array{Float64}(undef,ns)
λ = Array{Float64}(undef,ns)
d = Array{Float64}(undef,np)

μ[ns] = 1
λ = 1
d = 6 * ( (y-y)/h - conval ) / h
d[np] = 6 * ( conval- (y[np]-y[ns])/h[ns]) / h[ns]

for j in 1:ns-1
λ[j+1] = h[j+1] / ( h[j] + h[j+1] )
μ[j] = 1- λ[j+1]
d[j+1] = 6 * ( (y[j+2]-y[j+1])/h[j+1] - (y[j+1]-y[j])/h[j] ) / ( h[j] + h[j+1] )
end

M = Array{Float64}(undef,np)

# solving tridiagonal matrix algorithm
α = Array{Float64}(undef,np)
β = Array{Float64}(undef,ns)
β = λ/2
for i in 2:ns
β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
end
y_ch = Array{Float64}(undef,np)
y_ch = d/2
for i in 2:np
y_ch[i] = ( d[i] - μ[i-1] * y_ch[i-1] )/( 2 - μ[i-1] * β[i-1] )
end
M[np] = y_ch[np]
for i in ns
M[np-i] = y_ch[np-i] - β[np-i] * M[np-i+1]
end

j_x_i = 1
while x[j_x_i+1] < x_i
j_x_i = j_x_i + 1
end
return M[j_x_i] * (x[j_x_i+1]-x_i)^3/(6*h[j_x_i]) + M[j_x_i+1] * (x_i-x[j_x_i])^3/(6*h[j_x_i]) + (y[j_x_i]-M[j_x_i]*h[j_x_i]^2/6)*(x[j_x_i+1]-x_i)/h[j_x_i] + (y[j_x_i+1]-M[j_x_i+1]*h[j_x_i]^2/6)*(x_i-x[j_x_i])/h[j_x_i]
end

# input: data point (x,y)
x = [ 0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ]
y = [ 0.0, 0.79, 1.53, 2.19, 2.71, 3.03, 3.27, 2.89, 3.06, 3.19, 3.29]
# boundary conditions of clamped spline y'_0, y'_n
conval = [ 0.8, 0.2 ]

t = LinRange(0,10,100)
p_res = Array{Float64}(undef,100)
for i in 1:100
p_res[i] = myspline( x,y,conval,t[i] )
end
plot(t,p_res,label="Interp",legend=:bottomright )
scatter!(x,y,label="Date")
``````

I guess the instability might be caused by error, but I don’t know how to fix it. Also, I wonder that why it’s no smooth.
Thanks for the help.

Lots of code here, but my guess is that you’re reading from one of your `undef` array elements before writing to it. You can force the arrays to be initialized to zero by changing `Array{Float64}(undef, N)` to `zeros(Float64, N)`. If that fixes the problem, then you can go back to `undef` (if you want) and figure out which element is being read before being written.

1 Like

Changing `Array{Float64}(undef, N)` to `zeros(Float64, N)` fixed the instability problem. Thank you for your advice.

The curve is still not smooth. I compared this code with a Matlab code that using the same method to build a function for clamped spline, the curve plotted by that code is smooth. So I’m very confused. I supposed that the code should either plot curve that don’t fit the data at all, or should work fine and plots smooth curve as I expected for spline expect. Going one step further than what @rdeits correctly intuited, you can try to see exactly which vector element gets accessed before being initialized. One way of doing so consists in building arrays of type `Vector{Any}`, initialized with a value that is not a number. For example `nothing`.

For example:

``````h = Any[nothing for _ in 1:ns]   # instead of Array{Float64}(undef,ns)
``````

This way, if you try to compute something using an element of `h` without having initialized it, you’re likely to get a `MethodError`. (At the price of a longer computing time, because the code is now type-unstable; don’t forget to revert to `Vector{Float64}` after debugging).

Below is the code “instrumented” in this way, and indeed we get a `MethodError`:

Broken, "instrumented" code
``````using Plots

function init_array(n)
# Array{Float64}(undef, n)
Any[nothing for _ in 1:n]
end

function myspline( x::Array,y::Array,conval::Array, x_i::Float64 )
np = length(x) # the number of data points
ns = np - 1

h=init_array(ns)
for i in 1:ns
h[i]=x[i+1]-x[i]
end

μ = init_array(ns)
λ = init_array(ns)
d = init_array(np)

μ[ns] = 1
λ = 1
d = 6 * ( (y-y)/h - conval ) / h
d[np] = 6 * ( conval- (y[np]-y[ns])/h[ns]) / h[ns]

for j in 1:ns-1
λ[j+1] = h[j+1] / ( h[j] + h[j+1] )
μ[j] = 1- λ[j+1]
d[j+1] = 6 * ( (y[j+2]-y[j+1])/h[j+1] - (y[j+1]-y[j])/h[j] ) / ( h[j] + h[j+1] )
end

M = init_array(np)

# solving tridiagonal matrix algorithm
α = init_array(np)
β = init_array(ns)
β = λ/2
for i in 2:ns
β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
end
y_ch = init_array(np)
y_ch = d/2
for i in 2:np
y_ch[i] = ( d[i] - μ[i-1] * y_ch[i-1] )/( 2 - μ[i-1] * β[i-1] )
end
M[np] = y_ch[np]
for i in ns
M[np-i] = y_ch[np-i] - β[np-i] * M[np-i+1]  # this is where the problem occurs
end

j_x_i = 1
while x[j_x_i+1] < x_i
j_x_i = j_x_i + 1
end
return M[j_x_i] * (x[j_x_i+1]-x_i)^3/(6*h[j_x_i]) + M[j_x_i+1] * (x_i-x[j_x_i])^3/(6*h[j_x_i]) + (y[j_x_i]-M[j_x_i]*h[j_x_i]^2/6)*(x[j_x_i+1]-x_i)/h[j_x_i] + (y[j_x_i+1]-M[j_x_i+1]*h[j_x_i]^2/6)*(x_i-x[j_x_i])/h[j_x_i]
end

# input: data point (x,y)
x = [ 0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ]
y = [ 0.0, 0.79, 1.53, 2.19, 2.71, 3.03, 3.27, 2.89, 3.06, 3.19, 3.29]
# boundary conditions of clamped spline y'_0, y'_n
conval = [ 0.8, 0.2 ]

t = LinRange(0,10,100)
p_res = init_array(100)
for i in 1:100
p_res[i] = myspline( x,y,conval,t[i] )
end
plot(t,p_res,label="Interp",legend=:bottomright )
scatter!(x,y,label="Date")
``````
``````julia> include("splines.jl")
ERROR: LoadError: MethodError: no method matching *(::Float64, ::Nothing)
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...) at operators.jl:560
*(::Union{Float16, Float32, Float64}, ::BigFloat) at mpfr.jl:392
*(::AbstractFloat, ::Bool) at bool.jl:108
...
Stacktrace:
 myspline(x::Vector{Float64}, y::Vector{Float64}, conval::Vector{Float64}, x_i::Float64)
@ Main /tmp/discourse/splines.jl:48
 top-level scope
@ /tmp/discourse/splines.jl:67
 include(fname::String)
@ Base.MainInclude ./client.jl:444
 top-level scope
@ REPL:1
in expression starting at /tmp/discourse/splines.jl:66
``````

Line 48 in your code is the body of the following loop:

``````    for i in ns
M[np-i] = y_ch[np-i] - β[np-i] * M[np-i+1]   # this is where the problem occurs
end
``````

which is a classical mistake of iterating over `ns` instead of `1:ns`.

The corrected code looks like this, and produces the following plot:

Fixed code
``````using Plots

function init_array(n)
Array{Float64}(undef, n)
# Any[nothing for _ in 1:n]
end

function myspline( x::Array,y::Array,conval::Array, x_i::Float64 )
np = length(x) # the number of data points
ns = np - 1

h=init_array(ns)
for i in 1:ns
h[i]=x[i+1]-x[i]
end

μ = init_array(ns)
λ = init_array(ns)
d = init_array(np)

μ[ns] = 1
λ = 1
d = 6 * ( (y-y)/h - conval ) / h
d[np] = 6 * ( conval- (y[np]-y[ns])/h[ns]) / h[ns]

for j in 1:ns-1
λ[j+1] = h[j+1] / ( h[j] + h[j+1] )
μ[j] = 1- λ[j+1]
d[j+1] = 6 * ( (y[j+2]-y[j+1])/h[j+1] - (y[j+1]-y[j])/h[j] ) / ( h[j] + h[j+1] )
end

M = init_array(np)

# solving tridiagonal matrix algorithm
α = init_array(np)
β = init_array(ns)
β = λ/2
for i in 2:ns
β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
end
y_ch = init_array(np)
y_ch = d/2
for i in 2:np
y_ch[i] = ( d[i] - μ[i-1] * y_ch[i-1] )/( 2 - μ[i-1] * β[i-1] )
end
M[np] = y_ch[np]
for i in 1:ns
M[np-i] = y_ch[np-i] - β[np-i] * M[np-i+1]  # this is where the problem occurs
end

j_x_i = 1
while x[j_x_i+1] < x_i
j_x_i = j_x_i + 1
end
return M[j_x_i] * (x[j_x_i+1]-x_i)^3/(6*h[j_x_i]) + M[j_x_i+1] * (x_i-x[j_x_i])^3/(6*h[j_x_i]) + (y[j_x_i]-M[j_x_i]*h[j_x_i]^2/6)*(x[j_x_i+1]-x_i)/h[j_x_i] + (y[j_x_i+1]-M[j_x_i+1]*h[j_x_i]^2/6)*(x_i-x[j_x_i])/h[j_x_i]
end

# input: data point (x,y)
x = [ 0.0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ]
y = [ 0.0, 0.79, 1.53, 2.19, 2.71, 3.03, 3.27, 2.89, 3.06, 3.19, 3.29]
# boundary conditions of clamped spline y'_0, y'_n
conval = [ 0.8, 0.2 ]

t = LinRange(0,10,100)
p_res = init_array(100)
for i in 1:100
p_res[i] = myspline( x,y,conval,t[i] )
end
plot(t,p_res,label="Interp",legend=:bottomright )
scatter!(x,y,label="Date")
`````` 7 Likes

Since iterating over an integer is identified as a classical mistake, should one consider to display a warning on the specialized iterate method ?