Numerical Analysis: Computational Methods

I prepared my numerical-analysis notes during undergraduate Computer Engineering courses in the 2013-2014 period. This revision preserves the original progression from error analysis and matrix operations to linear-system solvers, nonlinear root finding, and interpolation, while making convergence and stability statements more explicit. Later references were used to review terminology and numerical qualifications without changing the historical course structure.

Unit 1: Introduction to Numerical Analysis

Subject of the course

Numerical analysis studies algorithms that approximate mathematical quantities when an exact symbolic solution is unavailable, impractical, or unnecessarily expensive. The result of a numerical method is meaningful only together with its approximation error, convergence behavior, and computational cost.

Power series

A power series centered at a has the form

Σ c_n (x-a)^n.

Inside its interval or radius of convergence it provides a representation that can be truncated to build polynomial approximations.

Maclaurin series

A Maclaurin series is a Taylor series centered at zero:

f(x) = f(0) + f'(0)x + f''(0)x^2/2! + ...

Examples such as the exponential, sine, and cosine series are useful because a finite number of terms can approximate transcendental functions.

Taylor series

Around x=a,

f(x) = f(a) + f'(a)(x-a) + f''(a)(x-a)^2/2! + ...

The truncated polynomial is a local approximation. Increasing the order can improve the approximation near the expansion point, but numerical effects and the convergence region still matter.

Error term in series

A finite Taylor polynomial leaves a remainder. The Lagrange form is commonly written as

R_n(x) = f^(n+1)(ξ) (x-a)^(n+1)/(n+1)!

for a suitable ξ between a and x. The remainder gives a theoretical way to bound truncation error when derivative bounds are available.

Polynomial approximation

Polynomials are computationally convenient because they require only addition and multiplication. Numerical methods therefore frequently replace a difficult function locally by a polynomial while controlling the approximation domain and error.

Types of error

Two broad error sources are important:

• truncation error, caused by replacing an infinite or exact mathematical process with a finite approximation;
• rounding error, caused by finite-precision representation and arithmetic.

Input uncertainty and modeling error are separate issues and should not be confused with floating-point error.

Approximate error

When the exact value is known, absolute and relative errors can be computed directly. In iterative methods the exact answer is normally unknown, so the change between successive approximations is used as an observable stopping measure.

Error control in iterations

A common approximate relative error is

ε_a = |(x_i - x_(i-1))/x_i|.

The iteration may stop when this quantity or another residual/error criterion falls below a specified tolerance. A small step does not always imply a small true error, so the criterion should match the method.

Rounding errors

Floating-point arithmetic represents only a finite subset of real numbers. Repeated operations can accumulate error, and subtraction of nearly equal values can cause cancellation. Numerical algorithms should therefore be judged not only by algebraic correctness but also by stability.

Unit 2: Matrices

Matrix concept

A matrix is a rectangular array of scalars. In numerical analysis, matrices encode systems of equations, transformations, discretized models, and many data-processing problems.

Addition and scalar multiplication

Matrices of the same size are added element by element. Scalar multiplication multiplies each entry by the same scalar. These operations preserve the matrix dimensions.

Matrix multiplication

If A is m x n and B is n x p, then AB is m x p, with

(AB)_ij = Σ_k a_ik b_kj.

Matrix multiplication is generally not commutative.

Special matrices

Important forms include zero, identity, diagonal, triangular, symmetric, and orthogonal matrices. Recognizing structure can reduce computation and improve numerical methods.

Elementary row operations

The three elementary row operations are row exchange, multiplication of a row by a nonzero scalar, and adding a multiple of one row to another. They form the basis of elimination.

Gaussian elimination

Gaussian elimination transforms a linear system to upper-triangular form and then uses back substitution. A pivot strategy is important when a prospective pivot is zero or numerically small.

LU factorization

For a suitable matrix,

A = LU,

where L is lower triangular and U is upper triangular. Once the factorization is available, multiple systems with the same coefficient matrix and different right-hand sides can be solved efficiently.

Cholesky factorization

For a real symmetric positive-definite matrix,

A = LL^T.

Cholesky uses the symmetry and positive-definite structure and is more efficient than a general factorization for this class.

Matrix norms

A matrix norm measures matrix magnitude and supports error and conditioning analysis. Common examples include the 1-norm, infinity norm, and Frobenius norm. Different norms emphasize different aspects of the matrix.

Determinants

The determinant is a scalar associated with a square matrix. A zero determinant indicates singularity. Determinant computation by cofactor expansion becomes expensive as dimension increases.

Determinant by Gaussian elimination

After triangularization, the determinant is related to the product of diagonal entries, with sign changes for row exchanges and scaling adjustments for row multiplication. This is normally more practical than large cofactor expansions.

Chio method

