Start notebook with --profile=sympy flag.

Following the 1968 paper by Spath. Consider special case with equally spaced abscissae

In [1]:
%qtconsole
In [2]:
(A, B, C, D)=symbols('A:D') # coefficients
p = symbols('p') # tension parameter
y1p = symbols('y1p') # derivative a left endpoint
ynp = symbols('ynp') # derivative a right endpoint
xcp = symbols('xcp') # x control point
ycp = symbols('ycp') # y control point

This is the $k^{th}$ piece of the $n$-piece interpolant,

$$f_k(x) = A_k+B_k (x-x_k) + C_k \exp(p_k (x-x_k)) +D_k \exp(-p_k (x-x_k))$$

given the derivative of the target function at $x(0)$ and at the other end $x(n-1)$.

In [3]:
f = A(k)+B(k)*(x-x(k)) + C(k)*exp(p(k)*(x-x(k))) +D(k)*exp(-p(k)*(x-x(k)))

Sample data to interpolate

In [4]:
X =[0,xcp,1]
Y =[0,ycp,1]

Set up each piece of interpolant

In [5]:
c=[f.subs(x(k),X[i]).subs(k,i) for i in range(2)]

Left-side Interpolation conditions

In [6]:
cond_i=[(Y[i]-f.subs(x,x(k)).subs(k,i)) for i in range(2)] # conditions for interpolation
cond_i+= [ Y[2]- c[1].subs(x,X[2])]

Match end-point 1st derivatives from givens

In [7]:
cond_end=[ diff(f,x).subs(k,0).subs(x(0),X[0]).subs(x,X[0]) - y1p,
           diff(f,x).subs(k,1).subs(x(1),X[1]).subs(x,X[2]) - ynp,
 ]
cond_end
Out[7]:
[-y1p + B(0) + C(0)*p(0) - D(0)*p(0),
 -ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))]

Inner continuity conditions

In [8]:
cond_cont=[] # continuity conditions
for i,j,v in zip(c[:-2],c[2:],range(1,3)):
    cond_cont.append( (i-j).subs(x,v) )
cond_cont

cond_cont=[(c[0]-c[1]).subs(x,X[1])]

Inner second derivatives must match

In [9]:
d2=[i.diff(x,2) for i in c]
cond2_cont=[] # 2nd derivative continuity conditions
for i,j,v in zip(d2[:-2],d2[2:],range(3)):
    cond2_cont.append( (i-j).subs(x,v) )
    
cond2_cont= [diff(c[0]-c[1],x,2).subs(x,X[1])]

Inner first derivatives must match

In [10]:
d=[i.diff(x) for i in c]
cond1_cont=[] # 2nd derivative continuity conditions
for i,j,v in zip(d[:-2],d[2:],range(3)):
    cond1_cont.append( (i-j).subs(x,v) )
    
cond1_cont=[diff(c[0]-c[1],x).subs(x,X[1])]
In [11]:
len(cond_i)+len(cond_cont)+len(cond_end)+len(cond2_cont)+len(cond1_cont)
Out[11]:
8
In [12]:
for i in (cond_i+cond_cont+cond_end+cond2_cont):
    print i.subs(p(k),0)

-A(0) - C(0) - D(0)
ycp - A(1) - C(1) - D(1)
-(-xcp + 1)*B(1) - A(1) - C(1)*exp((-xcp + 1)*p(1)) - D(1)*exp(-(-xcp + 1)*p(1)) + 1
xcp*B(0) + A(0) - A(1) + C(0)*exp(xcp*p(0)) - C(1) + D(0)*exp(-xcp*p(0)) - D(1)
-y1p + B(0) + C(0)*p(0) - D(0)*p(0)
-ynp + B(1) + C(1)*p(1)*exp((-xcp + 1)*p(1)) - D(1)*p(1)*exp(-(-xcp + 1)*p(1))
C(0)*p(0)**2*exp(xcp*p(0)) - C(1)*p(1)**2 + D(0)*p(0)**2*exp(-xcp*p(0)) - D(1)*p(1)**2
In [13]:
linsys=(cond_i+cond_cont+cond_end+cond2_cont+cond1_cont)
M=Matrix([[ l.coeff(i) for i in flatten([[A(j),B(j),C(j),D(j)] for j in range(2)]) ] for l in linsys])
In [14]:
sum([abs(diff(i,x,2)) for i in c]) # curvature metric
Out[14]:
Abs(C(0)*p(0)**2*exp(x*p(0)) + D(0)*p(0)**2*exp(-x*p(0))) + Abs(C(1)*p(1)**2*exp((x - xcp)*p(1)) + D(1)*p(1)**2*exp(-(x - xcp)*p(1)))
In [15]:
print c[0]
print c[1]

x*B(0) + A(0) + C(0)*exp(x*p(0)) + D(0)*exp(-x*p(0))
(x - xcp)*B(1) + A(1) + C(1)*exp((x - xcp)*p(1)) + D(1)*exp(-(x - xcp)*p(1))
In [15]: