Quick guide to Mathematica

To get help, place the cursor on the function and press F1
To run all commands in one cell, place the cursor there and press Shift + Enter or Enter on the numerical keyboard.
Output of many commands is suppressed by ; at the end of the command. If you want to see the output, delete ;
Argument of functions must be in square brackets [...], braces {...} are for arrays, ranges etc., use double brackets [[...]] to access elements of arrays
expr /. {x -> a} means substitute a for x in expr, expr // function means apply function to expr, i.e. it is equivalent to function[expr].
D[f[x],x] = derivative of f[x] with respect to x.
To define a function of x one can use f[x_] := 1+(x^2), notice the underscore and the colon, or # nad & like in f := 1+((#)^2) &
N[expr] evaluates expr with machine precision, if you use the decimal point in the expression it will be also evaluated with machine precision
Special characters can be inserted by pressing Esc ... Esc, e.g. Esc p Esc gives π
Clear[...] is used to unset any variables which could have been assigned previously.
Some other useful commands: Simplify, Expand, Factor; Integrate, Series, Sum

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Numerical solution of a system of linear equations

Run Basic.Functions.for.TMF057.m first !

Failure of LU decomposition without pivoting - simple example

Incorrect result due to round-off errors, solution should be {-1,1}:

A={{1.0 10^(-20),1.0},{1.0,1.0}};
b={1.0,0.0};
LU=MyLUDecomposition[A];
LU//MatrixForm
MyLUBackSubstitution[LU,b]

Correct solution by the built-in functions with pivoting

LUDecomposition[A]
LUDecomposition[A][[1]]//MatrixForm
LinearSolve[A,b]

Exact arithmetic:

A={{10^(-20),1},{1,1}};
b={1,0};
LU=MyLUDecomposition[A];
LU//MatrixForm
sol=MyLUBackSubstitution[LU,b]
N[sol]

QR factorization:

A={{1.0 10^(-20),1.0},{1.0,1.0}};
b={1.0,0.0};
{Q,R}=MyQRHouseholder[A];
Print[ConjugateTranspose[Q]//MatrixForm,R//MatrixForm]
MySolveUpperTriangular[R,Q.b]

LU decomposition - growth factor - the worst case

Detailed discussion in Trefethen, Bau: Numerical Linear Algebra, Lecture 22
First, only 5x5 matrix, everything is OK.

Clear[A];
n=5;
A=LowerTriangularize[ConstantArray[-1,{n,n}]]+2IdentityMatrix[n];
A[[All,n]]=1;
A//MatrixForm
LU=MyLUDecomposition[A];
LU//MatrixForm
Print["Growth factor: ",Max[LU]]

Second, the same matrix but this time 60x60 matrix. Note that Mathematica returns an incorrect solution without notice.

Clear[A];
(* Construct the matrix and determine the growth factor *)
n=60;
A=N[LowerTriangularize[ConstantArray[-1,{n,n}]]+2IdentityMatrix[n]];
A[[All,n]]=1.0;
LU=MyLUDecomposition[A];
Print["Growth factor: ",Max[N[LU]]]
(* Solve a problem with the known solution v *)
Print["Exact solution: "]
v=ConstantArray[1.0,{n}]
Print["Right-hand side: "]
b =A.v
Print["My numerical solution: "]
MyLUBackSubstitution[LU,b]
Print["Mathematica numerical solution: "]
x=LinearSolve[A,b]
Print["Checking the solution, Ax-b should be zero:"]
A.x-b
(* Determine the matrix condition number *)
{U,S,V}=SingularValueDecomposition[A];
σ=Diagonal[S];
κA=Max[σ]/Min[σ];
Print["Condition number: ",κA]

Solution by QR factorization:

Clear[A];
(* Construct the matrix and determine QR factorization *)
n=60;
A=N[LowerTriangularize[ConstantArray[-1,{n,n}]]+2IdentityMatrix[n]];
A[[All,n]]=1.0;
{Q,R}=MyQRHouseholder[A];
(* Print["A = ",A//MatrixForm,",  Q.R = ",ConjugateTranspose[Q].R//MatrixForm] *)
(* Solve a problem with the known solution v *)
Print["Exact solution: "]
v=ConstantArray[1.0,{n}]
Print["Right-hand side: "]
b =A.v
Print["Numerical solution: "]
MySolveUpperTriangular[R,Q.b]

LU decomposition - growth factor - random case

Detailed discussion in Trefethen, Bau: Numerical Linear Algebra, Lecture 22

Clear[A];
(* Construct the matrix and determine the growth factor *)
n=2500;
A=RandomReal[{-1.0,1.0},{n,n}];
{LU,P,c}=LUDecomposition[A];
Print["Growth factor: ",Max[LU]/Max[A]]
(* Solve a problem with the known solution v *)
(* Print["Exact solution: "] *)
v=RandomReal[{0.0,1.0},n];
(* Print["Right-hand side: "] *)
b =A.v;
x=LinearSolve[A,b];
Print["Error of Mathematica numerical solution: ", Norm[x-v]]

{U,S,V}=SingularValueDecomposition[A];
σ=Diagonal[S];
κA=Max[σ]/Min[σ];
Print["Condition number: ",κA]