程序代写代做代考 C algorithm Numerical Methods & Scientific Computing: lecture notes

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=n1..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 suciently 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 n1⇥ b(2) = b(1) li1b(1) 1⇥ ii1 ! (n1)(n+1)⇥/÷ for 1st stage (inc. n1 for b) Numerical Methods & Scientific Computing: lecture notes Linear Systems 2nd stage : n 7! n 1 ! (n2)(n)⇥/÷ for 2nd stage (inc. n2 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. Pn1 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÷ xn1 = [gn1 Un1,nxn]/Un1,n1 . 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 A1, solve takes ⇡ 1n3 +n·n2 ⇡ 4n3 operations x = A1b 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