∇ Multivariable Calculus

§01 Vectors & Vector Spaces

1.1 — Vectors in ℝⁿ

A vector v ∈ ℝⁿ is an ordered n-tuple v = (v₁, v₂, …, vₙ). In ℝ³: v = v₁i + v₂j + v₃k. Vectors have magnitude and direction — they are the atoms of multivariable calculus.

1.2 — Dot Product

1.3 — Cross Product (ℝ³)

1.4 — Triple Products

1.5 — Parametric Curves & Arc Length

A curve r(t) = (x(t), y(t), z(t)) for t ∈ [a,b]. The arc length is:

1.6 — Curvature

The Frenet-Serret frame {T, N, B} with T = r′/‖r′‖, N = T′/‖T′‖, B = T × N satisfies:

where τ is the torsion measuring how the curve twists out of its osculating plane.

§02 Partial Derivatives

2.1 — Definition

For f : ℝⁿ → ℝ, the partial derivative with respect to xᵢ is:

Hold all other variables constant and differentiate with respect to one — this is the fundamental operation of multivariable calculus.

2.2 — Clairaut's Theorem (Symmetry of Mixed Partials)

2.3 — Higher-Order Partials

For f : ℝ² → ℝ, the second-order partials form the Hessian matrix:

2.4 — Linearization & Differentiability

f is differentiable at a if there exists a linear map L : ℝⁿ → ℝ such that:

The matrix of L is the Jacobian [∂fᵢ/∂xⱼ]. For f : ℝⁿ → ℝ, this reduces to the gradient row vector.

2.5 — The Chain Rule (Multivariable)

§03 Gradient & Directional Derivatives

3.1 — The Gradient Vector

For f : ℝⁿ → ℝ, the gradient is the vector of partial derivatives:

The gradient points in the direction of steepest ascent and is perpendicular to level sets.

3.2 — Directional Derivative

where u is a unit vector. The rate of change of f in direction u.

3.3 — Gradient is Normal to Level Surfaces

3.4 — Interactive: Gradient Field Explorer

3.5 — Tangent Planes

For a surface F(x,y,z) = 0, the tangent plane at point (a,b,c) is:

For z = f(x,y), the tangent plane is z − f(a,b) = f_x(a,b)(x−a) + f_y(a,b)(y−b).

§04 Optimization & Lagrange Multipliers

4.1 — Critical Points

A critical point of f : ℝⁿ → ℝ is where ∇f = 0 (or ∇f undefined). These are candidates for local extrema.

4.2 — Second Derivative Test (ℝ²)

4.3 — Lagrange Multipliers

4.4 — Multiple Constraints

With constraints g₁(x) = 0, …, gₖ(x) = 0:

This gives n + k equations in n + k unknowns (x₁,…,xₙ,λ₁,…,λₖ).

4.5 — Interactive: Gradient Descent on a Surface

§05 Multiple Integration

5.1 — Double Integrals

Computes signed volume under the surface z = f(x,y) over region R.

5.2 — Fubini's Theorem

5.3 — Double Integrals in Polar Coordinates

5.4 — Triple Integrals

Applications: mass = ∭ ρ dV, center of mass = (1/M)∭ x·ρ dV, moment of inertia = ∭ r²·ρ dV.

5.5 — Cylindrical & Spherical Coordinates

5.6 — The Gaussian Integral via Polar Coordinates

§06 Change of Variables & Jacobians

6.1 — The Jacobian Matrix

For a transformation T : ℝⁿ → ℝⁿ, the Jacobian matrix is:

6.2 — Change of Variables Theorem

6.3 — Interactive: Jacobian Transformation Visualizer

6.4 — Inverse Function Theorem

6.5 — Implicit Function Theorem

Cryptographic application: The implicit function theorem guarantees that elliptic curves y² = x³ + ax + b have well-defined tangent lines (and thus a group law) wherever 2y ≠ 0 — this is the foundation of ECC point doubling.

