gram01.mw

Gram-Schmidt Examples.

Ingnore this stuff:

> restart; with(LinearAlgebra):

> ip := proc(u, v) DotProduct(u,v, conjugate=false): end;

ip := proc (u, v) LinearAlgebra:-LinearAlgebra(u, v, conjugate = false) end proc

> nrm := proc(v) sqrt(simplify(ip(v,v))); end;

nrm := proc (v) sqrt(simplify(ip(v, v))) end proc

> nrm2 := proc(v) simplify(ip(v,v)); end;

nrm2 := proc (v) simplify(ip(v, v)) end proc

The inner product will be denote by "ip(u,v)", where u and v are vectors.  For example:

> X := Vector(4, symbol=x);  Y := Vector(4, symbol=y);

X := Vector[column]([[x[1]], [x[2]], [x[3]], [x[4]]])

Y := Vector[column]([[y[1]], [y[2]], [y[3]], [y[4]]])

> ip(X,Y);

x[1]*y[1]+x[2]*y[2]+x[3]*y[3]+x[4]*y[4]

We'll use "nrm2(v)" for the norm squared of v.

> nrm2(X);

x[1]^2+x[2]^2+x[3]^2+x[4]^2

and "nrm(v)" will be the norm of v:

> nrm(X);

(x[1]^2+x[2]^2+x[3]^2+x[4]^2)^(1/2)

Example: Gram-Schmidit in 4 dimensional space.

Problem: apply the gram-schmidt process to the following 4 vectors.

> v1 := <1, 1, 0, 2>; v2 := <1, 0, 1, -1>; v3:=<1, 0, 0, 1>; v4 := <1, 1, 2, 1>;

v1 := Vector[column]([[1], [1], [0], [2]])

v2 := Vector[column]([[1], [0], [1], [-1]])

v3 := Vector[column]([[1], [0], [0], [1]])

v4 := Vector[column]([[1], [1], [2], [1]])

Check that these vectors form a basis:

> M := <v1|v2|v3|v4>;

M := Matrix([[1, 1, 1, 1], [1, 0, 0, 1], [0, 1, 0, 2], [2, -1, 1, 1]])

> Determinant(M);

3

Begin the Gram-Schmidt Process:

> u1 := v1/nrm(v1);

u1 := Vector[column]([[1/6*6^(1/2)], [1/6*6^(1/2)], [0], [1/3*6^(1/2)]])

> nrm(v1);

6^(1/2)

> u2p := v2 - ip(v2, u1)*u1;

u2p := Vector[column]([[7/6], [1/6], [1], [(-2)/3]])

> nrm(u2p);

1/6*102^(1/2)

> u2 := u2p/nrm(u2p);

u2 := Vector[column]([[7/102*102^(1/2)], [1/102*102^(1/2)], [1/17*102^(1/2)], [-2/51*102^(1/2)]])

> u3p := v3 - ip(v3, u1)*u1 - ip(v3, u2)*u2;

u3p := Vector[column]([[5/17], [(-9)/17], [(-3)/17], [2/17]])

> nrm(u3p);

1/17*119^(1/2)

> u3 := u3p/nrm(u3p);

u3 := Vector[column]([[5/119*119^(1/2)], [-9/119*119^(1/2)], [-3/119*119^(1/2)], [2/119*119^(1/2)]])

> u4p := v4 - ip(v4, u1)*u1 -ip(v4, u2)*u2 -ip(v4,u3)*u3;

u4p := Vector[column]([[(-3)/7], [(-3)/7], [6/7], [3/7]])

> nrm(u4p);

3/7*7^(1/2)

> u4 := u4p/nrm(u4p);

u4 := Vector[column]([[-1/7*7^(1/2)], [-1/7*7^(1/2)], [2/7*7^(1/2)], [1/7*7^(1/2)]])

Check that these are really orthogonal to each other:

> ip(u1, u2);

0

> ip(u1, u3);

0

