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

New to Julia
,
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