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

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