§07 Vector Fields

7.1 — Definition

A vector field F : ℝⁿ → ℝⁿ assigns a vector to each point. In ℝ³: F(x,y,z) = Pi + Qj + Rk.

7.2 — Divergence

Measures the rate of "expansion" of the field at a point. Positive = source, negative = sink.

7.3 — Curl

Measures the "rotation" of the field. If curl F = 0, the field is irrotational.

7.4 — Fundamental Identities

7.5 — The Laplacian

Laplace's equation ∇²f = 0 governs steady-state heat, electrostatics, and gravitational potentials. Solutions are called harmonic functions.

7.6 — Conservative Fields & Potential Functions

7.7 — Helmholtz Decomposition

Any sufficiently smooth vector field F can be decomposed:

into a curl-free (conservative) part and a divergence-free (solenoidal) part. This is fundamental in electrodynamics (E and B fields) and fluid mechanics.

§08 Line Integrals

8.1 — Scalar Line Integral

Integrates a scalar function along a curve weighted by arc length.

8.2 — Vector Line Integral (Work)

Computes the work done by force F along curve C.

8.3 — Fundamental Theorem for Line Integrals

8.4 — Circulation & Flux in 2D

8.5 — Complex Line Integrals

For f(z) analytic, ∮ f(z)dz = ∮ (u+iv)(dx+idy) connects complex analysis (residues) to real line integrals in ℝ². The Cauchy integral formula is a line integral that recovers all derivatives of f — see the cryptographic applications in §14.

§09 Surface Integrals

9.1 — Parametric Surfaces

A surface S is parameterized by r(u,v) = (x(u,v), y(u,v), z(u,v)) for (u,v) ∈ D.

The normal vector is n = r_u × r_v.

9.2 — Surface Area

For z = f(x,y): dS = √(1 + f_x² + f_y²) dx dy.

9.3 — Scalar Surface Integral

9.4 — Flux Through a Surface

Measures the flow of F through S. Orientation matters — the sign of r_u × r_v determines which side is "positive."

9.5 — Flux of Gravitational/Electric Fields

For the inverse-square field F = q r̂/‖r‖², the flux through a sphere of radius R centered at the origin:

This is Gauss's Law — the flux depends only on the enclosed charge, not the surface shape (proved generally via the Divergence Theorem §12).

9.6 — First & Second Fundamental Forms

The first fundamental form I = E du² + 2F du dv + G dv² where E = r_u·r_u, F = r_u·r_v, G = r_v·r_v encodes the intrinsic metric of the surface. The Gaussian curvature K = (LN−M²)/(EG−F²) from the second fundamental form is invariant under isometry (Gauss's Theorema Egregium).

§10 Green's Theorem

10.1 — Statement

Relates a line integral around a simple closed curve C (oriented counterclockwise) to a double integral over the region D it encloses.

10.2 — Area via Green's Theorem

This is exactly how GPS systems and GIS software compute polygon areas — the Shoelace formula is the discrete version.

10.3 — Green's Identities

First Green's Identity:

Second Green's Identity:

10.4 — Green's Functions

The Green's function G(x;x₀) satisfies ∇²G = δ(x−x₀). In 2D: G = (1/2π)ln‖x−x₀‖. Any Poisson equation ∇²u = f can be solved via u(x) = ∬ G(x;x₀)f(x₀) dA₀. This is how potential fields in cryptographic lattice problems are analyzed.

§11 Stokes' Theorem

11.1 — Statement

The line integral of F around the boundary ∂S of an oriented surface S equals the flux of curl F through S. Green's theorem is the special case where S is flat in the xy-plane.

11.2 — Consequences

If curl F = 0 on a simply connected domain, then ∮ F·dr = 0 for every closed curve, so F is conservative (F = ∇φ). This is the 3D extension of the path-independence theorem.

11.3 — Physical Interpretation

In electrodynamics, Faraday's law is Stokes' theorem applied to the electric field:

Ampère's law (with Maxwell's correction) is Stokes' for the magnetic field:

§12 The Divergence Theorem (Gauss's Theorem)

12.1 — Statement

The outward flux of F through the closed surface ∂V equals the total divergence of F inside V.

12.2 — Gauss's Law (Physics)

The total electric flux through any closed surface equals the enclosed charge divided by ε₀. This is the Divergence Theorem applied to Maxwell's first equation ∇·E = ρ/ε₀.

12.3 — The Grand Unified View

All are instances of the generalized Stokes' theorem: ∫_∂Ω ω = ∫_Ω dω (see §13).

§13 Differential Forms & the Generalized Stokes' Theorem

13.1 — 0-Forms, 1-Forms, 2-Forms

Differential forms unify all integration theorems into one framework.

13.2 — The Wedge Product

The wedge product ∧ is antisymmetric: dx∧dy = −dy∧dx, and dx∧dx = 0.

13.3 — The Exterior Derivative

The exterior derivative d maps k-forms to (k+1)-forms:

13.4 — Recovering the Classical Operations

13.5 — The Generalized Stokes' Theorem

This single equation encapsulates: the Fundamental Theorem of Calculus (Ω = [a,b]), Green's theorem (Ω = region in ℝ²), classical Stokes' (Ω = surface in ℝ³), and the Divergence theorem (Ω = volume in ℝ³).

13.6 — De Rham Cohomology (Preview)

The quotient space H^k = {closed k-forms} / {exact k-forms} captures the topology of the domain. A nonzero class in H¹ means there exists an irrotational field that is not conservative — because the domain has a "hole." This is why ∇ × F = 0 implies F = ∇φ only on simply connected domains.

§14 Cryptography — Multivariable Calculus in Security

14.1 — Elliptic Curve Cryptography: The Group Law via Calculus

An elliptic curve E: y² = x³ + ax + b defines a group where point addition uses implicit differentiation (a multivariable calculus operation).

14.2 — Gradient Descent in Machine Learning & Adversarial Cryptanalysis

Neural networks are trained by minimizing a loss function L(θ) over parameter space θ ∈ ℝⁿ using gradient descent: θ_t+1 = θ_t − η∇L(θ_t). This is a direct application of §3 (gradient = steepest descent direction).

Adversarial attacks: To fool a neural network, compute ∇_xL and perturb the input: x_adv = x + ε · sign(∇_xL) (FGSM attack). This is gradient ascent on the loss — maximizing misclassification via calculus.

14.3 — Lattice Cryptography & Gaussian Measures in ℝⁿ

Post-quantum schemes (CRYSTALS-Kyber, Dilithium) sample errors from discrete Gaussians over lattices Λ ⊂ ℝⁿ:

The normalization ρ_σ(Λ) = Σ_v∈Λ exp(−π‖v‖²/σ²) is computed via the Poisson summation formula:

where f̂ is the Fourier transform (a multivariable integral) and Λ* is the dual lattice.

14.4 — Differential Privacy via Multivariable Sensitivity

For a vector-valued query f : D → ℝⁿ, the ℓ₂-sensitivity is:

The Gaussian mechanism adds noise: output = f(D) + N(0, σ²Iₙ) where σ ≥ Δ₂f · √(2 ln(1.25/δ)) / ε for (ε,δ)-differential privacy. The multivariable Gaussian (§5.6) is essential.

14.5 — Homomorphic Encryption & Polynomial Rings

Fully Homomorphic Encryption (FHE) operates on polynomial rings ℤ[x]/(xⁿ+1). The noise growth under multiplication is analyzed via the ℓ₂ norm on coefficient vectors — a multivariable problem. The Kannan embedding reduces decryption to finding shortest vectors in high-dimensional lattices, where gradient-based methods in ℝⁿ provide heuristic attacks (the BKZ algorithm uses Gram-Schmidt orthogonalization, which is multivariable linear algebra).

14.6 — Zero-Knowledge Proofs & Polynomial Commitments

In a ZK-SNARK, a computation is encoded as a Quadratic Arithmetic Program (QAP) — a system of polynomial equations in multiple variables. The prover evaluates these polynomials over an elliptic curve group, using the multivariable chain rule to handle composed polynomial evaluations efficiently. The KZG commitment scheme commits to a polynomial p(x) by evaluating it at a secret point s hidden in an elliptic curve scalar multiplication: C = [p(s)]G, which combines ECC (§14.1) with polynomial evaluation.

14.7 — Information Geometry & Side-Channel Analysis

The space of probability distributions forms a Riemannian manifold with the Fisher information metric:

This is a multivariable inner product on the tangent space of the statistical manifold. In side-channel cryptanalysis, the Fisher metric measures how efficiently power traces reveal a secret key — the Cramér-Rao bound (the inverse of the Fisher information matrix) gives the minimum variance of any unbiased estimator of the key.

% --- 14A: ECC Point Operations on y²=x³+7 (mod 97) (secp-like) ---
p = 97; a = 0; b = 7;
function r = mod_inv(x, p)
  % Extended Euclidean via Fermat's little theorem
  r = 1; base = mod(x, p); exp = p - 2;
  while exp > 0
    if mod(exp, 2) == 1, r = mod(r * base, p); end
    exp = floor(exp / 2); base = mod(base * base, p);
  end
end

P = [15, 13]; % verify: 13²=169≡72, 15³+7=3382≡72 mod 97 ✓
m = mod((3*P(1)^2 + a) * mod_inv(2*P(2), p), p); % tangent slope (implicit diff)
x3 = mod(m^2 - 2*P(1), p); y3 = mod(m*(P(1)-x3) - P(2), p);
printf("ECC: 2·(%d,%d) = (%d,%d) on y²=x³+7 mod %d\n", P(1),P(2),x3,y3,p);

% --- 14B: Gradient Descent (find minimum of Rosenbrock) ---
f  = @(v) (1-v(1))^2 + 100*(v(2)-v(1)^2)^2;
gf = @(v) [-2*(1-v(1))-400*v(1)*(v(2)-v(1)^2); 200*(v(2)-v(1)^2)];
x = [-1.5; 2]; % start point
for k = 1:50000
  g = gf(x); if norm(g) < 1e-8, break; end
  alpha = 0.001;
  while f(x-alpha*g) > f(x) - 1e-4*alpha*(g'*g), alpha = alpha*0.5; end
  x = x - alpha*g;
end
printf("Gradient descent: min at (%.6f,%.6f), f=%.2e, iters=%d\n", x(1),x(2),f(x),k);

% --- 14C: Fisher Information Matrix for Normal Distribution ---
% p(x|μ,σ) = N(μ,σ²), θ=(μ,σ). Fisher = [[1/σ², 0]; [0, 2/σ²]]
sigma = 2.0; N = 100000;
x = sigma*randn(N,1); % sample from N(0,σ²)
% Score functions: ∂log p/∂μ = x/σ², ∂log p/∂σ = -1/σ + x²/σ³
s_mu = x / sigma^2;
s_sigma = -1/sigma + x.^2/sigma^3;
S = [s_mu, s_sigma];
I_hat = (S' * S) / N;  % empirical Fisher
I_exact = [1/sigma^2, 0; 0, 2/sigma^2];
printf("Fisher info (σ=%.1f):\n  Empirical: [[%.4f,%.4f],[%.4f,%.4f]]\n  Exact:     [[%.4f,%.4f],[%.4f,%.4f]]\n", ...
  sigma, I_hat(1,1),I_hat(1,2),I_hat(2,1),I_hat(2,2), I_exact(1,1),I_exact(1,2),I_exact(2,1),I_exact(2,2));

§15 GNU Octave — Multivariable Calculus Laboratory

15.1 — Numerical Gradient & Hessian

function g = num_gradient(f, x, h)
  if nargin < 3, h = 1e-20; end
  n = length(x); g = zeros(n,1);
  for k = 1:n
    e = zeros(n,1); e(k) = h;
    g(k) = imag(f(x + 1i*e)) / h;
  end
end

function H = num_hessian(f, x)
  h = 1e-5; n = length(x); H = zeros(n);
  for i = 1:n
    for j = i:n
      ei = zeros(n,1); ei(i)=h; ej = zeros(n,1); ej(j)=h;
      H(i,j) = (f(x+ei+ej)-f(x+ei-ej)-f(x-ei+ej)+f(x-ei-ej))/(4*h^2);
      H(j,i) = H(i,j); % symmetry (Clairaut)
    end
  end
end

% Test: f(x,y) = x²y + sin(xy)
f = @(v) v(1)^2*v(2) + sin(v(1)*v(2));
x0 = [1; 2];
g = num_gradient(f, x0);
H = num_hessian(f, x0);
printf("∇f(1,2) = [%.6f, %.6f]\n", g(1), g(2));
printf("Hessian:\n  [%.4f  %.4f]\n  [%.4f  %.4f]\n", H(1,1),H(1,2),H(2,1),H(2,2));

15.2 — Jacobian & Change of Variables

function J = num_jacobian(F, x, h)
  if nargin < 3, h = 1e-7; end
  n = length(x); f0 = F(x); m = length(f0);
  J = zeros(m, n);
  for k = 1:n
    e = zeros(n,1); e(k) = h;
    J(:,k) = (F(x+e) - F(x-e)) / (2*h);
  end
end

% Polar → Cartesian: T(r,θ) = (r cos θ, r sin θ)
T = @(u) [u(1)*cos(u(2)); u(1)*sin(u(2))];
pt = [3; pi/4];
J = num_jacobian(T, pt);
printf("Jacobian at (r=3,θ=π/4):\n  [%.4f  %.4f]\n  [%.4f  %.4f]\n", J(1,1),J(1,2),J(2,1),J(2,2));
printf("det(J) = %.4f (exact: r = %.4f)\n", det(J), pt(1));

15.3 — Double & Triple Integration

% Double integral: ∬ exp(-(x²+y²)) dA over [-3,3]²
f2 = @(x,y) exp(-(x.^2+y.^2));
I2 = dblquad(f2, -3, 3, -3, 3);
printf("∬ exp(-(x²+y²)) dA = %.8f (exact π = %.8f)\n", I2, pi);

% Triple integral: volume of unit sphere via ∭ 1 dV
% In spherical: ∫₀²π ∫₀^π ∫₀¹ ρ² sin φ dρ dφ dθ
f3 = @(rho,phi,theta) rho.^2 .* sin(phi);
% Use nested quad
I_theta = @(theta) quad(@(phi) arrayfun(@(p) quad(@(rho) rho.^2.*sin(p), 0, 1), phi), 0, pi);
V = quad(I_theta, 0, 2*pi);
printf("Volume of unit sphere: %.6f (exact 4π/3 = %.6f)\n", V, 4*pi/3);

15.4 — Divergence, Curl & Laplacian

% Numerical div, curl, and Laplacian on a 3D grid
N = 30; h = 2*pi/(N-1);
[X,Y,Z] = meshgrid(linspace(0,2*pi,N));

% Vector field F = (sin(y), sin(z), sin(x))
Fx = sin(Y); Fy = sin(Z); Fz = sin(X);

% Divergence: ∂Fx/∂x + ∂Fy/∂y + ∂Fz/∂z
% For this field: 0 + 0 + 0 = 0 (solenoidal)
[dFxdx] = gradient(Fx, h); % gradient gives partials along each dim
[~,dFydy] = gradient(Fy, h);
[~,~,dFzdz] = gradient(Fz, h);
divF = dFxdx + dFydy + dFzdz;
printf("div F max (should be ≈0): %.6f\n", max(abs(divF(:))));

% Curl: (∂Fz/∂y-∂Fy/∂z, ∂Fx/∂z-∂Fz/∂x, ∂Fy/∂x-∂Fx/∂y)
[~,dFzdy] = gradient(Fz, h); [~,~,dFydz] = gradient(Fy, h);
[~,~,dFxdz] = gradient(Fx, h); [dFzdx] = gradient(Fz, h);
[dFydx] = gradient(Fy, h); [~,dFxdy] = gradient(Fx, h);
curlx = dFzdy - dFydz; curly = dFxdz - dFzdx; curlz = dFydx - dFxdy;
printf("curl F at center: (%.4f, %.4f, %.4f)\n", ...
  curlx(N/2,N/2,N/2), curly(N/2,N/2,N/2), curlz(N/2,N/2,N/2));

% Laplacian of f = sin(x)sin(y)sin(z): ∇²f = -3f
f = sin(X).*sin(Y).*sin(Z);
[fx,fy,fz] = gradient(f, h);
[fxx] = gradient(fx, h); [~,fyy] = gradient(fy, h); [~,~,fzz] = gradient(fz, h);
lap = fxx + fyy + fzz;
exact_lap = -3*f;
printf("Laplacian max error: %.4e\n", max(abs(lap(:)-exact_lap(:))));

15.5 — Line & Surface Integrals

% Verify Green's theorem: ∮(y²dx + x²dy) = ∬(2x-2y)dA over unit disk

% Line integral: parameterize C as (cos t, sin t)
t = linspace(0, 2*pi, 10001); dt = t(2)-t(1);
x = cos(t); y = sin(t);
dx = -sin(t)*dt; dy = cos(t)*dt;
line_int = sum(y.^2.*dx + x.^2.*dy);

% Double integral: ∬(2x-2y)dA in polar
area_int = dblquad(@(r,th) (2*r.*cos(th) - 2*r.*sin(th)).*r, 0, 1, 0, 2*pi);

printf("Green's theorem verification:\n");
printf("  Line integral:   %.8f\n", line_int);
printf("  Double integral:  %.8f\n", area_int);
printf("  Difference:       %.2e\n", abs(line_int - area_int));

15.6 — Laplace's Equation (2D PDE Solver)

% Solve ∇²u = 0 on [0,1]² with u=sin(πx) on top, u=0 on other sides
N = 50; h = 1/(N-1);
u = zeros(N,N);
x = linspace(0,1,N);
u(1,:) = sin(pi*x); % top boundary

for iter = 1:5000
  u_old = u;
  u(2:end-1,2:end-1) = 0.25*(u(1:end-2,2:end-1) + u(3:end,2:end-1) + ...
                                 u(2:end-1,1:end-2) + u(2:end-1,3:end));
  if max(abs(u(:)-u_old(:))) < 1e-8, break; end
end

% Exact: u(x,y) = sin(πx)sinh(π(1-y))/sinh(π)
[X,Y] = meshgrid(x, x);
u_exact = sin(pi*X) .* sinh(pi*(1-Y)) / sinh(pi);
printf("Laplace solver: %d iters, max error = %.3e\n", iter, max(abs(u(:)-u_exact(:))));

15.7 — Monte Carlo Integration in Higher Dimensions

% Volume of the unit ball in n dimensions by Monte Carlo
% V_n = π^(n/2) / Γ(n/2 + 1)
for n = [2, 3, 5, 10, 20]
  N_samples = 1e6;
  X = 2*rand(N_samples, n) - 1;  % uniform in [-1,1]^n
  inside = sum(X.^2, 2) <= 1;
  V_mc = 2^n * sum(inside) / N_samples;
  V_exact = pi^(n/2) / gamma(n/2 + 1);
  printf("n=%2d: V_mc=%.4f, V_exact=%.4f, ratio=%.4f\n", n, V_mc, V_exact, V_mc/V_exact);
end

15.8 — Stokes' Theorem Numerical Verification

% F = (y, -x, z²), curl F = (0, 0, -2)
% Surface: z = 0 disk of radius 1, boundary is unit circle in xy-plane

% Line integral: ∮ F·dr around unit circle (cos t, sin t, 0)
t = linspace(0, 2*pi, 100001); dt = t(2)-t(1);
x = cos(t); y = sin(t);
% F·r' = y·(-sin t) + (-x)·(cos t) + 0
integrand = y.*(-sin(t)) + (-x).*(cos(t));
line_val = sum(integrand)*dt;

% Surface integral: ∬ curl F · dS = ∬ (0,0,-2)·(0,0,1) dA = -2·Area(disk) = -2π
surf_val = -2*pi;

printf("Stokes' theorem:\n  ∮ F·dr = %.8f\n  ∬ curl F·dS = %.8f\n  Error: %.2e\n", ...
  line_val, surf_val, abs(line_val - surf_val));

15.9 — Constrained Optimization (Lagrange Multipliers)

% Maximize f(x,y,z)=xyz subject to x²+y²+z²=3
% ∇f = λ∇g → (yz,xz,xy) = λ(2x,2y,2z)
% Solve: yz=2λx, xz=2λy, xy=2λz, x²+y²+z²=3
% By symmetry: x=y=z=1, λ=1/2, f=1

% Numerical verification via penalty method
penalty = @(v) -(v(1)*v(2)*v(3)) + 1000*(v(1)^2+v(2)^2+v(3)^2-3)^2;
grad_pen = @(v) [
  -v(2)*v(3) + 4000*v(1)*(v(1)^2+v(2)^2+v(3)^2-3);
  -v(1)*v(3) + 4000*v(2)*(v(1)^2+v(2)^2+v(3)^2-3);
  -v(1)*v(2) + 4000*v(3)*(v(1)^2+v(2)^2+v(3)^2-3)];
x = [0.5; 0.5; 0.5];
for k = 1:10000
  g = grad_pen(x); if norm(g)<1e-10, break; end
  x = x - 0.0001*g;
end
printf("Lagrange: max xyz on sphere = %.6f at (%.4f,%.4f,%.4f)\n", ...
  x(1)*x(2)*x(3), x(1), x(2), x(3));
printf("Constraint: x²+y²+z² = %.6f (should be 3)\n", sum(x.^2));

15.10 — Divergence Theorem Verification

% F = (x³, y³, z³), div F = 3x²+3y²+3z² = 3r²
% ∭ div F dV over unit sphere in spherical coords
% = ∫₀²π ∫₀^π ∫₀¹ 3ρ² · ρ² sin φ dρ dφ dθ = 3 · (2π)(2)(1/5) = 12π/5

vol_int = quad(@(rho) 3*rho.^4, 0, 1) * ...
          quad(@(phi) sin(phi), 0, pi) * ...
          quad(@(theta) 1, 0, 2*pi);

% Surface integral: ∯ F·n̂ dS on unit sphere
% n̂ = (sin φ cos θ, sin φ sin θ, cos φ), |dS| = sin φ dφ dθ
% F on sphere: (cos³θ sin³φ, sin³θ sin³φ, cos³φ)
% F·n̂ = cos⁴θ sin⁴φ + sin⁴θ sin⁴φ + cos⁴φ sin φ (... simplified)
surf_int = dblquad(@(phi,theta) ...
  (cos(theta).^4).*(sin(phi).^4) + (sin(theta).^4).*(sin(phi).^4) + ...
  (cos(phi).^4).*sin(phi), 0, pi, 0, 2*pi);

printf("Divergence theorem:\n");
printf("  ∭ div F dV  = %.8f\n", vol_int);
printf("  ∯ F·dS      = %.8f\n", surf_int);
printf("  Exact 12π/5  = %.8f\n", 12*pi/5);

⚡ Quick Reference Card

Core Operations

Vector Identities

Integration Theorems

Coordinate Systems

Second Derivative Test (ℝ²)