Chio condensation reduces the order of a determinant through a sequence of algebraic transformations. It is useful as a classical determinant technique but is not the main tool for modern large numerical systems.

Inverse matrix

For nonsingular A, the inverse satisfies

AA^-1 = A^-1A = I.

In practice, explicitly forming A^-1 is usually unnecessary when the actual task is to solve Ax=b. Direct solution with suitable pivoting or factorization normally uses less work and has better numerical behavior. The inverse should be computed when the inverse itself is required.

Unit 3: Systems of Linear Equations

Structure of the system

A linear system is written compactly as

Ax = b.

Its behavior depends on the rank and conditioning of A, not merely on the number of equations.

Solution rules

A system can have a unique solution, no solution, or infinitely many solutions. For a square nonsingular matrix, a unique solution exists. For general systems, rank conditions give the correct classification.

Cramer's rule

For a nonsingular square system,

x_i = det(A_i)/det(A).

Cramer's rule is valuable theoretically and for very small systems, but determinant-based evaluation is not efficient for large numerical computation.

Gaussian elimination

Elimination is a standard direct method. It reduces the coefficient matrix to triangular form and recovers the solution by substitution.

Pivoting

Pivoting exchanges rows, and in some variants columns, to avoid zero or very small pivots. Partial pivoting is a common practical strategy and substantially improves the robustness of Gaussian elimination for many problems.

Gauss-Jordan method

Gauss-Jordan elimination continues until the coefficient matrix is reduced toward the identity matrix. It is convenient for reduced-row-echelon form and small inverse calculations, but it performs more operations than basic Gaussian elimination for a single right-hand side.

Solution with LU factorization

With A=LU, solve

Ly = b
Ux = y.

The factorization is especially useful when many right-hand sides share the same matrix.

Solution with Cholesky factorization

If A is symmetric positive definite, Cholesky gives a specialized and efficient direct solver. The structural assumptions must be checked rather than inferred from appearance alone.

Iterative methods

Iterative methods start from an initial approximation and repeatedly improve it. They are especially attractive for large sparse systems, while direct methods are often suitable for smaller dense systems.

Convergence condition

Convergence depends on the iteration matrix. Conditions such as strict diagonal dominance can provide convenient sufficient criteria, but they are not the only possible convergence conditions.

Jacobi method

Jacobi computes every new component from the previous iteration:

x_i^(k+1) = (b_i - Σ_(j≠i) a_ij x_j^k)/a_ii.

Because all components use the old vector, the basic method exposes parallelism naturally.

Gauss-Seidel method

Gauss-Seidel immediately reuses newly computed components within the same iteration. On suitable problems it may converge faster than Jacobi. Its data dependency makes naïve direct parallelization harder than Jacobi, although block, coloring, and other parallel variants are possible.

Unit 4: Nonlinear Equations

Problem and method classes

Root finding seeks x such that

f(x) = 0.

Methods can be grouped broadly into bracketing methods, which maintain an interval containing a sign change under suitable continuity conditions, and open methods, which generate approximations without maintaining such a bracket.

Bracketing methods provide strong convergence behavior when their assumptions are satisfied. Open methods can be faster but do not have a general convergence guarantee.

Bisection method

If f is continuous on [a,b] and

f(a)f(b) < 0,

then at least one root lies inside. Bisection repeatedly halves the interval and retains the half with the sign change.

After n steps the interval width is reduced by 2^n, so an explicit upper bound on the number of steps needed for a target interval tolerance can be computed from the initial bracket.

A standard sign-change bisection test generally does not detect an even-multiplicity root where the function only touches the axis without changing sign.

Stopping iterations

Possible stopping rules include interval width, successive-iterate change, residual magnitude |f(x)|, or a combination. The criterion should reflect the required accuracy and the scaling of the problem.

Regula falsi

False position also preserves a sign-changing bracket but uses the secant line through the bracket endpoints to propose the next point. It can outperform bisection on some functions, although one endpoint may remain fixed for many iterations.

Fixed-point iteration

Rewrite the equation as

x = g(x)

and iterate

x_(k+1) = g(x_k).

Local convergence is related to the behavior of g' near the fixed point; a convenient sufficient condition is |g'(x)| < 1 on an appropriate neighborhood.

Newton-Raphson method

Newton's method uses the tangent line:

x_(k+1) = x_k - f(x_k)/f'(x_k).

Near a simple root and under appropriate smoothness and derivative conditions, convergence is locally quadratic. A poor initial guess or a derivative near zero can produce large steps, loss of stability, or failure to converge.

Secant method

The secant method approximates the derivative with the slope through the two most recent points:

x_(k+1) = x_k - f(x_k)(x_k-x_(k-1)) / [f(x_k)-f(x_(k-1))].

It avoids explicit derivative evaluation. It is an open method and, although it can converge rapidly from suitable starting points, it has no general convergence guarantee.

