Processing math: 100%

Consider solving the following augmented system

[22c2c112].

Assume we are working in 5 digit aritmetic and c106. We perform GE with partial pivoting. Initially, we do no row swaps. After one step, R21/2R1R2 we have

[22c2c0fl(1c)fl(2c)].

In 5-digit arithmetic fl(1c)=c and fl(2c)=c.

[22c2c0cc].

Performing backward subsitution, we get x2=c/(c)=1 x1=(2c2cx2)/2=0

But

[22c11][01]=[2c1][2c2]

If we instead interchange rows because of the large elements in the first row, we'd perform

[11222c2c].R22R1R2[1120fl(2c2)fl(2c4)].

In 5-digit arithmetic fl(2c2)=2c and fl(2c4)=2c.

[11202c2c]x2=2c/2c=1x1=2x2=1

In 5-digit arithmetic

[22c11][11]=[fl(2+2c)2]=[2c2]

Scaled partial pivoting

Let A be an n×n matrix. After k steps of Gaussian elimination, A has zeros below the first k diagonal entries. With partial pivoting, we look at the entries that are in the (k+1)th column and rows k+1 through n:

A=[a11a12a1n0a22a2n00ak+1,k+1ak+1,k+2ak+1,n00ak+2,k+1ak+2,k+1ak+2,n00an,kan,k+1an,n]

With partial pivoting, we swap rows so as to replace ak+1,k+1 with the largest element (in absolute value) among the set {ak+1,k+1,ak+2,k+1,,an,k+1}.

Often, what is more important is the relative size of the first non-zero entry in the row to the largest entry (in absolute value). Before any row operations define

si=max1jn|ai,j|.

The number si is computed once, is attached to row i, and will need to be swapped when row i is swapped.

One step of scaled partial pivoting is as follows:

Let A be an n×(n+1) matrix.

Compute α=max1in|Ai,1|/si.

Then choose the smallest value of p such that |Ap,1|/sp=α.

Perform R1Rp and s1sp

Then row reduce, below A11.

Full Gaussian elimination with scaled partial pivoting

In [ ]:
function A = GEspp(A)

INPUT: A is an n x m matrix

OUTPUT:  A an n x m upper-triangular matrix, or Inf if the method failed

STEP 1: Define s = max(abs(A(1:n,1:n)')' % the scale factors
STEP 2: For i = 1,2,...,n-1 do STEPS 3-8
  STEP 3: Set p = argmax(A(i:n,i)./s(i:n)) %This is a concise way to find the entry with the largest relative size
  STEP 4: If p = 0 then
    DISPLAY('Method failed: matrix is rank deficient')
    OUTPUT(A);
    STOP.
  STEP 5: Set p = p + i - 1  %Adjust p to correspond to the right row of A
  STEP 6: Do Ri <-> Rp on A  %if you swap rows
  STEP 7: Do si <-> sp on s  %you need to swap scale factors
  STEP 8: For j = i+1,i+2,...,n do STEP 9
    STEP 9: Do R_j - A(j,i)/A(i,i) R_i --> R_j on A
STEP 10: If A(n,n) = 0
  DISPLAY('Method failed: matrix is rank deficient')
  OUTPUT(A)
STEP 11: OUTPUT(A); STOP.
In [23]:
A = [5,2,3; 1,1,1; 4,3,1];
A
max(abs(A'))'
A =

     5     2     3
     1     1     1
     4     3     1


ans =

     5
     1
     4

Scaled partial pivoting is slightly better:

In [41]:
n = 100; A = rand(n); b = rand(n,1);
A(1,:) = (1:n).^4.*A(2,:) + A(1,:); % create a "difficult" matrix
U = GEpp([A,b]);
x = Backsub(U);
max(abs(x-A\b))

U = GEspp([A,b]);
x = Backsub(U);
max(abs(x-A\b))
ans =

   2.5175e-08


ans =

   1.0271e-08

Full (complete) pivoting

Complete pivoting involves both row and column interchanges. Let A be an n×n matrix. After k steps of Gaussian elimination, A has zeros below the first k diagonal entries:

A=[a11a12a1n0a22a2n00ak+1,k+1ak+1,k+2ak+1,n00ak+2,k+1ak+2,k+1ak+2,n00an,kan,k+1an,n]

We swap rows and columns to move the largest element (in absolute value) in the entire bottom right block to (k+1,k+1) entry.

Example

Perform Gaussian elimination with partial pivoting to solve Ax=b where

A=[111221030],b=[111].

[111122110301]unknown = [x1x2x3]
R1R3
[030122111111]unknown = [x1x2x3]
C1C2
[300122111111]unknown = [x2x1x3]
R2+2/3R1R2R3+1/3R1R3
[30010215/30114/3]unknown = [x2x1x3]
R31/2R2R3
[30010215/3001/21/2]unknown = [x2x1x3]
x3=1,x1=(5/3x3)/2=1/3,x2=1/3

Example (where partial pivoting failed)

In 5-digit arithmetic with c>106

[22c2c112]unknown = [x1x2]C1C2
[2c22c112]unknown = [x2x1]R2R1/(2c)R2
[2c22c0fl(1c1)1]unknown = [x2x1]fl(1c1)=1[2c22c011]unknown = [x2x1]x1=1,x2=fl((2c2)/2c)=1

In 5-digit arithmetic, this gives the solution of the system

[22c11][11]=[fl(2+2c)2]=[2c2]

Example (where standard Gaussian elimination failed)

In 5-digit arithmetic with c>106

[c111121]unknown = [x1x2]R1R2C1C2
[2111c11]unknown = [x2x1]R21/2R1R2
[2110fl(c11/2)1/2]unknown = [x2x1]fl(c11/2)=1/2
[21101/21/2]unknown = [x2x1]x1=1,x2=fl((1x1)/2)=1

This is the solution (in 5-digit arithmetic)

[c1112][11]=[fl(1c1)1]=[11]
In [ ]: