Numerical Methods For Computational Science and Engineering: Lecture 1, Sept 19, 2013: Introduction

Numerical Methods for Computational Science and Engineering
Numerical Methods for Computational Science

and Engineering
Lecture 1, Sept 19, 2013: Introduction
Peter Arbenz
Computer Science Department, ETH Zürich
E-mail: [email protected]
NumCSE, Lecture 1, Sept 19, 2013 1/40

Introduction
Outline of today’s lecture

I What is numerical methods for CSE
I Survey of the lecture
I Organization of the lecture (exercises/examination)
I References
I Start of the lecture

Introduction
Scientific Computing

Introduction
Survey on lecture
1. Introduction
2. Roundoff errors
3. Nonlinear equations in one variable (2 lectures)
4. Linear algebra review
5. Direct methods for linear system (2)
6. Linear least squares problems (2)
7. Iterative methods for linear system (2)
8. Eigenvalues and singular values (2)
9. Nonlinear systems and optimization (3)
10. (Piecewise) polynomial interpolation (3)
11. Best approximation

Introduction
Survey on lecture (cont.)

12. Filtering algorithms, Fourier transform
13. Numerical differentiation
14. Numerical integration (2)
15. Ordinary differential equations, initial value problems (3)

Introduction
About this course

Focus
I on algorithms (principles, scope, and limitations),
I on (efficient, stable) implementations in Matlab,
I on numerical experiments (design and interpretation).
No emphasis on
I theory and proofs (unless essential for understanding of
algorithms)
I hardware-related issues (e.g. parallelization, vectorization,
memory access)
(These aspects will be covered in the course
“High Performance Computing for Science and Engineering”
offered by D-INFK)
Introduction
Goals
• Knowledge of the fundamental algorithms in numerical
mathematics
• Knowledge of the essential terms in numerical mathematics
and the techniques used for the analysis of numerical
algorithms
• Ability to choose the appropriate numerical method for
concrete problems
• Ability to interpret numerical results
• Ability to implement numerical algorithms efficiently in
Matlab
Indispensable: Learning by doing (Ô exercises)

Introduction
Literature
Uri Ascher & Chen Greif: A First Course in Numerical
Methods. SIAM, 2011.
http://www.siam.org/books/cs07/
I Excellent reference.
I Main reference for large parts of this
course.
I Target audience: undergraduate students
in computer science.
I I will follow this book quite closely.

Introduction
Literature (cont.)
W. Dahmen & A. Reusken: Numerik für Ingenieure und
Naturwissenschaftler, Springer, 2006.
Good reference for large parts of this course; a lot of simple examples and lucid explanations, but also
rigorous mathematical treatment; Target audience: undergraduate students in science and engineering.
Available through Nebis.
H.-R. Schwarz & N. Köckler: Numerische Mathematik.

Teubner, 2006. 6. Auflage.
Main reference for large parts of this course; Target audience: undergraduate students in science and
engineering. Available through Nebis.
C. Moler: Numerical Computing with Matlab. SIAM 2004.

Good reference for some parts of this course; Target audience: Matlab users and programmers. See
http://www.mathworks.ch/moler/.

Introduction
Prerequisites
Essential prerequisite for this course is a solid knowledge in linear
algebra and calculus. Familiarity with the topics covered in the
first semester courses is taken for granted, see
• K. Nipp and D. Stoffer, Lineare Algebra, vdf
Hochschulverlag, Zürich, 5 ed., 2002.
• M. Gutknecht, Lineare algebra, lecture notes, SAM, ETH
Zürich, 2009.
http://www.sam.math.ethz.ch/~mhg/unt/LA/HS07/.
• M. Struwe, Analysis für Informatiker. Lecture notes, ETH
Zürich, 2009.

Organization
Organization
Lecturer:
Prof. Peter Arbenz [email protected]
Assistants:
Stefan Pauli [email protected]
Daniel Hupp [email protected]
Christian Schüller [email protected]
Robert Carnecky [email protected]
Laura Scarabosio [email protected]
Jonas Šukys [email protected]
Cecilia Pagliantini [email protected]
Alexander Lobbe [email protected]
Sebastian Grandis [email protected]
Sharan Jagathrakashakan [email protected]

Organization
Venue
Classes: Mon 10.15-12.00 (CAB G11); Thu 10.15-12.00 (HG G5)
Tutorials: Mon 13.15-15.00
Thu 8.15-10.00
Please register (on course website) for tutorial groups until
September 23th:
http://www.math.ethz.ch/education/bachelor/lectures/
hs2013/math/nummath_cse
Consulting hours: if needed, see the course website.

Organization
Assignments
• The assignment sheets will be uploaded on the course
webpage on Monday every week the latest.
• The exercise should be solved until the following tutorial class.
(Hand them in to the assistant or grade yourself.)