Müller method

Müller's method fits a quadratic through three recent points and uses a root of that quadratic as the next approximation. It can naturally produce complex iterates and is useful in polynomial and complex-root searches, but it does not guarantee that all roots will be found.

Comparing methods

No single method dominates every problem. Bisection emphasizes robustness under its bracketing assumptions. Newton emphasizes fast local convergence when a reliable derivative and suitable initial estimate are available. Secant avoids derivatives. Regula falsi keeps a bracket. Müller uses quadratic interpolation and can enter the complex plane.

Unit 5: Interpolation

Problem

Given data

(x0,y0), (x1,y1), ..., (xn,yn),

interpolation constructs a function that passes through the supplied data points. Polynomial interpolation is a classical approach.

Simple interpolation

For two points, linear interpolation uses the straight line joining them. It is local and inexpensive but cannot represent curvature between widely separated points.

Lagrange interpolation

The Lagrange polynomial is

P_n(x) = Σ_i y_i L_i(x),

where

L_i(x) = Π_(j≠i) (x-x_j)/(x_i-x_j).

It gives the unique polynomial of degree at most n through n+1 distinct nodes.

Newton interpolation polynomial

Newton's divided-difference form writes the same interpolating polynomial incrementally:

P_n(x) = f[x0]
       + f[x0,x1](x-x0)
       + ...

Its structure is convenient when new nodes are added because previous coefficients need not all be discarded.

Danger of high degree

Increasing polynomial degree does not guarantee increasing accuracy. With some node distributions, high-degree interpolation can oscillate strongly, especially near interval endpoints. This is the classical Runge phenomenon.

Interpolation versus curve fitting

Interpolation forces the model through every supplied point. Curve fitting or regression instead estimates a lower-dimensional relationship and need not pass exactly through noisy observations. Extrapolation evaluates outside the data range and is generally more risky because it depends more strongly on assumptions about behavior beyond the observed interval.

General Conceptual Framework

Numerical work can be organized around a few recurring questions:

Mathematical problem
      ↓
Approximation or factorization
      ↓
Algorithm
      ↓
Finite-precision computation
      ↓
Error / residual / convergence check
      ↓
Interpretation of the numerical result

An algebraically valid formula is not automatically a good numerical algorithm. Cost, conditioning, stability, stopping criteria, and problem structure all matter.

Direct and iterative linear solvers address different computational settings. Root-finding methods trade robustness, derivative requirements, bracketing, and local speed. Interpolation reproduces data points but does not automatically provide a reliable model outside their range.

Conceptual Distinctions

Truncation error ≠ rounding error. Truncation replaces an exact process with a finite approximation; rounding arises from finite-precision arithmetic.

True error ≠ approximate iterative error. The exact error requires the exact value; successive-iterate differences are only an observable estimate.

Linear-system solution ≠ explicit matrix inversion. Forming the inverse is usually unnecessary. Direct solution with suitable pivoting or factorization is commonly cheaper and numerically preferable.

LU ≠ Cholesky. LU is a general triangular factorization under suitable conditions; Cholesky specifically exploits symmetric positive-definite structure.

Jacobi ≠ Gauss-Seidel. Jacobi uses only the previous iterate; Gauss-Seidel immediately reuses updated components.

Bracketing method ≠ unconditional guarantee. Bisection and false position rely on continuity and a valid sign-changing bracket. Under these conditions the bracket provides controlled convergence.

Bisection gives an a priori step bound. The initial interval and target tolerance provide a direct upper bound on the number of halvings required.

Newton-Raphson is locally fast, not globally guaranteed. Near a simple root under appropriate conditions it is quadratically convergent; poor starting values or small derivatives can cause large steps or failure.

Interpolation ≠ extrapolation. Extrapolation leaves the observed interval and is generally more assumption-sensitive.

Higher polynomial degree ≠ automatically higher accuracy. Node placement, conditioning, noise, and oscillation can make a higher-degree interpolant worse.

References

• Ahmet Yesevi Üniversitesi Bilgisayar Mühendisliği Bölümü. Sayısal Çözümleme (TBIL301) ders materyalleri.
• Richard L. Burden, J. Douglas Faires, Annette M. Burden. Numerical Analysis, 10th Edition. Cengage, 2015.
• Steven C. Chapra, Raymond P. Canale. Numerical Methods for Engineers, 7th Edition. McGraw-Hill, 2015.
• Gene H. Golub, Charles F. Van Loan. Matrix Computations, 4th Edition. Johns Hopkins University Press, 2013.
• Lloyd N. Trefethen, David Bau III. Numerical Linear Algebra. SIAM, 1997.
• Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms, 2nd Edition. SIAM, 2002.
• William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery. Numerical Recipes: The Art of Scientific Computing, 3rd Edition. Cambridge University Press, 2007.