Engineering Mathematics 4: Vector Calculus - From Green's Theorem to Stokes' Theorem

Engineering Mathematics 4: Vector Calculus

Vector calculus is the mathematical backbone of virtually every engineering discipline — from classical physics and electromagnetism to fluid dynamics and modern machine learning. This guide builds systematically from vector functions and fields through the three great integral theorems: Green's theorem, the Divergence theorem, and Stokes' theorem.

1. Vector Functions and Vector Fields

Vector-Valued Functions

A vector-valued function maps a scalar parameter $t$ to a vector in space.

$\mathbf{r}(t) = x(t)\,\mathbf{i} + y(t)\,\mathbf{j} + z(t)\,\mathbf{k}$

A classic example is the helix:

$\mathbf{r}(t) = \cos t\,\mathbf{i} + \sin t\,\mathbf{j} + t\,\mathbf{k}$

Tangent and Normal Vectors

The tangent vector of a curve $\mathbf{r}(t)$ is obtained by differentiation:

$\mathbf{r}'(t) = \frac{d\mathbf{r}}{dt} = \left(\frac{dx}{dt},\, \frac{dy}{dt},\, \frac{dz}{dt}\right)$

The unit tangent vector is $\mathbf{T}(t) = \dfrac{\mathbf{r}'(t)}{|\mathbf{r}'(t)|}$, and the principal unit normal is $\mathbf{N}(t) = \dfrac{\mathbf{T}'(t)}{|\mathbf{T}'(t)|}$.

Vector Fields

A vector field assigns a vector to every point in space:

$\mathbf{F}(x, y, z) = P(x,y,z)\,\mathbf{i} + Q(x,y,z)\,\mathbf{j} + R(x,y,z)\,\mathbf{k}$

Canonical examples:

• Gravitational field: $\mathbf{g} = -g\,\mathbf{k}$
• Electric field: $\mathbf{E} = k_e \dfrac{q}{r^2}\hat{\mathbf{r}}$
• Fluid velocity field: $\mathbf{v}(x,y,z)$

import numpy as np
import matplotlib.pyplot as plt

# Visualize a 2D rotational vector field
x = np.linspace(-2, 2, 15)
y = np.linspace(-2, 2, 15)
X, Y = np.meshgrid(x, y)

# F = (-y, x): pure rotation
U = -Y
V = X

fig, ax = plt.subplots(figsize=(6, 6))
ax.quiver(X, Y, U, V, color='steelblue', alpha=0.8)
ax.set_title('Rotational Vector Field F = (-y, x)')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.tight_layout()
plt.show()

2. Scalar Fields and the Gradient

Directional Derivative

For a scalar field $f(x,y,z)$, the directional derivative in the direction of unit vector $\hat{\mathbf{u}}$ is:

$D_{\hat{u}} f = \nabla f \cdot \hat{\mathbf{u}}$

The Gradient Vector

The gradient captures both the direction of steepest ascent and the rate of increase:

$\nabla f = \frac{\partial f}{\partial x}\,\mathbf{i} + \frac{\partial f}{\partial y}\,\mathbf{j} + \frac{\partial f}{\partial z}\,\mathbf{k}$

Geometric and physical meaning:

• Direction of $\nabla f$: the direction of maximum rate of increase of $f$
• Magnitude $|\nabla f|$: the maximum rate of increase
• $\nabla f$ is always perpendicular to the level surfaces $f = c$

Example: For $f(x,y) = x^2 + y^2$, the gradient is $\nabla f = (2x, 2y)$.

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# Symbolic gradient computation
x_s, y_s = sp.symbols('x y')
f = x_s**2 + y_s**2 + x_s * y_s

grad_x = sp.diff(f, x_s)
grad_y = sp.diff(f, y_s)
print("grad_x =", grad_x)  # 2x + y
print("grad_y =", grad_y)  # x + 2y

