Docsity
Docsity
Prepare for your exams
Prepare for your exams

Study with the several resources on Docsity

Earn points to download
Earn points to download

Earn points by helping other students or get them with a premium plan

Guidelines and tips
Guidelines and tips
Sell on Docsity
Sell on Docsity
Docsity AI
Docsity AI
Log inSign up

Prepare for your exams

Study with the several resources on Docsity

Find documents

Prepare for your exams with the study notes shared by other students like you on Docsity

Search for your university

Find the specific documents for your university's exams

Docsity AINEW

Summarize your documents, ask them questions, convert them into quizzes and concept maps

Explore questions

Clear up your doubts by reading the answers to questions asked by your fellow students

Earn points to download

Earn points by helping other students or get them with a premium plan

Share documents
download point icon20 Points

For each uploaded document

Answer questions
download point icon5 Points

For each given answer (max 1 per day)

All the ways to get free points
Get points immediately

Choose a premium plan with all the points you need

Study Opportunities

Choose your next study program

Get in touch with the best universities in the world. Search through thousands of universities and official partners

Community

Ask the community

Ask the community for help and clear up your study doubts

Free resources

Our save-the-student-ebooks!

Download our free guides on studying techniques, anxiety management strategies, and thesis advice from Docsity tutors

Scientific Computing, Lecture Notes - Computer Science - 1, Study notes of Computer Numerical Control

University of London (Lond)Computer Numerical Control

Numerical Linear Algebra, Geometrical Interpretation, Forward elimination, Back substitution, Gaussian elimination, LU factorisation, LDU factorisation, Gauss-Seidel Method, Jacobi Method

Typology: Study notes

2010/2011

Uploaded on 09/09/2011

jennyfromtheblock
jennyfromtheblock 🇬🇧

2.3

(3)

225 documents

1 / 8

Toggle sidebar

This page cannot be seen from the preview

Don't miss anything!

bg1
1
CM0368
Scientific Computing
Spring 2009
Professor David Walker
[email protected]
Schedule
• Weeks 1 & 2 (6): Numerical linear algebra (DWW)
– Solving Ax= b by Gaussian elimination and Gauss-Seidel.
– Algebraic eigenvalue problem. Power method and QR method.
• Weeks 3 & 4 (6): Numerical solution of differential
equations (BMB)
– Ordinary differential equations (finite difference and Runge-Kutta
methods).
– Partial differential equations (finite difference method)
• Weeks 5 & 6 (5): Applications in Physics (BMB)
– Laplace and Helmholtz equations, Schrodinger problem for the
hydrogen atom, etc.
• Weeks 6 & 7 (4): Numerical optimisation (DWW)
• Weeks 2 – 11: Tutorials on Wednesdays at 1pm in
T2.07. These will be given by Kieran Robert.
Text Book
• “Numerical Computing with MATLAB” by
Cleve B. Moler, SIAM Press, 2004. ISBN
0898715601.
•http://www.readinglists.co.uk/rsl/student/sv
iewlist.dfp?id=20248
• Web edition at
http://www.mathworks.com/moler/
Web Site
• Lecture notes and other module material
can be found at:
http://users.cs.cf.ac.uk/David.W.Walker/CM0368
Numerical Linear Algebra
• A system of n simultaneous linear equations can
be represented in matrix notation as:
Ax= b
where A is an nμn matrix, and xand bare
vectors of length n.
• Can write solution as x= A-1b where A-1 is the
inverse of A.
• A square matrix is said to be non-singular if its
inverse exists. If the inverse does not exist then
the matrix is singular.
Geometrical Interpretation
•If A is a 2μ2 matrix, for example
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
−
=
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛−
1
3
43
12
2
1
x
x
•If A is a 3μ3 matrix each of the three equations
represents a plane, and the solution is the point lying at
the intersection of the three planes.
then 2x1–x
2= 3 and 3x1+ 4x2= -1. Each
represents a straight line and the solution of the
above is given by their intersection.
pf3
pf4
pf5
pf8

Partial preview of the text

Download Scientific Computing, Lecture Notes - Computer Science - 1 and more Study notes Computer Numerical Control in PDF only on Docsity!

CM

Scientific Computing

Spring 2009

Professor David Walker

[email protected]

