Mathematics for Engineering

J. Gordon Leishman

doi:https://doi.org/10.15394/eaglepub.2022.1066.n8

4 Mathematics for Engineering

Introduction

Omnia de mathematica agitur!^[1] Many will say that mathematics is the engineer’s language, so the first thing students need to do before they can study engineering is to learn (or relearn) the essential mathematics. Learning mathematics is a lifelong endeavor for most engineers and is essential to the profession. All rigorous engineering analysis ultimately rests on mathematical structure, even when the mathematics is implicit. Mathematics enables engineers to understand and predict the behavior of physical systems, deal with uncertainties, and find optimal solutions to complex problems.

Readers of this eBook will likely be already familiar with the essential mathematics required to enter various engineering and scientific disciplines. However, they will still need to review fundamental mathematical concepts to progress successfully. Geometry, algebra, trigonometry, calculus, and vectors provide engineers with essential mathematical tools for tracking processes and solving problems. Differential equations, for example, also arise in many aspects of engineering problem-solving. To this end, all engineers must become well-versed in solving various types of differential equations, analytically and numerically. In particular, when performing arithmetical and algebraic calculations, engineers must do so accurately and without error. Errores in engineeringo plerumque sunt catastrophes.^[2]

In engineering, quantities with both magnitude and direction are often required, such as forces, velocities, and accelerations, as well as their manipulations involving scalar and vector products. Vector quantities are most conveniently expressed using vector notation, a concise shorthand for the corresponding scalar equations. For example, the equations of motion for a fluid are often written as vector equations. Therefore, engineers must become comfortable with vector quantities, including shorthand notation for operators such as the gradient, the Laplace operator, and the substantial derivative. Exercitatio perfectos efficit.^[3]

Learning Objectives

Review the essential mathematical tools needed for aerospace engineering analysis.
Apply algebraic manipulation to rearrange and simplify engineering equations.
Review the basic rules of differentiation and integration and their use in engineering problems.
Use vectors to represent physical quantities with both magnitude and direction.
Interpret and apply primary vector operations, including the scalar product, cross-product, gradient, divergence, and curl.
Explain the meanings of line integrals, surface integrals, and volume integrals.
Describe the substantial derivative and the Laplace operator and recognize where they arise in engineering analysis.
Use complex numbers and recognize how quaternions extend similar ideas to three-dimensional rotations.
Recognize why ordinary and partial differential equations are needed to model engineering systems.

Algebra

Algebra is used extensively in engineering. It is the branch of mathematics that uses symbols to represent quantities and their relationships. These symbols are manipulated using the rules of arithmetic to solve equations, simplify expressions, and rearrange formulas into more useful forms. In engineering applications, algebra is not merely an abstract mathematical exercise; it is the means by which physical relationships are converted into useful numerical results.

Some of the most important algebraic operations include collecting like terms, expanding and factoring expressions, simplifying fractions, working correctly with exponents and radicals, and rearranging equations to solve for a specified variable. Care is needed when applying the same operation to both sides of an equation and when preserving signs, units, and dimensions throughout the calculation. For example, suppose an engineering relation has the form

(1) $\begin{equation*} F = C \left( \frac{V^2}{L+d} \right) + F_0 \end{equation*}$

where $F$ , $C$ , $L$ , $d$ , and $F_0$ are known quantities, and $V$ is a velocity to be determined. First, subtract $F_0$ from both sides to isolate the term containing $V$ , i.e.,

(2) $\begin{equation*} F - F_0 = C \left( \frac{V^2}{L+d} \right) \end{equation*}$

Next, divide by $C$ and multiply by $L+d$ , giving

(3) $\begin{equation*} V^2 = \frac{(F-F_0)(L+d)}{C} \end{equation*}$

Taking the square root then gives

(4) $\begin{equation*} V = \sqrt{\frac{(F-F_0)(L+d)}{C}} \end{equation*}$

provided that the quantity inside the square root is non-negative. This example also shows why units and signs must be tracked carefully. The combination $(F-F_0)(L+d)/C$ must have the dimensions of velocity squared, and the selected square-root branch must be consistent with the physical meaning of $V$ . Such rearrangements are common in engineering because equations are often derived first from physical principles, and only afterward rearranged to solve for the quantity of interest.

As a more demanding example, suppose another engineering relation has the form

(5) $\begin{equation*} T = T_\infty + (T_0-T_\infty)e^{-k\,t} \end{equation*}$

where ${T}$ , ${T_\infty}$ , and ${T_0}$ are temperatures, $k$ is a known positive constant, and $t$ is to be determined. First, subtract $T_\infty$ from both sides, i.e.,

(6) $\begin{equation*} T-T_\infty = (T_0-T_\infty)e^{-k\,t} \end{equation*}$

Next, divide by $T_0-T_\infty$ to isolate the exponential term to get

(7) $\begin{equation*} \frac{T-T_\infty}{T_0-T_\infty} = e^{-k\,t} \end{equation*}$

Taking the natural logarithm of both sides gives

(8) $\begin{equation*} \ln\left(\frac{T-T_\infty}{T_0-T_\infty}\right) = -k\,t \end{equation*}$

and solving for $t$ gives

(9) $\begin{equation*} t = -\frac{1}{k}\ln\left(\frac{T-T_\infty}{T_0-T_\infty}\right) \end{equation*}$

This result is valid only when the argument of the logarithm is positive. In addition, the result must be checked against the problem’s physical meaning. Because $k>0$ , the value of $t$ is non-negative only when the temperature difference has decreased in magnitude from its initial value.

A sound working knowledge of algebra is essential before progressing to more advanced engineering mathematics. Many errors in engineering calculations arise from incorrect algebraic manipulations, sign errors, inconsistent units, and the substitution of numerical values too early in the solution process. For this reason, algebraic fluency remains one of the most important foundations for reliable engineering problem-solving.

Check Your Understanding #1 – Review problems in algebra

Simplify the expression given by $2x^2 - 3x + 1 - (x^2 + 2x - 5)$ .
Solve for $x$ in the equation $2(x + 3) = 8 - x$ .
Simplify the expression $5(2x + 3y) - 2(3x - y)$ .
Solve for $x$ in the equation $\dfrac{x^2 + 3x - 10}{x + 2} = 0$ .
Solve for $U$ in the equation $P = P_0 + K U^3$ .
Solve for $h$ in the equation $p = p_0 e^{-h/H}$ .

Show solution/hide solution.

Removing parentheses and combining like terms gives
$x^2 - 5x + 6$
Distributing the 2 gives
$2x + 6 = 8 - x$

Combining like terms gives

$3x + 6 = 8$

Subtracting 6 from both sides gives

$3x = 2$

Finally, dividing both sides by 3 gives

$x = \frac{2}{3}$
Distributing the coefficients gives
$10x + 15y - 6x + 2y$

Combining like terms gives

$4x + 17y$
The expression is undefined when ${x = -2}$ because it would result in division by zero. For this reason, ${x = -2}$ is excluded from the domain and cannot be a solution. Factoring the numerator gives
$\frac{(x + 5)(x - 2)}{x + 2} = 0$

The fraction is zero only when the numerator is zero and the denominator is not zero. Setting the numerator factors equal to zero gives

$x + 5 = 0 \quad \mbox{or} \quad x - 2 = 0$

Solving each equation separately gives

$x = -5 \quad \mbox{or} \quad x = 2$
Subtracting $P_0$ from both sides gives
$P-P_0 = K U^3$

Dividing by $K$ gives

$U^3 = \frac{P-P_0}{K}$

Taking the cube root gives

$U = \left(\frac{P-P_0}{K}\right)^{1/3}$

For a physical speed, the result must be consistent with the allowable range of $U$ .
Dividing both sides by $p_0$ gives
$\frac{p}{p_0} = e^{-h/H}$

Taking the natural logarithm of both sides gives

$\ln\left(\frac{p}{p_0}\right) = -\frac{h}{H}$

Multiplying both sides by $-H$ gives

$h = -H \ln\left(\frac{p}{p_0}\right)$

This result is valid when the argument of the logarithm is positive.

Calculus

Engineering requires calculus, and to use it, one must know how to differentiate and integrate. Both processes have rules. Success with calculus means first learning the rules, i.e., the mechanical processes usually taught to engineers by the mathematics department. But to succeed with calculus in engineering means much more, in that the process of doing the mathematics must also be physically interpreted. For example, differentiation is akin to finding a rate of change, a slope, or a gradient, while integration involves finding a sum, such as an area, a volume, or the contribution of something that acts or is distributed along a particular path or over a particular area. The tangible application of mathematics within the engineering framework is an acquired skill, but one essential for success as an engineer.

Ordinary Differentiation

Finding the change in one thing with respect to another is the essence of differentiation. For a function of one variable, such as $f(x)$ , the derivative measures how rapidly the function changes as $x$ changes. It is defined as the limiting value of a ratio of small changes, i.e.,

(10) $\begin{equation*} \frac{df}{dx} \equiv f^{\prime} = \lim_{\Delta x \to 0} \frac{f(x+\Delta x)-f(x)}{\Delta x} = \lim_{\Delta x \to 0} \frac{\Delta f}{\Delta x} \end{equation*}$

This derivative is the local slope of the function at a point, as shown in Figure 1. It tells how much $\scriptstyle{f}$ changes per unit change in $x$ , in the limit as the change in $x$ approaches zero. In engineering, this idea appears constantly. For example, velocity is the derivative of position with respect to time, acceleration is the derivative of velocity with respect to time, and a pressure gradient is the derivative of pressure with respect to distance.

The derivative is the local slope of the function at a point in the limit as the change in $x$ approaches zero.

The fundamental rules of differentiation are:

1. Power Rule:

(11) $\begin{equation*} \frac{d}{dx}(x^n) = n x^{n-1}~\mbox{~where $n$ is a constant.} \end{equation*}$

2. Constant Rule:

(12) $\begin{equation*} \frac{d}{dx}(c) = 0~\mbox{~where ${c}$ is a constant.} \end{equation*}$

3. Sum/Difference Rule:

(13) $\begin{equation*} { \frac{d}{dx}(f(x) \pm g(x)) = \frac{d}{dx} (f(x)) \pm \frac{d}{dx} (g(x))~\mbox{~where ${\scriptstyle{f}}$ and $g$ are functions.} } \end{equation*}$

4. Product Rule:

(14) $\begin{equation*} \frac{d}{dx}(f(x)g(x)) = f'(x) \, g(x) + f(x) \, g'(x)~\mbox{~where ${\scriptstyle{f}}$ and $g$ are functions.} \end{equation*}$

5. Quotient Rule:

(15) $\begin{equation*} \frac{d}{dx}\left(\frac{f(x)}{g(x)}\right) = \frac{f'(x) \, g(x) - f(x) \, g'(x)}{(g(x))^2}~\mbox{~where ${\scriptstyle{f}}$ and $g$ are functions.} \end{equation*}$

6. Chain Rule:

(16) $\begin{equation*} \frac{d}{dx}(f(g(x))) = f'(g(x)) \, g'(x)~\mbox{~where ${\scriptstyle{f}}$ and $g$ are functions.} \end{equation*}$

Having laid down these basic rules, which are the first things to learn, they are necessary but insufficient for engineering. Interpreting and understanding the meaning of these processes is more challenging but essential to achieving success in engineering problem-solving. Experienced engineers take these “rules” for granted and do the mathematics valor facialis.^[4]

Check Your Understanding #2 – Review problems using differentiation

Find the derivative of the function $f(x) = 3x^2 + 2x - 1$
Find the derivative of the function $g(x) = \sin 2x + \cos x$
Find the derivative of the function $h(x) = \dfrac{2}{x^3}$
Find the derivative of the function $y = \ln(3x^2 + 1)$
Find the derivative of the function $f(x) = \sqrt{x^3 + 2x}$

Show solution/hide solution.

By taking the derivative term by term, then

$f'(x) = \frac{df}{dx} = 6x + 2$

2. Applying the chain rule leads to

$g'(x) = \frac{dg}{dx} = 2\cos(2x) - \sin x$

3. Using the power rule gives

$h'(x) = \frac{dh}{dx} = -6x^{-4}$

4. Applying the chain rule gives

$y'(x) = \frac{dy}{dx} = \frac{6x}{3x^2 + 1}$

5. Using the power rule and the chain rule leads to

$f'(x) = \frac{df}{dx} = \frac{3x^2 + 2}{2\sqrt{x^3 + 2x}}$

Partial Differentiation

Many engineering quantities depend on more than one independent variable. For example, the temperature in a room may vary with position and time so that it can be written as $T(x,y,z,t)$ . Similarly, the pressure in a fluid may vary with all three spatial coordinates and time so that it can be written as $p(x,y,z,t)$ . In such cases, the question “how fast does the quantity change?” is incomplete unless the variable with respect to which the change is being measured is specified. A partial derivative answers this question by measuring the rate of change with respect to one independent variable while all other independent variables are held fixed.

For a function of two variables, $f(x,y)$ , the partial derivative of $\scriptstyle{f}$ with respect to $x$ is defined as

(17) $\begin{equation*} \frac{\partial f}{\partial x} = \lim_{\Delta x \to 0} \frac{f(x+\Delta x,y)-f(x,y)}{\Delta x} \end{equation*}$

In this operation, $y$ is treated as a constant. Therefore, $\partial f / \partial x$ measures the local rate of change of $\scriptstyle{f}$ in the $x$ direction only. Similarly, the partial derivative of $\scriptstyle{f}$ with respect to $y$ is defined as

(18) $\begin{equation*} \frac{\partial f}{\partial y} = \lim_{\Delta y \to 0} \frac{f(x,y+\Delta y)-f(x,y)}{\Delta y} \end{equation*}$

where $x$ is held constant. Therefore, $\partial f / \partial y$ measures the local rate of change of $\scriptstyle{f}$ in the $y$ direction only.

In practice, partial derivatives are computed by differentiating with respect to the variable of interest while treating all other independent variables as constants. For example, if

(19) $\begin{equation*} f(x,y) = x^2 y + 3y^2 \end{equation*}$

then the partial derivative with respect to ${x}$ is found by treating ${y}$ as a constant. This gives

(20) $\begin{equation*} \frac{\partial f}{\partial x} = \frac{\partial}{\partial x}\left(x^2 y + 3y^2\right) = y\frac{\partial}{\partial x}\left(x^2\right) + 3y^2\frac{\partial}{\partial x}\left(1\right) \end{equation*}$

Because ${y}$ is being held fixed, both ${y}$ and ${3y^2}$ are constants with respect to ${x}$ . It follows that

(21) $\begin{equation*} \frac{\partial f}{\partial x} = 2xy \end{equation*}$

The partial derivative with respect to ${y}$ is found by treating ${x}$ as a constant. This gives

(22) $\begin{equation*} \frac{\partial f}{\partial y} = \frac{\partial}{\partial y}\left(x^2 y + 3y^2\right) = x^2\frac{\partial}{\partial y}\left(y\right) + 3\frac{\partial}{\partial y}\left(y^2\right) \end{equation*}$

Because ${x}$ is being held fixed, ${x^2}$ is a constant with respect to ${y}$ . Hence,

(23) $\begin{equation*} \frac{\partial f}{\partial y} = x^2 + 6y \end{equation*}$

The same function can have different partial derivatives, depending on which independent variable is being varied.

Partial differentiation is a fundamental tool in engineering because most physical quantities depend on several variables. As will become apparent in later chapters of this eBook, the use of partial derivatives is the norm rather than the exception. They are used to describe spatial gradients, time rates of change, conservation laws, fluid motion, heat transfer, structural deformation, electric and magnetic fields, and many other physical systems.

Check Your Understanding #3 – Review problems using partial derivatives

For each of the following functions, compute $\dfrac{\partial f}{\partial x}$ and $\dfrac{\partial f}{\partial y}$ :

$f(x,y)=2x^{2}+3xy+y^{2}$
$f(x,y)=4x^{3}y-y^{4}+5x$
$f(x,y)=\ln \left(x^{2}+y^{2}\right)$
$f(x,y)=e^{xy}\sin y$
$f(x,y)=\dfrac{x}{y}+y\cos x$

Show solution/hide solution.

$f(x,y)=2x^{2}+3xy+y^{2}$ . First, differentiate the function with respect to $x$ while treating $y$ as a constant, giving
$\frac{\partial f}{\partial x}=4x+3y$

Next, differentiate the function with respect to $y$ while treating $x$ as a constant, giving

$\frac{\partial f}{\partial y}=3x+2y$
$f(x,y)=4x^{3}y-y^{4}+5x$ . First, differentiate with respect to $x$ (treat $y$ as constant), i.e.,
$\frac{\partial f}{\partial x}=12x^{2}y+5$

Next, differentiate with respect to $y$ (treat $x$ as constant), i.e.,

$\frac{\partial f}{\partial y}=4x^{3}-4y^{3}$
$f(x,y)=\ln\left(x^{2}+y^{2}\right)$ . First, differentiate with respect to $x$ (use the chain rule), i.e.,
$\frac{\partial f}{\partial x} =\frac{1}{x^{2}+y^{2}}\, 2x =\frac{2x}{x^{2}+y^{2}}$

Next, differentiate with respect to $y$ (again using the chain rule), i.e.,

$\frac{\partial f}{\partial y} =\frac{1}{x^{2}+y^{2}}\, 2y =\frac{2y}{x^{2}+y^{2}}$
$f(x,y)=e^{xy}\sin y$ . First, differentiate with respect to $x$ (treat $y$ as constant and apply the chain rule to $e^{xy}$ ), i.e.,
$\frac{\partial f}{\partial x} =\left(y\,e^{xy}\right)\sin y =y\,e^{xy}\sin y$

Next, differentiate with respect to $y$ (product‐and chain‐rule on $e^{xy}\sin y$ ), i.e.,

$\frac{\partial f}{\partial y} =\left(xe^{xy}\right)\sin y + e^{xy}\cos y =xe^{xy}\sin y + e^{xy}\cos y$
$f(x,y)=\dfrac{x}{y}+y\cos x$ . First, differentiate with respect to $x$ (treat $y$ as constant), i.e.,
$\frac{\partial f}{\partial x} =\frac{1}{y} - y\sin x$

Next, differentiate with respect to $y$ (treat $x$ as constant), i.e.,

$\frac{\partial f}{\partial y} =-\frac{x}{y^{2}} + \cos x$

Integration

Integration is a process often described as the inverse of differentiation. Integration is more complex than differentiation, but it comes up in some form in most engineering problem-solving. Integration techniques involve using the power rule, trigonometric and logarithmic integrals, and the substitution method. Most engineers only recall some of the rules of integration; after all, they aspire to be more than just pure mathematicians. Still, they can look them up in books as standard integrals or use computational tools such as Wolfram Alpha or ChatGPT, provided the result is also checked against the mathematical and physical context. The successful use of calculus should never depend on guesswork or AI. Labora callide, non fortiter.^[5]

The fundamental rules of integration are:

1. Power Rule:

(24) $\begin{equation*} \int x^n \, dx = \frac{1}{n+1}x^{n+1} + c~\mbox{~where $n$ is $\ne$ $-1$, and ${c}$ is a constant.} \end{equation*}$

2. Constant Multiple Rule:

(25) $\begin{equation*} \int cf(x) \, dx = c\int f(x) \, dx~\mbox{~where ${c}$ is a constant.} \end{equation*}$

3. Sum/Difference Rule:

(26) $\begin{equation*} \int (f(x) \pm g(x)) \, dx = \int f(x) \, dx \pm \int g(x) \, dx~\mbox{~where ${f}$ and ${g}$ are functions.} \end{equation*}$

4. Integration by Parts:

(27) $\begin{equation*} \int u \, dv = uv - \int v \, du~\mbox{~where ${u}$ and ${v}$ are functions.} \end{equation*}$

5. Substitution Rule:

(28) $\begin{equation*} \int f( g(x))\, g'(x) \, dx = \int f(u) \, du~\mbox{~where $u = g(x)$.} \end{equation*}$

Again, the rules are the first things to learn, but the meaning of the integration processes must also be understood in engineering. For example, integration also has an important geometric interpretation. For a function $f(x)$ , the definite integral

(29) $\begin{equation*} \int_a^b f(x)\,dx \end{equation*}$

