I have been asked to write the Mathematica code to solve a 25x25 Hilbert matrix. The built-in function LinearSolve
would not work.
I started my solution by coding a classical Gaussian elimination:
m = 25;
HM1[n_] := Table[1/(i + j - 1), {i, n}, {j, n}];
a = HM[m];
b = Table[Random[],{m}];
A = MapThread[Append,{a,b}];
Table[
f = A[[i, j]]/A[[j, j]]; A[[i]] = A[[i]] - A[[j]]*f,
{j, 1, m - 1}, {i, j + 1, m}];
S = Array[0 &, m]
Table[
A[[i, m + 1]] = A[[i, m + 1]] - Sum[A[[i, j]]*S[[j]], {j, i, m}];
S[[i]] = A[[i, m + 1]]/A[[i, i]],
{i, m, 1, -1}];
However, I am totally lost in doing the partial pivoting for the matrix. Any help?