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

New to Julia
,
1

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.

4

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] = 1        
    d[1] = 6 * ( (y[2]-y[1])/h[1] - conval[1] ) / h[1]
    d[np] = 6 * ( conval[2]- (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)
    β[1] = λ[1]/2
    for i in 2:ns
       β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
    end
    y_ch = Array{Float64}(undef,np)
    y_ch[1] = d[1]/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.

2

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
3

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.
new

4

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] = 1
    d[1] = 6 * ( (y[2]-y[1])/h[1] - conval[1] ) / h[1]
    d[np] = 6 * ( conval[2]- (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)
    β[1] = λ[1]/2
    for i in 2:ns
       β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
    end
    y_ch = init_array(np)
    y_ch[1] = d[1]/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:
 [1] myspline(x::Vector{Float64}, y::Vector{Float64}, conval::Vector{Float64}, x_i::Float64)
   @ Main /tmp/discourse/splines.jl:48
 [2] top-level scope
   @ /tmp/discourse/splines.jl:67
 [3] include(fname::String)
   @ Base.MainInclude ./client.jl:444
 [4] top-level scope
   @ REPL[134]: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] = 1
    d[1] = 6 * ( (y[2]-y[1])/h[1] - conval[1] ) / h[1]
    d[np] = 6 * ( conval[2]- (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)
    β[1] = λ[1]/2
    for i in 2:ns
       β[i] = λ[i]/( 2 - μ[i-1] * β[i-1] )
    end
    y_ch = init_array(np)
    y_ch[1] = d[1]/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")

splines

7 Likes
5

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