# Programmierpraktikum

Claudius Gros, SS2012

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

# 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)$
• $$\left( \begin{array}{rrrr|r} 2&-3&-1& 2&4\\ 4&-4&-1& 4&4\\ 2&-5&-3& 3&9\\ \end{array} \right)$$
• 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
• $$\left( \begin{array}{ccccc} *&0&0&0&0\\ 0&*&0&0&0\\ 0&0&*&0&0\\ 0&0&0&*&0\\ 0&0&0&0&*\\ \end{array} \right)$$
• tri-diagonal
• $$\left( \begin{array}{ccccc} *&*&0&0&0\\ *&*&*&0&0\\ 0&*&*&*&0\\ 0&0&*&*&*\\ 0&0&0&*&*\\ \end{array} \right)$$
• block-diagonal (blocks are independent
• $$\left( \begin{array}{ccccc} *&*&0&0&0\\ *&*&0&0&0\\ 0&0&*&*&0\\ 0&0&*&*&0\\ 0&0&0&0&*\\ \end{array} \right)$$
• L (lower left) and R (upper right) triangular matrices
• $$L = \left( \begin{array}{ccccc} *&0&0&0&0\\ *&*&0&0&0\\ *&*&*&0&0\\ *&*&*&*&0\\ *&*&*&*&*\\ \end{array} \right) \qquad\quad R = \left( \begin{array}{ccccc} *&*&*&*&*\\ 0&*&*&*&*\\ 0&0&*&*&*\\ 0&0&0&*&*\\ 0&0&0&0&*\\ \end{array} \right)$$
• sparse (few non-zero entries, there are special algorithms)
• $$\left( \begin{array}{ccccc} 0&*&0&0&0\\ *&0&0&*&0\\ 0&*&0&0&0\\ *&0&0&0&*\\ 0&0&*&0&0\\ \end{array} \right)$$

# 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)$$