INM2012a.pdf

(182 KB) Pobierz
Introduction to Numerical Methods
Marcin Studniarski
Faculty of Mathematics and Computer Science
University of × z
ód´
1
Rounding of numbers
x
= 0:
1 2
:::
r
1
r
;
We consider a positive decimal number of the following form:
(1)
which uses
r
digits to the right of the decimal point.
De…nition 1
We say that
x
is
rounded
to
s
decimal places (where
s < r)
if the
following operartions are performed, depending on the values of
s
and
s+1
:
(a)
For
s+1
2 f0;
1; 2; 3; 4g,
the digit
s
remains unchanged and the digits
s+1
:::
r
are discarded.
(b)
For
s+1
2 f5;
6; 7; 8; 9g
and
s
6
= 9,
the digit
s
is increased by
1
and
the digits
s+1
:::
r
are discarded.
(c)
For
s+1
2 f5;
6; 7; 8; 9g
and
s
= 9,
the digit
s
is replaced by
0,
the
digit
s
1
is increased by
1
(except for the case where
s
1
= 9)
and the digits
s+1
:::
r
are discarded.
Remark 2
The case where
s+1
= 5
and
s+2
=
:::
=
r
= 0
can be handled
in di¤ erent ways. For example, we may choose to increase
s
by
1
only when
s
is even.
De…nition 3
We say that
x
is
chopped
or
truncated
to
s
decimal places (where
s < r)
if the digits
s+1
:::
r
are discarded (no matter what the value of
s+1
is).
Proposition 4
(a)
If
rd(x)
is obtained by rounding
x
to
s
decimal places, then
jx
rd(x)j
1
:
2 10
s
(2)
(b)
If
ch(x)
is obtained by chopping
x
to
s
decimal places, then
jx
ch(x)j
<
1
:
10
s
(3)
s
Proof.
(a) If
s+1
2 f0;
1; 2; 3; 4g,
then
x
= rd(x) +
",
where
" <
(1=2) 10
hence (2) follows. If
s+1
2 f5;
6; 7; 8; 9g,
then
rd(x) = ch(x) + 10
s
,
and
x
= ch(x) +
10
s
;
) 10
s
(4)
. Since
where
1=2
<
1.
Equalities (4) imply
rd(x)
0
<
1
1=2,
we obtain (2).
(b) We have
x
= ch(x) +
10
s
with
0
10
s
<
10
s
, and (3) holds.
1
x
= (1
<
1.
Then
jx
ch(x)j =
2
Floating-point number system
De…nition 5
(a)
A
‡
oating-point number system
F
is a subset of the set of
real numbers
R
such that each element
x
2
F
has the form
x
=
m
e t
;
(5)
where
m; ; e; t
are integer parameters, and
>
0
is called the
base
(most computers use
t
is called the
precision,
m
2
[0;
t
= 2),
1]
is called the
mantissa
or the
signi…cand,
e
2
[e
min
; e
max
]
is called the
exponent range.
(b)
The system
F
is called
normalized
if
m
t
1
for every
x
6
= 0
(6)
(the number
0
is an exceptional number which does not have a normalized rep-
resentation).
(c)
The elements of
F
are called
‡
oating-point numbers
or
machine numbers.
Proposition 6
The range of nonzero ‡oating-point numbers in a normalized
system
F
is given by
e
min
1
jxj
e
max
(1
t
):
(7)
Proof.
Let
x
2
F.
It follows from (5) and (6) that
jxj
=
m
e t
t
e t
t
1
e t
=
t
e
1
e
min
1
:
On the other hand, using the condition
m
jxj
=
m
(
1)
e
max
t
1,
we get
e
max
t
=(
t
1)
=
e
max
(1
t
):
De…nition 7
The
positional
representation of a number
x
2
F
is given by
x
=
e
d
1
+
d
2
2
+
:::
+
d
t
t
=
e
0:d
1
d
2
:::d
t
;
(8)
where
d
i
2
[0;
1]
are digits, and
d
1
6
= 0
for
x
6
= 0.
Here
d
1
is called the
most
signi…cant digit
and
d
t
is called the
least signi…cant digit.
De…nition 8
Let
G R
be the set of all real numbers of the form
(5)
with no
restriction on the exponent
e.
We denote by
:
R
!
G
any mapping such that
(x)
is an element of
G
nearest to
x
(since there may be two such elements for
some numbers
x,
the mapping may be de…ned in di¤ erent ways). We say that
(a)
(x)
over‡
ows
if
j
(x)j
>
maxfjyj :
y
2
Fg,
(b)
(x)
under‡
ows
if
j
(x)j
<
minfjyj : 0
6
=
y
2
Fg.
2
For any
2
R,
we will us the following notation:
b c
is the largest integer less than or equal to ,
d e
is the smallest integer greater than or equal to .
Theorem 9
[1]
Suppose that
x
2
R
lies in the range of
F
(that is, it satis…es
(7)).
Then there exists
2
R
such that
(x) =
x(1
+ )
and
j j
< u
R
:=
1
2
1
t
:
(9)
Proof.
We give the proof for the case of
x >
0
only. We may represent
x
in
the form
e t
x
=
;
for some
2
[
t
1
;
t
)
and
e
2
[e
min
; e
max
].
The number
x
lies between the
adjacent ‡
oating point numbers
y
1
=
b c
j
(x)
Hence, we can estimate
j
(x)
The conditions imposed on
x
=
xj
jy
2
2
y
1
j
e t
e t
and
y
2
=
d e
e t
:
Then either
(x) =
y
1
or
(x) =
y
2
. Moreover,
(x)
is chosen so that
xj
= minfjy
i
xj
:
i
= 1; 2g:
2
=
:
(10)
give that
e t
t
1
e t
e
1
:
(11)
Taking into account (10) and (11), we obtain
(x)
x
x
1
2
e t
e t
1
2
e t
e
1
=
1
2
1
t
:
(12)
The second inequality in (12) is strict, except for the case of
=
t
1
. But
in this case we have
(x) =
x,
and the …rst inequality is strict. Therefore, by
de…ning
:=
(x)
x
, we get (9).
x
De…nition 10
(a)
The number
u
R
appearing in
(9)
is called the
unit roundo¤
error.
(b)
The distance
"
M
from
1
to the next larger machine number is called the
machine epsilon.
Proposition 11
"
M
= 2u
R
=
1
t
:
Proof.
Using the positional representation (8), we see that
1=
0:1
1
0
2
:::0
t
1
0
t
;
where the lower indices denote the positions after the decimal point. The next
machine number
y
larger than
1
is given by
y
=
Hence, we have
"
M
=
y
1=
0:0
1
0
2
:::0
t
1
1
t
0:1
1
0
2
:::0
t
1
1
t
:
=
1
t
:
3
3
Error analysis
De…nition 12
Let a number
x
be approximated by another number
x
(where
^
x; x
2
R).
Then
^
(a)
x x
is called the
error;
^
(b)
jx
x
j
is called the
absolute error;
^
x x
^
(c)
x
is called the
relative error
(we assume that
x
6
= 0).
In particular, it follows from (12) that the relative error of approximating a
number
x
lying in the range of
F
by the nearest machine number is less than
u
R
.
We will now analyse roundo¤ errors appearing in the process of performing
basic arithmetic operations on a digital computer. Let
x; y
2
R
and let
2
f+;
; ; =g
Since each of these operations is carried out on machine numbers, the
result
x y
is de…ned for
x; y
2
F
only. We assume that
x y
is …rst computed
correctly, then normalized and rounded o¤ to
(x
y).
In real computers,
this assumption is not satis…ed exactly, but it is very close to the truth. In
practice, many computers carry out arithmetic operations in special registers
that have more bits than the usual machine numbers. Therefore, although the
value
x y
is not computed with full accuracy, it is initially determined with
greater precision that the result
(x
y)
later stored as a machine number.
We consider two possible cases:
(a)
x
and
y
are machine numbers (that is,
x; y
2
F)
and
x y
lies in the
range of
F.
Then Theorem 9 gives
(x
y)
= (x
y)(1
+ );
j j
< u
R
:
(13)
is
(b)
x
and
y
are not necessarily machine numbers. Then the operation
carried out on
(x)
and
(y),
and we get
( (x)
(y)) = (x(1 +
1
)
y(1
+
2
))(1
+
3
);
j
i
j
< u
R
,
i
= 1; 2; 3:
(14)
4
The fundamental operator
In many numerical problems (e.g. interpolation, numerical di¤erentiation and
integration, numerical solving of di¤erential equations) one uses approximations
of a given function by polynomials. In order to present a uni…ed approach to
such problems, we now introduce a linear operator which is a basis of all these
polynomial approximation methods.
De…nition 13
[3]The fundamental operator
L
:
C
m
(R)
L(f; x)
:=
k
f
(k)
R
!
R
is de…ned by
(15)
(x) +
m n
XX
i=0 j=1
q
ij
(x)f
(i)
(a
ij
);
where
f
(i)
is the
i-th
order derivative of
f
for
i
1,
f
(0)
=
f
,
a
ij
are constants,
k
2 f0;
1g,
and
q
ij
are polynomials of the variable
x.
By making di¤erent choices of parameters
k; m; n; a
ij
;
k
; q
ij
;
we can con-
struct various numerical methods for solving di¤erent approximation problems.
A general procedure is described below.
4
1. The part of the operator
L
which we want to approximate is formally
replaced by its approximation, and the resulting operator is denoted by
^
L.
2. Then we solve the given approximation problem by solving the equation
^
L(f; x)
= 0:
(16)
This gives a formula for the approximated quantity (e.g. function value,
derivative or integral).
3. It follows from (16) that the error of any approximation de…ned by steps
1– is equal to
2
E(x)
:=
L(f; x)
^
L(f; x)
=
L(f; x):
(17)
Since we are seeking polynomial approximations only, we assume that
the coe¢ cient (number or polynomial) of the approximated quantity in
(15) is constant. We now describe some special forms of the fundamental
operator which are used to solve speci…c problems.
4.1
Interpolation and extrapolation
Let
f
:
R
!
R
bed a function whose values (possibly with certain derivatives)
are known on a …nite set of points called
nodes.
The aim of
interpolation
is to
determine approximate values of
f
at points which are not nodes but lie between
nodes. The aim of
extrapolation
is to approximate
f
at points lying outside the
range of nodes (but su¢ ciently close to the smallest or to the largest node).
Suppose that in (15) we have
k
= 0,
0
= 1
and at least one coe¢ cient
q
0j
(x)
is di¤erent from
0.
Then the operator
L
de…nes an interpolation of
f
for
x
2
[min
a
ij
;
max
a
ij
]
and an extrapolation of
f
for
x
2
[min
a
ij
;
max
a
ij
].
A
=
general form of the interpolation/extrapolation operator is thus the following:
L(f; x)
:=
f
(x) +
m n
XX
i=0 j=1
q
ij
(x)f
(i)
(a
ij
):
(18)
Although there exist some interpolation methods that use the derivatives of
f
(e.g. Hermite interpolation), in this lecture we deal with the case of
m
= 0
only.
In this case, using the traditional symbols, we denote
q
0j
(x)
by
l
j
(x).
We also
write
a
j
instead of
a
0j
.Then (18) can be rewritten as
L(f; x)
=
f
(x)
n
X
j=1
l
j
(x)f (a
j
):
(19)
According to the general procedure, we now replace the function
f
in (19) by
^
its polynomial approximation
p.
The resulting operator
L
should be identically
zero; therefore
n
X
p(x)
=
l
j
(x)f (a
j
)
for all
x
2
R:
(20)
j=1
5

Plik z chomika:

Inne pliki z tego folderu:

Inne foldery tego chomika:

Zgłoś jeśli naruszono regulamin