represents the signed area under the curve between $x=a$ and $x=b$ , as shown in Figure 2. In engineering applications, this idea is often more important than the symbolic rule used to evaluate the integral. For example, integrating a velocity-time curve gives displacement, integrating a force-displacement curve gives work, and integrating a distributed load gives the resultant force.

A definite integral represents the signed area under a curve between two limits. In engineering, this area often corresponds to a physical quantity such as displacement, work, impulse, or resultant force.

Check Your Understanding #4 – Review problems using integration

Evaluate the integral $\displaystyle{ \int (3x^2 + 2x - 1)} \, dx$
Evaluate the integral $\displaystyle{\int (2\cos(2x) - \sin x) } \, dx$
Evaluate the integral $\displaystyle{\int -6x^{-4}} \, dx$
Evaluate the integral $\displaystyle{ \int \dfrac{6x}{3x^2 + 1} }\, dx$
Evaluate the integral $\displaystyle{ \int \dfrac{3x^2 + 2}{2\sqrt{x^3 + 2x}} } \, dx$

Show solution/hide solution.

1. Integrating term by term gives

$\int (3x^2 + 2x - 1)\,dx = \frac{3}{3}x^3 + \frac{2}{2}x^2 - x + C = x^3 + x^2 - x + C$

where $C$ is the constant of integration.

2. Integrating each term separately leads to

$\int (2\cos(2x)-\sin x)\,dx = 2\int \cos(2x)\,dx - \int \sin x\,dx$

$= 2\left(\frac{1}{2}\sin(2x)\right) - \left(-\cos x\right) + C$

$= \sin(2x) + \cos x + C$

where $C$ is the constant of integration.

3. Applying the power rule of integration gives

$\int -6x^{-4}\,dx = 2x^{-3} + C$

where $C$ is the constant of integration.

4. Using the substitution $u = 3x^2 + 1$ , then $du = 6x\,dx$ and the integral becomes

$\int \frac{6x}{3x^2 + 1}\,dx = \int \frac{du}{u} = \ln|u| + C$

By substituting back $u = 3x^2 + 1$ , the final result is

$\int \frac{6x}{3x^2 + 1}\,dx = \ln\!\bigl|\,3x^2 + 1\,\bigr| + C$

5. The integral to be solved is

$\int \frac{3x^{2}+2}{2\sqrt{x^{3}+2x}} \, dx$

Let $u = x^{3}+2x$ so $du = (3x^{2}+ 2)\,dx$ . Therefore,

$\int \frac{3x^{2}+2}{2\sqrt{x^{3}+2x}} \, dx = \frac12 \int u^{-1/2}\, du = u^{1/2} + C = \sqrt{x^{3}+2x} + C$

where $C$ is the constant of integration.

Working with Vectors

A vector is a quantity that possesses both magnitude and direction, such as force, velocity, or acceleration. Vectors come up everywhere in engineering. Si vectores nescis, tunc errabis. Vectors are represented by arrows with a length proportional to the magnitude, and the orientation denotes the direction. Vectors are also used to specify a position in space; such vectors are called position vectors. The location of an arbitrary point in space, say P, can be defined by specifying the values of three coordinates, i.e., in terms of $(x, y, z)$ in a standard Cartesian coordinate system. As shown in Figure 3, the point P can be located by the position vector, $\vec{r}$ , where

(30) $\begin{equation*} \vec{r} = x \, \vec{i} + y \, \vec{j} + z \, \vec{k} \end{equation*}$

Definition of a position vector $\vec{r}$ in a Cartesian coordinate system.

In general, if $\vec{A}$ is a given vector in a Cartesian coordinate system and $A_x, A_y$ , and $A_z$ are the scalar components of $\vec{A}$ in the $x$ , $y$ , and $z$ directions, then

(31) $\begin{equation*} \vec{A} = A_{x}\, \vec{i} + A_{y}\, \vec{j} + A_{z}\, \vec{k} \end{equation*}$

as shown in Figure 4 below.

The position of a point in space is defined by its scalar components in a Cartesian coordinate system.

Vectors can represent velocities (a speed in a given direction), forces, accelerations, or other quantities. Vectors can be added to find the total effect (such as total velocity, force, or acceleration) by adding components component-by-component, nose-to-tail, as shown in Figure 5. If $\vec{A}$ is directed from point $P_1$ to point $P_2$ and a second vector $\vec{B}$ is defined as

(32) $\begin{equation*} \vec{B} = B_{x}\, \vec{i} + B_{y}\, \vec{j} + B_{z}\, \vec{k} \end{equation*}$

that points from $P_2$ to $P_3$ , then the resultant vector $\vec{C}$ that points from $P_1$ to $P_3$ is

(33) $\begin{equation*} \vec{C} = (A_{x} + B_{x}) \, \vec{i} + (A_{y} + B_{y}) \, \vec{j} + (A_{z} + B_{z})\, \vec{k} \end{equation*}$

The process of adding vectors involves connecting them “nose-to-tail.”

Alternative coordinate systems can be used to describe any problem; for example, cylindrical or spherical coordinates can be used instead of Cartesian coordinates. Typically, selecting a suitable coordinate system is a matter of convenience. Still, it is also essential to determine whether the mathematics can be appropriately simplified in an alternative coordinate system. For example, in introductory engineering problems, Cartesian coordinates are used primarily, but polar coordinates are sometimes used when convenient. It is best to start with Cartesian coordinates to understand mathematical rules.

Scalar & Vector Fields

A field is a quantity that is defined at every point in space, and in general varies continuously from point to point. If physical systems are described in terms of such fields, the natural mathematical problem is to determine how these quantities change with position and with time. This problem leads directly to calculus. A simple example is a temperature distribution in a room: each point in space is associated with a specific temperature. A flow provides a second example: at each point, the air is characterized by a velocity with both magnitude and direction. In this description, attention is focused not on the motion of individual fluid elements, but on the value of the field at a specified point in space.

A scalar quantity given as a function of coordinate space $x, y, z$ (and perhaps time $t$ ) is called a scalar field. For example, fluid properties such as pressure $p$ , density $\varrho$ , and temperature $T$ are all scalar quantities, i.e.,

(34) $\begin{eqnarray*} p & = & p\left( x, y,z \right) \\[6pt] \varrho & = & \varrho\left( x, y, z \right) \\[6pt] T & = & T\left( x, y, z \right) \end{eqnarray*}$

Similarly, a vector quantity defined as a function of position and time is called a vector field. For example, velocity is a vector quantity. A velocity can be written in terms of its scalar components, i.e.,

(35) $\begin{equation*} \vec{V} = V_{x} \, \vec{i} + V_{y} \, \vec{j} + V_{z} \, \vec{k} \end{equation*}$

where the components are

(36) $\begin{eqnarray*} V_{x} & = & V_{x} \left( x, y, z \right) \\[6pt] V_{y} & = & V_{y} \left( x, y, z \right) \\[6pt] V_{z} & = & V_{z} \left( x, y, z \right) \end{eqnarray*}$

More often than not, the velocity vector is written as

(37) $\begin{equation*} \vec{V} = u\, \vec{i} + v \, \vec{j} + w \, \vec{k} \end{equation*}$

where ${u}$ , ${v}$ and ${w}$ are the components in the ${x}$ , ${y}$ and ${z}$ direction, respectively.

Similar expressions can be written for cylindrical and spherical coordinates. In many theoretical aerodynamic and other engineering problems, various scalar and vector fields are unknowns that must be determined as part of a solution subject to prescribed initial and boundary conditions.

Scalar Products & Vector Projections

Let the vectors $\vec{A}$ and $\vec{B}$ be given by

(38) $\begin{eqnarray*} \vec{A} & = & A_{x}\, \vec{i} + A_{y}\, \vec{j} + A_{z} \, \vec{k} \\[6pt] \vec{B} & = & B_{x}\, \vec{i} + B_{y}\, \vec{j} + B_{z} \, \vec{k} \end{eqnarray*}$

The scalar product, or dot product, of $\vec{A}$ and $\vec{B}$ is written as $\vec{A}\bigcdot\vec{B}$ and is defined by

(39) $\begin{equation*} \vec{A}\bigcdot\vec{B} = A_x\,B_x + A_y\,B_y + A_z\,B_z \end{equation*}$

where the centered dot ${\bigcdot}$ represents the scalar product operator. The result of this operation is a scalar, not a vector.

As shown in Figure 6, the scalar product has a useful geometric interpretation. It can also be written as

(40) $\begin{equation*} \vec{A}\bigcdot\vec{B} = | \vec{A} |\,| \vec{B} | \cos \theta \end{equation*}$

where $\theta$ is the angle between the vectors ${\vec{A}}$ and ${\vec{B}}$ . The factor $|\vec{A}|\cos\theta$ is the scalar projection of $\vec{A}$ in the direction of $\vec{B}$ , while $|\vec{B}|\cos\theta$ is the scalar projection of $\vec{B}$ in the direction of $\vec{A}$ . For this reason, the dot product can be interpreted as the magnitude of one vector multiplied by the scalar projection of the other vector onto it.

The physical interpretation of a scalar product is the projection of one vector onto another.

Vector projections isolate the part of one vector that lies along another. This idea recurs whenever a quantity must be decomposed into components parallel and perpendicular to a specified direction. For example, in calculating mechanical work, only the component of the force parallel to the displacement contributes. Similar decompositions also appear in aerospace engineering applications such as structural analysis, where loads may be split into axial and shear components.

The scalar projection, or scalar component, of $\vec{A}$ onto $\vec{B}$ is

(41) $\begin{equation*} \operatorname{comp}_{\vec{B}}\vec{A} =\frac{\vec{A}\bigcdot\vec{B}}{|\vec{B}|} =|\vec{A}|\cos\theta \end{equation*}$

This quantity is a signed scalar. It is positive if $\vec{A}$ has a component in the same general direction as $\vec{B}$ , negative if it has a component in the opposite direction, and zero if the two vectors are perpendicular.

The corresponding vector projection of $\vec{A}$ onto $\vec{B}$ is

(42) $\begin{equation*} \operatorname{proj}_{\vec{B}}\vec{A} = \left( \frac{\vec{A}\bigcdot\vec{B}}{|\vec{B}|^{2}} \right) \, \vec{B} \end{equation*}$

This expression gives the actual vector component of $\vec{A}$ that lies in the direction of $\vec{B}$ . Dividing by $|\vec{B}|$ gives the scalar projection, while dividing by $|\vec{B}|^{2}$ and multiplying by $\vec{B}$ gives the vector projection.

Vector Products

The vector product $\vec{A} \times \vec{B}$ , also called the cross product, produces a vector that is perpendicular to the plane containing $\vec{A}$ and $\vec{B}$ , as shown in Figure 7. The magnitude of the vector product is proportional to the area of the parallelogram formed by the two vectors, i.e.,

(43) $\begin{equation*} \vec{A} \times \vec{B} = | \vec{A} | \, | \vec{B} | \, \sin \theta \, \vec{n} \end{equation*}$

where $\theta$ is the angle between the vectors and $\vec{n}$ is a unit vector normal to the plane containing $\vec{A}$ and $\vec{B}$ , as determined by the right-hand rule. Notice that

(44) $\begin{equation*} |\vec{A} \times \vec{B}| = | \vec{A} | \, | \vec{B} | \, \sin \theta \end{equation*}$

which is the area of the parallelogram spanned by the two vectors.

The physical interpretation of a vector product is another vector that is perpendicular to the other two.

In general, the vector product can also be written as

(45) $\begin{equation*} \vec{A} \times \vec{B} = \left| \begin{array}{ccc} \, \vec{i} & \, \vec{j} & \, \vec{k} \\[6pt] A_x & A_y & A_z \\[6pt] B_x & B_y & B_z \end{array} \right| \end{equation*}$

Expanding this determinant gives

(46) $\begin{equation*} \vec{A}\times \vec{B} = \left( A_y B_z - A_z B_y \right) \, \vec{i} - \left( A_x B_z - A_z B_x \right) \, \vec{j} + \left( A_x B_y - A_y B_x \right) \, \vec{k} \end{equation*}$

Check Your Understanding #5 – Review problems with vector operations

Given two vectors $\vec{A} = -3 \, \vec{i} + 2 \, \vec{j}$ and $\vec{B} = 5 \, \vec{i} - 2 \, \vec{j}$ , then evaluate:

1. $\mbox{~~~}\vec{A} - 3\vec{B}$
2. $\mbox{~~~}\vec{A} \bigcdot \vec{B}$
3. ${\mbox{~~~}\vec{A} \times \vec{B}}$
4. The angle $\theta$ between $\vec{A}$ and $\vec{B}$
5. The projection of $\vec{A}$ onto $\vec{B}$

Show solution/hide solution.

Multiply $\vec{B}$ by 3 and subtract from $\vec{A}$ , i.e.,
${ \vec{A}-3\vec{B} =(-3\,\vec{i}+2\,\vec{j})-3(5\,\vec{i}-2\,\vec{j}) =(-3\,\vec{i}+2\,\vec{j})-(15\,\vec{i}-6\,\vec{j}) =-18\,\vec{i}+8\,\vec{j}}$

Notice that the result is a vector.
Take the dot product of matching components, i.e.,
$\vec{A} \bigcdot \vec{B} = (-3 \, \vec{i} + 2 \, \vec{j} ) \bigcdot ( 5 \, \vec{i} - 2 \, \vec{j} ) = (-3)(5) + (2)(-2) = -15 - 4 = -19$

Notice that the result is a scalar.
Form the determinant (treating the $k$ ‐components as zero), i.e.,
$\vec{A}\times\vec{B} =\begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\[4pt] -3 & 2 & 0 \\[4pt] 5 & -2 & 0 \end{vmatrix} = (-3)(-2) - (5)(2) = 6 - 10 = -4\,\vec{k}$
Find the angle $\theta$ between the vectors using $\displaystyle \cos\theta=\frac{\vec{A} \bigcdot \vec{B}}{|\vec{A}|\,|\vec{B}|}$ . The magnitudes of the vectors are
$|\vec{A}|=\sqrt{(-3)^2+2^2}=\sqrt{13} \quad \text{and} \quad |\vec{B}|=\sqrt{5^2+(-2)^2}=\sqrt{29}$

so that

$\cos\theta =\frac{-19}{\sqrt{13}\,\sqrt{29}} =\frac{-19}{\sqrt{377}} \quad \text{and} \quad \theta=\arccos \left(-\frac{19}{\sqrt{377}}\right) = 168.1^{\circ}$
Project $\vec{A}$ onto $\vec{B}$ using
$\displaystyle \operatorname{proj}_{\vec{B}}\vec{A} =\frac{\vec{A} \bigcdot \vec{B}}{|\vec{B}|^{2}}\, \vec{B}$

Therefore,

$|\vec{B}|^{2} = 29 \quad \text{and} \quad \operatorname{proj}_{\vec{B}}\vec{A} =\frac{-19}{29}\, \left( 5\,\vec{i} -2\,\vec{j} \right) = -\frac{95}{29}\,\vec{i} + \frac{38}{29}\,\vec{j}$

Differentiating a Vector

A vector can be differentiated with respect to any independent variable on which it depends. For example, if a vector $\vec{A}$ depends on a variable $s$ , then its derivative with respect to $s$ is

(47) $\begin{equation*} \frac{d\vec{A}}{ds} = \left( \frac{dA_x}{ds} \right) \, \vec{i} + \left( \frac{dA_y}{ds} \right) \, \vec{j} + \left( \frac{dA_z}{ds} \right) \, \vec{k} \end{equation*}$

where the Cartesian unit vectors $\vec{i}$ , $\vec{j}$ , and $\vec{k}$ are assumed to be fixed. The derivative describes how the vector changes with respect to the chosen independent variable. The change may be in the magnitude of the vector, its direction, or both. The derivative of a vector is also a vector.

In many engineering applications, the independent variable is time. For example, let the vector be $\vec{A}(t)$ , a position vector in space. The time derivative of $\vec{A}$ is a velocity vector, i.e.,

(48) $\begin{equation*} \vec{V} = \frac{d\vec{A}}{dt} = \left( \frac{dA_x}{dt} \right) \, \vec{i} + \left( \frac{dA_y}{dt} \right) \, \vec{j} + \left( \frac{dA_z}{dt} \right) \, \vec{k} = v_x \, \vec{i} + v_y \, \vec{j} + v_z \, \vec{k} = \left( v_x, v_y, v_z \right) \end{equation*}$

In other engineering examples, if $\vec{A}(t)$ were the momentum of something, then the time rate of change of momentum would be a force. If the vector were velocity, then its time derivative would be acceleration. Vector differentiation can yield valuable outcomes in engineering applications. Scalar quantities can also be differentiated with respect to time; for example, the time derivative of work or energy is power.

Check Your Understanding #6 – Differentiating a vector

Differentiate each of the following vector functions:

1. $\vec{A}(t) = \left( 3t^{2},\cos(t),\ln(t) \right)$
2. $\vec{B}(t) = \left( e^{2t},\,t^{3},\,\sqrt{t} \right)$
3. $\vec{C}(t) = \left( \sin 3t,\,1/t,\,\arctan(t) \right)$

Show solution/hide solution.

$\dfrac{d\vec{A}}{dt} =\dfrac{d}{dt}\left( 3t^{2},\cos t,\ln t \right) =\left( 6t,\,-\sin t,\,1/t \right)$
$\dfrac{d\vec{B}}{dt} =\dfrac{d}{dt}\left( e^{2t},\,t^{3},\,\sqrt{t} \right) =\left( 2e^{2t},\,3t^{2},\,\dfrac{1}{2\sqrt{t}} \right)$
$\dfrac{d\vec{C}}{dt} =\dfrac{d}{dt}\left( \sin 3t,\,1/t,\,\arctan t \right) =\left( 3\cos 3t,\,-1/t^{2},\,\dfrac{1}{1+t^{2}} \right)$

Line Integrals

Consider a vector field $\vec{A} = \vec{A} \left( x, y, z \right)$ . Also, consider a curve in space connecting two points $a$ and $b$ , as shown in Figure 8. Let ${ds}$ be an elemental length of the curve.

The line integral of a vector field along the curve from $a$ to $b$ is

(49) $\begin{equation*} \int_{a}^{b}\vec{A}\bigcdot d\vec{s} \end{equation*}$

which represents the accumulated component of the vector field tangent to the curve. If the curve is closed, i.e., points $a$ and $b$ are coincident, then the line integral is

(50) $\begin{equation*} \oint_{C}\vec{A}\bigcdot d\vec{s} \end{equation*}$

where, for a closed curve in a plane, the counterclockwise direction is usually taken as positive. More generally, the positive direction around a closed curve is defined relative to the chosen surface normal by the right-hand rule, as shown in Figure 9. Line integrals can also be based on scalar functions. For example, the scalar field $p(x,y)$ integrated along a curve $C$ is written as

(51) $\begin{equation*} \int_C p(x,y)\,ds \end{equation*}$

where ${ds}$ represents an infinitesimal length along the curve $C$ .

A closed-loop line integral sums the effects along a closed path and is frequently used in engineering.

Check Your Understanding #7 – Line Integrals

Evaluate each of the following line integrals:

$p(x,y)=x^{2}+y$ along the straight-line segment from (0, 0) to (2, 4).
$p(x,y)=x+y^{2}$ along the quarter-circle arc of radius 1 from (1, 0) to (0, 1) in the first quadrant.
$p(x,y)=x^{2}y$ along the straight-line segment from (0, 0) to (1, 1).

Show solution/hide solution.

First, parameterise the line by $x(t) = 2t,\, y(t) = 4t$ with $0 \le t\le 1$ . Then
$\left( \frac{dx}{dt}\right)^{2}+ \left( \frac{dy}{dt} \right)^{2} = 2^{2}+4^{2} = 20$

so $ds=2\sqrt5 \,dt$ , and

${ \int_{C}(x^{2}+y)\,ds=\int_{0}^{1}\left( 4t^{2} + 4t \right)2\sqrt5\,dt =2\sqrt5\int_{0}^{1}(4t^{2}+4t)\,dt =2\sqrt5 \left( \dfrac{4}{3} t^{3}+2t^{2}\right)_{0}^{1} =\dfrac{20}{3} \sqrt5 }$
Parameterise the quarter-circle using $x(\theta)=\cos\theta$ and $y(\theta)=\sin\theta$ with $0\le\theta\le\pi/2$ . Notice that $ds = d\theta$ because $\sqrt{(-\sin\theta)^{2}+(\cos\theta)^{2}}=1$ and so
$\int_{C}(x+y^{2})\,ds=\int_{0}^{\pi/2} \left( \cos\theta+\sin^{2}\theta \right)\,d\theta=\left[ \sin\theta+\dfrac{1}{2}\theta-\dfrac{1}{4}\sin2\theta \right]_{0}^{\pi/2}=1+\frac{\pi}{4}$
Parameterise the line using $x(t)=t,\, y(t)=t$ with $0\le t\le1$ . Then
${ ds=\sqrt{(dx/dt)^{2} + (dy/dt)^{2}}\ dt=\sqrt{1 + 1}\,dt=\sqrt2\,dt }$

The integrand becomes $x^{2} y =t^{3}$ , so that

$\int_{C}x^{2}y\,ds=\sqrt2\int_{0}^{1}t^{3}\,dt=\sqrt2 \left[ \dfrac{1}{4} t^{4}\right]_{0}^{1}=\frac{\sqrt2}{4}$

Surface Integrals

Consider now a surface $S$ bounded by a curve $C$ , as shown in Figure 10. Let $dS$ be an elemental surface area and let $\vec{{n}}$ be a unit normal vector, i.e., one perpendicular to the surface. Let $p$ be a scalar field in this space and $\vec{A}$ a vector field. A surface integral over the surface can be defined as

(52) $\begin{equation*} \iint_{S} p \, dS \quad \mbox{or} \quad \iint_{S} \vec{A}\bigcdot d\vec{S} \end{equation*}$

where $d\vec{S} = \vec{n}\, dS$ is the elemental vector area. If the surface is closed, then the integral is written as

(53) $\begin{equation*} \oiint_{S} p \, dS \quad \mbox{or} \quad \oiint_{S} \vec{A}\bigcdot d\vec{S} \end{equation*}$

Unit Normal Vector

The meaning of a unit normal vector, also known as a normal unit vector or simply a normal vector, should be understood because it comes up in many engineering problems. The unit normal is a vector perpendicular (or orthogonal) to a surface at a particular point, the idea being shown in Figure 11. It has a magnitude of 1, hence the term “unit.” The unit normal vector is perpendicular to the surface; for closed surfaces, the outward normal is used by convention.

The unit normal vector is often denoted by the symbol $\vec{{n}}$ . For a surface defined by a function $\phi (x, y, z) = 0$ , where $(x, y, z)$ are the coordinates of a point on the surface, the unit normal vector at that point can be obtained by taking the gradient of the function and normalizing it, as given by

(54) $\begin{equation*} \vec{n} = \cfrac{ \left( \displaystyle{\cfrac{\partial \phi}{\partial x}} , \, \displaystyle{\cfrac{\partial \phi}{ \partial y}} , \, \displaystyle{\cfrac{\partial \phi}{ \partial z}} \right) }{\left| \left( \displaystyle{\cfrac{\partial \phi}{ \partial x}} , \, \displaystyle{\cfrac{\partial \phi}{ \partial y}} , \, \displaystyle{\cfrac{\partial \phi}{ \partial z}} \right) \right| } \end{equation*}$

where $|~~|$ denotes the Euclidean L² norm or magnitude of the vector. By normalizing the gradient vector, the resulting unit normal vector will have unit magnitude. This normalization is achieved by dividing the gradient vector by its magnitude.

Various methods exist for computing the unit normal at discrete surface points, such as the cross-product approach and the least-squares fitting procedure. The cross-product approach is the easiest. The steps are:

Select three neighboring points on the surface defined by the vectors ${\vec{P}_1 = (x_1, y_1, z_1)}$ , ${\vec{P}_2 = (x_2, y_2, z_2)}$ , and ${\vec{P}_3 = (x_3, y_3, z_3)}$ .
Calculate the vectors between these points, i.e., $\vec{A}_1 = \vec{P}_2 - \vec{P}_1$ and $\vec{A}_2 = \vec{P}_3 - \vec{P}_1$ .
Compute the cross product of these vectors giving $\vec{N} = \vec{A}_1 \times \vec{A}_2$ .
Normalize the resulting vector to obtain the unit normal vector, i.e., $\vec{n} = \dfrac{\vec{N}}{|| \vec{N} ||}$ .

Remember that the unit normal vector is a crucial concept in vector calculus, used in various mathematical and physical applications, such as determining the direction of a force on a surface and solving problems involving surface integrals or differential equations.

Check Your Understanding #8 – Calculating unit normal vectors

In a three-dimensional Cartesian coordinate system, what is the unit vector in the positive $x$ -direction?
Consider the three-dimensional plane defined by $x + y = 1$ . What is one unit normal to this plane?
Consider a three-dimensional plane defined by $2x - 3y + z = 4$ . What is one unit normal to this plane?

Show solution/hide solution.

This vector is parallel to the ${x}$ -axis and points in the positive ${x}$ -direction, so $\vec{n} = (1, 0, 0)$ .
In this case, $\phi(x,y,z) = x + y - 1 = 0$ . Therefore, at any point on this plane, one unit normal vector is
$\vec{n} = \cfrac{ \left( \displaystyle{\cfrac{\partial \phi}{\partial x}}, \, \displaystyle{\cfrac{\partial \phi}{ \partial y}}, \, \displaystyle{\cfrac{\partial \phi}{ \partial z}} \right) }{ \left| \left( \displaystyle{\cfrac{\partial \phi}{ \partial x}}, \, \displaystyle{\cfrac{\partial \phi}{ \partial y}}, \, \displaystyle{\cfrac{\partial \phi}{ \partial z}} \right) \right| } = \cfrac{ (1, 1, 0)}{\sqrt{2}} = \left( \cfrac{1}{\sqrt{2}}, \, \cfrac{1}{\sqrt{2}}, \, 0 \right)$
In this final case, $\phi(x, y, z) = 2x - 3y + z - 4 = 0$ . Therefore, at any point on this plane, one unit normal vector is
${ \vec{n} = \cfrac{(2, -3, 1)}{\sqrt{14}} = \left(\cfrac{2}{\sqrt{14}}, \, \cfrac{-3}{\sqrt{14}}, \, \cfrac{1}{\sqrt{14}} \right) }$

Volume Integrals

Consider a volume $\cal{V}$ in space, as shown in Figure 12, which contains a small elemental volume $d \cal{V}$ . Let $\varrho$ be a scalar field in this space. The volume integral over the volume $\cal{V}$ of the quantity $\varrho$ is written as

(55) $\begin{equation*} \iiint_{\cal{V}} \varrho \, d\cal{V} \end{equation*}$

The result of the integration is a scalar.

If ${\vec{V}}$ is a vector field in space, then the volume integral is written as

(56) $\begin{equation*} \iiint_{\cal{V}} \vec{V} \, d\cal{V} \end{equation*}$

And in the latter case, the result will be a vector. If $\vec{V} = (V_1, V_2, V_3)$ , then the volume integral of ${\vec{V}}$ over a region $\cal{V}$ is evaluated using

(57) $\begin{equation*} \iiint_{\cal{V}} \vec{V} \, d{\cal{V}} = \left( \iiint_{\cal{V}} V_1 \, d{\cal{V}}, \ \iiint_{\cal{V}} V_2 \, d{\cal{V}}, \ \iiint_{\cal{V}} V_3 \, d{\cal{V}} \right) \end{equation*}$

which is a vector representing the accumulated effect of the vector field throughout the entire volume.

Check Your Understanding #9 – Volume integrals

Evaluate the volume integral that gives the volume of a solid sphere of radius $R$ , i.e.,
${\cal{V}} = \iiint_{\rm sphere} d{\cal{V}}$

where the sphere’s domain is ${x^{2}+y^{2}+z^{2}\le R^{2}}$ .
A right circular cone of height $H$ and base radius $R$ has a density that varies linearly with height according to
$\varrho(z) = \varrho_{0}\,\frac{z}{H}, \qquad 0\le z\le H$

Set up and evaluate the volume integral for the total mass, i.e.,

$M = \iiint_{\text{cone}} \varrho(z)\,d{\cal{V}}$

Show solution/hide solution.

Using spherical coordinates $(r,\theta,\phi)$ with $r\in[0,R]$ , $\theta\in[0,\pi]$ , and $\phi\in[0,2\pi]$ , then
${\cal{V}} = \iiint_{\rm sphere} d{\cal{V}} = \int_{0}^{2\pi}\!\!\int_{0}^{\pi}\!\!\int_{0}^{R} r^{2}\sin\theta \, dr\,d\theta\,d\phi = \frac{4}{3}\pi R^{3}$
In cylindrical coordinates $(r,\phi,z)$ the cone is given by $0\le z\le H$ and $0\le r\le (R/H)z$ . The mass integral is
$M = \iiint_{\text{cone}} \varrho_{0}\,\frac{z}{H}\,d{\cal{V}} = \varrho_{0}\!\!\int_{0}^{H}\!\!\int_{0}^{2\pi}\!\!\int_{0}^{(R/H)z} \frac{z}{H}\, r \, dr\,d\phi\,dz$

Carrying out the integrations gives

$M = \varrho_{0}\! \int_{0}^{H}\!\! \left( \frac{z}{H}\, \frac{1}{2}\Bigl(\tfrac{R}{H}z\Bigr)^{2} \right) 2\pi\, dz = \varrho_{0}\,\frac{\pi R^{2} H}{4}$

Integral Relations

The relationships among line, surface, and volume integrals arise naturally in many engineering problems. Consider an area $S$ bounded by a closed curve $C$ . Let $\vec{A}$ be a vector field. The line integral of $\vec{A}$ over $C$ is related to the surface integral of $\vec{A}$ over $S$ by using Stokes’ Theorem, i.e.,

(58) $\begin{equation*} \oint_{C}\vec{A}\bigcdot d\vec{s} = \iint_S \left( \nabla \times \vec{A} \right)\bigcdot d \vec{S} \end{equation*}$

This relation shows that the circulation of a vector field around a closed curve is determined by the curl of the field over the enclosed surface, i.e., by the local rotational behavior of the field.

Again, consider the volume ${\cal{V}}$ enclosed by the closed surface $S$ . Then, the surface and volume integrals of the vector field $\vec{A}$ are related by using Gauss’s divergence theorem, i.e.,

(59) $\begin{equation*} \oiint_{S} \vec{A}\bigcdot d\vec{S} = \iiint_{\cal{V}} \left( \nabla \bigcdot \vec{A} \right) \, d\cal{V} \end{equation*}$

which connects a volume integral to a surface integral. This relation shows that the net flux of a vector field through a closed surface is determined by the divergence of the field within the enclosed volume. Finally, if $\varrho$ represents a scalar field, then for a conservative field the line integral of its gradient between two points depends only on the end points, i.e.,

(60) $\begin{equation*} \int_{a}^{b} \nabla \varrho \bigcdot d\vec{s} = \varrho(b) - \varrho(a) \end{equation*}$

Gradient of a Scalar Field

If $p = p\left(x, y, z\right)$ , then the gradient $\nabla p$ (or grad $p$ ) is defined as

(61) $\begin{equation*} \mbox{grad } p \equiv \nabla p = \frac{\partial p}{\partial x}\, \vec{i}+ \frac{\partial p}{\partial y}\, \vec{j}+\frac{\partial p}{\partial z}\, \vec{k} \end{equation*}$

where the operator $\nabla$ (which is called nabla or del) is defined as

(62) $\begin{equation*} \nabla = \frac{\partial }{\partial x}\, \vec{i}+ \frac{\partial }{\partial y}\, \vec{j}+\frac{\partial }{\partial z}\, \vec{k} \end{equation*}$

Physically, $\nabla p$ is a vector normal to the surface defined by $p(x,y,z)$ = constant, and its magnitude represents the rate of change of $p$ in the direction of greatest increase. In fluid mechanics, this means that fluid accelerates in the direction of decreasing pressure because the force per unit volume is proportional to $-\nabla p$ .

A gradient of a scalar field is a vector, i.e., a vector that points in the direction of the slope of the scalar field.

Related to the gradient is the concept of a directional derivative. Consider some scalar quantity $p\left(x,y\right)$ defined at some point over a surface. Choose some arbitrary direction ${s}$ away from the point and let $\vec{{n}}$ be a unit vector in that direction. The rate of change of $p$ per unit length in the ${s}$ direction is

(63) $\begin{equation*} \frac{dp}{ds} = \nabla p \bigcdot \vec{n} \end{equation*}$

The directional derivative is defined as the rate of change of the scalar field in a specified direction, i.e., $\dfrac{dp}{ds} = \nabla p \bigcdot \vec{n}$ . Therefore, it can be appreciated that the rate of change of $p$ in any arbitrary direction is simply the component of $\nabla p$ in that particular direction.

Check Your Understanding #10 – Gradient of a scalar field

A two-dimensional pressure field is given by $p(x,y) = -2x - y + 10$ . If an isobar is a constant pressure line, plot the isobars for $p$ = 0, 2, 4, 6, and 8.

Show solution/hide solution.

The equations of the isobars are

$\begin{eqnarray*} 0 & = & -2x - y + 10 \implies y = -2x + 10 \\[6pt] 2 & = & -2x - y + 10 \implies y = -2x + 8 \\[6pt] 4 & = & -2x - y + 10 \implies y = -2x + 6 \\[6pt] 6 & = &-2x - y + 10 \implies y = -2x + 4 \\[6pt] 8 & = & -2x - y + 10 \implies y = -2x + 2 \end{eqnarray*}$

which are all straight lines, as shown in the figure below.
The gradient is given by

${ \nabla p = \frac{\partial p}{\partial x} \vec{i} + \frac{\partial p}{\partial y} \vec{j} }$

so that

$\nabla (-2x - y + 10) = \frac{\partial (-2x - y + 10)}{\partial x} \vec{i} + \frac{\partial (-2x - y + 10)}{\partial y} \vec{j} = -2 \vec{i} - 1 \vec{j}$

and the magnitude of the pressure gradient is

${ | \nabla p | = \sqrt{(-2)^2 + (-1)^2} = \sqrt{5} = 2.24 }$

The pressure gradient vector represents the magnitude and direction of the pressure change at a specific point in space. It points in the direction of the steepest increase in pressure, and its magnitude indicates the rate of change of pressure in that direction. In this case, the pressure gradient vector is given by $(-2, -1)$ , which means that the pressure increases most rapidly in the direction $(-2, -1)$ .

Pressure gradient vectors are typically plotted as arrows with their tails placed at the corresponding points in space. The direction and length of the arrows indicate the pressure gradient’s direction and magnitude at those points. Some arbitrary points can be selected, and then the corresponding pressure gradient vectors are plotted, e.g.,

${ \mbox{\small Point A: } (x, y) = (2, 4) \implies \mbox{\small $\nabla p$ at A: } (-2, -1) }$

${ \mbox{\small Point B: } (x, y) = (0, 6) \implies \mbox{\small $\nabla p$ at B: } (-2, -1) }$

${ \mbox{\small Point C: } (x, y) = (3, -2) \implies \mbox{\small $\nabla p$ at C: } (-2, -1) }$

Notice that the gradient vector runs perpendicular to the isobars. The pressure gradient vector is the same at points A, B, and C.

Check Your Understanding #11 – Gradient and directional change

A scalar field is given by $p(x,y) = x^2 + 2y^2$ . Answer the following:

1. Determine the gradient $\nabla p$ .

2. Evaluate $\nabla p$ at the point $(1,1)$ .

3. A unit vector in the direction $(3,4)$ is given by $\vec{n} = \dfrac{3}{5}\,\vec{i} + \dfrac{4}{5}\,\vec{j}$ . Determine the directional derivative of $p$ at $(1,1)$ in this direction.

4. Explain physically what this directional derivative represents.

Show solution/hide solution.

1. The gradient is

$\nabla p = \frac{\partial p}{\partial x}\,\vec{i} + \frac{\partial p}{\partial y}\,\vec{j} = 2 x\,\vec{i} + 4 y\,\vec{j}$

2. At $(1,1)$

$\nabla p = 2\,\vec{i} + 4\,\vec{j}$

3. The directional derivative is

$\frac{dp}{ds} = \nabla p \bigcdot \vec{n} = (2\,\vec{i} + 4\,\vec{j}) \bigcdot \left( \frac{3}{5}\,\vec{i} + \frac{4}{5}\,\vec{j} \right) = \frac{6}{5} + \frac{16}{5} = \frac{22}{5}$

4. This result represents the rate of change of the scalar field $p$ at the point $(1,1)$ in the specified direction. It gives the slope of the surface in that direction and indicates how rapidly $p$ increases along that path.

Divergence of a Vector Field

Consider a vector field $\vec{V} = \vec{V} \left(x, y, z \right)$ , then

(64) $\begin{eqnarray*} \mbox{div } \vec{V} = \nabla \bigcdot \vec{V} & = & \left(\frac{\partial }{\partial x}\, \vec{i}+\frac{\partial }{\partial y} \, \vec{j}+\frac{\partial }{\partial z}\, \vec{k}\right) \bigcdot \left(V_{x}\, \vec{i}+V_{y}\, \vec{j}+V_{z}\, \vec{k}\right) \nonumber \\[6pt] & = & \frac{\partial V_{x}}{\partial x}+\frac{\partial V_{y}}{\partial y}+\frac{\partial V_{z}}{\partial z} \end{eqnarray*}$

which is a scalar quantity. In aerodynamics and fluid dynamics, divergence can be interpreted as a flux density, i.e., the rate at which flux enters or leaves a point in space. Mass fluxes such as $\varrho u$ , $\varrho v$ , etc., are often used. Divergence is the rate of flux expansion (positive divergence) or contraction (negative divergence). By convention, the flux is defined as positive when it flows out of a closed surface. So, the divergence is just the net flux per unit volume or flux density of the contents of a given region of space.

Curl of a Vector Field

Again, consider the vector field defined by

(65) $\begin{equation*} \vec{V} = V_{x}\, \vec{i}+V_{y}\, \vec{j}+V_{z}\, \vec{k} \end{equation*}$

Then the curl of ${\vec{V}}$ is defined as

(66) $\begin{equation*} \mbox{curl } \vec{V} \equiv \nabla \times \vec{V} = \left| \begin{array}{ccc} \, \vec{i} &\, \vec{j} & \, \vec{k} \\[6pt] \displaystyle{\frac{\partial}{\partial x}} & \displaystyle{\frac{\partial}{\partial y}} & \displaystyle{\frac{\partial}{\partial z}} \\[10pt] V_x & V_y & V_z \end{array} \right| \end{equation*}$

and expanding this determinant gives

(67) $\begin{equation*} \nabla \times \vec{V} = \left( \frac{\partial V_{z}}{\partial y} - \frac{\partial V_{y}}{\partial z} \right) \, \vec{i} + \left( \frac{\partial V_{x}}{\partial z} - \frac{\partial V_{z}}{\partial x} \right) \, \vec{j} + \left( \frac{\partial V_{y}}{\partial x} - \frac{\partial V_{x}}{\partial y} \right) \, \vec{k} \end{equation*}$

A physical interpretation of the curl is that it measures the local rotation or “spin” of the contents of a given region of space. In aerodynamics, the curl of a velocity field is called vorticity; for a velocity field, the vorticity is twice the local angular velocity of a fluid element. If the vorticity is zero, the flow is irrotational; if it is nonzero, the flow is rotational.

Check Your Understanding #12 – Finding the curl of a vector field

In a solution to a specific engineering problem, the solution of the vector field $\vec{A}(x, y, z)$ is

$\begin{equation*} \vec{A} = xz \, \vec{i} - y^2 \, \vec{j} + 2 x^2y \, \vec{k} \end{equation*}$

Find the curl of this field.

Show solution/hide solution.

The curl of a vector requires the evaluation of a determinant, i.e.,

$\nabla \times \vec{A} = \left| \begin{array}{ccc} \, \vec{i} &\, \vec{j} & \, \vec{k} \\[6pt] \displaystyle{\frac{\partial}{\partial x}} & \displaystyle{\frac{\partial}{\partial y}} & \displaystyle{\frac{\partial}{\partial z}} \\[6pt] A_x & A_y & A_z \end{array} \right|$

