Numerical Methods & Scientific Computing: lecture notes
Root-finding
Week 6: aim to cover
Numerical linear algebra: Gauss Elimination with Partial Pivoting (GEPP), operations count (Lecture 11)
Newton’s method, 2D arrays (matrices) in MATLAB (Lab 6) LU factorization, special matrices (Lecture 12)
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Trefethen’s Maxims
In principle, the Taylor series of a function of n variables involves an n-vector, an n⇥n matrix, an n⇥n⇥n tensor, and so on. Actual use of orders higher than two, however, is so rare that the manipulation of matrices is a hundred times better supported in our brains and in our software tools than that of tensors.
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
The problem
Given A: n ⇥ n square matrix
b: n ⇥ 1 matrix (column vector) find x (column vector) where
Ax = b
We assume A is nonsingular so a unique solution exists.
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Direct vs. iterative methods
A
dense
sparse
(many entries = 0)
direct methods (only roundoff)
iterative methods truncation error also
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
In MATLAB
it’s very easy to solve a linear system of equations
1 define the matrix A
e.g. A = rand(100,100)
2 define the RHS column vector
e.g. b = rand(100,1)
3 then solve with backslash e.g. x = A\b
MAGIC
did it solve correctly?
compute b – A*x
We now explore how this magic is performed …
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Gauss elimination
Recall from Linear Algebra:
Use row operations to reduce augmented system to upper triangular form then solve by back-substitution.
Note: no need for Gauss-Jordan elimination (reduction to reduced row echelon form) — it takes more work!
As an algorithm:
Denote original A by A(1) = [A(1)]
original b by b(1) = [b(1)]T
ij i
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Step 1:
Assume A(1) 6= 0 11
A(1) i1 ,
Define multipliers l Put zeroes below A11:
i, j = 2..n b(2)=b(1) li1b(1), i,=2..n
0A(1) A(1) ··· A(1) 1 0b(1) 1 B 11 12 1n C B 1 C
B0A(2)···A(2)C Bb(2)C B . 22 2n Cx=B 2. C
@ . A @ . A 0 A(2) ··· A(2) b(2)
n2 nn n Repeat this step.
We now have
i1
A(1) 11
i = 2..n
=
A(2) = A(1) li1A(1),
ij ij 1j ii1
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
After k steps
A(k)x = b(k)
0 A(1) A(1) ··· ··· A(1) 1
B1112. 1nC B0A(2) ..···A(2)C
B B@ . . . . . . A ( k ) · · · A ( k ) C CA
Assume A(k) 6= 0 kk
B22 2nC A(k) =B . 0 … ··· . C
kk
kn
0 0 A(k) ··· A(k) nk nn
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Define multipliers l Put zeroes below Akk :
A(k )
= ik , i = k +1..n
A(k+1)=A(k) likA(k), i,j=k+1..n ij ij kj
b(k+1) = b(k) likb(k), i = k + 1..n iik
ik
A(k) kk
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
After n steps
we have
0(1) (1)1 0b(1)1 BA11 ···A1nC B1.C
@0 Ax=@.A 0 ··· A(n) b(n)
nn n
where A(n) is upper triangular; let’s call it U for ‘upper’
Ux=b(n) =g
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Triangular system
Solve this by back-substitution
xn = gn/Unn
1 Xn
xk=U [gk Ukjxj], k=n 1..1
kk j=k+1
These formulae ! algorithm for Gauss elimination (ready for coding) BUT
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Pivoting
1 what if A(k) = 0? kk
A(k) is the pivot kk
Remedy: swap rows
2 in order to put zeroes under pivot, we subtract rows
=) we can amplify roundo↵ error if we subtract 2 large numbers
(most common if the multiplier is large)
Gauss elimination is potentially vulnerable to subtractive cancellation!
an extreme 2 ⇥ 2 case
Example
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Partial pivoting
Remedy: choose pivot so that multipliers are never large — called a pivoting strategy
Simplest pivoting strategy (and usually enough)
Partial pivoting
At step k,
look at elements A(k) and below
choose the largest of these in magnitude to be the new pivot e.g.
A(k ) lk
swap rows l and k
=) multipliers formed from new pivot satisfy | lik |< 1
kk
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
This also handles zero pivots
Note: don’t actually swap rows; swap pointers or row indices
Example
Back to our extreme case
Other pivoting strategies exist but partial pivoting is regarded as usually su ciently stable.
Gauss Elimination with Partial Pivoting (GEPP) is the default algorithm for solving linear systems.
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Operations Count for Gaussian Elimination
Measure by number of multiply/divides.
Assume no pivoting required.
1st stage: for each i = 2..n A(1)
forml = i1 1÷ i1 A(1)
(2) (1) 11 (1)
Aij =Aij li1A1j , i,j=2..n n 1⇥
b(2) = b(1) li1b(1) 1⇥ ii1
! (n 1)(n+1)⇥/÷ for 1st stage (inc. n 1 for b)
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
2nd stage
:
n 7! n 1
! (n 2)(n)⇥/÷ for 2nd stage (inc. n 2 for b)
.
! (1)(3) ⇥ /÷ for last stage (inc. 1 for b)
last stage: n = 2
Total so far: Xn 1 1 5
(k2 1) = 3n3 + 2n2 6n k=1
inc. Pn 1 k = 1 n2 1 n for b k=1 2 2
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Solve triangular system
Now back-substitute
xn =gn/Unn 1÷
xn 1 = [gn 1 Un 1,nxn]/Un 1,n 1 .
x1 =···
2 ⇥ /÷ n⇥/÷
! total for back-substitution Pn k = 1 n2 + 1 n k=1 2 2
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Total work
1n3 + 1n2 5n + 1n2 + 1n 32622
= 1n3 + n2 1n ⇡ 1n3 333
of which
comes from processing RHS vector b
1n2 1n + 1n2 + 1n = n2 2222
Important: Work for Gauss elimination: ⇡ 1 n3 operations. Work for each RHS vector = n2 operations3
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
Many RHS vectors
So to solve
requires ⇡ 1 n3 + kn2 operations.
Ax = (b1,b2,··· ,bk)
By comparison: suppose we use inverse of matrix
3
To find A 1, solve
takes ⇡ 1n3 +n·n2 ⇡ 4n3 operations
x = A 1b
Ax = I = (e1,e2,··· ,en)
333
actually can be done in ⇡ n operations
Moral: DON’T FORM MATRIX INVERSE — 3 times more work than solving the system!
Numerical Methods & Scientific Computing: lecture notes
Linear Systems
End of Lecture 11