Revision history [back]
 | 1 | initial version |
In sage, there is a method gram_schnidt associated to a matrix.
After the matrix is declared, the method can be immediately called.
The following defines a matrix A , computes the Gram-Schnidt matrices G.M as in the examples from
A.gram_schmidt?
and does a check for the relation $A=MG$, it is a possible starting point:
var( 'x' )
def ip(f, g):
return integrate( f * g, x, -1, 1 )
f_list = [1, x, x^2, x^3, x^4]
A = matrix( QQ, len(f_list),
[ ip( f, g ) for f in f_list for g in f_list ] )
print "A is\n%s\n" % A
G, M = A.gram_schmidt()
print "G is\n%s\n" % G
print "M is\n%s\n" % M
print "Is MG = A? %s" % bool( M*G == A )
Results:
A is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[2/3 0 2/5 0 2/7]
[ 0 2/5 0 2/7 0]
[2/5 0 2/7 0 2/9]
G is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[ -424/5439 0 1376/9065 0 248/1813]
[ 0 -12/595 0 4/119 0]
[ 64/112035 0 -128/22407 0 64/9603]
M is
[ 1 0 0 0 0]
[ 0 1 0 0 0]
[675/1813 0 1 0 0]
[ 0 75/119 0 1 0]
[425/1813 0 950/1067 0 1]
Is MG = A? True
| | 2 | No.2 Revision |
In sage, there is a method gram_schnidt associated to a matrix.
After the matrix is declared, the method can be immediately called.
The following defines a matrix A , computes the Gram-Schnidt matrices G.M as in the examples from
A.gram_schmidt?
and does a check for the relation $A=MG$, it is a possible starting point:
var( 'x' )
def ip(f, g):
return integrate( f * g, x, -1, 1 )
f_list = [1, x, x^2, x^3, x^4]
A = matrix( QQ, len(f_list),
[ ip( f, g ) for f in f_list for g in f_list ] )
print "A is\n%s\n" % A
G, M = A.gram_schmidt()
print "G is\n%s\n" % G
print "M is\n%s\n" % M
print "Is MG = A? %s" % bool( M*G == A )
Results:
A is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[2/3 0 2/5 0 2/7]
[ 0 2/5 0 2/7 0]
[2/5 0 2/7 0 2/9]
G is
[ 2 0 2/3 0 2/5]
[ 0 2/3 0 2/5 0]
[ -424/5439 0 1376/9065 0 248/1813]
[ 0 -12/595 0 4/119 0]
[ 64/112035 0 -128/22407 0 64/9603]
M is
[ 1 0 0 0 0]
[ 0 1 0 0 0]
[675/1813 0 1 0 0]
[ 0 75/119 0 1 0]
[425/1813 0 950/1067 0 1]
Is MG = A? True
Later edit:
I try now hard to reproduce in sagemath the following result from maxima, obtained by copy+pasting into the maxima interpreter the code, thanks go to danielvolinski for this sample,
load(eigen);
ip(f,g) := integrate(f*g, u, 0, 1);
y : gramschmidt([1, u, u^2, u^3, u^4], ip);
which delivers
2 2
2 u - 1 6 u - 6 u + 1 (2 u - 1) (10 u - 10 u + 1)
(%o14) [1, -------, --------------, ----------------------------,
2 6 20
4 3 2
70 u - 140 u + 90 u - 20 u + 1
---------------------------------]
70
The situation is different, because there is no a priori knowledge of an orthonormal basis of the vector space generated by $1,u,u^2,u^3,u^4$ with the given scalar product $\langle f,g\rangle =\int_0^1f\;g$ .
To get a similar result i prefer to program from the scratch, there is no implemented facility...:
var( 'u' )
def ip(f, g):
return integrate( f * g, u, 0, 1 )
n = 5
B = vector( [u^k for k in range(n)] )
T = identity_matrix( QQ, n ) # in T we will have finally the transfromation matrix
C = vector( [u^k for k in range(n)] ) # in C we will have finally the Gram-Schmidt basis
for k in range(n-1):
S = identity_matrix( QQ, n )
for j in range(k+1, n):
S[j, k] = -ip( C[j], C[k] ) / ip( C[k], C[k] )
T = S*T
C = S*C
print T
print C
This gives the transformation matrix, and the new basis
[ 1 0 0 0 0]
[ -1/2 1 0 0 0]
[ 1/6 -1 1 0 0]
[-1/20 3/5 -3/2 1 0]
[ 1/70 -2/7 9/7 -2 1]
(1, u - 1/2, u^2 - u + 1/6, u^3 - 3/2*u^2 + 3/5*u - 1/20, u^4 - 2*u^3 + 9/7*u^2 - 2/7*u + 1/70)
in concordance with the maxima result.