Schedule

  • Weeks 1 & 2 (6): Numerical linear algebra (DWW)
    • Solving A x = b by Gaussian elimination and Gauss-Seidel.
    • Algebraic eigenvalue problem. Power method and QR method.
  • Weeks 3 & 4 (6): Numerical solution of differentialequations (BMB)
    • Ordinary differential equations (finite difference and Runge-Kuttamethods).
    • Partial differential equations (finite difference method)
  • Weeks 5 & 6 (5): Applications in Physics (BMB)
    • Laplace and Helmholtz equations, Schrodinger problem for thehydrogen atom, etc.
  • Weeks 6 & 7 (4): Numerical optimisation (DWW)
  • Weeks 2 – 11: Tutorials on Wednesdays at 1pm in

T2.07. These will be given by Kieran Robert.

Text Book

• “Numerical Computing with MATLAB” by

Cleve B. Moler, SIAM Press, 2004. ISBN

• http://www.readinglists.co.uk/rsl/student/sv

iewlist.dfp?id=

• Web edition at

http://www.mathworks.com/moler/

Web Site

• Lecture notes and other module material

can be found at:

http://users.cs.cf.ac.uk/David.W.Walker/CM

Numerical Linear Algebra

• A system of n simultaneous linear equations can

be represented in matrix notation as:

A x = b

where A is an nμn matrix, and x and b are

vectors of length n.

• Can write solution as x = A-1 b where A-1^ is the

inverse of A.

• A square matrix is said to be non-singular if its

inverse exists. If the inverse does not exist then

the matrix is singular.

Geometrical Interpretation

  • If A is a 2μ2 matrix, for example

x

x

  • If A is a 3μ3 matrix each of the three equations

represents a plane, and the solution is the point lying at

the intersection of the three planes.

then 2x 1 – x 2 = 3 and 3x 1 + 4x 2 = -1. Each

represents a straight line and the solution of the

above is given by their intersection.

MATLAB Solution

• In MATLAB we can find the solution to A x = b by

writing:

x = A\b

• Write: A = [10 -7 0;-3 2 6;5 -1 5]

• Then: b = [7;4;6]

• Then: x = A\b

Gaussian Elimination

• Eliminate x 1 from all the equations after the first.

• Then eliminate x 2 from all the equations after the

second.

• Then eliminate x 3 from all the equations after the

third.

• And so on, until after n-1 steps we have

eliminated xj from all the equations after the jth,

for j = 1, 2, …, n-1.

• These steps are referred to as the forward

elimination stage of Gaussian elimination.

Example

⎟⎟⎟^ =

3

2

1 x

x

x

Subtract -3/10 times equation 1 from equation 2, and 5/10 timesequation 1 from equation 3.

3

2

1 x

x

x

Next we swap equations 2 and 3. This is calledto get the largest absolute value on or below the diagonal in column 2 onto partial pivoting. It is done the diagonal. This makes the algorithm more stable with respect to round-off errors (see later).

Example (continued)

3

2

1 x

x

x

Now subtract -0.1/2.5 times equation 2 from equation 3.

⎟⎟⎟^ =

3

2

1 x

x

x

This completes the forward elimination stage of the

Gaussian elimination algorithm.

Pseudocode for Forward

Elimination

make b column n+1 of A

for k=1 to n-

find pivot A(m,k) in column k and row m¥k

swap rows k and m

for i=k+1 to n

mult = A(i,k)/A(k,k)

for j=i+1 to n+

A(i,j) = A(i,j) – mult*A(k,j)

end

end

end

Back Substitution

• After the forward elimination phase, the matrix

has been transformed into upper triangular form.

• Equation n just involves xn and so can now be

solved immediately.

• Equation n-1 just involves xn-1 and xn , and since

we already know xn we can find x n-.

• Working our way backwards through the

equations we can find xn , xn-1 ,…, x 1.

• This is called the back substitution phase of the

Gaussian elimination algorithm.

Pseudocode for LU Factorisation

make b column n+1 of A

for k=1 to n-

find pivot A(m,k) in column k and row m¥k

swap rows k and m

for i=k+1 to n

A(i,k) = A(i,k)/A(k,k)

for j=i+1 to n+

A(i,j) = A(i,j) – A(i,k)*A(k,j)

end

end

end

Explanation of LU

• At stage i of forward elimination we do

pivoting to find the largest absolute value

in column i on or below the diagonal, and

then exchange rows to bring it onto the

diagonal.

• Then we subtract multiples of row i from

rows i+1,…,n.

• Each of these operations can be

represented by a matrix multiplication.

  • An elementary matrix, M, is one that has 1’s along the

main diagonal and 0’s everywhere else, except for onenon-zero value (-m, say) in row i and column j.

  • Multiplying A by M has the effect of subtracting m timesrow j of matrix A from row i.
  • Ignoring pivoting, the GE algorithm applies a series of

elementary matrices to A to get U:

