Sympy funnies & Cubic hermite

hermite.py

from sympy import *

# https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf

a0, a1, a2, a3, t = symbols('a0 a1 a2 a3 t')

polynomial = a0 + a1*t + a2*t**2 + a3*t**3
polynomial_d = diff(polynomial, t)

p0, v0, p1, v1 = symbols('p0 v0 p1 v1')

spline = [
    Eq(p0, polynomial.subs(t, 0)),
    Eq(v0, polynomial_d.subs(t, 0)),
    Eq(p1, polynomial.subs(t, 1)),
    Eq(v1, polynomial_d.subs(t, 1))]

print '\n'.join(map(str, spline))
print 'solution:'
solution = solve(spline, (a0, a1, a2, a3))
for lhs, rhs in solution.items():
    print Eq(lhs, rhs)
cf = polynomial.subs(solution)
print 'c(t) =', cf

for label, var in [('h0', p0), ('h1', v0), ('h2', v1), ('h3', p1)]:
    cf = apart(cf, var)
    coeff = cf.coeff(var)
    cf = simplify(cf - coeff * var)
    print label, '=', coeff
assert cf == 0

# These would also "work", and be quicker. r
# But I wanted to extract coefficients safe rather than sorry.
#print 'h0 =', apart(cf, p0).coeff(p0)
#print 'h1 =', apart(cf, v0).coeff(v0)
#print 'h2 =', apart(cf, v1).coeff(v1)
#print 'h3 =', apart(cf, p1).coeff(p1)
#print 'h0 =', cf.subs({p0:1, v0:0, v1:0, p1:0})
#print 'h1 =', cf.subs({p0:0, v0:1, v1:0, p1:0})
#print 'h2 =', cf.subs({p0:0, v0:0, v1:1, p1:0})
#print 'h3 =', cf.subs({p0:0, v0:0, v1:0, p1:1})

# stdout
# p0 == a0
# v0 == a1
# p1 == a0 + a1 + a2 + a3
# v1 == a1 + 2*a2 + 3*a3
# solution:
# a3 == 2*p0 - 2*p1 + v0 + v1
# a0 == p0
# a2 == -3*p0 + 3*p1 - 2*v0 - v1
# a1 == v0
# c(t) = p0 + t**3*(2*p0 - 2*p1 + v0 + v1) + t**2*(-3*p0 + 3*p1 - 2*v0 - v1) + t*v0
# h0 = 2*t**3 - 3*t**2 + 1
# h1 = t**3 - 2*t**2 + t
# h2 = -2*t**3 + 3*t**2
# h3 = t**3 - t**2

simple.py

import sympy
from sympy import *

t = symbols('t')
function = (t, t**2)

print "C(t) =", function
print "sample t=0.0", [axis.subs(t, 0.0) for axis in function]
print "sample t=0.5", [axis.subs(t, 0.5) for axis in function]
print "sample t=1.0", [axis.subs(t, 1.0) for axis in function]

def curvature_2d((xf, yf)):
    x_d1 = diff(xf, t)
    x_d2 = diff(x_d1, t)
    y_d1 = diff(yf, t)
    y_d2 = diff(y_d1, t)
    return (x_d1*y_d2 - y_d1*x_d2) / sqrt(x_d1**2 + y_d1**2)**3.0

kf = curvature_2d(function)
print "k(t) =", kf
print "sample k(0.0) =", kf.subs(t, 0.0)
print "sample k(0.5) =", kf.subs(t, 0.5)
print "sample k(1.0) =", kf.subs(t, 1.0)

mec = Integral(kf, (t, -1.0, +1.0)).evalf()
print "mec =", mec

mvcc = Integral(diff(kf, t), (t, -1.0, +1.0)).evalf()
print "mvcc =", mvcc

x_d = diff(function[0], t)
y_d = diff(function[1], t)
length = Integral(sqrt(x_d**2 + y_d**2), (t, -1.0, +1.0)).evalf()
print "length =", length

print ccode(kf)
print ccode(diff(kf, t))
print ccode(sqrt(x_d**2 + y_d**2))