Organization
Examination
• Three-hour written examination involving coding problems to
be done at the computer on
TBA
• Dry-run for computer based examination:
Does not exist anymore.
Try out a computer in the student labs in HG.
• Pre-exam question session:
TBA

Organization
Examination (cont.)
• Topics of examination:
I All topics, that have been addressed in class or in a
homework assignment.
I One exam question will be one of the homework
assignment.
• Lecture slides will be available as (a single) PDF file during
the examination.
• The Ascher-Greif book will be made available, too.
• The exam questions will be asked both in German and in
English. You can chose among the two.

Computing environment: Matlab
Problem solving environment: Matlab

I We use Matlab for the exercises.
I Although most of the algorithm we are dealing with have been
implemented in Matlab, it is useful when you program them
yourselves.
I These (little) programs will be building blocks when you will
solve more complex problems in your future.
I Matlab help
I Matlab commands help/doc
I Matlab online documentation, e.g.,
http://www.mathworks.nl/help/pdf_doc/allpdf.html
I Numerous introductory textbooks / user guides / primers
I My own Matlab introduction:
http://people.inf.ethz.ch/arbenz/MatlabKurs/
matlabintro.pdf
Numerical algorithms and errors

I The most fundamental feature of numerical computing is the
inevitable presence of errors.
I The result of any interesting computation (and of many
uninteresting ones) is typically only approximate, and our goal
is to ensure that the resulting error is tolerably small.
I Example: How many loop iterations are there in this little
Matlab program?
x = 0; h = 1/10;
while x<1,
x=x+h;
% do something depending on x
end

How to measure errors

I Can measure errors as absolute or relative, or a combination
of both.
I The absolute error in v approximating given scalar quantity u
is |u − v |.
|u − v |
I The relative error (assuming u 6= 0) is .
|u|
u v Absolute Relative
Error Error
1 0.99 0.01 0.01
1 1.01 0.01 0.01
–1.5 –1.2 0.3 0.2
100 99.99 0.01 0.0001
100 99 1 0.01
Approximation example
The Stirling approximation
√ n n
v = Sn = 2πn ·
e
is used to approximate u = n! = 1 · 2 · · · n for large n.
The formula involves the constant e = exp(1) = 2.7182818 . . ..
We use the Matlab program Example1 1.m from
http://www.siam.org/books/cs07/programs.zip
to compute the values of u and v for n = 1, . . . , 10.
Notice that the quality of v approximating u increases although
the absolute error increases rapidly!

Types of errors
I Errors in the formulation of the problem to be solved.
I Errors in the mathematical model. Simplifications.
I Error in input data. Measurements.
I Approximation errors
I Discretization error.
I Convergence error in iterative methods.
I Discretization/convergence errors may be assessed by an
analysis of the method used.
I Roundoff errors
I Roundoff errors arise everywhere in numerical computation
because of the finite precision arithmetic.
I Roundoff errors behave quite erratic.

Discretization errors in action

Problem: want to approximate the derivative f 0 (x0 ) of a given
smooth function f (x) at the point x = x0 .
Example: Let f (x) = sin(x), −∞ < x < ∞, and set x0 = 1.2.
Thus, f (x0 ) = sin(1.2) ≈ 0.932 . . ..
Discretization: Function values f (x) are available only at a discrete
number of points, e.g. at grid points xj = x0 + jh, j ∈ Z.
Want to approximate f 0 (x0 ) by values f (xj ).

Discretization errors in action (cont.)

Taylor’s series gives us an algorithm to approximate f 0 (x0 ):
f (x0 + h) − f (x0 )
f 0 (x0 ) ≈ Dx0 ,h (f ) =
h

Discretization errors in action (cont.)

Expanding f (x) by a Taylor series around x = x0 gives
f (x0 + h) − f (x0 ) h
= f 0 (x0 ) − f 00 (ξ), x0 < ξ < x0 + h.
h 2
So, we expect the error to decrease linearly with h.

0
f (x0 ) − f (x0 + h) − f (x0 ) h 00 h 00
= f (ξ) ≈ f (x0 )
h 2 2
Or, using the big-O notation:
0
f (x0 ) − Dx ,h (f ) = O(h).
0

Results
Try for f (x) = sin(x) at x0 = 1.2.
(So, we are approximating cos(1.2) = 0.362357754476674 . . .)
h Absolute error
0.1 4.71667 · 10−2
0.01 4.666196 · 10−3
0.001 4.660799 · 10−4
10−4 4.660256 · 10−5
10−7 4.619326 · 10−8
These results reflect the discretization error as expected.

Note that f 00 (x0 )/2 = − sin(1.2)/2 ≈ −0.466.

Results for smaller h

