Riassunto esame Elementi di calcolo numerico, Prof. Morini Benedetta, libro consigliato Elementi di calcolo numerico:metodi ed algoritmi, R.Morandi

Introduction to numerical analysis

Numerical analysis: find a numerical solution to a mathematical problem. This process divides into five steps:

• Real problem
• Model
• Mathematical model
• Solve the mathematical model
• Interpret the solution

Sometimes the exact solution might not be available, and we need an approximation by using a numerical method. Other times we could have a system of equations or an equation of a high grade, both of which require high computation to solve.

Errors in numerical methods

In the formulation of a real problem into a numerical method, we introduce an error, which can have different sources/natures:

• Modelling errors
• Physical measurement errors
• Machine representation and arithmetic errors (computers are finite machines and can't represent infinite numbers)
• Mathematical approximation errors

The point is to guarantee that the error is under control and the different errors are of the same order. For example, if one error is 10^-16 and another is 10^-4, the effort required for the first one is wasted because the second one is huge compared to it.

Notation for errors

Absolute error on a scalar: \( \epsilon = | x - x_{\text{approx}} | \)

This does not take into account the order of magnitude of the value under examination. For example, compare 1.000 - 1.001 = 0.001 and 2 × 10^-3 - 1 × 10^-3 = 0.001; in the first case, the error is negligible, while in the second one it is not.

True relative error on a scalar: \( \epsilon = \left| \frac{x - x_{\text{approx}}}{x} \right| \)

This, compared to the absolute error, takes into account the order of magnitude of the value under examination.

Mathematical approximation errors

In most cases, the true value is not known, so the error is not computable. We approximate the error as: \( \epsilon = \left| \frac{x_{\text{approx}} - x}{x} \right| \)

If we apply a numerical method, it will produce better and better approximations:

And the iterative numerical method will stop once the approximation goes under a certain threshold: \( \epsilon < \text{tolerance} \).

Computer number representation

The fundamental unit is a word which consists of bits. How many digits are available defines a representation system. Only a limited range of integers can be represented. For example, in a 16-bit system, the range of integers that can be stored is -32,768 to 32,767.

In general, in an n-bit system, the range is \([-2^{n-1}, 2^{n-1} - 1]\).

Integers are stored by specifying the word size (how many bits are available). The first bit indicates the sign:

• 0 = positive
• 1 = negative

The remaining bits are used to store the absolute value of the number. For example, \( (10101101)_2 = 2^7 + 2^5 + 2^3 + 2^2 + 2^0 = 128 + 32 + 8 + 4 + 1 = 173_{10} \).

Floating point representation

\(\pm m \cdot b^e\) where:

• \(m\) is the mantissa or significant. It may be positive or negative.
• \(b\) is the base (2, 8, 10, 16).
• \(e\) is the integer exponent of the base used.

Note: This representation is normalized with only one digit at the left of the dot: \(0.00345 \cdot 10^{-3} = 3.45 \cdot 10^{-6}\).

The floating-point representation in computers is usually used with a binary base. Each digit (bit) is either 0 or 1. Since the representation is normalized, the first bit in \(m\) is always 1: \(\pm (1 + f) \cdot 2^e\), where \(f\) is the fractional part of the mantissa and the first bit of \(f\) is always 0.

For example, normalize the binary number \( (1101.1)_2 \cdot 2^{-5} = 1.1011 \cdot 2^{-2} \). We only have to store \(f\) and \(e\).

Most computers adopt IEEE double-precision; eight bytes (64 bits) are used to represent floating-point numbers:

• 1 bit for the sign
• 52 bits for the fractional part of the mantissa
• 11 bits for the exponent as an integer (ranging from -1023 to 1024)

Example

\(\epsilon\) = 10. Assume you have a normalized floating-point system with base 10 using 5 digits (1 for the sign, 2 for the exponent, and 2 for the mantissa). For simplicity, assume that one of the exponent digits is used for its sign. A general representation of the number following normalization would be \( \pm 0. m \cdot 10^e \).

What’s the largest possible positive number that can be represented? \(99.9 \cdot 10^{10}\) with both mantissa and exponent positive.

What’s the smallest possible positive number that can be represented? \(1.0 \cdot 10^{-9}\) with mantissa positive and exponent negative.

To obtain the respective negatives, change the sign in the first case to obtain the smallest possible negative number \(-9.9 \cdot 10^{19}\) and in the second case to obtain the largest possible negative number \(-1.0 \cdot 10^{-1}\).

Rounding and chopping

For any number, we denote its representation in the floating-point system. Assume that digits are available for the fractional part of the mantissa. There are two options to store the number:

• Chopping: Stop after we are out of digits available.
• Rounding: If the next digit is less than 5, keep the current digits. If it is 5 or more, round up.

The upper bound of the relative error is called machine precision:

• \(\epsilon_{\text{chopping}}\) = \(|\frac{x - fl(x)}{x}|\)
• \(\epsilon_{\text{rounding}}\) = \(|\frac{x - fl(x)}{x}|\)

This means that rounding has slightly higher machine precision, but rounding is overall better. For example, keeping 4 digits: 1.5636 → 1.563 (chopping) and 1.5636 → 1.564 (rounding).

Arithmetic operations

Sum of two real numbers:

• Express the numbers so they have the same exponent.
• Add the mantissas together.
• Normalize the number.
• Compute the floating-point result.

Multiplication and division of two real numbers:

• Perform the product or division of the mantissas of the two numbers.
• Sum or subtract the exponents.
• Normalize the number.
• Compute the floating-point result.

Main sources of errors in arithmetic operations

• Subtracting cancellation: Subtracting two nearly equal numbers can rise errors.
• Large computations: An individual round-off error could be small, but the cumulative effect over many operations can be significant.
• Loss of significant digits: Adding a very large number to a very small one can lead to a loss of precision.

Example: Add \(1.53 \times 10^8\) and \(2.01 \times 10^1\):

• \((1.53 \times 10^8) + (0.000000201 \times 10^8) = 1.530000201 \times 10^8\)

Conditioning and stability

A problem is said to be well-conditioned if small variations in the data correspond to small variations in the results. On the contrary, if small variations in the data lead to big variations in the results, the problem is ill-conditioned.

Example:

• \(4 \times (x - 1) = 0 \Rightarrow x = 1\)
• \(4 \times (x - 1) = 10^{-8} \Rightarrow x = 1 + 10^{-8} \)

In the second case, a perturbation on the right-hand side is introduced. This generates a perturbation also on the result. There are 6 orders of magnitude in difference between the perturbation introduced and the perturbation on the result, indicating that the problem is ill-conditioned.

An algorithm is said to be stable if it amplifies little the round-off errors introduced during computation.

Example: Consider calculating the series \( \{ \phi^n \} \), \( n = 0, 1, 2, \ldots \) of integer powers of the golden mean \(\phi = \frac{1 + \sqrt{5}}{2}\).

This can be done in two ways:

1. Set \( \phi^0 = 1 \) and calculate \( \phi^n = \phi \cdot \phi^{n-1} \) for \( n = 1, 2, \ldots \)
2. Set \(\phi^0 = 1\), \(\phi^1 = \phi\), and compute \( \phi^n = 2\phi^{n-1} - \phi^{n-2} \) for \( n = 2, 3, \ldots \)

Both methods are mathematically equivalent, but the second one is unstable since it gives entirely wrong results for large values of \( n \). This source of error in the second method is introduced from the loss of significant digits when adding/subtracting large and small numbers together.

Computational complexity

The speed of a computer is usually measured in FLOPS (Floating-point Operations Per Second).

The computational complexity is the rate at which the total number of FLOPS required to solve a problem of size \( n \) grows as \( n \) grows. Denoting the number of operations required to solve a problem by \( T(n) \), we have:

• The more efficient the algorithm, the slower \( T(n) \) grows as \( n \) grows.
• The more inefficient the algorithm, the faster \( T(n) \) grows as \( n \) grows.

We are always interested in the asymptotic bounds of the behavior of \( T(n) \) as \( n \to \infty \).

\( T(n) \sim O(f(n)) \Leftrightarrow |T(n)| \leq c \cdot f(n) \), for some positive \( c \).

This means that if \( T(n) \) is proportionate to \( f(n) \), it’s said to be superiorly limited by some scalar multiplied by the function \( f(n) \).

Linear systems and factorizations

Matrices are the fundamental element in computations involving linear algebra.

Matrix operations

• The sum of two matrices \( A \) and \( B \) of size \( m \times n \) is a matrix \( C = A + B \) of size \( m \times n \), where every element of \( C \) is \( c_{ij} = a_{ij} + b_{ij} \), for \( i = 1, \ldots, m \) and \( j = 1, \ldots, n \).
• The product of two matrices \( A \) of size \( m \times n \) and \( B \) of size \( n \times p \) is a matrix \( C = A \times B \) of size \( m \times p \), where every element of \( C \) is \( c_{ij} = \sum_{k=1}^{n} a_{ik}b_{kj} \), for \( i = 1, \ldots, m \) and \( j = 1, \ldots, p \).
• The transpose of a matrix \( A \), called \( A^T \), is obtained by turning every row into a column and vice versa.
• The determinant of a square matrix \( A \) is recursively defined.

If a square matrix \( A \) has \(\det(A) \neq 0\), it is said to be non-singular. In this case, the inverse \( A^{-1} \) of \( A \) is defined.

Matrices allow a compact notation for representing systems of linear equations:

\( a_{11}x_1 + a_{12}x_2 + \ldots + a_{1n}x_n = b_1 \)

\( a_{21}x_1 + a_{22}x_2 + \ldots + a_{2n}x_n = b_2 \)

\(\ldots\)

\( a_{m1}x_1 + a_{m2}x_2 + \ldots + a_{mn}x_n = b_m \)

Can be represented as \( A\mathbf{x} = \mathbf{b} \), where:

• \( A \) is the matrix containing the coefficients
• \( \mathbf{x} \) is the vector of unknowns
• \( \mathbf{b} \) is the vector of solutions

If the system has more equations (rows) than unknowns (columns), the system is said to be overdetermined, while if the system has fewer equations than unknowns, it is said to be underdetermined.

In principle, if \( A \) is non-singular, we could obtain the solution by computing the inverse:

\( A\mathbf{x} = \mathbf{b} \Rightarrow \mathbf{x} = A^{-1}\mathbf{b} \)

But this is NOT an efficient way to solve a system of equations since computing \( A^{-1} \) is very expensive and unstable.

Solving linear systems: Cramer’s Rule

Each unknown variable in a system of linear algebraic equations can be expressed as a fraction of two determinants:

• The denominator is the determinant of the matrix \( A \).
• The numerator is the determinant of the matrix where the column of the coefficients of the considered unknown variable is substituted with the known values \( \mathbf{b} \).

Note that Cramer’s Rule is inefficient as the number of equations grows.

Strategy: Elimination of the unknowns

1. Multiply the equations by constants.
2. Combine two equations so that one of the unknowns is eliminated.
3. Solve the resulting single equation for the remaining unknown.
4. Substitute this value into either of the original equations to compute the other variable.

Example

Apply the elimination method to the following equations:

\( a_{11}x_1 + a_{12}x_2 = b_1 \)

\( a_{21}x_1 + a_{22}x_2 = b_2 \)

1. Multiply eq.1 by \( a_{21} \) and eq.2 by \( a_{11} \):
• \( a_{21}a_{11}x_1 + a_{21}a_{12}x_2 = a_{21}b_1 \)
• \( a_{11}a_{21}x_1 + a_{11}a_{22}x_2 = a_{11}b_2 \)
2. Subtract eq.3 from eq.4 eliminating \( x_1 \):
• \( a_{11}a_{22}x_2 - a_{21}a_{12}x_2 = a_{11}b_2 - a_{21}b_1 \)
3. Solve the eq. for \( x_2 \):
• \( x_2 = \frac{a_{11}b_2 - a_{21}b_1}{a_{11}a_{22} - a_{21}a_{12}} \)
4. Substitute the value into the original eq. to get the value of the other unknown \( x_1 \):
• \( x_1 = \frac{b_1a_{22} - b_2a_{12}}{a_{11}a_{22} - a_{21}a_{12}} \)

This becomes challenging with many equations and unknowns.

General rule: Naive Gaussian Elimination

The elimination of the unknowns technique can be extended to large sets of equations by developing a systematic scheme or algorithm to eliminate unknowns and back-substitute: Gaussian Elimination. We need some modifications to obtain a reliable algorithm, particularly avoiding divisions by zero. The “naive” version does not address this problem. Nonetheless, the approach is designed to solve a general set of equations:

\( a_{11}x_1 + a_{12}x_2 + \ldots + a_{1n}x_n = b_1 \)

\( a_{21}x_1 + a_{22}x_2 + \ldots + a_{2n}x_n = b_2 \)

\(\ldots\)

\( a_{m1}x_1 + a_{m2}x_2 + \ldots + a_{mn}x_n = b_m \)

Steps for Gaussian Elimination

Step 1: Forward elimination of unknowns

The goal is to get an upper triangular matrix. We need to eliminate the corresponding unknowns accordingly:

• Eliminate the first unknown from the second through the m-th equations.
• Substitute into the equation to get \( x_{n-1} \), then substitute into the n-2 equation, and so on, proceeding backwards.

Substep 1: Eliminate the first unknown from the second through the m-th equations:

\( a_{21}x_1 + \ldots + a_{2n}x_n = b_2 - a_{21}x_1 \)

\( \ldots \)

Perform similar operations for the subsequent equations to obtain an upper triangular matrix.