Programmierpraktikum

Claudius Gros, SS2012

Institut für theoretische Physik
Goethe-University Frankfurt a.M.

Elimination

linear equations

$$ \left.\begin{array}{cccccc} 2x&-3 y& -z & + 2w& =&4\\ 4x&-4 y& -z & + 4w& =&4\\ 2x&-5 y& -3z& + 3w& =&9\\ \end{array}\right\} \qquad\quad A\,\vec x=\vec b \qquad\quad \vec x = \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) $$

$3\times 4$ matrix $\ \ A$
augmented matrix $\ \ (A|\vec b)$

solution via iterative elimination

Gauss elimination

$$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 4&-4&-1& 4&4\\ 2&-5&-3& 3&9\\ \end{array} \right) $$

solve for $x$ in first equation and substitute into the other equations
substract (multiples of) first row

$$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&5\\ 0&-2&-2& 1&5\\ \end{array} \right) $$

solve for $7$ in second equation and substitute into third and forth equation
substract (multiples of) second row

$$ \left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 0& 2& 1& 0&-4\\ 0& 0&-1& 1&1\\ \end{array} \right) $$

in matrix form

$$ A\,\vec x = b \qquad\quad \left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right) $$

back substitution

$$ \left( \begin{array}{rrrr} 2&-3&-1& 2\\ 0& 2& 1& 0\\ 0& 0&-1& 1\\ \end{array} \right) \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} 4 \\ -4 \\ 1 \end{array}\right) $$

last equation

$$ w-z=1 $$

take $z$ as free parameter and substitute back

$$ \begin{array}{rcl} w&=&z+1\\ 2y+z=-4\qquad\quad y &=&-\frac{z}{2}-2 \\ 2x-3y-z+2w=4\qquad\quad x &=&-\frac{5z}{4}-2 \end{array} $$

with

$$ 2x = 4-\frac{3}{2}z-6+z-2z-2 = -4-\frac{5}{2}z $$

final (one-dimensional) solution

$$ \vec x= \left(\begin{array}{c} x \\ y \\ z \\ w \end{array}\right) = \left(\begin{array}{c} -5z/4-2 \\ -z/2-2 \\ z \\ z+1 \end{array}\right) \qquad\qquad \mbox{for any}\ \ z $$

pivoting

Gauss elimination

$$ \left( \begin{array}{cccc} *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ 0&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right) $$

pivoting
select largest element $[*]$ in row

$$ \left( \begin{array}{cccc} *&[*]&*&*\\ *&*&*&*\\ *&*&*&*\\ *&*&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&[*]\\ *&0&*&*\\ *&0&*&*\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&[*]&0\\ *&0&*&0\\ \end{array} \right) \ \Rightarrow\ \left( \begin{array}{cccc} *&*&*&*\\ *&0&*&*\\ *&0&*&0\\ *&0&0&0\\ \end{array} \right) $$

one needs to divide by the pivot selected,
pivot need to be non-vanishing
numerical stable if pivot is large
partial pivoting
only largest element in given row/column
complete pivoting
largest element in matrix

Gauss elimination with pivoting

with swaping of pivot rows

public class GaussElimination {

  private static final double smallNumber = 1e-10;  // a small number

/** Gauss elimination of Ax=b with pivoting.
  * On input quadratic matrix A and vector b.
  * On ouput solution x.
  */
  public static double[] gaussElimination(double[][] A, double[] b) 
    {
//--- missing: testing that A is a NxN matrix
    int N  = b.length;   

//--- loop over all columns A[][p]
    for (int p = 0; p < N; p++) 
      {          
//--- find pivot row and swap rows
      int max = p; 
      for (int i = p + 1; i < N; i++) 
        if (Math.abs(A[i][p]) > Math.abs(A[max][p])) max = i;

//--- swap row A[p][] with A[max][]
      double[] temp = A[p]; 
      A[p] = A[max]; 
      A[max] = temp;

//--- swap elements b[p] with b[max]
      double tt = b[p];  
      b[p] = b[max]; 
      b[max] = tt;

//--- throw a warning if matrix is singular or nearly singular
      if (Math.abs(A[p][p]) <= smallNumber) {  
        throw new RuntimeException("Matrix is singular or nearly singular");
        }

//--- Gauss with pivot within A and b
      for (int i = p + 1; i < N; i++) {
        double alpha = A[i][p] / A[p][p];  // divide by pivot
        b[i] -= alpha * b[p];
        for (int j = p; j < N; j++) 
          A[i][j] -= alpha * A[p][j];
        }
      }  // end of loop over all columns

//--- back substitution
    double[] x = new double[N];
    for (int i = N - 1; i >= 0; i--) {
      double sum = 0.0;
      for (int j = i + 1; j < N; j++) 
         sum += A[i][j] * x[j];
      x[i] = (b[i] - sum) / A[i][i];   // pivots on diagonal after swapping
      }
    return x;
    }  // end of GaussElimination.gaussElimination


  public static void main(String[] args) {
//--- testing
    double[][] A = { { 0, 1,  1 },
                     { 2, 4, -2 },
                     { 0, 3, 15 }
                   };
     double[] b = { 4, 2, 36 };
     double[] x = gaussElimination(A, b);

     for (int i = 0; i < b.length; i++) 
        System.out.println(x[i]);

    }  // end of GaussElimination.main()
}      // end of GaussElimination

add formatted output
add testing

matrices

diagonal
tri-diagonal
block-diagonal (blocks are independent
L (lower left) and R (upper right) triangular matrices
sparse (few non-zero entries, there are special algorithms)

LR decomposition

decompose matrix into product of $L$ and $R$ matrices

$$ A = L\,R \qquad\qquad A\,\vec x =\vec b $$

always possible (Gauss, Hausholder, ...)

$$ L\,R\,\vec x =\vec b \qquad\qquad L\,\vec y =\vec b \qquad\qquad \vec y = R\,\vec x $$

forward substitution: solve for $\ \vec y$

$$ L\,\vec y = \vec b \qquad\qquad L = \left( \begin{array}{cccc} *&0&0&0\\ *&*&0&0\\ *&*&*&0\\ *&*&*&*\\ \end{array} \right) $$

backward substitution: solve for $\ \vec x$

$$ R\,\vec x = \vec y \qquad\qquad R = \left( \begin{array}{cccc} *&*&*&*\\ 0&*&*&*\\ 0&0&*&*\\ 0&0&0&*\\ \end{array} \right) $$