Skip to content

Instantly share code, notes, and snippets.

@jtravs
Created April 4, 2014 18:52
Show Gist options
  • Save jtravs/9980928 to your computer and use it in GitHub Desktop.
Save jtravs/9980928 to your computer and use it in GitHub Desktop.
Quintic spline interpolation
type quint{T}
dx::T
x0::T
Y::Array{T,1}
B::Array{T,1}
C::Array{T,1}
D::Array{T,1}
E::Array{T,1}
F::Array{T,1}
function quint(X::Array{T,1})
Y = similar(X)
B = similar(X)
C = similar(X)
D = similar(X)
E = similar(X)
F = similar(X)
new(X[2]-X[1], X[1], Y, B, C, D, E, F)
end
end
quint{T}(X::Array{T,1}) = quint{T}(X)
function qeval{T}(q::quint{T}, x::T)
P = (x - q.x0)/q.dx
I = int(P)
P -= I
I += 1
@inbounds return ((((q.F[I]*P+q.E[I])*P+q.D[I])*P+q.C[I])*P+q.B[I])*P+q.Y[I]
end
function qeval{T}(q::quint{T}, X::Array{T,1}, Y::Array{T,1})
N = length(X)
for i=1:N
@inbounds Y[i] = qeval(q, X[i])
end
end
function qset{T}(q::quint{T}, Y::Array{T,1})
N = length(Y)
q.Y[:] = Y
M = N - 3
P = 0.0
Q = 0.0
R = 0.0
S = 0.0
TT = 0.0
q.D[M+1] = 0.0
q.D[M+2] = 0.0
for I=1:M
@inbounds begin
U = P*R
q.B[I] = 1.0/(66.0-U*R-Q)
R = 26.0 - U
q.C[I] = R
q.D[I] = q.Y[I+3] - 3.0*(q.Y[I+2]-q.Y[I+1]) - q.Y[I] - U*S - Q*TT
Q = P
P = q.B[I]
TT = S
S = q.D[I]
end
end
I = N - 2
for M=4:N
I = I - 1
@inbounds q.D[I] = (q.D[I]-q.D[I+2]-q.D[I+1]*q.C[I])*q.B[I]
end
M = N - 1
Q = 0.0
R = q.D[1]
TT = R
V = R
for I=2:M
@inbounds begin
P = Q
Q = R
R = q.D[I]
S = TT
TT = P - Q - Q + R
q.F[I] = TT
U = 5.0*(-P+Q)
q.E[I] = U
q.D[I] = 10.0*(P+Q)
q.B[I] = 0.5*(q.Y[I+1]-q.Y[I-1]-S-TT) - q.D[I]
q.C[I] = 0.5*(q.Y[I+1]+q.Y[I-1]+S-TT) - q.Y[I] - U
end
end
@inbounds begin
q.F[1] = V
q.E[1] = 0.0
q.D[1] = 0.0
q.C[1] = q.C[2] - 10.0*V
q.B[1] = q.Y[2] - q.Y[1] - q.C[1] - V
q.E[N] = 0.0
q.D[N] = 0.0
q.C[N] = q.C[N-1] + 10.0*TT
q.B[N] = q.Y[N] - q.Y[N-1] + q.C[N] - TT
end
end
include("quint.jl")
x = linspace(-10.0,10.0,200000)
x2 = x[1:end-1] .+ (x[2]-x[1])/2.0
y = sin(x)
y2 = similar(x2)
q = quint(x)
qset(q, y)
@time qset(q, y)
@time qset(q, y)
@time qset(q, y)
qeval(q, x2, y2)
@time qeval(q, x2, y2)
@time qeval(q, x2, y2)
@time qeval(q, x2, y2)
Copy link

ghost commented Apr 4, 2014

Copy link

ghost commented Apr 4, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment