Numerical Analysis: Computational Methods

Numerical analysis explains systematically why a number produced by a computer is often only an approximation. Error analysis, root finding, interpolation, numerical differentiation and integration, and linear-system methods are connected through convergence and conditioning.

In numerical computation, an algorithm completing successfully is not the same as the result being trustworthy. Rounding error, truncation error, ill-conditioning, and algorithmic instability can look similar at the output while requiring different remedies. Each method is therefore read together with its error behaviour, not only as a sequence of operations.

Unit 1: Introduction to Numerical Analysis

1.1. 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.

1.2. 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.

1.3. 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.

1.4. 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.

1.5. 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.

1.6. 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.

1.7. 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.

1.8. 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.

1.9. 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.

1.10. 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.

1.11. Forward Error, Backward Error, and Conditioning

Three different questions should be separated when assessing a numerical result. Forward error compares the computed solution with the exact solution. Backward error asks for the small perturbed problem that the computed result solves exactly. Conditioning measures how strongly the exact solution reacts to small perturbations in the input data.

A small backward error does not imply a small forward error for an ill-conditioned problem. Conversely, a well-conditioned problem can be damaged by an unstable algorithm. It is therefore not enough to report that double precision was used or that the residual is small.

For a linear system, the residual

r = b - A x_hat

may be small while x_hat is inaccurate; the severity depends on the conditioning of A.

1.12. Reading Error, Stability, and Convergence Together in Numerical Analysis

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.

1.13. Distinguishing Conditioning, Stability, Accuracy, and Convergence

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.

Choosing a numerical method is a concrete engineering trade-off among accuracy, cost, and stability. More iterations do not automatically produce a better answer; on an ill-conditioned problem or with an unstable method, a larger computational budget can simply produce the wrong result more expensively.

1.14. Error Budgets from Models to Floating Point

Numerical error has several sources: modeling assumptions, measurement uncertainty, discretization, incomplete iteration, and floating-point rounding.

Using double does not automatically make a computation accurate. Conditioning comes first, algorithmic stability second, and representation precision after that.

Stopping criteria should combine absolute and relative tolerances with residual checks. A small change between iterations does not always mean the underlying equation is well satisfied.

1.15. Building an Error Budget for Numerical Results

Numerical discrepancy is not caused only by rounding. Modeling error, input uncertainty, discretization, iteration tolerance, and floating-point representation are distinct sources. Increasing precision may not help if another source dominates.

In a convergence study, refining the step or mesh should produce the expected order of convergence. If it does not, investigate implementation errors, boundary conditions, or violated method assumptions.

When no reference solution exists, residuals, conserved quantities, and an independent second method provide useful cross-checks. Running the same code with two parameter values is not an independent validation.

1.16. Choose the Numerical Method Before Trusting the Number

Producing a number is not the same as producing a reliable approximation. Error can originate in the data or model, in truncation introduced by the method, and in rounding caused by finite-precision arithmetic. Separating these sources helps explain why a result changes.

Bisection and Newton's method illustrate different guarantees. Bisection narrows a sign-changing interval slowly but predictably when the continuity assumptions hold. Newton's method can converge very rapidly near a suitable solution, yet it can fail when the derivative is small, the initial point is poor, or the function geometry is unfavourable.

Conditioning belongs to the problem; stability belongs to the algorithm. A well-conditioned problem can be damaged by an unstable method. In an ill-conditioned problem, small perturbations in the input can create large output changes even when the numerical algorithm is implemented correctly.

Interpolation and regression also answer different questions. Interpolation constructs a function that passes through the supplied points; regression can model an underlying trend in noisy data. A high-degree polynomial passing through every sample does not guarantee reliable behaviour between samples.

For iterative methods, iteration count alone is a weak stopping criterion. Step size, residual, and the scale of the problem should be considered together. A small |x_(k+1)-x_k| does not necessarily mean that the governing equation is well satisfied. The stopping rule should be tied to the mathematical quantity the method is intended to approximate.

1.17. Separating Numerical Stability from Statistical Robustness

Robust statistics asks how much an estimator degrades under small departures from an idealized model. In numerical experiments, the median and MAD can summarize outlying runtimes, residuals, or sampled errors; influence functions and breakdown points then help characterize whether the estimator that aggregates repeated runs is being dominated by a small number of corrupted observations. In heavy-tailed or contaminated data, decisions based only on the mean and standard deviation can be misleading.

• Separate amplification of floating-point error in ill-conditioned problems from estimator contamination by outliers.
• Distinguish residual magnitude from true forward error.
• Measure iterative-solver cost separately when using robust objective functions.

1.18. Condition Number, Backward Error, and Solution Stability

A large numerical solution error does not necessarily imply a poor algorithm: some mathematical problems amplify small perturbations in their inputs. For \(Ax=b\), the condition number \(\kappa(A)=\|A\|\,\|A^{-1}\|\), under compatible matrix norms, characterizes sensitivity to relative perturbations. When \(\kappa(A)\) is large, even a well-designed solver can produce substantial forward error as a consequence of input rounding or measurement error. Backward error instead asks which nearby problem has been solved by the computed result.

For a computed solution \(\widehat{x}\), define the residual \(r=b-A\widehat{x}\). A small \(\|r\|\) is not on its own a guarantee of small \(\|x-\widehat{x}\|\); the relationship also depends on conditioning. In a matrix whose columns are nearly linearly dependent, a tiny residual can coexist with a substantial change in the solution coefficients. A useful validation report therefore separates residual norm, relative backward error, estimated condition number, and, where available, comparison with an independent reference solution.

Partial pivoting in Gaussian elimination, QR factorization, and suitable scaling can improve algorithmic stability; they do not magically make an ill-conditioned problem well conditioned. Solver benchmarks must measure numerical error alongside execution time and state the limits imposed by the floating-point representation.

Unit 2: Matrices

2.1. 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.

2.2. 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.

2.3. 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.

2.4. Special Matrices

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

2.5. 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.

2.6. 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.

2.7. 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.

2.8. 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.

2.9. 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.

2.10. 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.

2.11. 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.

2.12. 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.

2.13. 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.

2.14. How a Small Singular Value Changes a Regularized Solution

Let A = diag(1, 0.01) and b = [1, 0.015]ᵀ. Direct inversion gives x = A⁻¹b = [1, 1.5]ᵀ. Yet if the second observation contains a disturbance of only 0.005, its underlying second component could be 1. The matrix has a 2-norm condition number of 100, allowing small relative perturbations in unfavorable data directions to cause much larger relative solution errors.

Tikhonov regularization minimizes ||Ax-b||² + λ||x||², with the algebraic condition (AᵀA + λI)x = Aᵀb. For λ = 0.0001, the second component is 0.01·0.015/(0.01²+0.0001) = 0.75, while the first is approximately 0.9999. Noise amplification has been reduced, but the answer is biased if the true second component is 1. The objective is not merely to shrink the solution; it is to balance data evidence against uncertainty and resolution. Explicitly forming normal equations can worsen conditioning, so practical implementations generally prefer an SVD-based or otherwise numerically appropriate factorization.

Unit 3: Systems of Linear Equations

Method selection for a linear system begins with size, sparsity, symmetry, and conditioning. Direct factorization may suit small dense systems, while large sparse problems call for checking iterative convergence and possible preconditioning. Report the residual norm, but distinguish residual from forward error when the problem is ill-conditioned. Memory layout and numerical stability matter alongside execution time.

3.1. 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.

3.2. 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.

3.3. 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.

3.4. Gaussian Elimination

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

3.5. 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.

3.6. 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.

3.7. 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.

3.8. 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.

3.9. 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.

3.10. 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.

3.11. 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.

3.12. 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

Root finding begins by identifying a plausible domain and checking the function's continuity assumptions. Bracketing methods can rely on a sign change, whereas derivative-based Newton iteration is more sensitive to its initial guess. Track the function residual as well as the difference between successive estimates. If steps become small while the residual remains large, or derivatives approach zero, reconsider the method.

4.1. 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.

4.2. 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.

4.3. 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.

4.4. 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.

4.5. 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.

4.6. 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.

4.7. 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.

4.8. 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.

4.9. Constructing a Numerically Stable Müller Step

For successive points x0, x1, and x2, define h0 = x1-x0, h1 = x2-x1, d0 = (f(x1)-f(x0))/h0, and d1 = (f(x2)-f(x1))/h1. The interpolating quadratic centered at x2 has coefficients:

a = (d1-d0)/(h0+h1)
b = a*h1+d1
c = f(x2)
D = b*b-4*a*c
x3 = x2-2*c/(b+sign_choice*sqrt(D))

Select the denominator with larger magnitude from b+sqrt(D) and b-sqrt(D); use complex magnitudes when necessary. This avoids subtracting nearly equal numbers in the root formula. If either point spacing, their sum, or the chosen denominator is close to zero, the interpolation step is ill-conditioned; choose new points or fall back to an appropriate secant step. A negative real discriminant calls for complex arithmetic when complex roots are allowed, not for silently clipping it to zero. A stopping rule should consider |f(x3)| as well as step size and the iteration budget.

4.10. 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.

4.11. Conditioning and Algorithmic Stability Are Different Questions

Conditioning describes the problem: how much the exact mathematical answer changes when the input is perturbed. Stability describes the algorithm: how much extra error the computational procedure introduces.

A stable algorithm cannot repair an inherently ill-conditioned problem, and a well-conditioned problem can still be solved poorly by an unstable algorithm.

This distinction should be checked before interpreting a small residual as proof of an accurate solution.

4.12. Richardson Extrapolation and Romberg Integration

If a numerical approximation has a known leading error term, two step sizes can be combined to cancel part of that error. Richardson extrapolation uses this idea generically.

Romberg integration applies repeated Richardson extrapolation to trapezoidal-rule estimates with successively halved step sizes. It can converge quickly for sufficiently smooth integrands, but discontinuities, endpoint singularities, and highly oscillatory functions require a different strategy.

Deriving the stopping guarantee before solving for a root. For f(x) = x³ − 7, the values f(1) = −6 and f(2) = 1 bracket a root in [1,2]. The function is continuous and strictly increasing there, so the root is unique. Each bisection iteration retains the half-interval with opposite signs at its endpoints; after n iterations, the remaining width is at most 1/2ⁿ. Requiring width below 0.001 gives 2ⁿ > 1000, hence n ≥ 10 suffices. Reporting the final midpoint guarantees absolute root error at most 1/2048 < 0.0005. This bound differs from the usually faster local convergence of Newton's method, which depends on the starting point and derivative behavior. The termination rule should follow from the error requirement rather than an arbitrary iteration count.

Unit 5: Interpolation

5.1. 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.

5.2. 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.

5.3. 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.

5.4. 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.

5.5. 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.

5.6. 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.

5.7. Interpolation in Imaging, Graphics, and Signal Processing

Interpolation is not limited to fitting tabulated data. When an image is enlarged, new pixel centers fall between samples of the original grid. Under a geometric transform, a destination coordinate usually maps to a noninteger source location. The missing sample must therefore be estimated from neighboring data.

Nearest-neighbor selects the closest sample. Bilinear interpolation performs a linear interpolation along two axes. Bicubic interpolation uses a larger neighborhood and a smoother local model. Other kernels are also possible. The choice trades sharpness, smoothness, ringing, computational cost, and repeated-resampling behavior.

The mathematical principle belongs here. Pixel coordinates, inverse mapping, color channels, and visible artifacts belong to Image and Audio Processing; sampling theory and aliasing belong to Signals and Systems.

5.7.1. Interpolation as a Kernel and Its two-dimensional Use

Resampling by interpolation can be written as a weighted sum: each neighboring sample is multiplied by a kernel function evaluated at its distance from the new position, and the products are added. Nearest-neighbor corresponds to a box kernel and bilinear interpolation to a triangle kernel. Keys (1981) derives a piecewise-cubic convolution kernel that uses four neighboring samples per axis and reaches third-order accuracy in the sampling interval; its order of accuracy lies between that of linear interpolation and cubic splines. For images the two-dimensional version is separable: the same one-dimensional kernel is applied along rows and then along columns. Mitchell and Netravali (1988) generalize this into a two-parameter (B, C) family of piecewise cubic filters, show that blurring, ringing, and angle-dependent (anisotropic) artifacts must be traded against one another, and report a reasonable compromise in a small region of the parameter space around B = C = 1/3.

