Chapter 1: Introduction

What are numerical methods?

Numerical methods approximate solutions to mathematical problems that are analytically intractable or prohibitively expensive to solve exactly. They translate mathematical models into algorithms that can be implemented on a computer to generate reliable, repeatable results. Numerical methods are widely used in engineering, physics, finance, and computer science to address differential equations, optimization, integration, and more.

Learning objectives

Identify when numerical methods are required and why.
Understand the sources of error and how they propagate.
Recognize the governing equations that motivate the numerical models used in this course.

Importance of numerical methods

Modern engineering relies heavily on mathematical models to describe physical, chemical, and biological phenomena. These models are most often expressed in the form of algebraic equations, ordinary differential equations, or partial differential equations. In realistic settings, such equations are rarely solvable in closed form, either because of their nonlinearity, the complexity of the geometry, or the coupling of multiple physical processes.

Numerical methods provide a systematic way to approximate the solutions of these models by transforming continuous mathematical problems into finite-dimensional problems that can be solved on a computer. Rather than producing exact formulas, numerical methods aim to deliver approximate solutions whose accuracy, stability, and reliability can be quantified and controlled.

The importance of numerical methods stems from several key aspects:

Complexity of real-world problems

Engineering systems often involve multiple spatial and temporal scales, nonlinear material behavior, and complex boundary conditions. Numerical methods make it possible to address such complexity in a systematic and reproducible manner.
Predictive simulation

Numerical simulations allow engineers to predict the behavior of systems before they are built or tested, reducing development costs and enabling design optimization. In many applications, numerical simulation is the only feasible tool for exploring extreme or inaccessible conditions.
Control of approximation errors

Every numerical approximation introduces errors, originating from modeling assumptions, discretization, and finite-precision arithmetic. A central objective of numerical analysis is to understand how these errors arise, how they propagate, and how they can be reduced to acceptable levels.
Algorithmic efficiency and scalability

Practical engineering problems often lead to very large computational systems. Numerical methods must therefore be designed not only for accuracy, but also for efficiency and scalability on modern computing architectures.
Bridging theory and computation

Numerical methods form the link between mathematical theory and practical computation. They provide the tools required to turn governing equations into algorithms whose behavior can be analyzed, implemented, and validated.

Physical models

Numerical methods are most often developed and assessed in the context of mathematical models describing physical phenomena. In this course, we focus on models arising in fluid mechanics and heat transfer, as they provide representative examples of partial differential equations.

The purpose of this section is not to derive these models in detail, but to:

introduce the governing equations that motivate the numerical problems studied later,
highlight the mathematical structure of these equations,
and identify the features that make their numerical approximation challenging.

We consider incompressible, viscous, and heat-conducting fluids occupying a spatial domain \(\Omega \subset \mathbb R^d\), with \(d \in \{1,2,3\}\).

Notation

Throughout this chapter, the following conventions are used:

Underlined symbols (e.g. \(\underline v\)) denote vectors,
Double-underlined symbols (e.g. \(\underline{\underline{s}}\)) denote second-order tensors,
\(\operatorname{grad}\) and \(\operatorname{div}\) denote the gradient and divergence operators, respectively.

1.1.1 Continuity equation (mass conservation)

The conservation of mass is expressed by the continuity equation, \[ \partial_t \rho + \operatorname{div} \left ( \rho \underline v \right ) = 0, \] where \(\rho\) denotes the mass density and \(\underline v\) the velocity field.

In the incompressible limit, where the density is constant, this equation reduces to the solenoidal constraint \[ \operatorname{div} \underline v = 0. \]

This constraint is equivalent to the statement that the fluid volume is preserved during the flow. It plays a central role in the numerical approximation of incompressible flows, as it couples the velocity field to the pressure.

1.1.2 Navier–Stokes equations (momentum conservation)

The conservation of linear momentum expresses Newton’s second law at the continuum level, \[ \partial_t \left ( \rho \underline v \right ) + \operatorname{div} \left ( \rho \underline v \underline v \right ) = \operatorname{div} \underline{\underline{s}}, \] where \(\underline{\underline{s}}\) denotes the Cauchy stress tensor.

The stress tensor is decomposed into an isotropic and a deviatoric part, \[ \underline{\underline{s}} = -P \, \underline{\underline{\delta}} + \underline{\underline{s}}^\circ, \] where:

\(P\) is the mechanical pressure,
\(\underline{\underline{\delta}}\) is the identity tensor,
\(\underline{\underline{s}}^\circ\) is the deviatoric stress tensor.

Substituting this decomposition into the momentum equation yields \[ \partial_t \left ( \rho \underline v \right ) + \operatorname{div} \left ( \rho \underline v \underline v \right ) + \operatorname{grad} P = \operatorname{div} \underline{\underline{s}}^\circ. \]