# Numerical visualization
x = np.linspace(-2, 2, 20)
y = np.linspace(-2, 2, 20)
X, Y = np.meshgrid(x, y)

f_val = X**2 + Y**2
U = 2 * X   # df/dx
V = 2 * Y   # df/dy

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
c = axes[0].contourf(X, Y, f_val, 20, cmap='viridis')
axes[0].contour(X, Y, f_val, 10, colors='white', alpha=0.5)
plt.colorbar(c, ax=axes[0])
axes[0].set_title('f(x,y) = x² + y² — Level Curves')

axes[1].quiver(X, Y, U, V, color='red', alpha=0.7)
axes[1].set_title('Gradient ∇f = (2x, 2y)')
plt.tight_layout()
plt.show()

3. Divergence

Definition

The divergence of a vector field $\mathbf{F} = (P, Q, R)$ is a scalar:

$\text{div}\,\mathbf{F} = \nabla \cdot \mathbf{F} = \frac{\partial P}{\partial x} + \frac{\partial Q}{\partial y} + \frac{\partial R}{\partial z}$

Physical Interpretation

Divergence measures the net outward flux per unit volume at a point:

• $\nabla \cdot \mathbf{F} > 0$: the point is a source — fluid flows outward
• $\nabla \cdot \mathbf{F} < 0$: the point is a sink — fluid flows inward
• $\nabla \cdot \mathbf{F} = 0$: incompressible flow — no net creation or destruction of fluid

Example: For $\mathbf{F} = (x^2, y^2, z^2)$:

$\nabla \cdot \mathbf{F} = 2x + 2y + 2z$

The rotational field $\mathbf{F} = (-y, x, 0)$ satisfies $\nabla \cdot \mathbf{F} = 0$ (incompressible).

import sympy as sp

x, y, z = sp.symbols('x y z')

P = x**2
Q = y**2
R = z**2

divergence = sp.diff(P, x) + sp.diff(Q, y) + sp.diff(R, z)
print("div F =", divergence)  # 2*x + 2*y + 2*z

4. Curl

Definition

The curl of $\mathbf{F} = (P, Q, R)$ is a vector, computed via the symbolic determinant:

$\text{curl}\,\mathbf{F} = \nabla \times \mathbf{F} = \begin{vmatrix} \mathbf{i} & \mathbf{j} & \mathbf{k} \\ \partial/\partial x & \partial/\partial y & \partial/\partial z \\ P & Q & R \end{vmatrix}$

Expanding:

$\nabla \times \mathbf{F} = \left(\frac{\partial R}{\partial y} - \frac{\partial Q}{\partial z}\right)\mathbf{i} - \left(\frac{\partial R}{\partial x} - \frac{\partial P}{\partial z}\right)\mathbf{j} + \left(\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}\right)\mathbf{k}$

Physical Meaning and Conservative Fields

• The curl measures local rotation or vorticity of a fluid element.
• $\nabla \times \mathbf{F} = \mathbf{0}$ implies the field is irrotational (conservative).
• A conservative field admits a potential function $\phi$ such that $\mathbf{F} = \nabla \phi$.

Curl in Maxwell's Equations

Faraday's law of induction:

$\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t}$

Ampere-Maxwell law:

$\nabla \times \mathbf{B} = \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \frac{\partial \mathbf{E}}{\partial t}$

import sympy as sp

x, y, z = sp.symbols('x y z')

# F = (y*z, x*z, x*y) — a conservative field
P = y * z
Q = x * z
R = x * y

curl_x = sp.diff(R, y) - sp.diff(Q, z)
curl_y = sp.diff(P, z) - sp.diff(R, x)
curl_z = sp.diff(Q, x) - sp.diff(P, y)

print("curl F =", (curl_x, curl_y, curl_z))
# Result: (0, 0, 0) => conservative (irrotational)

5. Line Integrals

Scalar Line Integral

