Start notebook with --profile=sympy flag.

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

In :
%qtconsole
In :
(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 :
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 :
X =[0,xcp,1]
Y =[0,ycp,1]

Set up each piece of interpolant

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

Left-side Interpolation conditions

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

Match end-point 1st derivatives from givens

In :
cond_end=[ diff(f,x).subs(k,0).subs(x(0),X).subs(x,X) - y1p,
diff(f,x).subs(k,1).subs(x(1),X).subs(x,X) - ynp,
]
cond_end
Out:
[-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 :
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-c).subs(x,X)]

Inner second derivatives must match

In :
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-c,x,2).subs(x,X)]

Inner first derivatives must match

In :
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-c,x).subs(x,X)]
In :
len(cond_i)+len(cond_cont)+len(cond_end)+len(cond2_cont)+len(cond1_cont)
Out:
8
In :
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 :
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 :
sum([abs(diff(i,x,2)) for i in c]) # curvature metric
Out:
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 :
print c
print c
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 :