In [21]:
# Continuous definition of functions (using indicator functions)

def indicator_above(x,vmin):
    x1 = Abs(vmin - x)
    return (x + x1 - vmin) / (2*x1)

def indicator_below(x,vmax):
    x1 = Abs(vmax - x)
    return 1 - (x + x1 - vmax) / (2*x1)

def indicator_range(x,vmin,vmax):
    x0 = vmax - x
    x1 = Abs(vmin - x)
    x2 = Abs(x0)
    return (x0 + x2) * (x + x1 - vmin) / (4*x1*x2)

def BS(i,k,t):
    if (k==1):
        return indicator_range(t,i,i+1)
    else:
        return BS(i,k-1,t) * (t - i)/((i+k-1) - i) + BS(i+1,k-1,t) * (i+k - t) / (i+k - (i+1));

def CBS(i,k,t):
    ind = indicator_below(t,i+k)
    s = sum( [BS(j,k,t) for j in range(i,i+k)] )
    return ind * s + (1-ind)

#t = symbols('t')
#plot(B(0,4,t), B(1,4,t), B(2,4,t), B(3,4,t), (t,0,7))
#plot(CB(0,4,t), CB(1,4,t), CB(2,4,t), CB(3,4,t), (t,0,7))
In [20]:
# Straightforward definition

def BS(i,k,t):
    if (k==1):
        if( i <= t and t < i+1):
            return 1.0
        else:
            return 0.0
    else:
        return BS(i,k-1,t) * (t - i)/((i+k-1) - i) + BS(i+1,k-1,t) * (i+k - t) / (i+k - (i+1));

def CBS(i,k,t):
    if(t > i+k):
        return 1.0
    else:
        return sum( [BS(j,k,t) for j in range(i,i+k)] )

xs = linspace(0,7,100)
for i in range(0,4):
    pyplot.plot(xs, [BS(i,4,x) for x in xs])
    
pyplot.figure(2)
for i in range(0,4):
    pyplot.plot(xs, [CBS(i,4,x) for x in xs])
In [92]:
from sympy import *
u = symbols('u')
#cse for common subexpr

# Matrix from General matrix representations for B-splines, Quin
M = Matrix([
    [1, -3,  3, -1],
    [4,  0, -6,  3],
    [1,  3,  3, -3],
    [0,  0,  0,  1]
    ]) / 6.0

Sum = Matrix((
    [-1, 0, 0],
    [-1,-1, 0],
    [-1,-1,-1]
    ))

C = Sum * M[0:3,:]

def U(u):
    return Matrix(([[1],[u],[u**2], [u**3]]))

def B(u):
    return M * U(u)

def B2(u):
    return Matrix([
     [BS(-2,4,u)],
     [BS(-1,4,u)],
     [BS(0,4,u)] 
    ])

def CB(u):
    return Matrix(([[1],[1],[1]])) + C * U(u)

def CB2(u):
    return Matrix([
     [CBS(-2,4,u)],
     [CBS(-1,4,u)],
     [CBS(0,4,u)] 
    ])

def dCB(u):
    return CB(u).diff(u,1)

def d2CB(u):
    return CB(u).diff(u,2)

plot(B(u)[0], B(u)[1], B(u)[2], (u,0,1))
#plot(B2(u)[0], B2(u)[1], B2(u)[2], (u,0,1))
plot(CB(u)[0], CB(u)[1], CB(u)[2], (u,0,1))
#plot(CB2(u)[0], CB2(u)[1], CB2(u)[2], (u,0,1))
for f in [CB, dCB, d2CB]:
    pprint( f(u).subs(u,0.0).transpose() )
[0.833333333333333  0.166666666666667  1.11022302462516e-16]
[0.5  0.5  0]
[-1.0  1.0  0]