The above results indicate that we can compute the derivative as
accurate as we like, provided that we take h small enough.
If we wanted

cos(1.2) − sin(1.2 + h) − sin(1.2) < 10−10 .

h
We have to set h ≤ 10−10 /0.466.
The following numbers and plot are generated by

.../Greif/programs/chap01/Example1 3Figure1 3.m

Results for smaller h
h Absolute error
10−8 4.36105 · 10−10
10−9 5.594726 · 10−8
10−10 1.669696 · 10−7
10−11 7.938531 · 10−6
10−13 6.851746 · 10−4
10−15 8.173146 · 10−2
10−16 3.623578 · 10−1
These results reflect both discretization and roundoff errors.

Results for all h
The solid curve interpolates the computed values of

|f 0 (x0 ) − f (x0 +h)−f
h
(x0 )
| for f (x) = sin(x) at x0 = 1.2.
The dash-dotted straight line depicts the discretization error
without roundoff error.
Algorithm properties
Performance features that may be expected from a good numerical
algorithm.
I Accuracy
Relates to errors. How accurate is the result going to be when
a numerical algorithm is run with some particular input data.
I Efficiency
I How fast can we solve a certain problem?
Rate of convergence. Floating point operations (flops).
I How much memory space do we need?
I These issues may affect each other.
I Robustness
(Numerical) software should run under all circumstances.
Should yield correct results to within an acceptable error or
should fail gracefully if not successful.
Complexity I
Complexity/computational cost of an algorithm :⇐⇒ number of
elementary operators
Asymptotic complexity = ˆ “leading order term” of complexity w.r.t.
large problem size parameters
The usual choice of problem size parameters in numerical linear
algebra is the number of independent real variables needed to
describe the input data (vector length, matrix sizes).
operation description #mul/div #add/sub
inner product (x ∈ Rn , y ∈ Rn ) 7→ xH y n n−1 O(n)
outer product (x ∈ Rm , y ∈ Rn ) 7→ xyH nm 0 O(mn)
tensor product
matrix product (A ∈ Rm×n , B ∈ Rn×k ) 7→ AB mnk mk(n − 1) O(mnk)

Big-O and Θ notation

For an error e depending on h we denote
e = O(hq )
if there are two positive constants q and C such that
|e| ≤ C hq ∀h > 0 small enough.
Similarly, for w = w (n) the expression
w = O(n log n)
means that there is a constant C > 0 such that
|w | ≤ Cn log n as n → ∞.

Big-O and Θ notation

More abstract:
Class O(f ) of functions is defined as
O(f ) = {g | ∃c1 , c2 > 0 : ∀N ∈ Z+ : g (N) ≤ c1 f (N) + c2 }
The Θ notation signifies a stronger relation than the O notation:

a function φ(h) for small h (resp., φ(n) for large n) is Θ(ψ(h))
(resp., Θ(ψ(n))) if φ is asymptotically bounded both above and
below by ψ.
Example:
O(h2 ) means at least “quadratic convergence” (see later). Θ(h2 )
is exact quadratic convergence.

Complexity II
To a certain extent, the asymptotic complexity allows to predict
the dependence of the runtime of a particular implementation of
an algorithm on the problem size (for large problems). For
instance, an algorithm with asymptotic complexity O(n2 ) is likely
to take 4× as much time when the problem size is doubled.
One may argue that the memory accesses are more decisive for run
times than floating point operations. In general there is a linear
dependence among the two. So, there is no difference in the O
notation.

Elementary operations
Scaling
Scaling ≡ multiplication with diagonal matrices (with non-zero
diagonal entries) from left and/or right.
It is important to know the different effects of multiplying with a
diagonal matrix from left or right:
DA vs. AD with Rn×n 3 A, D = diag(d1 , . . . , dn )
Scaling with D = diag(d1 , . . . , dn )

in Matlab:
y = diag(d)*x;
or
y = d.*x;
Elementary matrices
Matrices of the form A = I + αuvT are called elementary.
Again we can apply A to a vector x in a straightforward and a
more clever way:
Ax = (I + αuvT )x
or
Ax = x + αu(vT x)
Cf. exercises.

Problem conditioning and algorithm stability

Qualitatively speaking:
I The problem is ill-conditioned if a small perturbation in the
data may produce a large difference in the result.
The problem is well-conditioned otherwise.
I The algorithm is stable if its output is the exact result of a
slightly perturbed input.

An unstable algorithm
Ill-conditioned problem of computing output values y from input

values x by y = g (x): when x is slightly perturbed to x̄, the result
ȳ = g (x̄) is far from y .
A stable algorithm
An instance of a stable algorithm for computing y = g (x): the