For a curve $C$ parameterized by $\mathbf{r}(t)$, $a \le t \le b$:

$\int_C f\,ds = \int_a^b f(\mathbf{r}(t))\,|\mathbf{r}'(t)|\,dt$

Vector Line Integral (Work)

The work done by a force field $\mathbf{F}$ moving a particle along $C$:

$W = \int_C \mathbf{F} \cdot d\mathbf{r} = \int_a^b \mathbf{F}(\mathbf{r}(t)) \cdot \mathbf{r}'(t)\,dt$

Path Independence

For a conservative field $\mathbf{F} = \nabla \phi$, the line integral depends only on the endpoints:

$\int_C \mathbf{F} \cdot d\mathbf{r} = \phi(\mathbf{r}(b)) - \phi(\mathbf{r}(a))$

import numpy as np
from scipy.integrate import quad

# F = (2x, 2y), path: r(t) = (t, t^2), t in [0,1]
def work_integrand(t):
    x_t = t
    y_t = t**2
    Fx = 2 * x_t   # 2x
    Fy = 2 * y_t   # 2y
    dx = 1         # dx/dt
    dy = 2 * t     # dy/dt
    return Fx * dx + Fy * dy

work, _ = quad(work_integrand, 0, 1)
print(f"Work = {work:.4f}")  # = 2.0; phi = x^2 + y^2, phi(1,1)-phi(0,0) = 2

6. Surface Integrals

Scalar Surface Integral

For a parametric surface $\mathbf{r}(u,v)$:

$\iint_S f\,dS = \iint_D f(\mathbf{r}(u,v))\,|\mathbf{r}_u \times \mathbf{r}_v|\,dA$

Flux Integral

The flux of $\mathbf{F}$ through a surface $S$ measures the volume of fluid crossing $S$ per unit time:

$\Phi = \iint_S \mathbf{F} \cdot d\mathbf{S} = \iint_D \mathbf{F}(\mathbf{r}(u,v)) \cdot (\mathbf{r}_u \times \mathbf{r}_v)\,dA$

7. Green's Theorem

Statement

Let $C$ be a positively oriented (counterclockwise), piecewise-smooth, simple closed curve bounding a region $D$. If $P$ and $Q$ have continuous first partial derivatives on an open set containing $D$, then:

$\oint_C P\,dx + Q\,dy = \iint_D \left(\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}\right) dA$

Area Calculation

Green's theorem yields an elegant formula for the area enclosed by a curve:

$A = \frac{1}{2} \oint_C (x\,dy - y\,dx)$

This is the basis of the shoelace formula used in computational geometry to find polygon areas from vertex coordinates.

Fluid Mechanics Application

In 2D fluid flow, Green's theorem relates the circulation around a closed path to the net rotation (vorticity) within the enclosed region.

import numpy as np
from scipy.integrate import dblquad

# Green's theorem verification
# P = -y, Q = x => dQ/dx - dP/dy = 1 + 1 = 2
# Double integral over unit disk: 2 * pi * 1^2 = 2*pi

result, _ = dblquad(
    lambda y, x: 2.0,
    -1, 1,
    lambda x: -np.sqrt(max(0.0, 1 - x**2)),
    lambda x:  np.sqrt(max(0.0, 1 - x**2))
)
print(f"Double integral = {result:.5f}")  # ≈ 6.28318 = 2*pi
print(f"2*pi            = {2*np.pi:.5f}")

8. Gauss's Divergence Theorem

Statement

Let $V$ be a solid region with boundary surface $S = \partial V$ (outward orientation), and let $\mathbf{F}$ have continuous first partial derivatives throughout $V$. Then:

$\oiint_S \mathbf{F} \cdot \hat{\mathbf{n}}\,dS = \iiint_V (\nabla \cdot \mathbf{F})\,dV$

Physical Interpretation

The total outward flux through the closed surface equals the total amount of source/sink activity inside. This transforms a difficult surface integral into a potentially easier volume integral (or vice versa).