U = Mn-1….M 2 M 1 A

  • If Li-1=Mi then

L 1 L 2 …L n-1 U = A

so taking L = L 1 L 2 …Ln-1 we have LU=A.

  • Li is the same as Mi except the sign of the non-zero value

is changed.

  • With pivoting U = M (^) n-1 P (^) n-1 ….M 2 P 2 M 1 P 1 A, and it can be

shown in a similar way that LU=PA, where P = Pn-1 …P 2 P 1.

Elementary Matrices The Need for Pivoting

• Suppose we change our problem slightly to:

⎟⎟⎟^ =

3

2

1 x

x

x

where the solution is still (0,-1,1) T.

• Now suppose we solve the problem on a

computer that does arithmetic to five decimal

places.

Pivoting (continued)

• The first step of the

elimination gives: ⎟⎟

⎟ ⎠

⎞ ⎜⎜

⎜ ⎝ =⎛ ⎟⎟

⎟ ⎠

⎞ ⎜⎜

⎜ ⎝

⎛ ⎟⎟

⎟ ⎠

⎞ ⎜⎜

⎜ ⎝

⎛ −−

  1. (^001). 5 7 00 02 .. 500156

10 7 0 32

1 xx

x

• The (2,2) element is quite small. Now we

continue without pivoting.

• The next step is to subtract -2.5μ 103 times

equation 2 from equation 3:

(5-(-2.5 μ 103 )(6))x 3 = (2.5-(-2.5 μ 103 )(6.001))

• The righthand side is (2.5+1.50025 μ 104 ). The

second term is rounded to 1.5002 μ 104. When

we add the 2.5 the result is rounded to 1.

Pivoting (continued)

  • So the equation for x1.5005 μ 10 4 x 3 becomes: which gives x 3 = 0.99993 (instead of 1).^3 = 1.5004^ μ^104
  • Then x-0.001 x 2 is found from: which gives:^2 + (6)(0.99993) = 6. so x -0.001 x^2 = 1.5^ μ^10 -
  • Finally, x^2 = -1.5 (instead of -1). 1 is found from: which gives x10 x^1 + (-7)(-1.5) = 7
  • The problem arises from using a small pivot which leads to a larger^1 = -0.35 (instead of 0).
  • multiplier.Partial pivoting ensures that the multipliers are always less than or equal to 1 in magnitude, and results in a satisfactory solution.

Measuring Errors

• The discrepancy due to rounding error in a

solution can be measured by the error :

e = x – x *

and by the residual :

r = b - A x *

where x is the exact solution and x * is the

computed solution.

• e = 0 if and only if r = 0 , however, e and r are

not necessarily both small.

Error Example

  • Consider an example inwhich the matrix is almost

singular.

  • If GE is used to computethe solution with low

precision we get a largeerror but small residual.

  • Geometrically, the linesrepresented by the

equation are almostparallel.

  • GE with partial pivotingwill always produce small

residuals, but notnecessarily small errors.

⎟⎟⎠^ =⎜⎜⎝⎛ ⎟⎟⎠⎞

2

1

x

x

computed solution

exact solution

Sensitivity

• We want to know how sensitive the solution to

A x = b is to perturbations in A and b.

• To do this we first have to come up with some

measure of how close to singular a matrix is.

• If A is singular, a solution will not exist for some

b ’s (inconsistency), while for other b ’s it is not

unique.

• So if A is nearly singular we would expect small

changes in A and b to cause large changes in x.

• If A is the identity matrix then x = b , so if A is close

to the identity matrix we expect small changes in A

and b to result in small changes in x.

Singularity and Pivots

• In GE all the pivots are non-zero if and

only if A is non-singular, provided exact

arithmetic is used.

• So if the pivots are small we expect the

matrix to be close to singular.

• However, with finite precision arithmetic, a

matrix might still be close to singular even

if none of the pivots are small.

Norms

  • Size of pivots is not adequate to decide how close amatrix is to being singular.
  • Define l (^) p family of norms (1≤p≤¶):

x 1 = ⎜⎝⎛∑ i =^ n 1 xi ⎟⎠⎞

  • Most common norms use p = 1, 2, and ¶.

n^ p i

xp xip^1 /

=

1 / 2 1

=

n

x i xi x^ ∞=max i xi

Properties of Norms

• All these norms have the following properties:

x > 0 ifx ≠ 0

cx = c x

x + y ≤ x + y

for all scalars c

(the triangle inequality)

• In MATLAB use norm(x,p) to find a norm:

norm1 = norm(x,1) norm2 = norm(x)

norminf = norm(x,inf)

Matrix Norms