output ȳ is the exact result, ȳ = g (x̄), for a slightly perturbed
input, i.e., x̄ which is close to the input x. Thus, if the algorithm is
stable and the problem is well-conditioned, then the computed
result ȳ is close to the exact y .
Unstable algorithm
Problem statement: evaluate the integrals
Z 1
xn
yn = dx, for n = 0, 1, 2, . . . , 30.
0 x + 10
Algorithm development: observe that analytically, for n > 0,

Z 1 n Z 1
x + 10x n−1 1
yn + 10yn−1 = dx = x n−1 dx = .
0 x + 10 0 n
Also, Z 1
1
y0 = dx = log(11) − log(10).
0 x + 10
Algorithm:
I Evaluate y0 = log(11) − log(10).
1
I For n = 1, 2, . . . , 30, evaluate yn =
n − 10yn−1 .
Unstable algorithm (cont.)

Run the Matlab program Example1 6.m by Ascher and Greif to
see the catastrophic amplification of roundoff errors.
This code is available from
http://www.siam.org/books/cs07/programs.zip.

Unstable algorithm (cont.)

Roundoff error accumulation
I In general, if En is error after n elementary operations, cannot
avoid linear roundoff error accumulation
En ' c0 nE0 .
I Will not tolerate an exponential error growth such as
En ' c1n E0 , for some constant c1 > 1.
This is an unstable algorithm.

Numerical Methods For Computational Science and Engineering: Lecture 1, Sept 19, 2013: Introduction

Uploaded by

Copyright:

Available Formats

Numerical Methods For Computational Science and Engineering: Lecture 1, Sept 19, 2013: Introduction

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Methods For Computational Science and Engineering: Lecture 1, Sept 19, 2013: Introduction

Uploaded by

Copyright:

Available Formats

Numerical Methods for Computational Science and Engineering

Numerical Methods for Computational Science

NumCSE, Lecture 1, Sept 19, 2013 1/40

Outline of today’s lecture

NumCSE, Lecture 1, Sept 19, 2013 2/40

NumCSE, Lecture 1, Sept 19, 2013 3/40

NumCSE, Lecture 1, Sept 19, 2013 4/40

Survey on lecture (cont.)

NumCSE, Lecture 1, Sept 19, 2013 5/40

About this course

NumCSE, Lecture 1, Sept 19, 2013 7/40

NumCSE, Lecture 1, Sept 19, 2013 8/40

Available through Nebis.

H.-R. Schwarz & N. Köckler: Numerische Mathematik.

engineering. Available through Nebis.

C. Moler: Numerical Computing with Matlab. SIAM 2004.

NumCSE, Lecture 1, Sept 19, 2013 9/40

NumCSE, Lecture 1, Sept 19, 2013 10/40

NumCSE, Lecture 1, Sept 19, 2013 11/40

NumCSE, Lecture 1, Sept 19, 2013 12/40

NumCSE, Lecture 1, Sept 19, 2013 13/40

NumCSE, Lecture 1, Sept 19, 2013 14/40

NumCSE, Lecture 1, Sept 19, 2013 15/40

Problem solving environment: Matlab

Numerical algorithms and errors

NumCSE, Lecture 1, Sept 19, 2013 17/40

How to measure errors

NumCSE, Lecture 1, Sept 19, 2013 19/40

NumCSE, Lecture 1, Sept 19, 2013 20/40

Discretization errors in action

NumCSE, Lecture 1, Sept 19, 2013 21/40

Discretization errors in action (cont.)

NumCSE, Lecture 1, Sept 19, 2013 22/40

Discretization errors in action (cont.)

NumCSE, Lecture 1, Sept 19, 2013 23/40

These results reflect the discretization error as expected.

NumCSE, Lecture 1, Sept 19, 2013 24/40

Results for smaller h

We have to set h ≤ 10−10 /0.466.

The following numbers and plot are generated by

NumCSE, Lecture 1, Sept 19, 2013 25/40

Results for smaller h

These results reflect both discretization and roundoff errors.

NumCSE, Lecture 1, Sept 19, 2013 26/40

Results for all h

The solid curve interpolates the computed values of

NumCSE, Lecture 1, Sept 19, 2013 29/40

Big-O and Θ notation

if there are two positive constants q and C such that

|e| ≤ C hq ∀h > 0 small enough.

Similarly, for w = w (n) the expression

means that there is a constant C > 0 such that

NumCSE, Lecture 1, Sept 19, 2013 30/40

Big-O and Θ notation

O(f ) = {g | ∃c1 , c2 > 0 : ∀N ∈ Z+ : g (N) ≤ c1 f (N) + c2 }

The Θ notation signifies a stronger relation than the O notation:

NumCSE, Lecture 1, Sept 19, 2013 31/40

NumCSE, Lecture 1, Sept 19, 2013 32/40

Scaling with D = diag(d1 , . . . , dn )