5.7.2. Enlarging and Reducing Are Different Problems

When an image is enlarged, the destination samples fall between source samples and the task is to estimate the missing ones; the kernel governs reconstruction quality. When an image is reduced, the sampling rate drops. The source can contain spatial frequencies that the destination grid cannot carry, and those components produce aliasing; the principle behind this condition is the Nyquist-Shannon sampling theorem, treated in Signals and Systems. Mitchell and Netravali separate aliasing into two sources: prealiasing, caused by insufficient filtering before sampling, and postaliasing, caused by poor reconstruction. The classical remedy is to filter before sampling. Reading a bilinear or bicubic value at the destination coordinates therefore does not suppress high-frequency components when an image is shrunk; a low-pass prefilter or area averaging is required. Heckbert (1989) treats, systematically, the filtering methods that prevent aliasing during resampling in image warping and texture mapping.

5.7.3. Resampling Adds No Information; learning-based super-resolution Is Different

Interpolation does not recover high-frequency detail that the source never contained; it only produces intermediate values under the smoothness assumption built into the chosen kernel. Learning-based methods for super-resolution instead learn a mapping from pairs of low- and high-resolution images (Dong et al., 2016) and estimate the detail from a learned prior. That detail is a generated estimate rather than a measurement, so the two approaches cannot be trusted equally in measurement and verification tasks.

5.7.4. Application Pipeline: Perceptual Hashing

Resampling is often not a goal in itself but part of the definition of a downstream feature. Perceptual hash methods reduce the input to a small fixed working area, derive frequency-domain features from the intensity values, and quantize them into a short binary signature. In such a pipeline the reduction step does two jobs: it decouples computational cost from the input resolution, and it removes fine detail so that coarse structure dominates. By construction, a signature bit changes only when a coefficient crosses its threshold; coefficients close to the threshold can therefore be sensitive to small differences in the kernel and prefilter. Signatures that will be compared within one collection should thus be produced by the same resampling pipeline, and the acceptance threshold should be chosen on a validation set built with that pipeline.

A worked example of this chain is documented on this site. The article Perceptual Image Similarity with DCT and SIMD describes reducing the image to a fixed 64×64 working area while preserving the aspect ratio, explains that this makes the computational cost independent of the input resolution, and notes that a large part of the high-frequency fine detail is removed during resampling; the corresponding C# implementation is archived as a persistent software record on Zenodo (doi:10.5281/zenodo.22117445; DCT Perceptual Hash). The article explains the role of the standardization step in cost and detail removal rather than a particular kernel choice. The general framework for images, color channels, and artifacts continues in Image and Audio Processing. A broader published provenance layer for this application domain is the co-authored Image and Audio Processing — Book Chapter; its published table of contents covers image recording and communication, image feature vectors, restoration/enhancement, segmentation, and AI applications. The published contents do not document a particular interpolation kernel, so the chapter is linked as verifiable author provenance in image processing rather than as the source of the interpolation theory used here.

Unit 6: Numerical Differentiation and Integration

When an analytical derivative or integral is unavailable, sampled function values can be used to construct numerical approximations. Truncation error and floating-point roundoff must be considered together.

6.1. Finite-difference Derivatives

The forward difference is:

f'(x) ≈ [f(x+h) - f(x)] / h

and the central difference is:

f'(x) ≈ [f(x+h) - f(x-h)] / (2h)

For sufficiently smooth functions, the central formula has a higher-order truncation error for the same step size.

Making h indefinitely smaller does not indefinitely improve the answer. When h is too small, subtracting nearly equal values amplifies cancellation and roundoff. A useful step lies in the region where truncation and floating-point error are balanced.

6.2. Trapezoidal and Simpson Rules

For an integral over [a,b], the trapezoidal rule approximates the function by a line:

∫ f(x) dx ≈ (b-a) [f(a)+f(b)] / 2

The composite rule applies the same idea over subintervals.

Simpson's rule uses a quadratic approximation and can provide higher accuracy for sufficiently smooth functions. Higher formal order is not automatically better for every dataset; smoothness and sampling geometry must satisfy the assumptions of the method.