Evaluating this determinant gives

$\nabla \times \vec{A} = \left( \frac{\partial A_{z}}{\partial y} - \frac{\partial A_{y}}{\partial z} \right) \, \vec{i} + \left( \frac{\partial A_{x}}{\partial z} - \frac{\partial A_{z}}{\partial x} \right) \, \vec{j} + \left( \frac{\partial A_{y}}{\partial x} - \frac{\partial A_{x}}{\partial y} \right) \, \vec{k}$

After differentiation then

$\nabla \times \vec{A} = \left( 2x^2 - 0 \right) \, \vec{i} + \left( x - 4xy \right) \, \vec{j} + \left( 0 - 0 \right) \, \vec{k}$

or

$\nabla \times \vec{A} = 2x^2 \, \, \vec{i} +( x - 4xy) \, \vec{j}$

Other Mathematical Operators

Mathematical operators are fundamental in fields such as calculus, differential equations, and physics, and are essential tools for describing and analyzing engineering problems. In addition to operators such as the gradient operator ${\nabla}$ , the divergence ${\nabla \bigcdot \,}$ , and the curl ${\nabla \times \,}$ , the Laplace operator and substantial derivative operator are frequently used.

Laplace Operator

The Laplace operator appears in several of the equations of fluid mechanics and is a second-order differential operator defined as the divergence (i.e., $\nabla \bigcdot$ ) of the gradient (i.e., $\nabla$ ). The Laplace operator can be applied to both scalar and vector fields.

In the case of a scalar field, say $\phi$ , the Laplace operator would be written in Cartesian coordinates as

(68) $\begin{equation*} \nabla^2 \,\phi = \nabla \bigcdot \nabla \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} \end{equation*}$

Of course, the Laplace operator appears in the familiar Laplace equation, i.e.,

(69) $\begin{equation*} \nabla^2 \,\phi = 0 \end{equation*}$

In aerodynamics, the Laplace equation is, in fact, the governing equation for an incompressible, irrotational flow. To this end, the linearity of the Laplace equation is helpful in that if $\phi_1$ and $\phi_2$ are both flow solutions to the Laplace equation, then $\phi_1 + \phi_2$ is also a solution, i.e., the idea that incompressible, irrotational flows can be combined using the principle of linear superposition, which is a significant advantage when dealing with complex flows.

The Laplace operator can also be applied to a vector field, i.e., $\vec{V} = u \, \vec{i} + v \, \vec{j} + w \, \vec{k}$ . In this case, then

(70) $\begin{equation*} \nabla^2 \,\vec{V} = \nabla^2 \,u \, \, \vec{i} + \nabla^2 \,v \, \, \vec{j} + \nabla^2 \,w \, \, \vec{k} \end{equation*}$

where

(71) $\begin{eqnarray*} \nabla^2 \,u & = & \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \\[6pt] \nabla^2 \,v & = & \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial z^2} \\[6pt] \nabla^2 \,w & = & \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} \end{eqnarray*}$

Check Your Understanding #13 – Finding the Laplacian of a scalar field

A two-dimensional function $\psi (x, y) = x^2 + y^2$ , where ${x}$ and ${y}$ are the Cartesian coordinates. Determine the Laplacian of this function.

Show solution/hide solution.

1. To find the Laplacian, the second partial derivatives must be computed and then summed. Calculate the first partial derivatives:

$\frac{\partial \psi}{\partial x} = 2x \mbox{~~and~~}\frac{\partial \psi}{\partial y} = 2y$

2. Calculate the second partial derivatives, i.e.,

$\frac{\partial^2 \psi}{\partial x^2} = \frac{\partial}{\partial x}\left(\frac{\partial \psi}{\partial x}\right) = \frac{\partial}{\partial x}(2x) = 2$

$\frac{\partial^2 \psi}{\partial y^2} = \frac{\partial}{\partial y}\left(\frac{\partial \psi}{\partial y}\right) = \frac{\partial}{\partial y}(2y) = 2$

3. Sum the second partial derivatives to obtain the Laplacian, i.e.,

$\nabla^2 \,\psi = \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} = 2 + 2 = 4$

So, in this example, the Laplacian of the function $\psi (x, y) = x^2 + y^2$ is 4.

Substantial Derivative Operator

The substantial derivative operator in Cartesian coordinates is written as

(72) $\begin{equation*} \frac{D}{Dt} = \frac{\partial}{\partial t} + \frac{dx}{dt}\frac{\partial}{\partial x} + \frac{dy}{dt}\frac{\partial}{\partial y} + \frac{dz}{dt}\frac{\partial}{\partial z} \end{equation*}$

or because the velocity components are given by

(73) $\begin{equation*} u = \frac{dx}{dt}; \quad v = \frac{dy}{dt}; \quad w = \frac{dz}{dt} \end{equation*}$

then

(74) $\begin{equation*} \frac{D}{Dt} = \frac{\partial}{\partial t} + u\frac{\partial}{\partial x} + v\frac{\partial}{\partial y} + w\frac{\partial}{\partial z} \end{equation*}$

Recall that the vector gradient operator $\nabla$ is defined as

(75) $\begin{equation*} \nabla = \frac{\partial}{\partial x} \, \vec{i} + \frac{\partial}{\partial y} \, \vec{j} + \frac{\partial}{\partial z} \, \vec{k} \end{equation*}$

Hence, the substantial derivative operator can be written as

(76) $\begin{equation*} \frac{D}{Dt} = \frac{\partial}{\partial t} + \vec{V} \bigcdot \nabla \end{equation*}$

The substantial derivative can be applied to any flow-field variable, scalar, or vector. For example, if ${u}$ is the component of the velocity field $\vec{V} = u\,\vec{i} + v\,\vec{j} + w\,\vec{k}$ in the $x$ direction, then ${Du/Dt}$ can be written as

(77) $\begin{equation*} \frac{Du}{Dt} = \underbrace{ \frac{\partial u}{\partial t}}_{\rm Local~derivative} \hspace*{-3mm} + \underbrace{ u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}}_{\rm Convective~derivative} \end{equation*}$

The term $\partial/\partial t$ is called the local derivative, which can be interpreted physically as the time rate of change of a given quantity at a fixed point. The term $\vec{V} \bigcdot \nabla$ is called the convective derivative, which can be interpreted as the contribution to the rate of change caused by motion through a spatially varying flow field. It will be apparent that terms like $u (\partial /\partial x)$ have units of time $^{-1}$ .

Notice that for a vector ${\vec{V}}$ then $D\vec{V}/Dt$ is

(78) $\begin{equation*} \frac{D\vec{V}}{Dt} = \frac{\partial \vec{V}}{\partial t} + u\frac{\partial \vec{V}}{\partial x} + v\frac{\partial \vec{V}}{\partial y} + w\frac{\partial \vec{V}}{\partial z} \end{equation*}$

which, in terms of the scalar components, is

(79) $\begin{equation*} \frac{D\vec{V}}{Dt} = \left( \frac{D u}{Dt}, \, \frac{D v}{Dt}, \, \frac{D w}{Dt} \right) \end{equation*}$

where

(80) $\begin{eqnarray*} \frac{D u}{Dt} & = & \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} \\[8pt] \frac{D v}{Dt} & = & \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} \\[8pt] \frac{D w}{Dt} & = & \frac{\partial w}{\partial t} + u\frac{\partial w}{\partial x} + v\frac{\partial w}{\partial y} + w\frac{\partial w}{\partial z} \end{eqnarray*}$

Check Your Understanding #14 – Calculating a substantial derivative

A fluid flows through a pipe, and the temperature distribution is given by $T(x, y, z, t)$ = $10x + 5y - 2z + 3t$ . What is the change in temperature at the point where the fluid velocity is $\vec{V} = (2, 3, -1)$ ?

Show solution/hide solution.

The temporal change and the convective terms involving the gradients and velocity components must be determined to calculate the substantial derivative of temperature at that point. The substantial derivative is

$\frac{D T}{Dt} = \frac{\partial T }{\partial t} + u \frac{\partial T }{\partial x} + v \frac{\partial T}{\partial y} + w\frac{\partial T}{\partial z}$

The contributing terms are

$\frac{\partial T }{\partial t} = 3 ; \quad \frac{\partial T }{\partial x} = 10; \quad \frac{\partial T }{\partial y} = 5; \quad \frac{\partial T }{\partial z} = -2$

Therefore,

$\frac{D T}{Dt} = 3 + 10 (2) + 5 (3) - 2 (-1) = 40~\mbox{units of temperature per unit time}$

So, the substantial temperature derivative at the point in this example is 40.

Notice that the substantial derivative incorporates both the local temporal change in temperature and the convective transport from fluid motion. In this case, the local temporal change in temperature alone contributes a rate of 3 units per unit of time. The remaining terms, which involve convective transport and account for temperature gradients and fluid velocity, contribute a total rate of 37 units per unit time.

Matrices & Systems of Linear Equations

A matrix is a rectangular array of numbers organized into rows and columns. It is denoted as having dimensions $m \times n$ , where $m$ is the number of rows and $n$ is the number of columns. Matrices are widely used to represent and solve systems of equations, perform transformations, and model phenomena in various engineering applications.

Types of Matrices

A square matrix has the same number of rows and columns, i.e., ${m} = {n}$ . For example,

(81) $\begin{equation*} \begin{bmatrix} 1 & 2 & 3 \\[6pt] 4 & 5 & 6 \\[6pt] 7 & 8 & 9 \end{bmatrix} \end{equation*}$

The identity matrix, $I$ , is a square matrix with 1s on the main diagonal and 0s elsewhere. For a square matrix $A$ of compatible size, multiplication by the identity matrix leaves $A$ unchanged, i.e.,

(82) $\begin{equation*} A I = I A = A \end{equation*}$

For example,

(83) $\begin{equation*} \begin{bmatrix} 1 & 2 & 3 \\[6pt] 4 & 5 & 6 \\[6pt] 7 & 8 & 9 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\[6pt] 0 & 1 & 0 \\[6pt] 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 2 & 3 \\[6pt] 4 & 5 & 6 \\[6pt] 7 & 8 & 9 \end{bmatrix} \end{equation*}$

A zero matrix has all elements equal to 0, i.e.,

(84) $\begin{equation*} \begin{bmatrix} 0 & 0 & 0\\[6pt] 0 & 0 & 0\\[6pt] 0 & 0 & 0 \end{bmatrix} \end{equation*}$

A diagonal matrix has zero elements everywhere off the main diagonal; the diagonal entries themselves may be zero or nonzero. For example,

(85) $\begin{equation*} \begin{bmatrix} 5 & 0 & 0\\[6pt] 0 & 7 & 0\\[6pt] 0 & 0 & 9 \end{bmatrix} \end{equation*}$

A symmetric matrix is a square matrix that is equal to its transpose, i.e., $A = A^T$ . This means that the element in the $i$ -th row and $j$ -th column, $A_{ij}$ , is equal to the element in the $j$ -th row and $i$ -th column, $A_{ji}$ , for all $i$ and $j$ . For example, a symmetric matrix is

(86) $\begin{equation*} A = \begin{bmatrix} 1 & 2 & 3 \\[6pt] 2 & 5 & 6 \\[6pt] 3 & 6 & 9 \end{bmatrix} \end{equation*}$

To verify that $A$ is symmetric, its transpose $A^T$ is obtained by flipping the matrix over its diagonal, i.e.,

(87) $\begin{equation*} A^T = \begin{bmatrix} 1 & 2 & 3 \\[6pt] 2 & 5 & 6 \\[6pt] 3 & 6 & 9 \end{bmatrix} \end{equation*}$

Because $A = A^T$ , this matrix is symmetric.

Matrix Operations

Addition and subtraction are performed element by element, and the matrices must have the same dimensions. For example, if

(88) $\begin{equation*} A = \begin{bmatrix} 1 & 2 \\[6pt] 3 & 4 \end{bmatrix} \quad \text{and} \quad B = \begin{bmatrix} 5 & 6 \\[6pt] 7 & 8 \end{bmatrix} \end{equation*}$

then

(89) $\begin{equation*} A + B = \begin{bmatrix} 6 & 8 \\[6pt] 10 & 12 \end{bmatrix} \end{equation*}$

Matrix multiplication is defined when the number of columns in the first matrix equals the number of rows in the second matrix. If $C=AB$ , then

(90) $\begin{equation*} C_{ij} = \sum_{k=1}^{n} A_{ik}B_{kj} \end{equation*}$

For example, if

(91) $\begin{equation*} A = \begin{bmatrix} 1 & 2 \\[6pt] 3 & 4 \end{bmatrix} \quad \text{and} \quad B = \begin{bmatrix} 5 & 6 \\[6pt] 7 & 8 \end{bmatrix} \end{equation*}$

then

(92) $\begin{equation*} A B = \begin{bmatrix} 19 & 22 \\[6pt] 43 & 50 \end{bmatrix} \end{equation*}$

The transpose of a matrix flips it over its main diagonal. For example,

(93) $\begin{equation*} A = \begin{bmatrix} 1 & 2 \\[6pt] 3 & 4 \end{bmatrix} \quad \text{so that} \quad A^T = \begin{bmatrix} 1 & 3 \\[6pt] 2 & 4 \end{bmatrix} \end{equation*}$

A square matrix is invertible if its determinant is nonzero, i.e., $\det(A) \neq 0$ . The inverse matrix, $A^{-1}$ , satisfies

(94) $\begin{equation*} A A^{-1} = A^{-1} A = I \end{equation*}$

The determinant is a scalar quantity associated with a square matrix. It is denoted by $\det(A)$ and provides important information about whether the matrix is invertible. If $\det(A)=0$ , then the matrix is called singular and no inverse exists. If $\det(A)\neq 0$ , then the matrix is nonsingular and has an inverse. Geometrically, the determinant can also be interpreted as a scale factor for area or volume under the linear transformation represented by the matrix. A determinant of zero means that the transformation collapses area or volume, so the original variables cannot be uniquely recovered.

Jacobian Matrix

A Jacobian matrix is a matrix of first partial derivatives. It is used when several dependent variables are each influenced by several independent variables. For example, suppose

(95) $\begin{equation*} \vec{f}(x,y,z) = \begin{bmatrix} f_1(x,y,z) \\[6pt] f_2(x,y,z) \\[6pt] f_3(x,y,z) \end{bmatrix} \end{equation*}$

Then the Jacobian matrix of $\vec{f}$ with respect to $x$ , $y$ , and $z$ is

(96) $\begin{equation*} J = \frac{\partial(f_1,f_2,f_3)}{\partial(x,y,z)} = \begin{bmatrix} \dfrac{\partial f_1}{\partial x} & \dfrac{\partial f_1}{\partial y} & \dfrac{\partial f_1}{\partial z} \\[12pt] \dfrac{\partial f_2}{\partial x} & \dfrac{\partial f_2}{\partial y} & \dfrac{\partial f_2}{\partial z} \\[12pt] \dfrac{\partial f_3}{\partial x} & \dfrac{\partial f_3}{\partial y} & \dfrac{\partial f_3}{\partial z} \end{bmatrix} \end{equation*}$

The Jacobian describes how small changes in the independent variables produce small changes in the dependent variables. In differential form, then

(97) $\begin{equation*} d\vec{f} = J\,d\vec{x} \end{equation*}$

where

(98) $\begin{equation*} d\vec{x} = \begin{bmatrix} dx \\[6pt] dy \\[6pt] dz \end{bmatrix} \end{equation*}$

This relationship is important in many engineering applications, including coordinate transformations, nonlinear equation solving, finite element methods, flight dynamics, robotics, and optimization.

For example, if

(99) $\begin{equation*} f_1=x^2+y, \qquad f_2=xy, \qquad f_3=z^2+x \end{equation*}$

then

(100) $\begin{equation*} J = \begin{bmatrix} 2x & 1 & 0 \\[6pt] y & x & 0 \\[6pt] 1 & 0 & 2z \end{bmatrix} \end{equation*}$

The determinant of the Jacobian is also important in coordinate transformations because it gives the local scale factor for area or volume elements.

Systems of Linear Equations

A system of linear equations can be written in matrix form as

(101) $\begin{equation*} A\vec{x} = \vec{b} \end{equation*}$

where $A$ is the coefficient matrix, $\vec{x}$ is the column vector of unknowns, and $\vec{b}$ is the right-hand-side vector.

Solving systems of linear equations is a fundamental task in engineering. The rank of a matrix is the number of linearly independent rows or columns. In solving systems of equations, rank helps determine whether the system has a unique solution, no solution, or infinitely many solutions. Row reduction is one practical way to determine the rank of a matrix. Different methods may be used depending on the size and properties of the coefficient matrix $A$ , as well as the purpose of the calculation. These methods include row reduction (also known as Gaussian elimination), matrix inversion, LU decomposition, and Cramer’s Rule. Each method has advantages and limitations.

Row Reduction (Gaussian Elimination)

Gaussian elimination solves a system of linear equations by applying elementary row operations to the augmented matrix. These operations do not change the solution of the system. The objective is to convert the system to row echelon form, or to reduced row echelon form, from which the unknowns can be found directly.

The basic procedure is to choose a nonzero pivot element, use row operations to eliminate the entries below the pivot, and then repeat the process for the remaining rows. Once an upper triangular form has been obtained, the variables are found by back substitution. In practical numerical computations, row exchanges may be used to avoid zero or very small pivot elements.

For example, to solve the system

(102) $\begin{equation*} 2x + y = 5 \quad \text{and} \quad 4x + 3y = 6 \end{equation*}$

the augmented matrix is

(103) $\begin{equation*} \begin{bmatrix} 2 & 1 & 5 \\[6pt] 4 & 3 & 6 \end{bmatrix} \end{equation*}$

Dividing the first row by 2 to make the pivot element in the first column equal to 1 gives

(104) $\begin{equation*} \begin{bmatrix} 1 & 0.5 & 2.5 \\[6pt] 4 & 3 & 6 \end{bmatrix} \end{equation*}$

Subtracting 4 times the first row from the second row gives

(105) $\begin{equation*} \begin{bmatrix} 1 & 0.5 & 2.5 \\[6pt] 0 & 1 & -4 \end{bmatrix} \end{equation*}$

The second row gives $y=-4$ . Substituting this value into the first row gives $x=4.5$ . Therefore,

(106) $\begin{equation*} x = 4.5 \quad \text{and} \quad y = -4 \end{equation*}$

Matrix Inversion

Inverse matrices arise in many engineering applications because they describe how one set of variables can be recovered from another. They appear in coordinate transformations, control systems, structural flexibility methods, least-squares estimation, covariance propagation, and many other problems where the inverse matrix itself has physical or computational significance. Matrix inversion can also be used to solve a system of linear algebraic equations.

A system of linear equations can be written in matrix form as

(107) $\begin{equation*} A\vec{x} = \vec{b} \end{equation*}$

where $A$ is the coefficient matrix, $\vec{x}$ is the vector of unknowns, and $\vec{b}$ is the right-hand-side vector. If the matrix $A$ is invertible, then $\det(A) \neq 0$ and the system has a unique solution. Multiplying both sides of $A\vec{x}=\vec{b}$ by $A^{-1}$ gives

(108) $\begin{equation*} A^{-1}A\vec{x} = A^{-1}\vec{b} \end{equation*}$

so that

(109) $\begin{equation*} \vec{x} = A^{-1} \, \vec{b} \end{equation*}$

If $\det(A)=0$ , then $A^{-1}$ does not exist, and the system does not have a unique solution.

For small matrices, the inverse may be found using the adjugate, or cofactor, method. In this approach,

(110) $\begin{equation*} A^{-1} = \frac{1}{\det(A)}\,\operatorname{adj}(A) \end{equation*}$

where $\operatorname{adj}(A)$ is the transpose of the cofactor matrix of $A$ . For a $2 \times 2$ matrix,

(111) $\begin{equation*} A = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \end{equation*}$

the determinant is

(112) $\begin{equation*} \det(A)=a \, d - b \, c \end{equation*}$

and the inverse is

(113) $\begin{equation*} A^{-1} = \frac{1}{a\,d-b\,c} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix} \end{equation*}$