Gauss's Law for Electric Fields

$\oiint_S \mathbf{E} \cdot d\mathbf{A} = \frac{Q_{enc}}{\varepsilon_0}$

This is a direct application of the divergence theorem combined with $\nabla \cdot \mathbf{E} = \rho/\varepsilon_0$. It allows us to compute electric fields for highly symmetric charge distributions (sphere, cylinder, infinite plane).

import sympy as sp

x, y, z = sp.symbols('x y z')

# F = (x^3, y^3, z^3), integrate over unit sphere
P, Q, R = x**3, y**3, z**3
div_F = sp.diff(P, x) + sp.diff(Q, y) + sp.diff(R, z)
print("div F =", div_F)  # 3*x**2 + 3*y**2 + 3*z**2

# In spherical coords: integral of 3*r^2 over unit sphere = 3*(4*pi/5)
# Total flux = 12*pi/5
print("Expected total flux =", sp.Rational(12, 1) * sp.pi / 5)

9. Stokes' Theorem

Statement

Let $S$ be an oriented piecewise-smooth surface bounded by a simple, closed, piecewise-smooth curve $C = \partial S$ (oriented consistently with $S$). Then:

$\oint_C \mathbf{F} \cdot d\mathbf{r} = \iint_S (\nabla \times \mathbf{F}) \cdot \hat{\mathbf{n}}\,dS$

• Left side: line integral of $\mathbf{F}$ around $C$ (circulation)
• Right side: flux of $\text{curl}\,\mathbf{F}$ through $S$

Green's Theorem as a Special Case

When $S$ is a flat region $D$ in the $xy$-plane with $\hat{\mathbf{n}} = \mathbf{k}$, Stokes' theorem reduces to Green's theorem:

$\oint_C \mathbf{F} \cdot d\mathbf{r} = \iint_D (\nabla \times \mathbf{F}) \cdot \mathbf{k}\,dA = \iint_D \left(\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}\right) dA$

Ampere's Law

$\oint_C \mathbf{B} \cdot d\mathbf{l} = \mu_0 I_{enc}$

Stokes' theorem converts this integral form to the differential form $\nabla \times \mathbf{B} = \mu_0 \mathbf{J}$, one of Maxwell's equations.

import numpy as np

# Numerical verification of Stokes' theorem
# F = (-y, x, 0), S: upper hemisphere of unit sphere, C: unit circle
# curl F = (0, 0, 2)
# Surface integral: 2 * projected area of S onto xy-plane = 2 * pi
# Line integral: circulation around unit circle = 2*pi

R = 1.0
t = np.linspace(0, 2 * np.pi, 10000)
x_c = R * np.cos(t)
y_c = R * np.sin(t)
dx = np.gradient(x_c, t)
dy = np.gradient(y_c, t)
F_dot_dr = (-y_c) * dx + x_c * dy
line_integral = np.trapz(F_dot_dr, t)
print(f"Line integral  ≈ {line_integral:.5f}")
print(f"2*pi           = {2*np.pi:.5f}")

10. Engineering and AI Applications

Computational Fluid Dynamics (CFD)

The incompressible Navier-Stokes equations — the cornerstone of CFD — rely heavily on divergence and curl:

$\rho\left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \mathbf{f}$

Incompressibility constraint: $\nabla \cdot \mathbf{v} = 0$

Maxwell's Equations

Computer Graphics

In 3D rendering, surface normal vectors for lighting calculations (Phong shading, bump mapping) are computed exactly as $\mathbf{r}_u \times \mathbf{r}_v$ from surface parameterization — a direct application of vector calculus.

Gradient Descent in Machine Learning

The gradient descent optimizer updates parameters in the direction opposite to the gradient of the loss:

$\theta_{n+1} = \theta_n - \eta\,\nabla_\theta \mathcal{L}(\theta_n)$

