kill(all);

maxima:gram.schmidt and upper triangle: solve linear equations,version 1.0
peter.vlasschaert@gmail.com,08/04/2022

youtube :

ref:https://www.youtube.com/watch?v=9tMDEEm-pJw&list=
PLOROtRhtegr4HyXIIHKHketOm2k3JWSpK&index=15

load packages , option
---------------------------

load(eigen)$

ratprint:false$

given: linear system: Mx = b
----------------------------------

→ 1e step : system of three linear equations

p1:1*x+3*y+3*z = 1;
p2:2*x+2*y-2*z = 2;
p3:-2*x+2*y+1*z =3;

unknowns :

p4:listofvars(p1);

p4[1];

p4a:transpose(p4);

→ text string, concat : connect two symbols (p,i).

p4aa:makelist(ev(concat(p,i),nouns),i,1,length(p4));

→ evaluation: text string.

p4bb:ev(p4aa, nouns);

→ 2e step : matrix ,system.

matrix : M

p5a:coefmatrix(p4bb,p4);

p5b:augcoefmatrix(p4bb,p4);

rem : fourth column , negative sign front,we 
want remove it.

rhs :system of equations.

p5b1:col((p5b)*(-1),4);

new : augmented matrix (use:p5a)

p5b2:addcol(p5a,p5b1);

→ last column , unknowns from :1e step.

p5c:addcol(p5b,p4a);

→ total matrix: with all column,corrected

use:p5b2 & p4a

p5c1:addcol(p5b2,p4a);

→ 3e step : matrix , linear system.

simp;           /* default : simp*/

simp:false$ /* automatic evaluation : stop*/

use:p5a,p4a,p5b1

p6:p5a.p4a =p5b1;

representation: M.x =b , see "above"

→ start : again evaluation.

simp;

simp:true$

→ 4e step : Gram-Schmidt wrong result : problem.

ref :https://maths.cnam.fr/Membres/wilk
/MathMax/help/Maxima/maxima_25.html

part a:  → Gram-Schmidt orthogonalization

use : p5a

p7: gramschmidt (p5a);

norm(xx) := sqrt(xx . xx);

p7cc:makelist(p7[i]/norm(p7[i]),i,1,3);

matrix: representation:"U"

p7a:transpose(apply(matrix,p7cc));

p7a1:length(p7a);/* number of rows*/
p7a2:length(transpose(p7a));/*number of columns */

float(combine(determinant(p7a)));

unitair matrix: GM

p7cc1:apply(matrix,p7cc);
p7cc2:float(combine(determinant(p7cc1)));

check : "dot product,innerproduct"

map (innerproduct, [p7a[1], p7a[2], p7a[3]], [p7a[2], p7a[3], p7a[1]]);

ev(p7a[1].p7a[2],nouns);

rem : vect package declares the dot operator . to be a commutative operator. 

part b : → unknown matrix "T"

use :{p7a1,p7a2}

use : Menu → Generate Matrix :
Array     : T
rows       : p7a1
columns  :p7a2
Name       :T

T: genmatrix(T, 3, 3);

p7b:listofvars(T);

→ step 5: implimentation :Gram-Schmidt

W[k]:=V[k]-sum((U[j].V[k])*U[j],j,1,k-1);

column space(p5a)

makelist(V[i]:col(p5a,i),i,1,3);

computation :{W[i],U[i]},i=1,2,3

m1:W[1]$
U[1]:m1/norm(m1);

m2:W[2]$
U[2]:m2/norm(m2);

m3:W[3]$
U[3]:m3/norm(m3);

result :{p5a},V[i],i=1,2,3

makelist(V[i]:col(p5a,i),i,1,3);

→ step 6e : T[i,j], upper triangle.

part a: lower triangle T[i,j]=0

[T[2,1]:0,T[3,1]:0,T[3,2]:0];

→ input values : T matrix, "see:above"

q1:T;

→ find T[i,j]:V[i].U[j]

S[i,j]:=V[i].U[j];

T[1,1]:S[1,1];

q1a:q1;

for i:1 thru 3 do(
 
for j:1 thru 3 do(
  if i<j then  T[i,j]:S[i,j]        
  else T[i,j]:S[i,j]      
    )
)$

q1b:transpose(q1);

→ determinant :q1b

q1b1:determinant(q1b);

      determinant : product of the trace "q1b"

q1b2:product(q1b[i,i],i,1,length(q1b));

rem : determinant = product of trace. 
1e) diagonaal        matrix
2e) upper triangle matrix (see:q1b)
3e) lower triangle matrix

info:

column:U[i] → list

f1:flatten(args(U[1]));

→ construct list 

flist:makelist(flatten(args(U[i])),i,1,length(f1));

fmatrix:transpose(apply(matrix,flist));

fmatrix : unitair matrix,det(fmatrix) = 1

f2:determinant(fmatrix);

check : p5a ? fmatrix.uppertriangle.
------------------------------------------

simp:false;

→ decomposition "p5a" = z1

z1:fmatrix.q1b;

simp:true;/* start : re-evaluation*/

z2:z1;

is(z2=p5a);

solution : system of equations " = linear system"
----------------------------------------------------------
x:unknowns " = variable"

M = matrix

M.x=b
decomposition:M
M = U.T, U =unitair matrix → U^(-1) = transpose(U)
U = Gram-Schmidt : orthogonalization,norm ,see code : above.
T  = upper triangle matrix
b   = rhs system of equations
x   = unknown variables
U^(-1).U = I unit matrix.
            U.T.x = b
         U^(-1). U.T.x  = U^(-1) b
                           I.x  = transpose(U).b
                           T.x = transpose(U).b  (**)

rem : video transpose(U) = U*,complex numbers

simp:false;

representation :(**)
T ="q1b"
U ="fmatrix"
b  = "p5b1"
x  = "p4"

c1:q1b.transpose(p4) = transpose(fmatrix).p5b1;

simp:true$

c2:c1;

principal : backsubstitution→[x,y,z]
--------------------------------------------
rem : first : z,last equation,then up .
here : we will use solve.

→ lhs : c2

c3a:flatten(args(lhs(c2)));

→ rhs : c2

c3b:flatten(args(rhs(c2)));

c3c:map("=",c3a,c3b);

solve:{x,y,z}

c3d:solve(c3c,p4)[1];

check : is good.

d41:subst(c3d,c3c);

rem : see video ( det =determinant)
------------------
GENERAL :det(transpose(A)) =det(A)

rem :M:[V[1],V[2],V[3]],see: example above,video.

→det(M) =det(U*T) = det(U)*det(T)
→U^(-1)=U*
→det(U*) = 1/det(U*)
rem :det(U*) = complexconjugate(U)
det(U)*det(complexconjugate(U)) =1
ref : complex number , i=sqrt(-1)
(a+ib)*(a-ib) =a^2+b^2
(abs(det(U)))^2 = 1
→ abs(det(U)) = 1, abs: no sign 
see:abs
abs(det(M)) =abs(product(T[i,i],i,1,n)),trace
abs(det(M)) =abs(product(V[i].U[i],i,1,n))
see :dot product,V[i].U[i] =T[i,i]

rem : V[i]→U[i],both columns, i=1,2,3; example.