provided that $a\,d-b\,c \neq 0$ . For larger matrices, the same idea can be extended using cofactors, but the arithmetic becomes lengthy.

The basic procedure is to check that $\det(A) \neq 0$ , find the inverse matrix $A^{-1}$ using cofactors, row reduction, or another suitable method, and then multiply $A^{-1}$ by $\vec{b}$ to obtain $\vec{x}$ . Although matrix inversion is straightforward in principle, explicitly forming $A^{-1}$ is not always the most efficient way to solve large systems of equations. If the only objective is to find $\vec{x}$ from $A\vec{x}=\vec{b}$ , then direct solution methods such as Gaussian elimination, LU decomposition, or other numerical algorithms are usually preferred. However, when the inverse matrix itself is needed, or when the same matrix must be applied repeatedly to many different right-hand-side vectors, forming or factoring the inverse relationship may be useful.

LU Decomposition

LU decomposition is a method for solving systems of linear algebraic equations by factoring the coefficient matrix into a lower-triangular matrix and an upper-triangular matrix. For a system written as

(114) $\begin{equation*} A\vec{x}=\vec{b} \end{equation*}$

the matrix $A$ is decomposed as

(115) $\begin{equation*} A = L \, U \end{equation*}$

where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix. The original system then becomes

(116) $\begin{equation*} L U\vec{x}=\vec{b} \end{equation*}$

If

(117) $\begin{equation*} U\vec{x}=\vec{y} \end{equation*}$

then the solution can be found in two steps. First, solve

(118) $\begin{equation*} L\vec{y}=\vec{b} \end{equation*}$

for $\vec{y}$ using forward substitution. Then solve

(119) $\begin{equation*} U\vec{x}=\vec{y} \end{equation*}$

for $\vec{x}$ using back substitution.

This method is especially useful in engineering computations when the same coefficient matrix $A$ must be used with several different right-hand-side vectors. In that case, the matrix is factored once, and the triangular systems are solved repeatedly for each new $\vec{b}$ . In practical numerical work, row pivoting is often used for stability so that the decomposition may be written as

(120) $\begin{equation*} P A = L U \end{equation*}$

where $P$ is a permutation matrix that accounts for row exchanges.

Consider, for example, the system

(121) $\begin{equation*} A = \begin{bmatrix} 2 & 3 \\[6pt] 4 & 7 \end{bmatrix} \quad \text{and} \quad \vec{b} = \begin{bmatrix} 5 \\[6pt] 10 \end{bmatrix} \end{equation*}$

One possible LU decomposition is

(122) $\begin{equation*} L = \begin{bmatrix} 1 & 0 \\[6pt] 2 & 1 \end{bmatrix} \quad \text{and} \quad U = \begin{bmatrix} 2 & 3 \\[6pt] 0 & 1 \end{bmatrix} \end{equation*}$

so that $A=L U$ . The first step is to solve $L\vec{y}=\vec{b}$ , i.e.,

(123) $\begin{equation*} \begin{bmatrix} 1 & 0 \\[6pt] 2 & 1 \end{bmatrix} \begin{bmatrix} y_1 \\[6pt] y_2 \end{bmatrix} = \begin{bmatrix} 5 \\[6pt] 10 \end{bmatrix} \end{equation*}$

The first equation gives $y_1=5$ . The second equation gives

(124) $\begin{equation*} 2y_1+y_2=10 \end{equation*}$

and so $y_2=0$ . Therefore,

(125) $\begin{equation*} \vec{y} = \begin{bmatrix} 5 \\[6pt] 0 \end{bmatrix} \end{equation*}$

The second step is to solve $U\vec{x}=\vec{y}$ , i.e.,

(126) $\begin{equation*} \begin{bmatrix} 2 & 3 \\[6pt] 0 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\[6pt] x_2 \end{bmatrix} = \begin{bmatrix} 5 \\[6pt] 0 \end{bmatrix} \end{equation*}$

The second equation gives $x_2=0$ . The first equation then gives

(127) $\begin{equation*} 2x_1+3x_2=5 \end{equation*}$

and so $x_1 =$ 2.5. Therefore, the solution is

(128) $\begin{equation*} \vec{x} = \begin{bmatrix} 2.5 \\[6pt] 0 \end{bmatrix} \end{equation*}$

Cramer’s Rule

Cramer’s Rule is a determinant-based method for solving a square system of linear equations. For the system

(129) $\begin{equation*} A\vec{x}=\vec{b} \end{equation*}$

each unknown $x_i$ is found from

(130) $\begin{equation*} x_i = \frac{\det(A_i)}{\det(A)} \end{equation*}$

where $A_i$ is obtained by replacing the $i$ -th column of $A$ with $\vec{b}$ . This method requires $\det(A)\neq 0$ ; otherwise, the system does not have a unique solution.

Cramer’s Rule is useful for small systems and for illustrating the relationship between determinants and linear systems. However, it is not usually used for large engineering systems because it requires computing several determinants, which becomes inefficient as the number of unknowns increases.

For example, if

(131) $\begin{equation*} A = \begin{bmatrix} 2 & 3 \\[4pt] 1 & 4 \end{bmatrix} \quad \text{and} \quad \vec{b} = \begin{bmatrix} 5 \\[4pt] 6 \end{bmatrix} \end{equation*}$

then

(132) $\begin{equation*} \det(A) = \begin{vmatrix} 2 & 3 \\[4pt] 1 & 4 \end{vmatrix} = 8 - 3 = 5 \end{equation*}$

The solution for $x_1$ is

(133) $\begin{equation*} x_1 = \frac{ \det \begin{bmatrix} 5 & 3 \\[4pt] 6 & 4 \end{bmatrix} } { \det \begin{bmatrix} 2 & 3 \\[4pt] 1 & 4 \end{bmatrix} } = \frac{20-18}{5} = \frac{2}{5} \end{equation*}$

and the solution for $x_2$ is

(134) $\begin{equation*} x_2 = \frac{ \det \begin{bmatrix} 2 & 5 \\[4pt] 1 & 6 \end{bmatrix} } { \det \begin{bmatrix} 2 & 3 \\[4pt] 1 & 4 \end{bmatrix} } = \frac{12-5}{5} = \frac{7}{5} \end{equation*}$

Therefore,

(135) $\begin{equation*} \vec{x} = \begin{bmatrix} 2/5 \\[6pt] 7/5 \end{bmatrix} \end{equation*}$

Check Your Understanding #15 – Solving a set of linear simultaneous equations

Solve the system of equations given by $2x + 3y = 5$ and $4x - y = 6$ .

Show solution/hide solution.

The system of equations is

$2x + 3y = 5 \quad \text{and} \quad 4x - y = 6$

In matrix form, these equations can be written as

$\begin{bmatrix} 2 & 3 \\[4pt] 4 & -1 \end{bmatrix} \begin{bmatrix} x \\[4pt] y \end{bmatrix} = \begin{bmatrix} 5 \\[5pt] 6 \end{bmatrix}$

The determinant of the coefficient matrix $A$ is

$\det(A) = \begin{vmatrix} 2 & 3 \\[4pt] 4 & -1 \end{vmatrix} = (2)(-1) - (4)(3) = -2 - 12 = -14$

The inverse of a $2 \times 2$ matrix is

$A^{-1} = \frac{1}{\det(A)} \begin{bmatrix} d & -b \\[4pt] -c & a \end{bmatrix}$

Substituting the values gives

$A^{-1} = \frac{1}{-14} \begin{bmatrix} -1 & -3 \\[4pt] -4 & 2 \end{bmatrix} = \begin{bmatrix} \dfrac{1}{14} & \dfrac{3}{14} \\[10pt] \dfrac{2}{7} & -\dfrac{1}{7} \end{bmatrix}$

To solve for the unknown vector $\vec{x}$ , use

$\vec{x} = A^{-1}\vec{b}$

Substituting gives

$\vec{x} = \begin{bmatrix} \dfrac{1}{14} & \dfrac{3}{14} \\[10pt] \dfrac{2}{7} & -\dfrac{1}{7} \end{bmatrix} \begin{bmatrix} 5 \\[12pt] 6 \end{bmatrix}$

Performing the matrix multiplication gives

$\vec{x} = \begin{bmatrix} \dfrac{23}{14} \\[12pt] \dfrac{4}{7} \end{bmatrix}$

Therefore, the solution to the system is

$x = \dfrac{23}{14} \quad \text{and} \quad y = \dfrac{4}{7}$

Complex Numbers

A complex number $z = x + i\,y$ combines a “real” part ${x}$ and an “imaginary” part ${y}$ . Engineering work can leverage complex number theory to simplify the mathematics. The use of complex numbers appears throughout core engineering classes because they combine magnitude and direction into a single quantity, simplifying the bookkeeping for differential equations. All of this reduces to ordinary algebra or trigonometry once the rules, such as $i^{2} = -1$ , are kept in mind.

The addition and subtraction of complex numbers combine the like parts, i.e.,

(136) $\begin{equation*} (x + i\,y) \pm (x' + i\,y') = (x \pm x') + i\,(y \pm y') \end{equation*}$

The multiplication of complex numbers follows the distributive law of algebra. Collecting real and imaginary terms gives

(137) $\begin{equation*} (x + i\,y) \, (x' + i\,y') = (x\,x' - y\,y') + i\,(x\,y' + y\,x') \end{equation*}$

Complex numbers are often expressed in their polar (exponential) form, i.e.,

(138) $\begin{equation*} z = r \,e^{i\,\theta} \end{equation*}$

where the modulus $r = |z|$ sets the length of the vector and the exponential factor $e^{i\,\theta}$ rotates the vector through an angle $\theta$ . If

