# My code for cubic spline interpolation plots different image everytime and not smooth
**URL:** https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564
**Category:** New to Julia
**Tags:** interpolations, splines
**Created:** [October 28, 2021, 3:49pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564 "2021-10-28T15:49:25Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@pessoalzy](https://discourse.julialang.org/u/pessoalzy)
#### Post date: [October 28, 2021, 3:49pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564/1 "2021-10-28T15:49:25Z")
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.
```julia
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.
---
### Author:  [@rdeits](https://discourse.julialang.org/u/rdeits)
#### Post date: [October 28, 2021, 3:59pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564/2 "2021-10-28T15:59:42Z")
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.
---
### Author:  [@pessoalzy](https://discourse.julialang.org/u/pessoalzy)
#### Post date: [October 28, 2021, 4:24pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564/3 "2021-10-28T16:24:58Z")
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.
---
### Author:  [@ffevotte](https://discourse.julialang.org/u/ffevotte)
#### Post date: [October 28, 2021, 4:49pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564/4 "2021-10-28T16:49:36Z")
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:
```julia
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**
>
> ```julia
> 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
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:
```julia
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**
>
> ```julia
> 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")
>
> ```
---
### Author:  [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)
#### Post date: [October 28, 2021, 8:12pm UTC](https://discourse.julialang.org/t/my-code-for-cubic-spline-interpolation-plots-different-image-everytime-and-not-smooth/70564/5 "2021-10-28T20:12:38Z")
Since iterating over an integer is identified as a classical mistake, should one consider to display a warning on the specialized iterate method ?