6.3. Gaussian Quadrature

Rather than evaluating at equally spaced points, Gaussian quadrature chooses specific nodes and weights so that a high polynomial degree can be integrated exactly with relatively few function evaluations.

The engineering choice is therefore not only "which formula has higher order?" Function-evaluation cost, singularities, fixed sample locations, and the required error bound all matter.

6.4. Adaptive Differentiation and Quadrature

For finite-difference derivatives, reducing the step decreases truncation error only until floating-point cancellation and rounding begin to dominate. The mathematical limit h -> 0 must not be interpreted computationally as "choose the smallest representable h."

Similarly, a fixed integration mesh may waste samples in smooth regions and undersample difficult ones. Adaptive quadrature subdivides intervals where a local error estimate is large and leaves easy regions coarse.

Richardson extrapolation combines results at different step sizes to cancel a leading error term. Romberg integration systematizes the idea with the trapezoidal rule. A sound stopping criterion is based on a scale-aware error estimate rather than merely on the number of subintervals.

Unit 7: Initial-Value Problems

An initial-value problem has the form:

y' = f(t, y)
y(t0) = y0

When no analytical solution is available, time is discretized and an approximate state is advanced step by step.

7.1. Euler's Method

The simplest explicit method is:

y_(n+1) = y_n + h f(t_n, y_n)

It advances along the local tangent. Euler's method is valuable pedagogically, but error may accumulate rapidly. Reducing the step generally improves accuracy while increasing computational cost.

7.2. Heun and Runge-Kutta Methods

Heun's method combines the slope at the beginning of the interval with a predicted slope near its end. Classical fourth-order Runge-Kutta (RK4) evaluates four slopes:

k1 = f(tn, yn)
k2 = f(tn+h/2, yn+h k1/2)
k3 = f(tn+h/2, yn+h k2/2)
k4 = f(tn+h, yn+h k3)

and updates:

y_(n+1) = y_n + h(k1 + 2k2 + 2k3 + k4)/6

7.3. Stability and Stiff Systems

A numerical method can make an analytically stable system appear unstable when the step size is inappropriate. Formal order of accuracy and numerical stability are therefore separate properties.

Stiff systems may force explicit methods to use very small steps. Implicit methods can be more appropriate even though each step is computationally more expensive.

A reliable numerical ODE workflow should:

1. reduce the step and check convergence,
2. verify physical or mathematical invariants,
3. compare with another method when practical,
4. report tolerances and failure criteria explicitly.

7.4. Adaptive Step Size and Local versus Global Error

For an ODE solver, local truncation error describes the error introduced by one idealized step starting from the exact state. Global error is the accumulated difference after many steps.

Adaptive methods estimate local error by comparing two approximations of different orders or step sizes. The step is accepted, rejected, enlarged, or reduced according to a tolerance contract.

A smaller step is not automatically better without limit. It increases the number of operations and can amplify roundoff, while stiff systems may require methods with a different stability region rather than simply a smaller explicit step.

Balancing error scale and step cost in an ODE solver. For y' = −y, y(0) = 1, the exact solution is y(t) = e^(−t), with reference value 0.367879 at t = 1. Forward Euler with fixed step h = 0.25 gives y₄ = (1−0.25)^4 = 0.31640625, an absolute error of about 0.05147. Halving the step to h = 0.125 requires eight steps and gives y₈ = (0.875)^8 ≈ 0.343609, reducing error to about 0.02427. An error reduction of approximately twofold is consistent with first-order global convergence, while the number of steps has doubled. Step refinement must also respect the solver's stability region and the system's fast dynamics. When the exact solution is unknown, the difference between refined solutions can support an appropriate error estimator, but a small difference is not universally a proof of accuracy.

Unit 8: Adaptive ODE Solvers, Boundary-Value Problems, PDEs, and Automatic Differentiation

8.1. Local and Global Error with Adaptive Step Size

The local truncation error of an ODE method measures one idealized step, while the global error reflects accumulation across the complete numerical trajectory. Methods of the same formal order can still behave differently because their stability regions and step-size policies differ.

Adaptive solvers commonly compare two approximations of different order to estimate local error and then increase or decrease the step size. A typical scale-aware target resembles

