Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

well, I had to tinker a little :code on sagecell

from sage.plot.plot3d.plot3d import axes
aL=[2,1,3]   # coordinates list of the vector a on which b will be projeted
bL=[3,7,5] # coordinates list vector b (vector example to be projeted)
lL=5/2 # factor expansion to display red line in vector a direction

# here I adopt the convention a vector is a column matrix
A=matrix(QQ,aL).transpose() # matrix(2x1)

print( "P is the B projected vector on vector A ")
print("below all Capital letters are matrices(vectors are matrices),tiny letters are scalars ")
print("E=B-P with P=k*A " ,"E=B-k*A as E perpendicular to A => dot product(A,E)) = 0 means A.transpose()*E=0 " )
print("then A^T*(B-k*A)=0 so k*A^T*A)=A^T*B and k*(A^T*A)^-1 *A^T*A=(A^T*A)^-1 * A^T*B " )
print("so k=(A^T*A)^-1 * A^T*B")
#try: Pr
#except NameError: Pr = None
#if Pr == None :
#    Pr = A
@interact
def _(xb=slider([1..10],default=3),yb=slider([1..10],default=4),zb=slider([1..10],default=7)):
    #show(xb,yb,zb)
    bL=[xb,yb,zb]
    B=matrix(QQ,bL).transpose() # matrix(3x1)
    # Projection Matrix , called Pr 
    #The general formula for the orthogonal projection onto the column space of a matrix A is
    # P = A (A^T A)^(−1) A^T
    Pr_num=A * (A.transpose()) # matrix 3x3
    Pr_den=(( (A.transpose()) * A ).det()) # scalar
    #Pr_den=1 # scalar
    #Pr=Pr_num/Pr_den # matrix(3x3) 
    # other compact form 
    Pr=A*((A.transpose()*A).inverse() * A.transpose())
    g1=line([[-lL*e for e in A.list()],[lL*e for e in A.list()]], legend_label=' a projection line', legend_color='red',color='red')
    g2=plot(vector(B),color="green",width=20,legend_label=' vector b', legend_color='green')
    g3=plot(vector(A),color="black",width=3,legend_label=' vector A', legend_color='black')
    aT = text3d("a",vector(A)*1.1,color='black',fontsize=20)
    bT = text3d("b",vector(B),color='green',fontsize=20)
    # projected vector P
    P=Pr*B # matrix 3x1

    k0=((A.transpose()*A).inverse() * A.transpose()*B).det()
    k1=sqrt(((P.transpose()*P).det()/(A.transpose()*A).det()))

    g4=plot(vector(P),color='blue',width=20,legend_label=' projected vector P', legend_color='blue',fontsize=20)
    pT = text3d("p",vector(P)*1.1,color='blue',fontsize=40)
    E=B-P # matrix 3x1
    g5=plot(arrow(P.coefficients(),B.coefficients(),color='pink',width=5,legend_label=' vector e', legend_color='pink'))
    eT = text3d("e",vector(E)/2+vector(P)*1.05,color='pink',fontsize=20)
    # unit vectors
    gx=plot(vector([1,0,0]),color='black',width=5,legend_label='x', legend_color='black')
    xT = text3d("x",vector([1,0.1,0.1])*1.1,color='black',fontsize=20)
    gy=plot(vector([0,1,0]),color='black',width=5,legend_label='x', legend_color='black')
    yT = text3d("y",vector([0.1,1,0.1])*1.1,color='black',fontsize=20)
    gz=plot(vector([0,0,1]),color='black',width=5,legend_label='x', legend_color='black')
    zT = text3d("z",vector([0.1,0.1,1])*1.1,color='black',fontsize=20)

    show(g1+g2+g3+g4+g5+aT+bT+pT+eT+gx+xT+gy+yT+gz+zT+axes(7, color='black') \
    ,figsize=(6, 6),xmin=-vector(B)[0],xmax=vector(B)[0]+3,ymax=vector(A)[1]+1,ymin=-1,zmax=vector(A)[2]+1,zmin=-1)

    show(" A : ",A, " B : ",B, " Pr_num : ",Pr_num, " Pr_den : ",Pr_den, " Pr : ",Pr)
    show(" A : ",A, " B : ",B, " projected vector   P : ",P  )
    print ("verify that k=(A^T*A)^(-1) * A^T*B = ||P|| / ||A||")
    print (" so P=A*k= A*(A^T*A)^-1 * A^T*B ")
    show ("P also equals to : A*(A^T*A)^(-1) * A^T*B : ",A*((A.transpose()*A).inverse() * A.transpose())*B )
    show( "P= Projection Matrix * B= Pr*B so PR=A*((A.transpose()*A).inverse() * A.transpose() :",A*((A.transpose()*A).inverse() * A.transpose()))
    show ( " sqrt(((P.transpose()*P).det()/(A.transpose()*A).det())) = ||P|| / ||A|| : ",k1)

well, I had to tinker a little :code on sagecell

from sage.plot.plot3d.plot3d import axes
aL=[2,1,3]   # coordinates list of the vector a on which b will be projeted
bL=[3,7,5] # coordinates list vector b (vector example to be projeted)
lL=5/2 # factor expansion to display red line in vector a direction

# here I adopt the convention a vector is a column matrix
A=matrix(QQ,aL).transpose() # matrix(2x1)