For a Newtonian fluid, linear non-equilibrium thermodynamics leads to the relations \[ P = p - \eta \, \operatorname{div} \underline v, \] and \[ \underline{\underline{s}}^\circ = 2 \mu \left [ \left ( \operatorname{grad} \underline v \right )_{\operatorname{sym}} - \frac{1}{d} \left ( \operatorname{div} \underline v \right ) \underline{\underline{\delta}} \right ], \] where \(\mu\) is the dynamic viscosity, \(\eta\) the bulk viscosity, and \(d\) the spatial dimension.

It is convenient to introduce the viscous stress tensor \[ \underline{\underline{\tau}} \equiv 2 \mu \left ( \operatorname{grad} \underline v \right )_{\operatorname{sym}} + \left ( \eta - \frac{2 \mu}{d} \right ) \left ( \operatorname{div} \underline v \right ) \underline{\underline{\delta}}, \] so that the momentum equation can be written in the compact form \[ \partial_t \left ( \rho \underline v \right ) + \operatorname{div} \left ( \rho \underline v \underline v \right ) + \operatorname{grad} p = \operatorname{div} \underline{\underline{\tau}}. \]

These equations form a nonlinear, time-dependent system whose numerical solution poses significant challenges, particularly in the incompressible limit.

1.1.3 Energy equation (heat transfer)

The conservation of total energy is written as \[ \partial_t \left ( \rho \varepsilon \right ) + \operatorname{div} \left [ \left ( \rho \varepsilon + P \right ) \underline v \right ] = \operatorname{div} \left ( \underline v \cdot \underline{\underline{s}}^\circ - \underline q \right ), \] where \(\varepsilon\) is the specific total energy and \(\underline q\) denotes the heat flux.

Using Fourier’s law, \[ \underline q = -\kappa \, \operatorname{grad} T, \] and separating kinetic and internal energy contributions, one obtains an equation for the internal energy \(e\), \[ \partial_t \left ( \rho e \right ) + \operatorname{div} \left ( \rho e \underline v + \underline q \right ) = \Phi - p \operatorname{div} \underline v, \tag{1}\] where the viscous dissipation term \(\Phi\) is given by \[ \Phi = \eta \left ( \operatorname{div} \underline v \right )^2 + 2 \mu \left [ \left ( \operatorname{grad} \underline v \right )_{\operatorname{sym}} - \frac{1}{d} \left ( \operatorname{div} \underline v \right ) \underline{\underline{\delta}} \right ] : \left [ \left ( \operatorname{grad} \underline v \right )_{\operatorname{sym}} - \frac{1}{d} \left ( \operatorname{div} \underline v \right ) \underline{\underline{\delta}} \right ]. \] By construction, \(\Phi\) is non-negative, reflecting irreversible dissipation.

1.1.4 Boundary conditions

To obtain a well-posed mathematical problem, the governing equations must be supplemented with appropriate boundary conditions.

For the velocity field, common choices include:

no-slip conditions at solid walls,
free-slip conditions,
inflow/outflow conditions at open boundaries.

For the temperature field, typical boundary conditions are:

prescribed temperature (Dirichlet),
prescribed heat flux (Neumann),
convective heat transfer (Robin).

The choice of boundary conditions has a strong influence on both the mathematical properties of the problem and the design of numerical methods.

Model scalar transport equation

The aforementioned systems of equations form a set of coupled, nonlinear partial differential equations whose direct numerical solution is challenging. To introduce and test numerical ideas in a controlled setting, we will frequently use a simpler model problem that retains the main mathematical features while remaining easier to analyze.

Starting from the internal energy equation (Equation 1), we neglect viscous dissipation and pressure work and assume the equation of state \[ \mathrm d _ t e = \rho c \, \mathrm d _ t T, \] valid for a perfect gas. \(\mathrm d _ t\) denotes the material derivative and \(c\) is the specific heat capacity. With Fourier’s law and constant coefficients, the temperature field satisfies the scalar advection-diffusion equation \[ \partial _ t T + \underline v \cdot \operatorname{grad} T = \alpha \, \Delta T + f, \] where \(f\) is a volumetric source term and \[ \alpha \equiv \frac{\kappa}{\rho c} \] is the thermal diffusivity. This equation combines hyperbolic (advection) and parabolic (diffusion) effects.

In this course we will focus on purely conductive phenomena, so we will set \(\underline v = \underline 0\) and study the resulting diffusion model in Cartesian coordinates (\(x\), \(y\), \(z\)). These simplified equations provide a clean test bed for discretization, stability analysis, and error estimation. By the end of the class, you will be able to numerically solve the following equations.

1.2.1 The heat equation

In one spatial dimension, the heat equation reads, \[ \partial _ t T = \alpha \, \partial _ {x x} T + f. \]

1.2.2 The Poisson equation

In steady-state conditions, the heat equation reduces to the Poisson equation, \[ 0 = \alpha \, \left ( \partial _ {xx} T + \partial _ {yy} T \right ) + f. \] (Written for two spatial dimensions.)