error <= atol + rtol * scale.

Explicit methods compute the new value directly from known information. Implicit methods require solving an equation containing the new state. For stiff systems, a more expensive implicit step can be more efficient than forcing an explicit scheme to use extremely small steps.

8.2. Boundary-value Problems

In an initial-value problem all required conditions are specified at one initial point. In a boundary-value problem, conditions can be distributed across the domain:

y'' = f(x,y,y')
y(a) = alpha
y(b) = beta.

A shooting method guesses the missing initial slope and tries to satisfy the terminal boundary. A finite-difference method instead approximates derivatives on a mesh and solves a coupled algebraic system over the complete domain.

A second-order boundary-value problem is defined by both its differential equation and its boundary conditions. For y''(x)=0 with y(0)=2 and y(1)=5, the solution is y(x)=2+3x. Setting the discrete second difference to zero identifies a family of linear profiles, not a unique solution until endpoint constraints are imposed. Grid refinement should therefore assess boundary residuals and relevant physical invariants alongside the residual of the differential equation. A small interior residual alone does not establish a correct numerical solution.

8.3. Introduction to PDE Discretization

In a PDE, the unknown depends on more than one independent variable. Heat, diffusion, wave, and potential equations share a common numerical workflow:

continuous model
-> spatial/temporal mesh
-> discrete derivative approximation
-> algebraic system
-> solution
-> mesh refinement and error control.

Finite differences approximate derivatives from neighboring nodes. Finite volume balances conserved fluxes across cell faces. Finite element methods use a weak form with piecewise basis functions. Particle methods such as SPH need not rely on a fixed mesh.

The course introduces their common numerical structure rather than developing full FEM, FVM, or particle-mechanics formulations. From a computing perspective, the shared transition is from a continuous model to a large structured or sparse algebraic system.

In a conservation equation, flux crossing a shared cell face is an outflow for one cell and an inflow for its neighbor. A finite-volume discretization constructs the mass or energy balance from face fluxes; interior contributions cancel with opposite signs when the cells are aggregated. Transferring 0.3 units between two cells with impermeable external boundaries changes their individual states but not the total conserved quantity. This invariant supplies a verification criterion distinct from pointwise derivative error. Computing inconsistent fluxes independently on each side of a shared face can break conservation.

8.4. Symbolic, Numerical, and Automatic Differentiation

Three derivative-computation strategies should be distinguished:

• symbolic differentiation transforms an expression into another expression,
• numerical differentiation estimates derivatives from function values,
• automatic differentiation applies the chain rule to elementary program operations and produces derivatives close to machine precision.

Forward mode propagates an input direction toward outputs; reverse mode propagates an output sensitivity back to many inputs. Rather than materializing a large Jacobian or Hessian, software can compute JVPs, VJPs, or Hessian-vector products.

Known sparsity in a Jacobian or Hessian provides another optimization by avoiding structurally zero blocks. Mathematics, computational graphs, and data structures therefore meet in the same engineering problem.

Automatic differentiation is not the same as finite-difference numerical differentiation. For f(x)=x²+sin(x), local derivatives of the square and sine operations can be propagated through the computation graph by the chain rule, subject to machine-arithmetic limitations. Finite differences instead trade truncation error against cancellation and rounding as the step size shrinks. Reverse mode is often efficient for a scalar output with many inputs, while forward mode can be attractive for one input and many outputs; the choice depends on the graph's dimensions and computational structure.

8.5. The Role of Numerical Analysis in AI Computation

AI models may be defined over real numbers, but computers execute finite-precision arithmetic and approximate algorithms. Numerical analysis therefore remains directly relevant. Larger models do not remove numerical error; they can amplify its visibility through long reduction chains and wide dynamic ranges.

8.6. Floating Point and Accumulated Error

Simple decimal values may not be represented exactly in binary floating point. Small rounding errors can accumulate across millions of operations or become amplified by ill-conditioned computations. Matrix products, reductions, and normalizations in training and inference are all affected by data type and operation order.

FP16, BF16, and integer quantization can improve performance, but reduce precision. Mixed-precision execution is a modern form of classical error management: some operations use narrow formats while critical accumulations use wider ones.

8.7. Conditioning and Stability

An ill-conditioned problem can turn small input perturbations into large output changes. A stable algorithm cannot remove that inherent sensitivity. Conversely, a well-conditioned problem can still be solved poorly by an unstable algorithm.

For least squares, forming normal equations:

(A^T A)x = A^T b

can square the condition number. QR or SVD methods can be safer. Even elementary linear regression therefore has a numerical-method dimension.

8.8. Optimization Is a Numerical Iteration

Training commonly has the form:

θ_(k+1) = θ_k + Δθ_k

Step size, stopping criteria, gradient scale, and tolerances affect the result. Overflow, underflow, or unstable updates can occur even while the high-level model is correctly specified.

8.9. Stable Evaluation of Nonlinear Expressions

Products of small probabilities can underflow, motivating log-likelihood:

log Π p_i = Σ log p_i

Softmax can overflow if exponentials are evaluated directly. Subtracting the maximum leaves the mathematical result unchanged while improving numerical stability:

softmax(z_i) = exp(z_i-m) / Σ exp(z_j-m)

with m=max(z).

Direct multiplication of many small likelihood factors can underflow in floating-point arithmetic. Summing log likelihoods avoids repeated tiny products; for log-sum-exp, choose m = max(aᵢ) and compute log(Σexp(aᵢ)) = m + log(Σexp(aᵢ − m)). This evaluates the same expression with substantially reduced overflow risk.

In a Laplace approximation, an ill-conditioned Hessian near the mode can make the local quadratic approximation unreliable. A converged variational optimizer may still underrepresent posterior uncertainty because of its restricted approximation family. For MCMC, strong autocorrelation reduces effective sample size; a large draw count alone does not demonstrate convergence. Numerical correctness and adequacy of the statistical model are distinct concerns.

8.10. Approximation Beyond Floating Point

Approximate nearest-neighbor search, low-rank approximation, quantization, and pruning intentionally trade accuracy for lower cost. “Approximate” should therefore be quantified with recall, error norms, task metrics, and latency rather than treated as a qualitative label.

Parallel reductions also affect reproducibility because floating-point addition is not associative at the bit level:

(a+b)+c ≠ a+(b+c)

Small differences across hardware or execution order are not automatically defects; acceptable tolerance must be defined by the task.

Numerical analysis does not decide what an AI model should learn. It decides how reliably the intended computation survives finite precision, approximation, and iteration.

8.11. Proximal Gradients and Split Optimization

For a composite objective F(x) = f(x) + g(x) with a differentiable loss and a nonsmooth regularizer, ordinary gradient descent need not be defined for g. If the gradient of f is L-Lipschitz continuous and g is an appropriate closed convex function, a proximal-gradient step is

x_(k+1) = prox_(αg)(x_k − α ∇f(x_k))
prox_(αg)(v) = argmin_z [g(z) + ||z−v||₂² / (2α)]

The fixed step α ≤ 1/L is a common conservative choice under the stated assumptions. For g(x) = λ‖x‖₁, the proximal map becomes componentwise soft-thresholding; sparsity in lasso then arises directly from the update. Convergence rates must be interpreted in relation to smoothness and, where applicable, strong convexity rather than iteration count alone.

ADMM splits min_(x,z) [f(x) + g(z)] subject to Ax + Bz = c into separate subproblems. With the scaled dual variable u, the method alternates updates of x and z and then updates u to reflect constraint mismatch. The primal residual r_k = Ax_k + Bz_k − c measures constraint violation; a dual residual measures changes associated with dual feasibility. Small residuals are meaningful only relative to suitable stopping tolerances and problem scaling. Convergence results available for convex formulations cannot automatically be extended to arbitrary nonconvex problems.

Alternating minimization and majorization–minimization (MM) replace difficult objectives with simpler subproblems during each iteration. Solving those subproblems well does not prove that the overall global optimum has been reached. A numerical report should state the objective value, constraint violations, a gradient or other stationarity measure, the step size and the computational budget.

Low-precision arithmetic introduces two distinct error mechanisms: rounding numbers that cannot be represented exactly and overflow or underflow during intermediate operations. Adding many tiny increments one at a time to a large sum can yield a different result from first accumulating the increments separately; floating-point addition is not generally associative. FP16 and BF16 must be compared by their exponent and significand layouts, not merely by storage size. Numerical verification should establish a reference-precision result, acceptable absolute or relative error and representative boundary inputs before choosing an optimized computation mode.

8.12. Gauss–Newton and Damping in Nonlinear Least Squares

For F(θ)=||r(θ)||²/2, the residual Jacobian J gives gradient Jᵀr and an approximate Hessian JᵀJ. A Gauss–Newton step solves (JᵀJ)Δ=-Jᵀr; QR or appropriate Cholesky factorization is preferable to explicitly computing an inverse. Large residuals and nearly dependent Jacobian columns may make normal equations poorly conditioned. A Levenberg–Marquardt system (JᵀJ+λD)Δ=-Jᵀr uses damping to control unreliable steps. This is a local curvature-based step, not another name for gradient computation.

For the scalar residual r(θ)=θ²-4 at θ=3, the residual equals 5, the Jacobian equals 6, and the Gauss–Newton step is Δ=-30/36=-5/6. The next value is 13/6, yielding residual approximately 0.6944. One step illustrates local linearization, not guaranteed global convergence.

References

• Ahmet Yesevi Üniversitesi Bilgisayar Mühendisliği Bölümü. Sayısal Çözümleme (TBIL301) ders materyalleri.

• Chao Dong, Chen Change Loy, Kaiming He, Xiaoou Tang. “Image Super-Resolution Using Deep Convolutional Networks.” IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2), 2016, 295–307. doi:10.1109/TPAMI.2015.2439281