> ip(u1, u4);

0

> ip(u2, u3);

0

> ip(u2, u4);

0

> ip(u3, u4);

0

>

>

Modified Gram-Schmidt on the same vectors:

>

> u1p := v1;

u1p := Vector[column]([[1], [1], [0], [2]])

> u2p := v2 - ip(v2, u1p)*u1p/nrm2(u1p);

u2p := Vector[column]([[7/6], [1/6], [1], [(-2)/3]])

> u3p := v3 - ip(v3,u1p)*u1p/nrm2(u1p) -ip(v3, u2p)*u2p/nrm2(u2p);

u3p := Vector[column]([[5/17], [(-9)/17], [(-3)/17], [2/17]])

> u4p := v4 - ip(v4, u1p)*u1p/nrm2(u1p)-ip(v4,u2p)*u2p/nrm2(u2p) -ip(v4, u3p)*u3p/nrm2(u3p);

u4p := Vector[column]([[(-3)/7], [(-3)/7], [6/7], [3/7]])

> u1 := u1p/nrm(u1p);

u1 := Vector[column]([[1/6*6^(1/2)], [1/6*6^(1/2)], [0], [1/3*6^(1/2)]])

> u2 := u2p/nrm(u2p);

u2 := Vector[column]([[7/102*102^(1/2)], [1/102*102^(1/2)], [1/17*102^(1/2)], [-2/51*102^(1/2)]])

> u3 := u3p/nrm(u3p);

u3 := Vector[column]([[5/119*119^(1/2)], [-9/119*119^(1/2)], [-3/119*119^(1/2)], [2/119*119^(1/2)]])

> u4 := u4p/nrm(u4p);

u4 := Vector[column]([[-1/7*7^(1/2)], [-1/7*7^(1/2)], [2/7*7^(1/2)], [1/7*7^(1/2)]])

>

Consider the following vector:

> w := <1, 2, 4, 1>;

w := Vector[column]([[1], [2], [4], [1]])

>

Calculate the projection of w onto the span of u1 and u2:

> p := ip(w,  u1)*u1 + ip(w, u2)*u2;

p := Vector[column]([[48/17], [19/17], [29/17], [9/17]])

Calculate the matrix of the projection onto the span of u1 and u2:

> U := <u1|u2>;

U := Matrix([[1/6*6^(1/2), 7/102*102^(1/2)], [1/6*6^(1/2), 1/102*102^(1/2)], [0, 1/17*102^(1/2)], [1/3*6^(1/2), -2/51*102^(1/2)]])

> Pr := U.Transpose(U);

Pr := Matrix([[11/17, 4/17, 7/17, 1/17], [4/17, 3/17, 1/17, 5/17], [7/17, 1/17, 6/17, (-4)/17], [1/17, 5/17, (-4)/17, 14/17]])

>

Check that this gives the same result for the projection of w:

> Pr.w;

Vector[column]([[48/17], [19/17], [29/17], [9/17]])

Example in 3 dimensions

> v1 := <1, 1, 0>; v2 := < 0, 1, 1>; v3 := <1, -1, 1>;

v1 := Vector[column]([[1], [1], [0]])

v2 := Vector[column]([[0], [1], [1]])

v3 := Vector[column]([[1], [-1], [1]])

> M := <v1|v2|v3>;

M := Matrix([[1, 0, 1], [1, 1, -1], [0, 1, 1]])

> Determinant(M);

3

> u1 := v1/nrm(v1);

u1 := Vector[column]([[1/2*2^(1/2)], [1/2*2^(1/2)], [0]])

> u2p := v2 - ip(v2, u1)*u1;

u2p := Vector[column]([[(-1)/2], [1/2], [1]])

> u2 := u2p/nrm(u2p);

u2 := Vector[column]([[-1/6*6^(1/2)], [1/6*6^(1/2)], [1/3*6^(1/2)]])

> u3p := v3 - ip(v3, u1)*u1 - ip(v3, u2)*u2;

