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](https://global.discourse-cdn.com/julialang/original/3X/6/a/6ab2470166ba5efce8cefa83c83816fe81bd02f3.png)