• Don P. Mitchell, Arun N. Netravali. “Reconstruction Filters in Computer Graphics.” ACM SIGGRAPH Computer Graphics, 22, 1988, 221–228. doi:10.1145/378456.378514

• Gene H. Golub, Charles F. Van Loan. Matrix Computations, 4th Edition. Johns Hopkins University Press, 2013.

• Göran Kauermann; Helmut Küchenhoff; Christian Heumann. Statistical Foundations, Reasoning and Inference: For Science and Data Science. Springer, 2021. https://doi.org/10.1007/978-3-030-69827-0
• Lloyd N. Trefethen, David Bau III. Numerical Linear Algebra. SIAM, 1997.

• Mor Harchol-Balter. Performance Modeling and Design of Computer Systems: Queueing Theory in Action. Cambridge University Press, 2013.

• Muhammet Ali Köker. dct-perceptual-hash-csharp (research software). Zenodo, 2026. doi:10.5281/zenodo.22117445

• Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms, 2nd ed. SIAM, 2002.

• Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms, 2nd Edition. SIAM, 2002.

• Paul S. Heckbert. Fundamentals of Texture Mapping and Image Warping. Technical Report UCB/CSD-89-516, University of California, Berkeley, 1989. eecs.berkeley.edu

• Paulius Micikevicius et al. “Mixed Precision Training.” International Conference on Learning Representations, 2018.

• Peter J. Huber; Elvezio M. Ronchetti. Robust Statistics. 2nd ed. Wiley, 2009.

• Richard L. Burden, J. Douglas Faires, Annette M. Burden. Numerical Analysis, 10th Edition. Cengage, 2015.

• Robert G. Keys. “Cubic Convolution Interpolation for Digital Image Processing.” IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(6), 1981, 1153–1160. doi:10.1109/TASSP.1981.1163711

• Steven C. Chapra, Raymond P. Canale. Numerical Methods for Engineers, 7th Edition. McGraw-Hill, 2015.

• Walter W. Piegorsch, Richard A. Levine, Hao Helen Zhang, Thomas C. M. Lee (eds.). Computational Statistics in Data Science. Wiley, 2022. ISBN 9781119561071.

• 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.
• Necsulescu, D. Advanced Mechatronics: Monitoring and Control of Spatially Distributed Systems. World Scientific, 2009. ISBN 978-981-277-181-0.
• Chunwei Zhang, Tianpeng Li, Ying Dai, Li Sun, and Ardashir Mohammadzadeh. Practical Neural Networks in Python and MATLAB. Springer Nature, 2026. DOI: https://doi.org/10.1007/978-3-032-14746-2.