(139) $\begin{equation*} z = r \,e^{i\,\theta} \quad \text{and} \quad z' = r' \,e^{i\,\theta'} \end{equation*}$

then their product becomes

(140) $\begin{equation*} z\,z' = r\,r' \,e^{i\,(\theta + \theta')} \end{equation*}$

Therefore, their magnitudes multiply, and their arguments add.

Division is equally simple in polar form, but in Cartesian form it is usually performed by multiplying the numerator and denominator by the complex conjugate, i.e., $z^{\!*} = x - i\,y$ , so that

(141) $\begin{equation*} \frac{z_{1}}{z_{2}} = \frac{z_{1} \, z_{2}^{\!*}}{|z_{2}|^{2}} \end{equation*}$

where the modulus is

(142) $\begin{equation*} |z| = \sqrt{x^{2} + y^{2}} \end{equation*}$

The argument (phase) is defined by

(143) $\begin{equation*} \theta = \arg z = \operatorname{atan2}(y,x) \end{equation*}$

so $z$ can be reconstructed from its modulus and argument as

(144) $\begin{equation*} z = |z| \,e^{i\,\theta} \quad \text{and} \quad e^{i\,\theta} = \cos \theta + i \,\sin \theta \end{equation*}$

which is the famous Euler formula, condensing the cosine and sine components into a single exponential function. The special case $e^{i\pi}+1=0$ is known as Euler’s identity.

Powers and roots of complex numbers follow directly. In polar form, then

(145) $\begin{equation*} (r\,e^{i\,\theta})^{n} = r^{n} \, e^{i\,n\theta} \end{equation*}$

so raising a complex number to the $n$ th power multiplies its modulus by itself $n$ times and adds its argument $n$ times. Correspondingly, the $n$ th roots are

(146) $\begin{equation*} z^{1/n} = |z|^{1/n} \, e^{i\,(\theta + 2\pi k)/n},\quad k = 0,1,\dots ,n-1 \end{equation*}$

thereby revealing the equally spaced root pattern commonly encountered in solving differential equations and various engineering problems. Aerospace engineers often first encounter complex number theory in the study of dynamics and control, as well as in two-dimensional potential flows. In advanced courses, they appear in descriptions of unsteady aerodynamics and aeroelasticity.

Mastering the use of complex numbers largely involves recognizing when to switch between Cartesian and polar forms. For example, a lumped spring-mass-damper system obeys the canonical form

(147) $\begin{equation*} m\overbigddot{x} + c\overbigdot{x} + k x = F_0 \, e^{i\omega \, t} \end{equation*}$

Assuming a complex displacement $x(t) =\hat{x} \, e^{i\omega \, t}$ converts the differential equation into the algebraic relation

(148) $\begin{equation*} \hat{x} = \dfrac{F_0}{\left( k - m\omega^{2} + i \, c\omega \right)} \end{equation*}$

The complex denominator reveals resonance when its magnitude is at its minimum, and its phase indicates how far the mass lags the forcing.

In aerodynamics, complex potentials $w(z) = \phi(x,y) + i \psi(x,y)$ can help solve two-dimensional, incompressible, irrotational “potential” flows because both the velocity potential $\phi$ and the stream function $\psi$ satisfy Laplace’s equation, i.e., $\nabla^2 \phi = 0$ and $\nabla^2 \psi = 0$ . This allows the entire flow field to be represented by a single analytic function of the complex variable $z = x + i \, y$ , thereby significantly simplifying the analysis of potential flow.

Quaternions

Quaternions are an extension of complex numbers that can be used to represent rotations in three-dimensional space. A complex number has one real part and one imaginary part, i.e., in the form $z = x + i\,y$ . A quaternion has one real part and three imaginary parts and is written as

(149) $\begin{equation*} q = q_0 + q_1\,i + q_2\,j + q_3\,k \end{equation*}$

where $q_0$ is the scalar part and $q_1\,i + q_2\,j + q_3\,k$ is the vector part. The quantities $i$ , $j$ , and $k$ obey special multiplication rules, including

(150) $\begin{equation*} i^2 = j^2 = k^2 = i\,j\,k = -1 \end{equation*}$

These rules make quaternion multiplication order-dependent, which is useful for describing rotations because finite rotations in three-dimensional space are also order-dependent. For example, rotating an object about one axis and then about another axis generally does not yield the same final orientation as applying the same two rotations in the opposite order.

A rotation through an angle $\theta$ about a unit vector $\vec{n}$ can be represented by the unit quaternion

(151) $\begin{equation*} q = \cos\frac{\theta}{2} + \vec{n}\,\sin\frac{\theta}{2} \end{equation*}$

where ${\vec{n}}$ specifies the axis of rotation. The appearance of the half-angle $\theta/2$ is a characteristic feature of the quaternion representation of rotations.

Although quaternions are less intuitive than ordinary angles, they are very useful in numerical simulations because they can represent arbitrary three-dimensional orientations compactly and robustly. In aerospace engineering, quaternions are used in spacecraft attitude dynamics, inertial navigation systems, autopilots, flight simulators, and full six-degree-of-freedom simulations of aircraft, rotorcraft, missiles, and uncrewed aerial vehicles. They are especially valuable when a vehicle may rotate through large angles or assume orientations for which ordinary angle-based descriptions become inconvenient or poorly conditioned. If the orientation must be expressed in terms of ordinary angles, the quaternion can be converted back to roll, pitch, and yaw angles using trigonometric relationships.

Eigenvalues, Eigenvectors, & Eigenfunctions

In many areas of aerospace engineering, particularly in dynamics, vibrations, stability analysis, and fluid mechanics, mathematical solutions naturally yield eigenvalue formulations. These arise whenever a linear system can be expressed in the form

(152) $\begin{equation*} A\,\vec{x} = \lambda\,\vec{x} \end{equation*}$

where $A$ is a square matrix, $\vec{x}$ is a nonzero vector, and $\lambda$ is a scalar. The scalar $\lambda$ is called an eigenvalue, and the corresponding vector $\vec{x}$ is the eigenvector.^[6] An “eigenvalue” is literally a value that is intrinsic to or belongs to a given transformation, and an “eigenvector” (or eigenfunction) is the shape or direction that the transformation leaves unchanged except for scaling. Physically, this relation means that the action of the linear operator $A$ on the vector $\vec{x}$ produces a scaled version of the same vector rather than changing its direction. The set of all eigenvalues of $A$ forms the spectrum of the system. For an $n\times n$ matrix, the eigenvalues are found from the characteristic equation given by

(153) $\begin{equation*} \det(A - \lambda I) = 0 \end{equation*}$

whose $n$ roots give the possible values of $\lambda$ . Each eigenvalue $\lambda_i$ has one or more associated eigenvectors $\vec{x}_i$ .

A classic example illustrating the appearance of eigenvalues and eigenvectors is the determination of the natural frequencies of a spring–mass system. For a system of $n$ discrete masses, Newton’s second law gives a set of coupled linear equations, i.e.,

(154) $\begin{equation*} M\,\overbigddot{\vec{x}} + K\,\vec{x} = \vec{0} \end{equation*}$

where $M$ is the $n\times n$ mass matrix, $K$ is the $n\times n$ stiffness matrix, and $\vec{x}(t)$ is the vector of displacements measured from an equilibrium position. If harmonic motion is assumed in the form of

(155) $\begin{equation*} \vec{x}(t) = \vec{\phi}\,e^{i\omega \, t} \end{equation*}$

then $\overbigddot{\vec{x}} = -\,\omega^2 \vec{\phi}\,e^{i\omega \, t}$ and substitution gives

(156) $\begin{equation*} (K - \omega^2 M)\,\vec{\phi} = \vec{0} \end{equation*}$

This equation represents an eigenvalue problem. The scalar $\omega^2$ acts as the eigenvalue, and the corresponding vector $\vec{\phi}$ is the eigenvector. In vibration problems, this eigenvector defines the relative motion of the masses and is called a mode shape. Each eigenvalue is the square of a natural frequency, i.e., $\lambda=\omega^2$ , and each eigenvector describes the shape of vibration associated with that frequency. Therefore, finding the normal modes of a mechanical system is mathematically equivalent to finding the eigenvalues and eigenvectors of its governing matrix equations.

In a continuous system, the motion of a system is described not by a finite set of coordinates $x_i(t)$ but by a field variable that varies with position and time, such as the displacement of a vibrating string or beam, or the velocity or pressure field in a fluid flow. Here, the concept of eigenvalues and eigenvectors extends naturally to eigenfunctions, which describe the spatial modes of the system. As an example, consider the one-dimensional wave equation given by

(157) $\begin{equation*} \frac{\partial^2 \phi}{\partial t^2} = a^2\,\frac{\partial^2 \phi}{\partial x^2} \end{equation*}$

where $\phi(x,t)$ represents a dependent variable such as displacement or pressure fluctuation, $x$ is the spatial coordinate, $t$ is time, and $a$ is the wave propagation speed, e.g., the speed of sound in air. A standard method of solution is to assume that the spatial and temporal dependences are separable, i.e.,

(158) $\begin{equation*} \phi(x,t) = X(x)\,T(t) \end{equation*}$

where $X(x)$ is a purely spatial function and $T(t)$ is a purely temporal function. Substituting this form into the governing equation gives

(159) $\begin{equation*} X(x)\,\frac{d^2 T}{dt^2} = a^2\,T(t)\,\frac{d^2 X}{dx^2} \end{equation*}$

Dividing both sides by $a^2\,X(x)\,T(t)$ gives

(160) $\begin{equation*} \frac{1}{a^2}\,\frac{1}{T}\,\frac{d^2 T}{dt^2} = \frac{1}{X}\,\frac{d^2 X}{dx^2} = -\,\lambda \end{equation*}$

where the common constant of separation, $-\lambda$ , must be introduced because the left-hand side depends only on $t$ and the right-hand side only on $x$ . This separation produces the two ordinary differential equations

(161) $\begin{equation*} \frac{d^2 X}{dx^2} + \lambda X = 0 \quad \text{and} \quad \frac{d^2 T}{dt^2} + \lambda a^2 T = 0 \end{equation*}$

Here, $\lambda$ is the eigenvalue, and the corresponding spatial function $X(x)$ is the eigenfunction, which are then determined by the boundary conditions of the problem.

For example, for a string of length $L$ fixed in displacement at both ends ( $X = 0$ at $x = 0$ and $x = L$ ), the eigenfunctions are

(162) $\begin{equation*} X_n(x) = \sin \left(\frac{n\pi x}{L}\right) \end{equation*}$

with eigenvalues $\lambda_n = \left(\dfrac{n\pi}{L}\right)^2$ . The corresponding temporal functions are harmonic, i.e.,

(163) $\begin{equation*} T_n(t) = A_n \cos(\omega_n \, t) + B_n \sin(\omega_n \, t) \end{equation*}$

where $\omega_n = a\sqrt{\lambda_n}$ is the natural angular frequency of the $n$ th mode. The mode shapes of a string fixed at both ends, corresponding to the first four eigenfunctions, $X_1(x)$ , $X_2(x)$ , $X_3(x)$ , and $X_4(x)$ , are shown in Figure 14. Each mode satisfies the boundary conditions that $X(0) = X(L) = 0$ , with standing wavelengths $\Lambda_n = 2L/n$ and natural frequencies $\omega_n = a \, n \, \pi / L$ .

Mode shapes of a fixed–fixed string corresponding to the first four modes.

Because the eigenfunctions $\{X_n(x)\}$ form an orthogonal set over the domain $0 \le x \le L$ , any physically allowable solution can be expressed as a linear combination (or series expansion) of these modes, i.e.,

(164) $\begin{equation*} \phi(x,t) = \sum_{n=1}^{\infty} a_n\,X_n(x)\,T_n(t) \end{equation*}$

These spatial eigenfunctions allow any vibration pattern to be expressed as a linear combination of the normal modes. This principle underpins modal analysis in fields such as structural dynamics, aeroelasticity, and acoustics, where the overall response of a system is represented as a superposition of independent modes, each associated with its own eigenvalue and eigenfunction.

Check Your Understanding #16 – Eigenvalues of a two-mass system

Two identical masses of $m_1 = m_2 = 1\,\text{kg}$ are connected by three springs, as shown in the figure below. Each outer spring has stiffness $k$ =100 N/m, and the center spring has a stiffness of $k_c$ =50 N/m. Let $\vec{x} = [\,x_1\;\;x_2\,]^T$ represent the small horizontal displacements of the masses. The equations of motion are given by

$M\,\overbigddot{\vec{x}} + K\,\vec{x} = 0 \qquad M = I, \quad K = \begin{bmatrix}150 & -50 \\[3pt] -50 & 150 \end{bmatrix}$

Find the natural frequencies and mode shapes of the system.

Show solution / hide solution.

Assuming harmonic motions of the form $\vec{x}=\vec{\phi}e^{i\omega \, t}$ gives

$\left( K - \omega^2 M \right) \,\vec{\phi} = 0$

which is an eigenvalue problem with $\lambda=\omega^2$ . The characteristic equation is

$\det(K-\lambda I) = \begin{vmatrix} 150-\lambda & -50\\[3pt] -50 & 150-\lambda \end{vmatrix} = 0$

giving

$(\lambda-100)(\lambda-200)=0$

Hence the eigenvalues are $\lambda_1$ =100 and $\lambda_2$ = 200, and the corresponding natural frequencies are $\omega_1 = \sqrt{100} = 10$ rad/s and $\omega_2 = \sqrt{200} = 14.14$ rad/s.

To find the eigenvectors, substitute each $\lambda$ into $(K-\lambda I)\vec{\phi} = 0$ . For $\lambda_1$ = 100 then $\vec{\phi}_1 = [\,1\;\;1\,]^T$ , i.e., the masses move in phase. For $\lambda_2$ =200 then $\vec{\phi}_2 = [\,1\;\;-1\,]^T$ , i.e., the masses move out of phase.

Notice that the two eigenvalues correspond to the squares of the natural frequencies, and the eigenvectors describe the relative motion in each mode. Mode 1 has both masses moving in phase, and Mode 2 has them moving out of phase.

Ordinary Differential Equations (ODEs)

An ordinary differential equation (ODE) is an equation that involves the derivatives of a function with respect to a variable. A common engineering problem is solving an ordinary differential equation (ODE), i.e., determining the function (or functions) that will satisfy a given differential equation or a set of such equations. For example, given the ODE

(165) $\begin{equation*} \frac{dx}{dt} = \sin t \end{equation*}$

then what will be the function $x(t)$ ? The antiderivative of $\sin t$ is $-\cos t$ , so

(166) $\begin{equation*} x(t) = -\cos t + c \end{equation*}$

where ${c}$ is some unknown constant. Therefore, solving an ODE is more complicated than just anti-differentiation because the value of the constant ${c}$ must also be determined. To this end, additional information, such as an initial or boundary condition, is required. Commonly, the value of $x(t)$ will be known at some specific value of $t$ , e.g., $x(t_0) = x_0$ . In this case, $x_0 = -\cos t_0 + c$ , so $c = x_0 + \cos t_0$ . Therefore, the particular solution is

(167) $\begin{equation*} x(t) = x_0 + \cos t_0 - \cos t \end{equation*}$

for all values of $t$ .

Check Your Understanding #17 – Solving ODEs

Use the method of separation of variables to solve
$\frac{dx}{dt} = 5x - 3$

Determine the function $x(t)$ . Take the particular case $x(2)=1$ to evaluate the integration constant.
Solve the first-order linear ordinary differential equation
$\frac{dy}{dt} + 3y = 6e^{-2t}$

subject to the initial condition $y(0)=1$ .
Find the general solution of the second-order constant-coefficient ODE
$\frac{d^{2}z}{dt^{2}} + 4z = 0$

Show solution / hide solution.

Separate the variables to get
$\frac{dx}{5x-3} = dt$

Integrating gives

$\frac{1}{5}\ln|5x-3| = t + c$

so that

$5x - 3 = K e^{5t}$

Hence,

$x(t) = \frac{3 + K e^{5t}}{5}$

Use the result that $x(2) = 1$ to find $K$ . Substituting into $5x - 3 = K e^{5t}$ gives $5(1)-3 = K e^{10}$ , so $K = 2 e^{-10}$ . Therefore

$x(t) = \frac{3}{5} + \frac{2}{5}\,e^{5(t-2)}$
The integrating factor is
$\mu(t) = e^{\int 3\,dt} = e^{3t}$

Multiply through the equation to obtain

$e^{3t} \frac{dy}{dt} + 3e^{3t} y = 6e^{(3-2)t}$

so that

$\frac{d}{dt} \bigl(e^{3t}y\bigr) = 6e^{t}$

Integrating gives

$e^{3t}y = 6e^{t} + c$

so that

$y(t) = 6e^{-2t} + c e^{-3t}$

Apply $y(0) = 1$ to get $1 = 6 + c$ , i.e., $c = -5$ . Therefore,

$y(t) = 6e^{-2t} - 5e^{-3t}$
The characteristic equation is
$r^{2} + 4 = 0$

so that $r = \pm 2i$ . Therefore, the general solution is

$z(t) = A \cos(2t) + B \sin(2t)$

where $A$ and $B$ are arbitrary constants.

State-Space Form

Many engineering systems are modeled by differential equations, often of second or higher order. To analyze such systems systematically, it is convenient to rewrite them in state-space form, which expresses the dynamics as a set of coupled first-order equations. A linear, time-invariant system can be written as

(168) $\begin{equation*} \overbigdot{\vec{x}}(t) = A\,\vec{x}(t) + B\,\vec{u}(t) \qquad \text{and} \qquad \vec{y}(t) = C\,\vec{x}(t) + D\,\vec{u}(t) \end{equation*}$

Here, $\vec{x}(t)$ is the state vector, $\vec{u}(t)$ is the input vector, and $\vec{y}(t)$ is the output vector. The matrices $A$ , $B$ , $C$ , and $D$ have dimensions that depend on the number of states, inputs, and outputs. For a system with $n$ states, $m$ inputs, and $p$ outputs, $A$ is $n \times n$ , $B$ is $n \times m$ , $C$ is $p \times n$ , and $D$ is $p \times m$ . The first equation, known as the state equation, determines how the system evolves in time, while the second equation, called the output equation, relates the states and inputs to the outputs.

Consider a linear differential equation of order $n$ , i.e.,

(169) $\begin{equation*} a_n y^{(n)} + a_{n-1} y^{(n-1)} + \cdots + a_1 \overbigdot y + a_0 y = b_0 u \end{equation*}$

where $a_n \ne 0$ . Defining the state variables as successive derivatives of the output gives

(170) $\begin{equation*} x_1 = y,\quad x_2 = \overbigdot y,\quad \ldots,\quad x_n = y^{(n-1)} \end{equation*}$

This definition leads to the chain of first-order equations

(171) $\begin{equation*} \overbigdot x_1 = x_2,\quad \overbigdot x_2 = x_3,\quad \ldots,\quad \overbigdot x_{n-1} = x_n \end{equation*}$

The last equation comes from solving the original differential equation for the highest derivative, i.e.,

(172) $\begin{equation*} \overbigdot x_n = -\frac{a_{n-1}}{a_n}x_n - \cdots - \frac{a_1}{a_n} x_2 - \frac{a_0}{a_n} x_1 + \frac{b_0}{a_n} u \end{equation*}$

Collecting these equations in matrix form gives

(173) $\begin{equation*} \overbigdot{\vec{x}} = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ -\dfrac{a_0}{a_n} & -\dfrac{a_1}{a_n} & -\dfrac{a_2}{a_n} & \cdots & -\dfrac{a_{n-1}}{a_n} \end{bmatrix} \vec{x} + \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \dfrac{b_0}{a_n} \end{bmatrix} u \end{equation*}$

For this particular choice of states, if the output is the original dependent variable $y=x_1$ , then

(174) $\begin{equation*} y = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \end{bmatrix} \vec{x} \end{equation*}$

so that $C = [\,1 \quad 0 \quad 0 \quad \cdots \quad 0\,]$ and $D=0$ .

The solution of the state-space equations follows directly. For the homogeneous system with zero input, i.e.,

(175) $\begin{equation*} \overbigdot{\vec{x}}(t) = A\,\vec{x}(t) \end{equation*}$

the solution is

(176) $\begin{equation*} \vec{x}(t) = e^{At}\,\vec{x}(0) \end{equation*}$

where $e^{At}$ is the state-transition matrix. In the presence of inputs, the complete solution is

(177) $\begin{equation*} \vec{x}(t) = e^{At}\,\vec{x}(0) + \int_0^t e^{A(t-\tau)} B\,\vec{u}(\tau)\, d\tau \end{equation*}$

This second term is known as Duhamel’s integral or the convolution integral, which represents the forced response obtained by superposing the effects of the input over time. The corresponding output is

(178) $\begin{equation*} \vec{y}(t) = C\,\vec{x}(t) + D\,\vec{u}(t) \end{equation*}$

Therefore, the eigenvalues of $A$ govern the natural modes of the system, while Duhamel’s integral describes the forced response to the input. This solution framework is widely used in both control theory and flight dynamics, where the state-space form offers a compact and powerful representation of system behavior.

In more general cases, the system equations may be nonlinear and so will take the form

(179) $\begin{equation*} \overbigdot{\vec{x}} = \vec{f}(\vec{x},\vec{u},t) \quad \text{and} \quad \vec{y} = \vec{g}(\vec{x},\vec{u},t) \end{equation*}$

Around an equilibrium or reference point $(\vec{x}_0,\vec{u}_0)$ , a first-order Taylor expansion gives the linearized state-space equations, i.e.,

(180) $\begin{equation*} \delta\overbigdot{\vec{x}} = A\,\delta\vec{x} + B\,\delta\vec{u},\qquad \delta\vec{y} = C\,\delta\vec{x} + D\,\delta\vec{u} \end{equation*}$

with the Jacobian matrices

(181) $\begin{equation*} A = \left.\frac{\partial \vec{f}}{\partial \vec{x}}\right|_{(\vec{x}_0,\vec{u}_0)},\qquad B = \left.\frac{\partial \vec{f}}{\partial \vec{u}}\right|_{(\vec{x}_0,\vec{u}_0)},\qquad C = \left.\frac{\partial \vec{g}}{\partial \vec{x}}\right|_{(\vec{x}_0,\vec{u}_0)},\quad \text{and} \quad D = \left.\frac{\partial \vec{g}}{\partial \vec{u}}\right|_{(\vec{x}_0,\vec{u}_0)} \end{equation*}$

Check Your Understanding #18 – Mass-spring-damper in state-space form

Consider a classic mass-spring-damper system governed by

$m \, \overbigddot q + c \, \overbigdot q + k \, q = u$

First, put this governing equation into state-space form by defining the appropriate state variables. Then find the eigenvalues of the state matrix and express them in terms of the natural frequency ${\omega_n}$ and damping ratio ${\zeta}$ . State the damped natural frequency for the underdamped case.

Show solution / hide solution.

Define the states as $x_1 = q$ and $x_2 = \overbigdot q$ . Then from the definitions $\overbigdot{x}_1 = x_2$ . From the ODE, then solve for the highest derivative, i.e.,

$\overbigddot q = \frac{1}{m}\big(u - c\,\overbigdot q - k\,q\big)$

Substituting $x_1, x_2$ gives

$\overbigdot{x}_2 = -\frac{c}{m}\,x_2 - \frac{k}{m}\,x_1 + \frac{1}{m}\,u$

With the states as

$\vec{x}=\begin{bmatrix}x_1\\[2pt]x_2\end{bmatrix}$

then

$\overbigdot{\vec{x}} = \begin{bmatrix} 0 & 1 \\[8pt] -\dfrac{k}{m} & -\dfrac{c}{m} \end{bmatrix}\vec{x} + \begin{bmatrix} 0 \\ \dfrac{1}{m} \end{bmatrix}u$

If the output is the displacement, $y=q$ , then

$y=\begin{bmatrix}1 & 0\end{bmatrix}\vec{x}$

Therefore, the system matrices are

$A = \begin{bmatrix} 0 & 1 \\[8pt] -\dfrac{k}{m} & -\dfrac{c}{m} \end{bmatrix}, \qquad B = \begin{bmatrix} 0 \\[8pt] \dfrac{1}{m} \end{bmatrix}, \qquad C = \begin{bmatrix} 1 & 0 \end{bmatrix}, \quad \text{and} \quad D = 0$

The characteristic polynomial of $A$ is

$\det(\lambda I-A)=\lambda^{2}+\frac{c}{m}\lambda+\frac{k}{m} = 0$

Now the natural frequency and damping ratio become

$\omega_n =\sqrt{\frac{k}{m}} \quad \text{and} \quad \zeta=\frac{c}{2\sqrt{k \, m}}=\frac{c}{2 \, m\,\omega_n}$

Then the eigenvalues are

$\lambda_{1,2}=\omega_n\!\left(-\zeta \pm \sqrt{\zeta^{2}-1}\right)$

For the underdamped case, i.e., $0 < \zeta < 1$ , the eigenvalues are a complex-conjugate pair,

$\lambda_{1,2} = -\zeta\,\omega_n \pm i\,\omega_n\sqrt{1-\zeta^2}$

The damped natural frequency is

$\omega_d = \omega_n\sqrt{1-\zeta^2}$

For the critically damped case, $\zeta=1$ , the eigenvalues reduce to a repeated real negative root.

Laplace Transforms

The Laplace transform is a mathematical technique that converts a function of time, $f(t)$ , into a function of a complex variable $s$ , thereby linking the time domain with the frequency domain. This transformation simplifies the analysis of dynamic systems by converting differential equations into algebraic equations that are easier to manipulate and solve. In the $s$ -domain, a system’s behavior can be expressed in terms of its natural frequencies, damping, and stability characteristics, providing clear insight into its response to external inputs. After solving the algebraic equation in the transformed domain, the inverse Laplace transform is applied to obtain the corresponding time-domain solution.

The Laplace transform of a time-dependent function $f(t)$ is defined as

(182) $\begin{equation*} \mathcal{L}\{f(t)\} = F(s) = \int_{0}^{\infty} f(t)\,e^{-st}\,dt \end{equation*}$

where $s$ is a complex variable written as $s = \sigma + i\omega$ . The exponential factor $e^{-st}$ acts as a weighting function that measures how strongly each part of $f(t)$ contributes to different rates of exponential decay and oscillation. The result, i.e., $F(s)$ , describes the behavior of $f(t)$ in the combined frequency-growth domain. For the transform to exist, the function $f(t)$ must not grow faster than an exponential function as $t \to \infty$ .

The behavior of many engineering systems is described by differential equations that express how some quantity changes with time. The Laplace transform allows these time-dependent equations to be expressed and solved in terms of the complex variable $s$ , which combines both frequency and damping. In the transformed domain, differentiation in time produces multiplication by $s$ together with terms involving the initial conditions, while integration in time corresponds to division by $s$ when the integral is taken from zero. Applying the Laplace transform to a differential equation such as

(183) $\begin{equation*} a_2\frac{d^2x}{dt^2} + a_1\frac{dx}{dt} + a_0x = f(t) \end{equation*}$

uses the properties

(184) $\begin{equation*} \mathcal{L}\!\left\{\frac{dx}{dt}\right\} = sX(s)-x(0) \qquad \text{and} \qquad \mathcal{L}\!\left\{\frac{d^2x}{dt^2}\right\} = s^2X(s)-s\,x(0)-\overbigdot{x}(0) \end{equation*}$

so that each time derivative of $x(t)$ produces powers of $s$ multiplying $X(s)$ , together with the corresponding initial-condition terms. Substituting these relations gives

(185) $\begin{equation*} (a_2 s^2 + a_1 s + a_0)X(s) = F(s) + a_2[sx(0)+\overbigdot{x}(0)] + a_1x(0) \end{equation*}$

which is an algebraic equation in $s$ . This simplification arises because the exponential kernel $e^{-st}$ in the Laplace transform converts differentiation into algebraic terms in $s$ upon integration by parts, with additional terms that account for the initial conditions. Once $X(s)$ is obtained algebraically, the original time-domain solution follows from the inverse transform $\mathcal{L}^{-1}$ , i.e.,

(186) $\begin{equation*} x(t) = \mathcal{L}^{-1}\{X(s)\} \end{equation*}$

In most practical cases, this is done by expanding $X(s)$ into partial fractions and consulting tables of standard transforms.

The Laplace transform is used widely in engineering analysis to describe the behavior of dynamic systems, including mechanical vibrations, electrical circuits, and control systems. In flight dynamics, for example, the linearized equations of motion are often written in the Laplace or $s$ -domain so that stability and transient response can be examined in terms of system poles and damping.

Check Your Understanding #19 – Solving a first-order ODE using Laplace transforms

Solve the first-order ordinary differential equation

$\frac{dx}{dt} + a\,x(t) = b\,u(t), \quad \text{with} \quad x(0) = x_0$

where $a > 0$ and $b$ are constants, and $u(t)$ is the unit step input.

Show solution / hide solution.

Taking Laplace transforms gives

$\mathcal{L}\!\left\{\frac{dx}{dt}\right\} + a\,\mathcal{L}\{x\} = \mathcal{L}\{b\,u(t)\}$

so that

$(sX(s) - x_0) + aX(s) = \frac{b}{s}$

Solving for $X(s)$ gives

$X(s) = \frac{x_0}{s+a} + \frac{b}{s(s+a)}.$

Using partial fractions, then

$\frac{b}{s(s+a)} = \frac{b}{a}\!\left(\frac{1}{s} - \frac{1}{s+a}\right)$

so that

$X(s) = \frac{b}{a}\frac{1}{s} + \left(x_0 - \frac{b}{a}\right)\frac{1}{s+a}$

Taking the inverse Laplace transform gives

$x(t) = \frac{b}{a} + \left(x_0 - \frac{b}{a}\right)e^{-a \, t}$

This result represents a step response that exponentially approaches the steady-state value $b/a$ with a time constant of $1/a$ , starting from the initial condition $x_0$ .

Fourier Series

While the Laplace transform is useful for analyzing the transient and general time-dependent behavior of a system, the Fourier series applies specifically to periodic functions, expressing them as sums of harmonic sinusoids in the frequency domain. This is a very important concept in engineering. A periodic function $f(t)$ with period $T$ can be represented as a sum of sine and cosine terms whose frequencies are integer multiples of the fundamental frequency $\omega_0 = 2\pi/T$ . In practical use, the Fourier series is usually approximated by retaining a finite number of harmonics, i.e.,

(187) $\begin{equation*} f_J(t) = a_0 + \sum_{j=1}^{J} \left( a_j \cos(j \, \omega_0 \, t) + b_j \sin(j \, \omega_0 \, t)\right) \end{equation*}$

where $J$ is the number of harmonics retained in the reconstruction. The coefficients determine the mean value and the contribution of each harmonic, and are given by

(188) $\begin{equation*} a_0 = \frac{1}{T}\int_{0}^{T} f(t)\,dt \end{equation*}$

and, for $j \ge 1$ , then

(189) $\begin{equation*} a_j = \frac{2}{T}\int_{0}^{T} f(t)\cos(j \, \omega_0 \, t)\,dt \qquad \text{and} \qquad b_j = \frac{2}{T}\int_{0}^{T} f(t)\sin(j \, \omega_0 \, t)\,dt \end{equation*}$

The quantity $f_J(t)$ is a finite approximation to the periodic function $f(t)$ . Increasing $J$ adds more harmonics and generally improves the reconstruction. In the limiting case, as more and more harmonics are retained, the Fourier series approaches the original periodic function, subject to the usual convergence behavior near discontinuities.

The Fourier series shows that a periodic signal can be expressed as a superposition of harmonically related sinusoids. Each term corresponds to a distinct frequency component, and the relative magnitudes and phases of these components describe the frequency content of the signal. The energy or power associated with each component is related to the squared magnitude of its coefficient, with the precise normalization depending on the Fourier-series convention used.

A useful way to interpret the Fourier series is to examine how a periodic waveform is reconstructed as more harmonics are retained. Figure 15 shows this process for a square wave. The first harmonic gives only a smooth sinusoidal approximation, but additional odd harmonics progressively sharpen the waveform and bring the reconstruction closer to the original periodic function. The remaining oscillations near the jump discontinuities are not plotting errors; they are a characteristic feature of Fourier reconstruction of discontinuous functions, known as the Gibbs phenomenon.

Fourier-series reconstruction of a square wave — Fourier-series reconstruction of a periodic square wave using increasing numbers of odd harmonics.

In practical applications, signals are often available only as discrete samples rather than as continuous functions of time. Examples include measured acoustic pressure, structural vibration, rotor loads, or unsteady aerodynamic data recorded at equally spaced time intervals. Suppose that a signal has been sampled at $N$ equally spaced instants, giving the finite sequence $x_0, x_1, x_2, \ldots, x_{N-1}$ . Here, $N$ is the number of samples and $n$ is the sample index.

The Discrete Fourier Transform (DFT) determines the harmonic content of this finite set of samples. The integer $n$ is the sample number, so $x_n$ denotes the value of the signal at sample number $n$ , where $n = 0, 1, 2, \ldots, N-1$ . For real-valued sampled data, the reconstruction can be interpreted as a sum of discrete sine and cosine components, i.e.,

(190) $\begin{equation*} x_n = A_0 + \sum_{m=1}^{M} \left[A_m \cos\left(\frac{2\pi m n}{N}\right) + B_m \sin\left(\frac{2\pi m n}{N}\right)\right] \end{equation*}$

where $N$ is the total number of samples, $m$ is the harmonic number, and $M$ is the highest harmonic number retained in the reconstruction. For a finite sampled sequence, only a finite number of independent harmonic modes are available; for real-valued data, this limit is set by the Nyquist harmonic. The coefficients $A_m$ and $B_m$ determine how much of each cosine and sine component is present in the sampled data. With a small value of $M$ , only the slowly varying components are included, so the reconstruction gives a smoothed approximation to the sampled data. Increasing $M$ adds higher-frequency content and progressively improves the reconstruction.

The interpretation of the DFT is illustrated in Figure 16. The symbols show the original sampled data. The curves show reconstructions obtained by increasing the retained harmonic order $M$ . With only the lowest harmonics retained, the reconstruction captures only the broad variation in the data. As additional harmonics are included, the reconstruction follows progressively finer features of the sampled sequence. When all independent DFT harmonics are retained, the reconstruction passes through the original sampled values, apart from numerical roundoff error.

This discrete Fourier representation is the basis of most modern signal-processing methods. In computational work, the corresponding harmonic coefficients are usually found efficiently using the Fast Fourier Transform (FFT) algorithm. This allows sampled data from acoustics, vibration, communications, and aerodynamics to be analyzed by frequency content.

Check Your Understanding #20 – Interpreting Fourier reconstruction of sampled data

Figure 16 (above) shows a finite set of sampled data and several reconstructions obtained by retaining increasing numbers of harmonic components.

Why is the original signal shown as discrete sample points with symbols rather than as a continuous curve?
What does it mean to retain only the lowest harmonics in the reconstruction?
Why does the reconstruction with only a few harmonics appear smoother than the sampled data?
What changes as more harmonics are included in the reconstruction?
Why does retaining enough harmonics allow the reconstruction to pass through the sampled values?

Show solution/hide solution.

The original signal is shown as discrete sample points because measured data are usually recorded at finite time intervals. The available information is the finite sequence of sampled values, given by discrete symbols, not the complete continuous signal between those samples.
Retaining only the lowest harmonics means that only the slowly varying sine and cosine components are included. These components capture the broad variation in the sampled data but not the finer details of the signal.
The reconstruction with only a few harmonics appears smoother because the higher-frequency components have been omitted. Without these higher harmonics, rapid changes and small-scale variations in the sampled sequence cannot be represented.
As more harmonics are included, the reconstruction more closely matches the sampled data. The added harmonics allow the reconstruction to represent progressively finer features of the original sequence.
Retaining all independent DFT harmonics allows the reconstruction to pass through the sampled values because these harmonic components provide a complete representation of the finite sampled sequence, apart from numerical roundoff error.

Partial Differential Equations (PDEs)

Partial differential equations (PDEs) are equations that involve multiple independent variables and their corresponding partial derivatives. They model various engineering phenomena, including heat conduction, fluid dynamics, structural dynamics, elasticity, aeroelasticity, and acoustics.

PDEs can be classified into several types based on their order and linearity. The highest order of the partial derivatives involved determines the order of a PDE. For example, a second-order PDE involves second-order partial derivatives. Linearity refers to whether the equation is linear or nonlinear in terms of the unknown function and its derivatives. Examples of PDEs include:

Heat Equation

$\frac{\partial u}{\partial t} = \alpha \left( \frac{ \partial^2 u}{ \partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{ \partial^2 u }{ \partial z^2} \right)$

This equation is second-order, linear in ${u}$ , and parabolic. It has a single-time derivative balanced by a spatial Laplacian, so disturbances diffuse without oscillations; in the classical heat equation, the mathematical propagation speed is infinite, although the response decays rapidly with distance.

Wave Equation

$\frac{ \partial^2 u }{ \partial t^2} = a^2 \left (\frac{ \partial^2 u}{ \partial x^2} + \frac{ \partial^2 u} {\partial y^2} + \frac{ \partial^2 u }{\partial z^2} \right)$

This equation is second-order, linear, and hyperbolic. The second time derivative pairs with a spatial Laplacian, giving rise to propagating waves that carry information at the finite speed $a$ .

Navier-Stokes Equations (incompressible, convective form)

$\frac{ \partial \vec{V} }{\partial t} + (\vec{V} \bigcdot \nabla) \vec{V} = -\frac{1}{\varrho} \nabla p + \nu \nabla^2 \,\vec{V}$

It is second-order and nonlinear because of the convective term $(\vec{V} \bigcdot \nabla)\vec{V}$ . With viscosity present ( $\nu> 0$ ), the momentum equation contains a parabolic diffusion term, while the incompressibility constraint $\nabla \bigcdot \vec{V} = 0$ couples the velocity field to the pressure. In many incompressible-flow formulations, this coupling leads to an elliptic pressure equation. For this reason, the system is commonly treated as a coupled, nonlinear, mixed parabolic-elliptic set of PDEs.

Solving PDEs can be challenging, and different techniques are employed depending on the type of equation and its properties. For example, in incompressible flow theory, PDEs are used to determine the velocity potential and the stream function. Analytical methods involve finding exact solutions using techniques such as the separation of variables, a Fourier series, or Laplace transforms.

In general, exact or analytical solutions are feasible when the PDE is linear, and the geometry or boundary data are simple. Typical techniques include separation of variables, Fourier series expansions, and Laplace transforms. When nonlinearities or complex boundaries preclude closed-form solutions, numerical approximations based on finite-difference, finite-element, or finite-volume discretizations are employed.

Examples for Potential Flows

For two-dimensional incompressible, irrotational flow, the velocity potential $\phi$ and the stream function $\psi$ both satisfy Laplace’s equation, i.e.,

(191) $\begin{equation*} \nabla^{2}\phi = 0 \quad\text{and}\quad \nabla^{2}\psi = 0 \end{equation*}$

For example, in a rectangular domain, a separated solution of Laplace’s equation may be written in the form

(192) $\begin{equation*} \phi(x,y)=X(x)Y(y) \end{equation*}$

where

(193) $\begin{equation*} \frac{X''}{X}=-\frac{Y''}{Y}=-\lambda \end{equation*}$

so that

(194) $\begin{equation*} X''+\lambda X=0 \quad \text{and} \quad Y''-\lambda Y=0 \end{equation*}$

For boundary conditions with $X(0)=X(L_x)=0$ , the admissible eigenvalues are

(195) $\begin{equation*} \lambda_m=\left(\frac{m\pi}{L_x}\right)^2 \end{equation*}$

and the corresponding separated terms have the form

(196) $\begin{equation*} \sin\left(\frac{m\pi x}{L_x}\right) \left[ A_m\sinh\left(\frac{m\pi y}{L_x}\right) + B_m\cosh\left(\frac{m\pi y}{L_x}\right) \right] \end{equation*}$

with the constants determined from the remaining boundary conditions.

Numerical Discretisation

When analytical techniques fail, discretization is applied. On a uniform Cartesian grid, the two-dimensional Laplacian is approximated by the five-point stencil

(197) $\begin{equation*} \nabla^{2}\phi\approx \frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^{2}} + \frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta y^{2}} \end{equation*}$

This second-order scheme produces a linear system of equations, i.e.,

(198) $\begin{equation*} A\,\vec{\phi}= \vec{b} \end{equation*}$

which can be solved efficiently using the Gauss-Seidel method, successive overrelaxation, or multigrid iteration.

Verification relies on grid-convergence studies and manufactured solutions, whereas validation compares results with experimental or benchmark data, such as potential flow over a cylinder or the lid-driven cavity. Combining analytical insight with robust numerical schemes enables engineers to address the wide variety of incompressible-flow PDEs encountered in practice.

Check Your Understanding #21 – Partial differential equations

In a specific incompressible potential flow, the velocity field $(u,v)$ in the $(x,y)$ -plane is given by
${ u = \frac{\partial \phi}{\partial x} = 2x, \qquad v = \frac{\partial \phi}{\partial y} = -2y }$

Find the velocity potential $\phi(x,y)$ .
Solve Laplace’s equation, i.e.,
$\frac{\partial^{2}u}{\partial x^{2}} +\frac{\partial^{2}u}{\partial y^{2}} = 0, \qquad 0 < x< L, 0<y<H$

subject to the boundary conditions

$u(0,y) = 0, u(L,y) = 0, u(x,0) = 0, u(x,H) = f(x)$

Determine the harmonic function $u(x,y)$ .
Solve the one-dimensional heat equation
$u_t = \alpha u_{xx},\qquad 0 < x < L, \quad t > 0$

with the boundary conditions $u(0,t) = u(L,t) = 0$ and initial condition $u(x,0)=g(x)$ . Find the temperature distribution $u(x,t)$ .

Show solutions / hide solutions.

Integrate the relation $u = \partial\phi/\partial x = 2x$ with respect to $x$ , giving
$\phi = x^{2} + f_1(y)$

Differentiating this result with respect to $y$ and equating to the given velocity component $v = -2y$ gives $f_1'(y) = -2y$ . Hence $f_1(y) = -y^{2} + c$ , and so

$\phi(x,y) = x^{2} - y^{2} + c$
Assume a separable form $u = X(x) Y(y)$ . Separation of variables gives
$\dfrac{X''}{X} = -\dfrac{Y''}{Y} = -\lambda$

Enforcing $u(0,y) = u(L,y) = 0$ gives the eigen-pairs

${ \lambda_n = \left(\frac{n\pi}{L}\right)^2, \qquad X_n(x) = \sin \left(\frac{n\pi x}{L}\right) }$

The corresponding $Y_n$ solutions are proportional to $Y_n(y) = \sinh \left(n\pi y/L\right)$ . Express the top-wall data $f(x)$ as a sine series, i.e.,

$f(x) = \sum_{n=1}^{\infty} B_n \sin \left(\frac{n\pi x}{L}\right) \qquad B_n = \frac{2}{L}\int_{0}^{L} f(\xi)\sin \left(\frac{n\pi \xi}{L}\right)\,d\xi$

Therefore, the harmonic solution is

$u(x,y) = \sum_{n=1}^{\infty} B_n \frac{\sinh \bigl(n\pi y/L\bigr)} {\sinh \bigl(n\pi H/L\bigr)} \sin \left(\frac{n\pi x}{L}\right)$
Let $u = X(x) \, T(t)$ . After separation, then
$\dfrac{T'}{\alpha\,T} = \dfrac{X''}{X} = -\lambda$

The boundary conditions $u(0,t) = u(L,t) = 0$ give

$\lambda_n = \left(\frac{n\,\pi}{L}\right)^2,\qquad X_n(x) = \sin \left(\frac{n\,\pi\,x}{L}\right)$

The temporal factors satisfy $T_n' = -\alpha\,\lambda_n\,T_n$ , giving

$T_n(t) = \exp \left[ -\alpha \left(\frac{n\,\pi}{L}\right)^2 t \right]$

Expand the initial profile $g(x)$ as

$g(x) = \sum_{n=1}^{\infty} A_n \sin \left(\frac{n\,\pi\,x}{L}\right), \qquad A_n = \frac{2}{L}\int_{0}^{L} g(\xi)\sin \left(\frac{n\,\pi\,\xi}{L}\right)\,d\xi$

Therefore, the temperature distribution can be expressed as

$u(x,t) = \sum_{n=1}^{\infty} A_n \exp \left[ -\alpha \left(\frac{n\,\pi}{L}\right)^2 t \right] \sin \left(\frac{n\,\pi\,x}{L}\right)$

Other Coordinate Systems

In engineering, several coordinate systems are used to describe and analyze spatial relationships and physical phenomena. Cartesian coordinates, using orthogonal axes ${x}$ , ${y}$ , and ${z}$ , are natural coordinates and are fundamental for geometric modeling and structural analysis in three-dimensional space. Polar coordinates, defined by radial distance $r$ and azimuthal angle $\theta$ , are valid for applications involving radial symmetry. Cylindrical coordinates expand upon polar coordinates with an added vertical dimension ${z}$ , making them ideal for radial problems with cylindrical and axial symmetry. Spherical coordinates arise in specific contexts and are defined by the radial distance $r$ , the polar angle $\theta$ , and the azimuthal angle $\phi$ . Lastly, curvilinear coordinates, which encompass systems such as elliptic and parabolic coordinates, also find applications in engineering.

Polar Coordinates

Polar coordinates often arise in engineering problem-solving. A polar coordinate system is a two-dimensional coordinate system in which each point in a plane is defined by its distance from a reference point and the angle it forms with a reference direction. This axis system is advantageous for problems involving radial symmetry, such as fluid dynamics and structural analysis.

In the polar system, the radial coordinate $r$ is the distance from the origin O to another point P. The corresponding angular coordinate, $\theta$ , is shown in Figure 17; $\theta$ is measured counterclockwise from the positive ${x}$ -axis to the line connecting the origin to the point.

The polar coordinate system is well-suited to specific problems.

Polar coordinates are easily converted into Cartesian coordinates by using the transformation

(199) $\begin{equation*} x = r \cos\theta \quad \text{and} \quad y = r \sin\theta \end{equation*}$

and from Cartesian coordinates to polar coordinates using

(200) $\begin{equation*} r = \sqrt{x^2 + y^2} \quad \text{and} \quad \theta = \operatorname{atan2}(y,x) \end{equation*}$

Notice that the differential elements in polar coordinates are $ds = r \, d\theta$ for the arc length at fixed $r$ and $dA = r \, dr \, d\theta$ for the area element, both of which frequently occur in engineering problem-solving.

The principal value of the inverse tangent function typically returns an angle in the range $-\dfrac{\pi}{2}$ to $\dfrac{\pi}{2}$ (or -90 $^{\circ}$ to 90 $^{\circ}$ ). The tangent function is periodic, and the inverse function needs a unique range for every output value. The inverse tangent function alone does not account for the correct quadrant of the angle because $\tan\theta$ is positive for both the first and third quadrants and negative for both the second and fourth quadrants. To handle this, the atan2 function is used, which is a function available on all scientific calculators, i.e.,

(201) $\begin{equation*} \theta = \text{atan2}(y, x) \end{equation*}$

This function returns the angle in the correct quadrant by considering the signs of both ${x}$ and ${y}$ . It returns the angle and its correct sign in the range $-\pi$ to $\pi$ (or -180 $^{\circ}$ to 180 $^{\circ}$ ).

In polar coordinates, the commonly used vector operators are the gradient of a scalar field, i.e.,

(202) $\begin{equation*} \nabla p = \frac{\partial p}{\partial r} \, \hat{e}_r + \frac{1}{r} \frac{\partial p}{\partial \theta} \, \hat{e}_\theta \end{equation*}$

where $p$ is a scalar field, and $\hat{e}_r$ and $\hat{e}_\theta$ are the unit vectors in the $r$ and $\theta$ directions, respectively. The divergence of a vector ${\vec{V}}$ is

(203) $\begin{equation*} { \nabla \bigcdot \vec{V} = \frac{1}{r} \frac{\partial}{\partial r}(r V_r) + \frac{1}{r} \frac{\partial V_\theta}{\partial \theta} } \end{equation*}$

Finally, the Laplacian operator is

(204) $\begin{equation*} \nabla^2 p = \frac{1}{r} \frac{\partial}{\partial r} \left( r \, \frac{\partial p}{\partial r} \right) + \frac{1}{r^2} \, \frac{\partial^2 p}{\partial \theta^2} \end{equation*}$

Cylindrical & Spherical Coordinates

Cylindrical coordinates are a three-dimensional extension of polar coordinates. They are helpful for modeling problems involving cylindrical objects, such as pipes, shafts, and columns. In cylindrical coordinates $(r, \theta, z)$ , each component has a specific role in defining the position of a point in three-dimensional space, as shown in Figure 18. Spherical coordinates $(r, \theta, \phi)$ define the position of a point in space using three values: radial distance, polar angle, and azimuthal angle.

Cylindrical and spherical coordinate systems are convenient for problems involving radial and axial symmetry.

The radial coordinate $r$ represents the radial distance from the origin to the point projected onto the $(x, y)$ plane. It is similar to the $r$ in polar coordinates but extends vertically into three dimensions. The azimuthal angle, $\theta$ , is the counterclockwise angle from the positive ${x}$ -axis to the point projection onto the $(x, y)$ -plane. The vertical coordinate ${z}$ denotes the vertical height above (or below) the $(x, y)$ plane.

Converting from Cartesian to cylindrical coordinates is done by using the transformation

(205) $\begin{equation*} r = \sqrt{x^2 + y^2}, \quad \theta = \operatorname{atan2}(y,x), \quad \text{and} \quad z = z \end{equation*}$

and from cylindrical to Cartesian coordinates by using

(206) $\begin{equation*} x = r \cos\theta, \quad y = r \sin \theta, \quad \text{and} \quad z = z \end{equation*}$

In spherical coordinates, the radial distance is denoted here by $r$ . Some texts use $\rho$ for this coordinate, but $r$ avoids conflict with density notation. The polar angle, $\theta$ , is the angle between the positive ${z}$ -axis and the line from the origin to the point. The azimuthal angle, $\phi$ , is the angle between the positive ${x}$ -axis and the projection onto the $(x, y)$ -plane.

Converting from Cartesian to spherical coordinates is performed using

(207) $\begin{equation*} r = \sqrt{x^2 + y^2 + z^2}, \quad \theta = \cos^{-1}\left(\frac{z}{r}\right), \quad \phi = \operatorname{atan2}(y,x) \end{equation*}$

and spherical to Cartesian coordinates using

(208) $\begin{equation*} x = r \sin\theta \cos\phi, \quad y = r \sin\theta \sin\phi, \quad z = r \cos\theta \end{equation*}$

Transformations

In mathematics, transformations are operations that map a set of points or coordinates defining a body or some equation(s) to another set of points or equations, altering the position, size, shape, or orientation of the initial values. Transformation operations are fundamental to a wide range of engineering applications.

Translation

One type of transformation is a translation, as shown in Figure 19 for two-dimensional space, which shifts every point of a body by the same distance in a given direction. For example, translating a point $\text{P}(x, y) = \text{P}(1, 2)$ by 3 units to the right and 4 units up, the new coordinates P’ will be

(209) $\begin{equation*} \text{P'}(x', y') = (x + 3, y + 4) = (4, 6) \end{equation*}$

Translation moves bodies or coordinates without altering their orientation or size.

Rotation

Rotation turns a set of points defining a body about a fixed point through a specified angle. For instance, by rotating P 90 $^{\circ}$ counterclockwise (e.g., about the origin), the new coordinates P’ are

(210) $\begin{equation*} \text{P'}(x', y') = (-y, x) = (-2, 1) \end{equation*}$

Rotation is a transformation used in solving many types of engineering problems.

Reflection

Reflection reverses the set of points that define a body’s shape, producing its mirror image. For example, reflecting P over the ${y}$ -axis, the new coordinates P’ are

(211) $\begin{equation*} \text{P'}(x', y') = (-x, y) = (-1, 2) \end{equation*}$

This type of transformation is fundamental to understanding and solving problems involving planar symmetry.

Scaling

Scaling, or dilation, enlarges or reduces a body by a scale factor relative to a fixed point. If point P is scaled by a factor of 2 with respect to the origin, the new coordinates P’ are

(212) $\begin{equation*} \text{P'}(x', y') = (2x, 2y) = (2, 4) \end{equation*}$

Scaling is essential in similarity transformations, which occur frequently in engineering.

Shear

Shear distorts the body’s shape by sliding one part of it parallel to a fixed line. For example, in a shear transformation where the ${x}$ -coordinates are unchanged but the ${y}$ -coordinates are shifted by twice the ${x}$ -coordinate value, the new coordinates P’ are

(213) $\begin{equation*} \text{P'}(x', y') = (x, y + 2x) = (1, 4) \end{equation*}$

Shear transformations occur in engineering problems involving stress and deformation.

Affine Transformations

Affine transformations include translation, rotation, scaling, and shear. An affine transformation maps points to points and preserves collinearity, i.e., points that lie on a straight line remain on a straight line after the transformation, and ratios of distances along the same line are preserved; for example, the midpoint of a line segment remains the midpoint after the transformation. Projective transformations extend affine transformations and include perspective projections. An affine transformation can be represented as

(214) $\begin{equation*} T(\vec{x}) = A\, \vec{x} + \vec{b} \end{equation*}$

where $A$ is a matrix representing the linear part of the transformation, $\vec{b}$ is a vector representing the translation, and $\vec{x}$ is the input vector.

Transformation in the Lifting Line Theory

Lifting line theory is an aerodynamic method for analyzing finite wings. The lifting line theory uses the $\theta$ transformation to relate the spanwise position along the wing to the angle $\theta$ . The spanwise position ${y}$ along the wing is transformed into an angular coordinate $\theta$ using

$y = -\frac{b}{2} \, \cos \theta = -s \, \cos \theta$

where $b$ is the wing span, $s=b/2$ is the semi-span, ${y}$ is the spanwise coordinate, ranging from $-\dfrac{b}{2}$ to $\dfrac{b}{2}$ or ${-s}$ to ${s}$ , and $\theta$ is the transformed angular coordinate ranging from $0$ to $\pi$ . This transformation maps the spanwise coordinate to the angular coordinate, thereby facilitating the application of Fourier-series solution techniques.

Nonlinear Transformations

Nonlinear transformations can take various forms. Examples include polynomial transformations such as

(215) $\begin{equation*} T(x) = a_0 + a_1 \, x + a_2 \, x^2 + \cdots + a_n \, x^n \end{equation*}$

There are also exponential transformations, such as

(216) $\begin{equation*} T(x) = a \, e^{b \, x} = a \, \exp ({b \, x}) \end{equation*}$

as well as logarithmic transformations, i.e.,

(217) $\begin{equation*} T(x) = a \, \log(b \, x) \end{equation*}$

Error Function

The error function is a standard special function that first arose during the early nineteenth century and was first studied systematically by Carl Gauss within the context of the method of least squares, where integrals of the density distribution ${\exp(-x^{2})}$ or ${e^{-x^{2}}}$ appeared naturally. The error function is defined as

(218) $\begin{equation*} \mathrm{erf}(z) = \frac{2}{\sqrt{\pi}} \int_{0}^{z} e^{-x^{2}} \, dx = \frac{2}{\sqrt{\pi}} \int_{0}^{z} \exp({-x^{2}}) \, dx \end{equation*}$

which is shown in Figure 20.

Plots of the error function and complementary error function.

The error function has the following properties:

(219) $\begin{eqnarray*} \mathrm{erf}(-z) & = & -\,\mathrm{erf}(z) \\[6pt] \mathrm{erf}(0) & = & 0 \\[6pt] \lim_{z \to \infty} \mathrm{erf}(z) & = & 1 \\[6pt] \lim_{z \to -\infty} \mathrm{erf}(z) & = & -1 \end{eqnarray*}$

Its derivative is

(220) $\begin{equation*} \frac{d}{dz}\,\mathrm{erf}(z) = \frac{2}{\sqrt{\pi}}\,e^{-z^{2}} \end{equation*}$

The complementary error function is defined as

(221) $\begin{equation*} \mathrm{erfc}(z) = 1 - \mathrm{erf}(z) = \frac{2}{\sqrt{\pi}} \int_{z}^{\infty} e^{-x^{2}} \, dx \end{equation*}$

as also shown in the plot.

The error function arises naturally in problems of heat conduction and diffusion, where it describes how thermal or concentration disturbances spread in time. Joseph Fourier showed that solutions to the heat conduction equation could be expressed in terms of the same integral. However, the term “error function” became more commonly used because of its role in describing the distribution of random errors in observational data. Today, $\mathrm{erf}(x)$ and $\mathrm{erfc}(x)$ are recognized as standard special functions and crop up from time to time in engineering work.

Tensors

Many engineering quantities are adequately described using scalars or vectors. Mass is a scalar, while velocity and force are vectors. In these cases, the quantity has a single value or a single direction at a point. However, some quantities, such as stress in a solid or a fluid, depend on direction in more than one way, and for these quantities, a vector description is no longer sufficient. The force acting at a point depends not only on the direction of the force but also on the orientation of the surface on which that force acts. A single force vector cannot capture this dependence.

To represent such quantities, a mathematical object is needed that accepts a direction as input and returns a force vector as output. The simplest object with this capability is a ${3\times 3}$ matrix. In engineering, such objects are called tensors and can be viewed as matrices that encode directional dependence. For example, in Cartesian coordinates, the stress at a point is written as the matrix

(222) $\begin{equation*} \boldsymbol{\sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{bmatrix} \end{equation*}$

Each entry in the matrix represents a force per unit area. The column indicates the direction normal to a surface, and the row indicates the direction of the force acting on that surface. This matrix is called the stress tensor. The force per unit area acting on a surface with a unit normal vector

(223) $\begin{equation*} \vec{n} = \begin{bmatrix} n_x \\ n_y \\ n_z \end{bmatrix} \end{equation*}$

is obtained by direct matrix multiplication, i.e.,

(224) $\begin{equation*} \vec{t} = \boldsymbol{\sigma}\,\vec{n} \end{equation*}$

Written in component form, this gives

(225) $\begin{equation*} t_x = \sigma_{xx} n_x + \sigma_{xy} n_y + \sigma_{xz} n_z \end{equation*}$

(226) $\begin{equation*} t_y = \sigma_{yx} n_x + \sigma_{yy} n_y + \sigma_{yz} n_z \end{equation*}$

(227) $\begin{equation*} t_z = \sigma_{zx} n_x + \sigma_{zy} n_y + \sigma_{zz} n_z \end{equation*}$

In a fluid, the stress consists of pressure and viscous effects and may be written as

(228) $\begin{equation*} \boldsymbol{\sigma} = -p\,\boldsymbol{I} + \boldsymbol{\tau} \end{equation*}$

where $p$ is the pressure, $\boldsymbol{I}$ is the identity matrix, and $\boldsymbol{\tau}$ represents viscous stresses. Consider an example of a one-dimensional steady shear flow with velocity

(229) $\begin{equation*} \vec{v} = (u(y),\,0,\,0) \end{equation*}$

In this case, the stress matrix becomes

(230) $\begin{equation*} \boldsymbol{\sigma} = \begin{bmatrix} -p & \mu \dfrac{du}{dy} & 0 \\ \mu \dfrac{du}{dy} & -p & 0 \\ 0 & 0 & -p \end{bmatrix} \end{equation*}$

where $\mu$ is the dynamic viscosity. For a surface whose normal points in the $y$ -direction, substitution into the traction relation gives

(231) $\begin{equation*} \vec{t} = \begin{bmatrix} \mu \dfrac{du}{dy} \\ -p \\ 0 \end{bmatrix} \end{equation*}$

This result directly shows that pressure produces a force normal to the surface, velocity gradients produce shear forces, and changes in surface orientation alter the force, even though the fluid state at the point remains unchanged. Therefore, stress cannot be reduced to a single number or vector, and a $3\times 3$ matrix is the minimal mathematical structure required to consistently relate surface orientation to force.

Optimization

Engineering problems rarely have a single unique solution. Instead, engineers are typically confronted with a range of feasible alternatives that satisfy the same physical and operational requirements. Mathematics is used not only to predict system behavior, but also to determine which of these alternatives is preferable according to a specified criterion. This process is known as optimization. In practice, optimization arises whenever system performance depends continuously on one or more variables. Typical examples include selecting geometric dimensions, operating conditions, or control parameters to minimize losses, reduce weight, limit power consumption, or improve efficiency. By expressing a chosen performance metric as a mathematical function of the relevant variables, calculus provides systematic tools for identifying conditions that improve or optimize performance, replacing ad hoc trial-and-error approaches with a rational and defensible methodology.

Gradient-Based Optimization

To explain the process of optimization, let an objective (or cost) function be defined as

(232) $\begin{equation*} J = f(\vec{x}) \end{equation*}$

where $\vec{x} = (x_1, x_2, \ldots, x_n)$ is a vector of independent variables. These variables may represent geometric parameters, operating conditions, or other controllable or selectable quantities. The fundamental mathematical quantity underlying most optimization methods is the gradient of the function, i.e.,

(233) $\begin{equation*} \nabla f(\vec{x}) = \left( \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \ldots, \frac{\partial f}{\partial x_n} \right) \end{equation*}$

The gradient describes the local sensitivity of the function to changes in each variable. For a small change $\Delta \vec{x}$ , the corresponding change in the function can be approximated using a first-order Taylor expansion, i.e.,

(234) $\begin{equation*} f(\vec{x} + \Delta \vec{x}) \approx f(\vec{x}) + \nabla f \bigcdot \Delta \vec{x} \end{equation*}$

This expression indicates that the gradient points in the direction of the steepest local increase of the function. The direction of steepest local decrease is the negative gradient, $-\nabla f$ . Therefore, for an unconstrained local maximum or minimum that occurs at an interior point of the design space, a necessary condition is

(235) $\begin{equation*} \nabla f = \vec{0} \end{equation*}$

which generalizes the single-variable condition $df/dx = 0$ .

In many practical problems, solving $\nabla f = \vec{0}$ analytically is not convenient or possible. Instead, the optimization process is carried out iteratively. Starting from an initial estimate $\vec{x}^0$ , the variables are updated according to

(236) $\begin{equation*} \vec{x}^{k+1} = \vec{x}^k - \alpha^k \, \nabla f(\vec{x}^k) \end{equation*}$

where $\alpha^k$ is a positive step size. The step size must be chosen judiciously because a value that is too large may overshoot the minimum or prevent convergence. For a suitably chosen step size, the update moves the solution toward a decrease in $\scriptstyle f$ , following the negative local gradient, as shown in Figure 21. The step size may be fixed or adjusted during the iteration to improve convergence. In simple applications, engineers often assess convergence by examining plots of the objective function or by monitoring changes in the variables.

In a gradient-based optimization problem, the optimal solution is often reached fairly quickly with sufficient information.

Imposition of Constraints

In many engineering problems, the design variables cannot be chosen freely. Instead, they are constrained by physical, geometric, or operational requirements. Examples include limits on stress, power, weight, size, or operating conditions. Optimization problems that include such restrictions are called constrained optimization problems. Therefore, optimization under constraints is not about finding a perfect solution, but about identifying the best feasible solution.

An equality constraint restricts the variables to satisfy a relationship of the form

(237) $\begin{equation*} g(\vec{x}) = 0 \end{equation*}$

Such a constraint can also reduce the number of independent variables. In practice, this means that one variable can often be expressed in terms of the others and substituted directly into the objective function. For example, if the constraint allows one variable to be written as

(238) $\begin{equation*} x_n = h(x_1, x_2, \ldots, x_{n-1}) \end{equation*}$

then the objective function becomes a function of fewer variables, i.e.,

(239) $\begin{equation*} J = f(x_1, x_2, \ldots, x_{n-1}) \end{equation*}$

and the optimization can proceed using the same gradient-based ideas developed for unconstrained problems. From a geometric standpoint, the constraint defines a surface in the design space, and optimization is performed along this surface rather than over the entire space.

Many engineering constraints appear in inequality form, such as

(240) $\begin{equation*} g(\vec{x}) \le 0 \end{equation*}$

or as simple bounds on variables, i.e.,

(241) $\begin{equation*} x_{\min} \le x \le x_{\max} \end{equation*}$

In these cases, the optimum may occur either at an interior point, where the unconstrained optimum satisfies all constraints, or on the boundary of the feasible region. A common and effective engineering approach is first to determine the unconstrained optimum, then verify that it satisfies all constraints. If it does not, the objective function is evaluated along the relevant constraint boundaries, and the best feasible solution is selected.

Check Your Understanding #22 – Gradient-based optimization

Consider the objective function

$J(x,y) = x^2 + 2y^2 + 2x - 4y$

Compute the gradient $\nabla J$ .
Determine the stationary point of $J$ .
Starting from the initial estimate $(x_0,y_0) = (1,0)$ , perform one step of a gradient-based minimization using a step size $\alpha = 0.25$ .

Show solution / hide solution.

1. The gradient of $J$ is

$\nabla J = \left( \frac{\partial J}{\partial x}, \frac{\partial J}{\partial y} \right) = \left( 2x + 2,\; 4y - 4 \right)$

2. At a stationary point, the gradient vanishes, so

$2x + 2 = 0, \qquad 4y - 4 = 0$

Solving gives

$x = -1, \qquad y = 1$

Because the Hessian matrix is

$H = \begin{bmatrix} 2 & 0 \\ 0 & 4 \end{bmatrix}$

which is positive definite, this stationary point corresponds to a minimum.

3. A single gradient-based update is given by

$x_{1} = x_0 - \alpha (2x_0 + 2), \qquad y_{1} = y_0 - \alpha (4y_0 - 4)$

Substituting $(x_0,y_0) = (1,0)$ and $\alpha = 0.25$ gives

$x_1 = 1 - 0.25(4) = 0, \qquad y_1 = 0 - 0.25(-4) = 1$

This update moves the estimate from $(1,0)$ to $(0,1)$ , which is closer to the minimum at $(-1,1)$ .

Summary & Closure

The effective use of mathematics is fundamental to engineering problem-solving. Mathematics is not an abstract exercise but a practical means of modeling physical systems, quantifying behavior, and guiding design decisions. Geometry, algebra, trigonometry, calculus, vector analysis, complex numbers, and differential equations provide a common language for expressing relationships among physical quantities and for translating physical intuition into predictive models.

Many engineering problems involve scalar, vector, and field quantities, including forces, velocities, accelerations, fluxes, potentials, and flow fields. These quantities require systematic use of mathematical operations such as dot and cross products, gradients, divergence, curl, and line, surface, and volume integrals. Calculus further enables engineers to describe change, accumulation, rates, and balances, forming the basis for conservation laws and governing equations used throughout aerospace engineering.

The value of mathematics in engineering lies not only in formal manipulation but also in interpretation. Equations must be read as statements about physical systems, not merely as symbolic expressions. Engineers must understand the assumptions underlying a model, recognize the meaning of each term, and assess whether the mathematical result is consistent with the behavior of the system under analysis. This combination of mathematical skill and physical interpretation is what turns equations into useful engineering tools.

5-Question Self-Assessment Quickquiz

For Further Thought or Discussion

Explain the fundamentals of a scalar quantity. List some familiar scalar quantities used in engineering.
Explain the fundamentals of a vector quantity. Make a list of some examples of vector quantities used in engineering.
Explain the physical meaning of the curl of a vector field. Is the curl of a uniform flow positive, negative, or zero?
Explain why line integrals, surface integrals, and volume integrals are useful in engineering analysis.
What is the difference between an ordinary differential equation and a partial differential equation? Give an aerospace example of each.
If the function $A(x,y) = K(x^2 + y^2)$ where $K$ is a constant, then what is $A(r, \theta)$ ? Your answer may depend on whether you are biased toward pure mathematics, physics, or engineering.

Other Useful Online Resources

There is more to explore with these mathematical concepts:

This video offers a more in-depth examination of vector fields, divergence, and curl.
If you want to learn more about the Laplace equation, this website provides a wealth of information on the topic, as well as other mathematical operations.
All The Math You Need For Engineering: The Ultimate Guide (Step-by-Step).
Maths Every Engineer Must Know To Be Successful & Why.
You can’t suck at maths and be an engineer – an informative and motivational video.

"It’s all about mathematics!" Whether it’s fluid dynamics, aerodynamics, or orbital mechanics, the foundation often rests on understanding and applying mathematical principles. ↵
"Errors in engineering are often catastrophes." Even minor miscalculations can lead to severe consequences, which is why double-checking, simulations, and testing are vital to avoiding catastrophic outcomes. ↵
"Practice makes perfect." Consistent practice is essential for mastery. Whether designing a flight control system or finding optimal airfoil shapes, practice refines understanding and long-term skills. ↵
At face value, i.e., on-demand, as needed. ↵
"Work smart, not hard." Efficiency and cleverness are crucial in accomplishing tasks, rather than relying solely on sheer effort. ↵
The word eigen is German for “own” or “self.” “Eigenvalue” means an intrinsic or characteristic value of a system; “eigenvector” or “eigenfunction” means the corresponding intrinsic direction or shape. ↵

License

Icon for the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License

Introduction to Aerospace Flight Vehicles Copyright © 2022–2026 by J. Gordon Leishman is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, except where otherwise noted.

Digital Object Identifier (DOI)

https://doi.org/https://doi.org/10.15394/eaglepub.2022.1066.n8

Introduction

Algebra

Calculus

Ordinary Differentiation

Partial Differentiation

Integration

Working with Vectors

Scalar & Vector Fields

Scalar Products & Vector Projections

Vector Products

Differentiating a Vector

Line Integrals

Surface Integrals

Unit Normal Vector

Volume Integrals

Integral Relations

Gradient of a Scalar Field

Divergence of a Vector Field

Curl of a Vector Field

Other Mathematical Operators

Laplace Operator

Substantial Derivative Operator

Matrices & Systems of Linear Equations

Types of Matrices

Matrix Operations

Jacobian Matrix

Systems of Linear Equations

Row Reduction (Gaussian Elimination)

Matrix Inversion

LU Decomposition

Cramer’s Rule

Complex Numbers

Quaternions

Eigenvalues, Eigenvectors, & Eigenfunctions

Ordinary Differential Equations (ODEs)

State-Space Form

Laplace Transforms

Fourier Series

Partial Differential Equations (PDEs)

Examples for Potential Flows

Numerical Discretisation

Other Coordinate Systems

Polar Coordinates

Cylindrical & Spherical Coordinates

Transformations

Translation

Rotation

Reflection

Scaling

Shear

Affine Transformations

Nonlinear Transformations

Error Function

Tensors

Optimization

Gradient-Based Optimization

Imposition of Constraints

Summary & Closure

License

Digital Object Identifier (DOI)

Share This Book