Tcl Library Source Code

Let "mt1" and "b" be the matrix and vector of a system of linear equations:

mt1*x = b (System 1)

where x = (x1,x2,x3,x4) are the unknows of the system

and "mt2", "b", the matrix and vector of the same system but with the order of the unknows changed:

mt2*y = b (System 2)

where y = (x2,x4,x3,x1)

for the values in the following example, the solutions of System 1 and System 2 differ for solveGauss but not for solvePGauss. Solutions of System 2 with solveGauss and solvePGauss are almost identical and coincide with the solution obtained with numpy





(HM13) 7 % package require math::linearalgebra
1.1.5
(HM13) 7 % set mt1 {{-0.0012585709756272178 -235219735842.05304 117609867921.10306 117609867921.10408} {-0.0007706083587115171 -0.0026517381941796457 42416766269.665436 -42416766269.66554} {1010999714627.5773 974275470171.4105 24625609209.839363 24625609209.839558} {450880107.9198128 1219061933.8972623 38265151.50552135 38265151.50552134}}
{-0.0012585709756272178 -235219735842.05304 117609867921.10306 117609867921.10408} {-0.0007706083587115171 -0.0026517381941796457 42416766269.665436 -42416766269.66554} {1010999714627.5773 974275470171.4105 24625609209.839363 24625609209.839558} {450880107.9198128 1219061933.8972623 38265151.50552135 38265151.50552134}
(HM13) 8 % set mt2 {{-235219735842.05304 117609867921.10408 117609867921.10306 -0.0012585709756272178} {-0.0026517381941796457 -42416766269.66554 42416766269.665436 -0.0007706083587115171} {974275470171.4105 24625609209.839558 24625609209.839363 1010999714627.5773} {1219061933.8972623 38265151.50552134 38265151.50552135 450880107.9198128}}
{-235219735842.05304 117609867921.10408 117609867921.10306 -0.0012585709756272178} {-0.0026517381941796457 -42416766269.66554 42416766269.665436 -0.0007706083587115171} {974275470171.4105 24625609209.839558 24625609209.839363 1010999714627.5773} {1219061933.8972623 38265151.50552134 38265151.50552135 450880107.9198128}
(HM13) 9 % set b {3848.0 0.0 20449.0 13.0}
3848.0 0.0 20449.0 13.0
(HM13) 10 % math::linearalgebra::solveGauss $mt1 $b
1.6982058148140882e-8 2.8514750715307717e-9 1.921064657651454e-8 1.921064657651401e-8
(HM13) 11 % math::linearalgebra::solveGauss $mt2 $b
3.5603980701693165e-9 1.9919569575152068e-8 1.991956957515263e-8 1.5825057309843255e-8
(HM13) 12 % math::linearalgebra::solvePGauss $mt1 $b
1.582505730984325e-8 3.560398070169323e-9 1.9919569575152637e-8 1.9919569575152074e-8
(HM13) 13 % math::linearalgebra::solvePGauss $mt2 $b
3.560398070169323e-9 1.9919569575152078e-8 1.9919569575152637e-8 1.582505730984325e-8

The described issue is actually not a problem as such. The given example is extreme, as the entries in the matrix differ by many orders of magnitude and the solution is very near to the null vector. The example suffers from the finite precision that is inherent in floating-point arithmetic.

Nonetheless, I have added a test case to demonstrate that consistent permutation of the rows does not lead to a different solution.

The difference between solveGauss and solvePGauss is that the first uses the Gauss-Seidel algorithm in a straightforward way to solve the system and the second uses partial pivoting to do this. The effect of pivoting is that the solution is slightly more stable. I do not know how ill-conditioned the matrix is, but looking at the solution (almost the null vector) and the wide variation in coefficients, I would hazard to say that the matrix is quite ill-conditioned.

I agree that you would expect a permutation of the solution's components, but there is no guarantee given the finite precision of floating-point operations.

Can you provide an example with a solution vector that is not near zero?