Task

Solve   Ax=b   using Gaussian elimination then backwards substitution.

A   being an   n by n   matrix.

Also,   x and b   are   n by 1   vectors.

To improve accuracy, please use partial pivoting and scaling.



See also

   the Wikipedia entry:   Gaussian elimination


with javascript_semantics
function gauss_eliminate(sequence a, b)
    {a, b} = deep_copy({a,b})
    integer n = length(b)
    atom tmp
    for col=1 to n do
        integer m = col
        atom mx = a[m][m]
        for i=col+1 to n do
            tmp = abs(a[i][col])
            if tmp>mx then
                {m,mx} = {i,tmp}
            end if
        end for
        if col!=m then
            {a[col],a[m]} = {a[m],a[col]}
            {b[col],b[m]} = {b[m],b[col]}
        end if
        for i=col+1 to n do
            tmp = a[i][col]/a[col][col]
            for j=col+1 to n do
                a[i][j] -= tmp*a[col][j]
            end for
            a[i][col] = 0
            b[i] -= tmp*b[col]
        end for
    end for
    sequence x = repeat(0,n)
    for col=n to 1 by -1 do
        tmp = b[col]
        for j=n to col+1 by -1 do
            tmp -= x[j]*a[col][j]
        end for
        x[col] = tmp/a[col][col]
    end for
    return x
end function
 
constant a = {{1.00, 0.00, 0.00,  0.00,  0.00,   0.00},
              {1.00, 0.63, 0.39,  0.25,  0.16,   0.10},
              {1.00, 1.26, 1.58,  1.98,  2.49,   3.13},
              {1.00, 1.88, 3.55,  6.70, 12.62,  23.80},
              {1.00, 2.51, 6.32, 15.88, 39.90, 100.28},
              {1.00, 3.14, 9.87, 31.01, 97.41, 306.02}},
         b = {-0.01, 0.61, 0.91,  0.99,  0.60,   0.02}
 
pp(gauss_eliminate(a, b))
