# Sage session illustrating similarity.

# Standard Basis for P_2.
t = var('t')
S = vector([1,t,t^2])

def SCoordVec(p):
    """Standard-basis coordinate vector of a quadratic polynomial."""
    v=vector([p.subs(t==0),diff(p,t).subs(t==0),(1/2)*diff(p,t,2).subs(t==0)])
    return(v)

# LaGrange Basis for P_2 based on 0,1,2.
# For each i from 0 to 2 G[i] is the
# quadratic satisfying G[i](i) = 1 and
# G[i](j) = 0 for j not equal to i.

G = vector([(t-1)*(t-2)/((0-1)*(0-2)),
            (t-0)*(t-2)/((1-0)*(1-2)),
            (t-0)*(t-1)/((2-0)*(2-1))])

# Here are the basis polynomials:

print(G)
 (1/2*(t - 2)*(t - 1), -(t - 2)*t, 1/2*(t - 1)*t)

# ... and here are each polynomial's values at 0,1 and 2:

for i in range(3):
    print([(G[i]).subs(t==j) for j in range(3)]),
 [1, 0, 0] [0, 1, 0] [0, 0, 1]

# The G-coordinate vector of a polynomial is just
# the vector of its values at 0,1,2:

def GCoordVec(p):
    """LaGrange-basis coordinate vector of a quadratic polynomial."""
    v = vector([p.subs(t==i) for i in range(3)])
    return(v)

# Here is a quadratic polynomial q.

q = t^2 - 3*t + 5

# Here are the standard and LaGrange coordinate vectors of q.

[SCoordVec(q), GCoordVec(q)]
 [(5, -3, 1), (5, 3, 3)]

# The LaGrange coordinate vector is the vector of coefficients
# expressing q as a linear combination of G.

l = GCoordVec(q); l
l.dot_product(G).expand()
 t^2 - 3*t + 5

# We define a linear transformation L: P_2 --> P_2
# by L = D + I (D=derivative, I=identity).
# We will find the matrices that represent this 
# linear transformation in the bases S and G, respectively.

def L(p):
    """Linear transformation D+I : P_2 --> P_2"""
    return(diff(p,t)+p)

# The matrix representing L in the basis S.
# (We have to transpose the result because SCoordVec produces a row vector.)

A = matrix(QQ, [SCoordVec(L(S[i])) for i in range(3)]).transpose(); A
 [1 1 0]
 [0 1 2]
 [0 0 1]

# The matrix representing L in the basis G.
# (We must transpose the result because GCoordVec produces a row vector.)

B = matrix(QQ, [GCoordVec(L(G[i])) for i in range(3)]).transpose(); B
 [-1/2    2 -1/2]
 [-1/2    1  1/2]
 [ 1/2   -2  5/2]

# Now we show that B = P^(-1) A P, matrix P being 
# the S <-- G transition matrix.

P = matrix(QQ,[SCoordVec(G[i]) for i in range(3)]).transpose(); P
 [   1    0    0]
 [-3/2    2 -1/2]
 [ 1/2   -1  1/2]

g,s = GCoordVec(q),SCoordVec(q); print(g,s)
 ((5, 3, 3), (5, -3, 1))

s == P*g
 True

B == P.inverse()*A*P
 True