(* This is the Gaussian Elimination solver made symbolic, and its check;
see the file gaussian_3.math for comments on its basic construction.
The data are the matrix and vector used in the lecture.  *)

a = { {  4.2 ,      2.2 + p , -3.9 ,      9.3 ,      0.1     } ,
      {  8.6 - p ,  r ,        0.7 ,     -2.3 ,     -0.3     } ,
      {  8.4 ,     -5.9 ,     -8.1 + q ,  9.6 ,      3.8     } ,
      { -0.8 ,     -9.4 ,     -9.9 ,      9.9 ,      5.0 - r } ,
      { -1.3 ,     -8.1 ,      0.6 ,     -9.2 + r , -7.3     } } ;

c = { p , 1.0 , p + q , 1.0 , q } ;

(* We apply t2 recursively every time we do a division, as described in
the lecture.  We use Expand and apply t1 to clobber high-order terms
following every multiplication and division+t2 - they obviously can't be
created by an addition or subtraction. *)

symsolver [ a_ , c_ , t1_ , t2_ ] :=
        Module [ { a1 , c1 , i, j, n, x } ,
    a1 = a ; c1 = c ;
    n = Dimensions [ a1 ] [[ 1 ]] ;
    p = p1 * e ; q = q1 * e ; r = t1 * e ;
    For [ j = 1 , j < n , ++ j ,
        For [ i = j + 1 , i <= n , ++ i ,
            x = Expand [ ( a1 [[ i , j ]] / a1 [[ j , j ]] )
                //. t2 ] /. t1 ;
            a1 [[ i , ;; ]] = a1 [[ i , ;; ]] -
                Expand [ x * a1 [[ j , ;; ]] ] /. t1 ;
            c1 [[ i ]] = c1 [[ i ]] -
                Expand [ x * c1 [[ j ]] ] /. t1 ;
        ]
    ] ;
    For [ j = n , j > 1 , -- j ,
        For [ i = 1 , i < j , ++ i ,
            x = Expand [ ( a1 [[ i , j ]] / a1 [[ j , j ]] )
                //. t2 ] /. t1 ;
            a1 [[ i , ;; ]] = a1 [[ i , ;; ]] -
                Expand [ x * a1 [[ j , ;; ]] ] /. t1 ;
            c1 [[ i ]] = c1 [[ i ]] -
                Expand [ x * c1 [[ j ]] ] /. t1 ;
        ]
    ] ;
    Expand [ ( c1 / Diagonal [ a1 ] ) //. t2 ] /. t1
] ;

(* These are the rules we shall use, assuming that the infinitesimal is
called e ; the first clobbers all second-order terms and above, and the
second converts a division into multiplicative form. *)

rule1 = e ^ k_ /; k > 1 -> 0 ;
rule2 = k1_ / ( k2_ + k3_ * e ) -> \
    k1 * ( 1 - ( k3 / k2 ) * e ) / k2 ;

(* We convert p, q and r to multiples of the infinitesimal, e, apply the
solver, cancel the conversions, and print the result converted back to
the original variables. *)

p = p1 * e ; q = q1 * e ; r = r1 * e ;
z = symsolver [ a , c , rule1 , rule2 ] ;
p =. ; q =. ; r =. ;
Print [ z /. e p1 -> p /. e q1 -> q /. e r1 -> r ] ;

(* We do exactly the same process, but print the solution multiplied by
the matrix to check that we get back to the original vector, rather than
solving the equations.  We need to use Expand and rule1 again, as we may
reintroduce some second-order terms by doing this. *)

p = p1 * e ; q = q1 * e ; r = r1 * e ;
s = Expand [ a . z ] /. rule1 ;
p =. ; q =. ; r =. ;
t = s /. e p1 -> p /. e q1 -> q /. e r1 -> r
Print [ t ] ;

(* We print the result with tiny floating-point values converted to
zero, because it clarifies the result; this is not something that should
be done if you are going to reuse the result. *)

Print [ t /. x_ /; Abs[x] < 1.0*^-11 -> 0 ] ;