This is gradient vector calculus applied directly to a high-dimensional parameter space. Variants such as Adam and RMSProp adaptively scale the gradient using its history.

Physics-Informed Neural Networks (PINNs)

PINNs embed physical laws (PDEs involving divergence, curl, gradient) directly into the loss function, using automatic differentiation to evaluate differential operators:

$\mathcal{L} = \mathcal{L}_{data} + \lambda\,\mathcal{L}_{physics}$

where $\mathcal{L}_{physics}$ might be $\|\nabla \cdot \mathbf{v}\|^2$ enforcing incompressibility.

11. Comprehensive Python Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad, quad
import sympy as sp

print("=== Vector Calculus Python Examples ===\n")

# --- 1. Gradient ---
x_s, y_s, z_s = sp.symbols('x y z')
f = x_s**3 + y_s**2 * z_s - x_s * y_s

grad = [sp.diff(f, v) for v in (x_s, y_s, z_s)]
print("1. Gradient")
print(f"   f = {f}")
print(f"   grad f = {grad}\n")

# --- 2. Divergence ---
P = x_s**2 * y_s
Q = y_s**2 * z_s
R = z_s**2 * x_s
div_F = sp.diff(P, x_s) + sp.diff(Q, y_s) + sp.diff(R, z_s)
print("2. Divergence")
print(f"   F = ({P}, {Q}, {R})")
print(f"   div F = {div_F}\n")

# --- 3. Curl ---
curl_x = sp.diff(R, y_s) - sp.diff(Q, z_s)
curl_y = sp.diff(P, z_s) - sp.diff(R, x_s)
curl_z = sp.diff(Q, x_s) - sp.diff(P, y_s)
print("3. Curl")
print(f"   curl F = ({sp.simplify(curl_x)}, {sp.simplify(curl_y)}, {sp.simplify(curl_z)})\n")

# --- 4. Green's theorem (numerical) ---
result, _ = dblquad(
    lambda y, x: 2.0,
    -1, 1,
    lambda x: -np.sqrt(max(0, 1 - x**2)),
    lambda x:  np.sqrt(max(0, 1 - x**2))
)
print("4. Green's Theorem verification")
print(f"   Double integral ≈ {result:.5f}, 2*pi ≈ {2*np.pi:.5f}\n")

# --- 5. Vector field visualization ---
x_arr = np.linspace(-2, 2, 20)
y_arr = np.linspace(-2, 2, 20)
X, Y = np.meshgrid(x_arr, y_arr)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

axes[0].quiver(X, Y, 2*X, 2*Y, color='crimson', alpha=0.7)
axes[0].set_title('Gradient of x²+y²')

axes[1].quiver(X, Y, -Y, X, color='steelblue', alpha=0.7)
axes[1].set_title('Rotational Field (-y, x)')

axes[2].quiver(X, Y, X, Y, color='forestgreen', alpha=0.7)
axes[2].set_title('Source Field (x, y)')

plt.tight_layout()
plt.savefig('vector_fields_en.png', dpi=120)
plt.show()

12. Quiz

References

• Kreyszig, E. — Advanced Engineering Mathematics (10th ed.), Wiley. The standard reference for engineering mathematics.
• MIT OCW 18.02 — Multivariable Calculus (Denis Auroux). ocw.mit.edu
• Stewart, J. — Calculus: Early Transcendentals (9th ed.), Cengage. Chapters on line integrals, surface integrals, and vector theorems.
• Griffiths, D.J. — Introduction to Electrodynamics (4th ed.), Cambridge University Press. Physical applications of all three vector theorems.
• SciPy Documentation — scipy.integrate for numerical multi-dimensional integration. scipy.org
• SymPy Documentation — Symbolic differentiation and integration library. sympy.org
• Raissi, M. et al. — "Physics-informed neural networks," Journal of Computational Physics, 2019. PINNs and vector calculus in deep learning.