(* intro to linear algebra with Mathematica *)
(* Used in several classes to illustrate the concepts *)
(* Cut and paste the commands to make it easy to follow. *) 

(* ============================================================= *)
(*
(* CLASS I. illustrate the importance of well-conditioning *)


(* define your matrix in  Ax = b *) 


A = {{ 9.7, 6.6}, {4.1, 2.8}}

(* define the right-hand side vector *)

b = {9.7, 4.1} 

(* Solve using Gaussian elimination *) 

x = LinearSolve[A, b] 

Print[x]

(* let's add a bit of noise into b *)

b1 = {9.7, 4.11} 

x1 = LinearSolve[A, b1]

print[x1] 

(* Observe how dramatically different the new solution is *) 
(* This is because the matrix is ill-conditioned. *) 
(* Generally, || Dx ||/|| x || <= Cond|A| * || Db ||/||b||, where *) 
(* Db = "delta"b = change in the input vector, Dx = change in solution *)
(* Note the use of relative changes. In general, if CondA ~ 10^k, you lose *) 
(* k digits of precision in the output relative to the input *)

(* Use the inverse to compute Cond = ||A||*||A^{-1}||. *)
(* Use any norm ||A||, such as ||A||_1 or ||A||_infinity  *)

Ainv = Inverse[A]

(* You will find that Cond(A) ~ 2000 for the matrix A in the above example *)
(* A matrix is well-conditions when Cond(A) ~ 1. *)
(* Note that Det(A) ~ 0 does not necessarily mean that a matrix is *)
(* ill-conditioned, e.g. A = diag(0.1 0.1... 0.1), (100x100).      *)


(*=================Class II============================================ *)
*)

(* Check LU-decomposition. Same as "factorization" *) 

(* input matrix we used in class *)
B = {{10, -7, 0}, {-3, 2, 6}, {5, -1, 5}}

(* show what it looks like in familiar notation *)
MatrixForm[B]

LU = LUDecomposition[B]

(* First row -- combined L and U. *)
(* Second row -- permutation vector used in pivoting *)
(* 3-rd row -- approximate Cond(B)   *)

MatrixForm[LU]

{lu,p,n}=LUDecomposition[B];
l = LowerTriangularize[lu, -1] + IdentityMatrix[3];
MatrixForm[l]

u = UpperTriangularize[lu];
MatrixForm[u]

(* ======================Calss III======================================= *)
(*

(* Now play with eigenvalues *)

M = {{0, 1}, {1, 0}}

MatrixForm[M]

Eigenvalues[M]

MatrixForm[Eigenvectors[M]]

p = CharacteristicPolynomial[M, x]

NSolve[p == 0.0, x]

(* define a random matrix *)

MBIG[n_] := Table[Random[Real, 0.5, 1.5], {i, 1, n}, {j, 1, n}]

MatrixForm[MBIG[3]]

(* Do the timing *) 

Do[Print [i, Timing[N[Eigenvalues[MBIG[2^i] ] ] ] ], {i, 1, 7, 1}] 

(* maximum memory (bytes) used in this session *)
MaxMemoryUsed[]
*)