• The norm of matrix A is ||A|| = M, i.e.,

x

A =max Ax

• Since ||A-1^ || = 1/m it follows that the

condition number can also be defined as:

k(A) = ||A|| ||A-1^ ||

• In MATLAB, cond(A,p) computes the condition

number of A relative to the lp -norm.

Iterative methods for A x = b

• Gaussian Elimination is a direct method

that computes the solution of A x = b is

O(n 3 /3) operations.

• If n is large we might want a faster, less

accurate, method.

• With an iterative method we can stop the

iteration whenever we think we have a

sufficiently accurate solution.

Iterative methods

• Suppose we split the matrix as A = S-T, then S x

= T x + b.

• We can turn this into an iteration:

S x k+1 = T x k + b

or

x k+1 = S-1^ T x k + S-1 b

• So if this sequence converges then we can start

the iteration with a guess at the solution x 0 and

get an approximate solution.

• We need S to be easily invertible

Common Iterative Methods

1. S = diagonal part of A (Jacobi’s method)

2. S = triangular part of A (Gauss-Seidel method)

3. S = combination of 1 and 2 (successive over-

relaxation or SOR)

S is called a pre-conditioner. The choice of

S affects the convergence properties of the

solution.

Jacobi Method

  • S = diag(A) so formula for iteration k+1 becomes: i

n k ji ij jk

i ii ik j ij j

a x = −∑ a x −∑ a x + b

=+

−

  • = 1

1 (^1 )

  • Expressed in matrices:

D x k+1 = - (L+U) x k + b

where D, L, and U are the diagonal, strictly lowertriangular, and strictly upper triangular parts of A,

respectively. LDU factorisation. Note: these are not the D, L, and U of the

Jacobi Example

  • Example: ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎟⎟^ =⎛ ⎠ ⎜⎜ ⎞ ⎝ ⎟⎟^ =⎛ ⎠ ⎜⎜ ⎞ ⎝ ⎟⎟^ =⎛ ⎠ ⎜⎜ ⎞ ⎝ = ⎛ − − − 12 02 ,^01 10 ,^01 02 ,^20 1 2 A^2 1 S T S^1 T

Then:

⎟⎟⎜⎜⎝⎛ ⎟⎟⎠⎞^ +

2

1 (^21)

(^1) bb x

x x

x k k

Pseudocode for Jacobi Method

choose an initial guess x 0

for k=0,1,2,….

for i=1 to n

sum = 0.

for j=1 to i-1 and i+1 to n

sum = sum + A(i,j)*xk (j)

end

xk+1(i) = (b(i)–sum)/A(i,i)

end

check convergence and continue if needed

end

Gauss-Seidel Method

• A problem with the Jacobi method is that we

have to store all of x k until we have finished

computing x k+.

• In the Gauss-Seidel method the x k+1 are used

as soon as they are computed, and replace the

corresponding x k on the righthand side

i

n k ji ij jk

i ii ik j ij j

a x = −∑ a x −∑ a x + b

  • =+

−

  • = (^11)

1 (^1 )

• This uses half the storage of the Jacobi method.

Gauss-Seidel in Matrix Notation

• Expressed in matrices:

(D+L) x k+1 = -U x k + b

where D+L is the lower triangular part of A,

and U is the strictly upper triangular part of

A.

Gauss-Seidel Example

• Example:

A = (^) ⎜⎜⎝⎛ −^21 − 21 ⎟⎟⎠⎞, S =⎜⎜⎝⎛ −^2120 ⎟⎟⎠⎞, T =⎜⎜⎝⎛ 0001 ⎟⎟⎠⎞, S −^1 T =⎜⎜⎝⎛ 0011 // 42 ⎟⎟⎠⎞

Then:

⎜⎜⎝⎛ −^2 120 ⎟⎟⎠⎞⎜⎜⎝⎛ xx 21 ⎟⎟⎠⎞ k + 1 =⎜⎜⎝⎛ 00 01 ⎟⎟⎠⎞⎜⎜⎝⎛ xx 21 ⎟⎟⎠⎞ k^ +⎜⎜⎝⎛ bb 21 ⎟⎟⎠⎞

• Gauss-Seidel is better than Jacobi because ituses half the storage and often converges

faster.

Pseudocode for Gauss-Seidel

choose an initial guess xfor k=0,1,2,…. 0 for i=1 to nsum = 0. for j=1 to i-1sum = sum + A(i,j)x end k+1(j) for j=i+1 to nsum = sum + A(i,j)x end k(j) endx^ k+1^ (i) = (b(i)–sum)/A(i,i) endcheck convergence and continue if needed