print( "P is the B projected vector on vector A ")
print("below all Capital letters are matrices(vectors are matrices),tiny letters are scalars ")
print("E=B-P with P=k*A " ,"E=B-k*A as E perpendicular to A => dot product(A,E)) = 0 means A.transpose()*E=0 " )
print("then A^T*(B-k*A)=0 so k*A^T*A)=A^T*B and k*(A^T*A)^-1 *A^T*A=(A^T*A)^-1 * A^T*B " )
print("so k=(A^T*A)^-1 * A^T*B")
#try: Pr
#except NameError: Pr = None
#if Pr == None :
#    Pr = A
@interact
def _(xb=slider([1..10],default=3),yb=slider([1..10],default=4),zb=slider([1..10],default=7)):
    #show(xb,yb,zb)
    bL=[xb,yb,zb]
    B=matrix(QQ,bL).transpose() # matrix(3x1)
    # Projection Matrix , called Pr 
    #The general formula for the orthogonal projection onto the column space of a matrix A is
    # P = A (A^T A)^(−1) A^T
    Pr_num=A * (A.transpose()) # matrix 3x3
    Pr_den=(( (A.transpose()) * A ).det()) # scalar
    #Pr_den=1 # scalar
    #Pr=Pr_num/Pr_den # matrix(3x3) 
    # other compact form 
    Pr=A*((A.transpose()*A).inverse() * A.transpose())
    g1=line([[-lL*e for e in A.list()],[lL*e for e in A.list()]], legend_label=' a projection line', legend_color='red',color='red')
    g2=plot(vector(B),color="green",width=20,legend_label=' vector b', legend_color='green')
    g3=plot(vector(A),color="black",width=3,legend_label=' vector A', legend_color='black')
    aT = text3d("a",vector(A)*1.1,color='black',fontsize=20)
    bT = text3d("b",vector(B),color='green',fontsize=20)
    # projected vector P
    P=Pr*B # matrix 3x1

    k0=((A.transpose()*A).inverse() * A.transpose()*B).det()
    k1=sqrt(((P.transpose()*P).det()/(A.transpose()*A).det()))

    g4=plot(vector(P),color='blue',width=20,legend_label=' projected vector P', legend_color='blue',fontsize=20)
    pT = text3d("p",vector(P)*1.1,color='blue',fontsize=40)
    E=B-P # matrix 3x1
    g5=plot(arrow(P.coefficients(),B.coefficients(),color='pink',width=5,legend_label=' vector e', legend_color='pink'))
    eT = text3d("e",vector(E)/2+vector(P)*1.05,color='pink',fontsize=20)
    # unit vectors
    gx=plot(vector([1,0,0]),color='black',width=5,legend_label='x', legend_color='black')
    xT = text3d("x",vector([1,0.1,0.1])*1.1,color='black',fontsize=20)
    gy=plot(vector([0,1,0]),color='black',width=5,legend_label='x', legend_color='black')
    yT = text3d("y",vector([0.1,1,0.1])*1.1,color='black',fontsize=20)
    gz=plot(vector([0,0,1]),color='black',width=5,legend_label='x', legend_color='black')
    zT = text3d("z",vector([0.1,0.1,1])*1.1,color='black',fontsize=20)

    show(g1+g2+g3+g4+g5+aT+bT+pT+eT+gx+xT+gy+yT+gz+zT+axes(7, color='black') \
    ,figsize=(6, 6),xmin=-vector(B)[0],xmax=vector(B)[0]+3,ymax=vector(A)[1]+1,ymin=-1,zmax=vector(A)[2]+1,zmin=-1)

    show(" A : ",A, " B : ",B, " Pr_num : ",Pr_num, " Pr_den : ",Pr_den, " Pr : ",Pr)
    show(" A : ",A, " B : ",B, " projected vector   P : ",P  )
    print ("verify that k=(A^T*A)^(-1) * A^T*B = ||P|| / ||A||")
    print (" so P=A*k= A*(A^T*A)^-1 * A^T*B ")
    show ("P also equals to : A*(A^T*A)^(-1) * A^T*B : ",A*((A.transpose()*A).inverse() * A.transpose())*B )
    show( "P= Projection Matrix * B= Pr*B so PR=A*((A.transpose()*A).inverse() * A.transpose() :",A*((A.transpose()*A).inverse() * A.transpose()))
    show ( " sqrt(((P.transpose()*P).det()/(A.transpose()*A).det())) = ||P|| / ||A|| : ",k1)

another simple example: Point -> Vector to a point (post)

E = EuclideanSpace(3)
O = E.point((0,0,0), name='Origin')
u= E.tangent_space(O)((1/2,1,0), name='u')  # u vector at O
v = E.tangent_space(O)((1,2,1), name='v')  # v vector at O
#vu=E.tangent_space(v.)((1,2,1), name='v')  # v vector at O
P=O.plot(size=20, color='red', opacity=0.5) + v.plot(color='green') + u.plot(color='red')
P+=plot(vector([1,0,0]),color='black',width=5,legend_label='x', legend_color='black')
P+= text3d("x",vector([1,0.1,0.1])*1.1,color='black',fontsize=20)
P+=plot(vector([0,1,0]),color='black',width=5,legend_label='x', legend_color='black')
P+= text3d("y",vector([0.1,1,0.1])*1.1,color='black',fontsize=20)
P+=plot(vector([0,0,1]),color='black',width=5,legend_label='x', legend_color='black')
P+= text3d("z",vector([0.1,0.1,1])*1.1,color='black',fontsize=20)
show(P, viewer="jmol")
#show(P, viewer="threejs")