u3p := Vector[column]([[1], [-1], [1]])

> u3 := u3p/nrm(u3p);

u3 := Vector[column]([[1/3*3^(1/2)], [-1/3*3^(1/2)], [1/3*3^(1/2)]])

> u1, u2, u3;

Vector[column]([[1/2*2^(1/2)], [1/2*2^(1/2)], [0]]), Vector[column]([[-1/6*6^(1/2)], [1/6*6^(1/2)], [1/3*6^(1/2)]]), Vector[column]([[1/3*3^(1/2)], [-1/3*3^(1/2)], [1/3*3^(1/2)]])

> ip(u1, u2);

0

> ip(u1, u3);

0

> ip(u2, u3);

0

>

>

>

Example in the space C[0,1]

>

We use the inner product: `<,>`(f, g) = Int(f(x)*g(x), x = 0 .. 1)

> ip := proc(f,g) int(f*g, x=0..1); end;

ip := proc (f, g) int(f*g, x = 0 .. 1) end proc

> ip(x,x);

1/3

Consider the subspaces spanned by the following polynomials.  Use Gram-Schmidt to find an ON basis of this subspace.

> v1 := 1; v2 := x; v3 := x^2; v4 := x^3;

v1 := 1

v2 := x

v3 := x^2

v4 := x^3

>

>

> u1 := v1/nrm(v1);

u1 := 1

> u2p := v2 - ip(v2,u1)*u1;

u2p := x-1/2

> u2 := u2p/nrm(u2p);

u2 := 2*(x-1/2)*3^(1/2)

> u3p := v3 - ip(v3, u1)*u1 -ip(v3, u2)*u2;

u3p := x^2+1/6-x

> u3 := u3p/nrm(u3p);

u3 := 6*(x^2+1/6-x)*5^(1/2)

> u4p := v4 - ip(v4,u1)*u1 - ip(v4, u2)*u2 -ip(v4,u3)*u3;

u4p := x^3-1/20+3/5*x-3/2*x^2

> u4 := u4p/nrm(u4p);

u4 := 20*(x^3-1/20+3/5*x-3/2*x^2)*7^(1/2)

>

>

Same problem, modified Gram-Schmidt.

>

> u1p := v1;

u1p := 1

> u2p := v2 - ip(v2, u1p)*u1p/nrm2(u1p);

u2p := x-1/2

> u3p := v3 - ip(v3, u1p)*u1p/nrm2(u1p) - ip(v3, u2p)*u2p/nrm2(u2p);

u3p := x^2+1/6-x

> u4p := v4 - ip(v4, u1p)*u1p/nrm2(u1p) - ip(v4, u2p)*u2p/nrm2(u2p) - ip(v4, u3p)*u3p/nrm2(u3p);

u4p := x^3-1/20+3/5*x-3/2*x^2

> u1 := u1p/nrm(u1p);

u1 := 1

> u2 := u2p/nrm(u2p);

u2 := 2*(x-1/2)*3^(1/2)

> u3 := u3p/nrm(u3p);

u3 := 6*(x^2+1/6-x)*5^(1/2)

> u4 := u4p/nrm(u4p);

u4 := 20*(x^3-1/20+3/5*x-3/2*x^2)*7^(1/2)

Plot these othogonal polynomials to see what they look like:

> plot([u1, u2, u3, u4], x=0..1);

[Plot]

Check that u3 and u4 are othogonal:

> Int(u3*u4,x=0..1);

Int(120*(x^2+1/6-x)*5^(1/2)*(x^3-1/20+3/5*x-3/2*x^2)*7^(1/2), x = 0 .. 1)

> value(%);

0

Check that u4 is a unit vector:

> plot(u4^2, x=0..1);

[Plot]

> Int(u4^2, x=0..1);

Int(2800*(x^3-1/20+3/5*x-3/2*x^2)^2, x = 0 .. 1)

> value(%);

1

>