Calculus of Variations

In standard calculus, we find a number $x$ that minimizes a function $f(x)$. In the Calculus of Variations, we find a function $y(x)$ that minimizes an integral $J[y]$, called a Functional.

$J[y] = \int_{x_1}^{x_2} L(x, y, y') dx$

This is the math behind “nature is lazy”: physics always chooses the path that minimizes “Action.”

The Euler-Lagrange Equation

To find the function $y(x)$ that makes $J[y]$ stationary (a minimum or maximum), we solve the Euler-Lagrange Equation:

$\frac{\partial L}{\partial y} - \frac{d}{dx} \left( \frac{\partial L}{\partial y'} \right) = 0$

Famous Problems

1. The Geodesic

What is the shortest path between two points? In flat space, the Euler-Lagrange equation tells us it is a straight line. On a curved surface (like Earth), it is a Great Circle.

2. The Brachistochrone

What shape of a wire allows a bead to slide from $A$ to $B$ in the shortest amount of time under gravity? Hint: It is not a straight line. It is a Cycloid (the path traced by a point on a rolling wheel).

3. Fermat’s Principle

Light travels between two points along the path that takes the least time. This principle alone allows us to derive Snell’s Law of refraction and the law of reflection.

import numpy as np

# Let's compare the time taken for a bead to slide down a straight line
# vs. a simple parabola (approximation of a cycloid)
g = 9.81
h = 10.0 # Height
L = 10.0 # Horizontal distance

# Straight line path
time_line = np.sqrt(2 * (L**2 + h**2) / (g * h))

# A curved path (e.g. y = x^2/10) actually allows the bead to pick up 
# speed faster at the start, potentially reducing total time.
print(f"Time for straight line: {time_line:.2f}s")
print("Curved paths allow the object to trade potential energy for kinetic energy")
print("earlier in the run, resulting in a faster average velocity.")

Use Cases

1. Classical Mechanics: Formulating the “Lagrangian” to find equations of motion.
2. Structural Engineering: Finding the shape of a bridge that minimizes stress.
3. Machine Learning: Variational Inference and optimizing neural network weights over a continuous space.

Exercises

Complex Analysis

Complex Analysis is the study of functions of a complex variable $z = x + iy$. It reveals a profound rigidity: if a complex function is differentiable once, it is differentiable infinitely many times.

Holomorphic Functions

A function $f(z) = u(x, y) + i v(x, y)$ is holomorphic if it is complex-differentiable. This requires the Cauchy-Riemann Equations: $\frac{\partial u}{\partial x} = \frac{\partial v}{\partial y}, \quad \frac{\partial u}{\partial y} = -\frac{\partial v}{\partial x}$

def check_cauchy_riemann(f, z, h=1e-7):
  # Numerical gradient check
  dz_x = f(z + h) - f(z)
  dz_y = f(z + h*1j) - f(z)
  
  df_dx = dz_x / h
  df_dy = dz_y / (h*1j)
  
  # In a holomorphic function, df/dz should be the same 
  # regardless of the direction of the limit.
  is_holomorphic = abs(df_dx - df_dy) < 1e-5
  return is_holomorphic, df_dx, df_dy

# Test f(z) = z^2 (Holomorphic)
f1 = lambda z: z**2
print(f"z^2 at 1+i: {check_cauchy_riemann(f1, 1+1j)}")

# Test f(z) = conj(z) (Non-holomorphic)
f2 = lambda z: z.conjugate()
print(f"conj(z) at 1+i: {check_cauchy_riemann(f2, 1+1j)}")

Complex Integration and Cauchy’s Theorem

If $f$ is holomorphic in a simply connected domain, the line integral along any closed path $\gamma$ is zero: $\oint_\gamma f(z) dz = 0$

However, if the domain has “holes” (singularities), we apply the Residue Theorem: $\oint_\gamma f(z) dz = 2\pi i \sum \text{Res}(f, z_k)$

Conformal Mappings

Holomorphic functions are conformal maps, meaning they preserve angles between curves. This property is vital in fluid dynamics and electromagnetism.

import numpy as np

def transform_grid(f, points):
  return [f(p) for p in points]

# Grid of points
grid = [x + y*1j for x in np.linspace(-1, 1, 5) for y in np.linspace(-1, 1, 5)]

# Exponential map: mapping lines to circles/rays
f_exp = lambda z: np.exp(z)
transformed = transform_grid(f_exp, grid)

print(f"First 5 Original: {grid[:5]}")
print(f"First 5 Transformed: {[f'{t:.2f}' for t in transformed[:5]]}")

Laurent Series and Singularities

Functions can be expanded into Laurent series near singularities: $f(z) = \sum_{n=-\infty}^{\infty} a_n (z-z_0)^n$

• Removable singularity: No negative power terms.
• Pole of order m: Finite negative terms up to $n = -m$.
• Essential singularity: Infinitely many negative power terms.

Derivatives: Rates of Change

The derivative is the mathematical tool for measuring change. If limits are about “approaching,” derivatives are about “moving.”

1. The Instantaneous Rate of Change

We define the derivative $f'(x)$ as the limit of the average rate of change as the interval goes to zero: $f'(x) = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}$

Instead of jumping to rules, let’s calculate the slope of $f(x) = x^2$ at $x=1$ by making $h$ smaller and smaller. This is exactly what your computer does when it performs Numerical Differentiation.

def f(x): return x**2

x_target = 1.0
h_values = [0.1, 0.01, 0.001, 0.0001]

print(f"{'h':<10} | {'Estimated Slope':<15}")
print("-" * 30)

for h in h_values:
  slope = (f(x_target + h) - f(x_target)) / h
  print(f"{h:<10} | {slope:<15.6f}")

print("-" * 30)
print("The limit seems to be 2.0, which matches the power rule (2x).")

2. Linear Approximation

Around a point $a$, we can approximate a complex function $f(x)$ with a simple line defined by the derivative: $f(x) \approx f(a) + f'(a)(x-a)$

This is why derivatives are so powerful: they turn complex, curvy problems into simple, linear ones.

import numpy as np
import matplotlib.pyplot as plt

def f(x): return np.sin(x)
def df(x): return np.cos(x) # Derivative of sin is cos

a = 1.0 # Point of approximation
x_range = np.linspace(0, 2, 100)

# Linear approx: L(x) = f(a) + f'(a)(x-a)
linear_approx = f(a) + df(a) * (x_range - a)

plt.plot(x_range, f(x_range), label='True sin(x)', linewidth=2)
plt.plot(x_range, linear_approx, '--', label='Linear Approx', color='r')
plt.scatter([a], [f(a)], color='black', zorder=5)
plt.title("Linear Approximation (Tangent Line) at x=1")
plt.legend()
plt.show()

3. The Rules of Calculus

To avoid numerical limits, we use rules. The Chain Rule is particularly vital for modern AI, as it powers backpropagation. $\frac{d}{dx} f(g(x)) = f'(g(x)) \cdot g'(x)$

4. Optimization: Finding Extremes

If a smooth function reaches a peak or valley, its slope must be zero. These are Critical Points.

import numpy as np
import matplotlib.pyplot as plt

# f(x) = x^4 - 4x^2
def f(x): return x**4 - 4*x**2
def df(x): return 4*x**3 - 8*x # 4x(x^2 - 2) = 0 => x=0, x=±sqrt(2)

x = np.linspace(-2.5, 2.5, 500)
plt.plot(x, f(x))
plt.axhline(0, color='k', alpha=0.3)
plt.title("Find the valleys: f'(x) = 0")
plt.show()

5. Summary Check

Differential Equations

In basic calculus, we learn how to find the derivative of a function. In the real world, we often encounter the reverse problem: we know how a system changes, and we want to find the function that describes the system’s state. These are Differential Equations (DEs).

A DE is an equation that relates a function $y$ to its derivatives. For example: $\frac{dy}{dx} = ky$ This says that the rate of change of $y$ is proportional to its current value. This is the model for population growth, radioactive decay, and continuously compounded interest.

Separation of Variables

The simplest way to solve a first-order DE is to “separate” the variables so that all $y$‘s are on one side and all $x$‘s are on the other. Example: Solve $\frac{dy}{dx} = xy$.

1. $\frac{1}{y} dy = x dx$
2. $\int \frac{1}{y} dy = \int x dx$
3. $\ln|y| = \frac{1}{2}x^2 + C$
4. $y = Ae^{x^2/2}$ (where $A = e^C$)

Modeling: Newton’s Law of Cooling

Newton’s Law of Cooling states that the rate of change of the temperature $T$ of an object is proportional to the difference between its temperature and the ambient temperature $T_a$: $\frac{dT}{dt} = -k(T - T_a)$

This equation tells us that a hot cup of coffee cools down quickly at first (when the difference is large) and then slows down as it approaches room temperature.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define the model: dT/dt = -k(T - Ta)
def model(T, t, k, Ta):
  dTdt = -k * (T - Ta)
  return dTdt

# Initial conditions
T0 = 90       # Coffee temperature at t=0
Ta = 20       # Room temperature
k = 0.1       # Cooling constant
t = np.linspace(0, 60, 100) # 60 minutes

# Solve ODE
T = odeint(model, T0, t, args=(k, Ta))

plt.plot(t, T)
plt.axhline(Ta, color='r', linestyle='--', label='Room Temp')
plt.xlabel('Time (min)')
plt.ylabel('Temperature (C)')
plt.title('Newton\'s Law of Cooling')
plt.legend()
plt.grid(True)
plt.show()

print(f"Temperature at 10 mins: {T[16][0]:.2f} C")
print(f"Temperature at 60 mins: {T[-1][0]:.2f} C")

Linear Systems and Matrices

When we have multiple interdependent variables (like a predator and its prey, or a set of connected water tanks), we use Systems of Differential Equations. $\mathbf{x}'(t) = A\mathbf{x}(t)$ The solution to this system is deeply connected to the Eigenvalues and Eigenvectors of the matrix $A$. In fact, the general solution is based on the matrix exponential $e^{At}$.

Exercises

Differential Forms and Integration

Differential forms provide a unified framework for the study of integration, Stokes’ theorem, and the geometry of manifolds.

Exterior Algebra and Wedge Products

At each point of a manifold, a differential $k$-form is a purely antisymmetric $(0, k)$-tensor. The wedge product $\wedge$ allows us to combine forms: $\omega \wedge \eta = (-1)^{p q} \eta \wedge \omega$ where $\omega$ is a $p$-form and $\eta$ is a $q$-form.

import numpy as np

def wedge_product_2d(v1, v2):
  """Calculates the 2-form (area) spanned by two vectors in R3."""
  # This is effectively the cross product for 1-forms in R3
  return np.cross(v1, v2)

a = np.array([1, 0, 0])
b = np.array([0, 1, 0])

area_element = wedge_product_2d(a, b)
print(f"Area element (wedge product): {area_element}")
print(f"Magnitude (Area): {np.linalg.norm(area_element)}")

The Exterior Derivative

The exterior derivative $d$ is a operator that transforms a $k$-form into a $(k+1)$-form. It generalizes the concepts of gradient, curl, and divergence from vector calculus.

• If $f$ is a 0-form (function), $df$ is its gradient.
• $d(d\omega) = 0$ for any form $\omega$.

Generalized Stokes’ Theorem

The fundamental theorem of calculus, Green’s theorem, and the divergence theorem are all special cases of the generalized Stokes’ theorem: $\int_{\Omega} d\omega = \int_{\partial \Omega} \omega$ where $\partial \Omega$ is the boundary of $\Omega$.

Cohomology

The study of closed forms ($d\omega = 0$) that are not exact ($\omega \neq d\eta$) leads to De Rham Cohomology, which reveals topological information about the underlying space (e.g., the presence of holes).

Fourier Analysis

Fourier Analysis allows us to decompose complex, periodic signals into a sum of simple sine and cosine waves.

Fourier Series

Any periodic function $f(x)$ with period $2L$ can be represented as: $f(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left( a_n \cos\frac{n\pi x}{L} + b_n \sin\frac{n\pi x}{L} \right)$

import numpy as np

def square_wave_fourier(x, terms):
  res = 0
  for n in range(1, terms + 1, 2):
      # f(x) = (4/pi) * summation sin(n*x)/n
      res += (4 / np.pi) * (np.sin(n * x) / n)
  return res

x_val = np.pi/2
print(f"Square wave approx (5 terms) at pi/2: {square_wave_fourier(x_val, 5)}")
print(f"Square wave approx (100 terms) at pi/2: {square_wave_fourier(x_val, 100)}")

The Fourier Transform

For non-periodic functions, the Fourier Transform $\hat{f}(\xi)$ maps a function from the time (or space) domain to the frequency domain: $\hat{f}(\xi) = \int_{-\infty}^{\infty} f(x) e^{-2\pi i x \xi} dx$

import numpy as np

# Fast Fourier Transform (FFT)
signal = np.array([1, 0, 1, 0, 1, 0, 1, 0])
fft_res = np.fft.fft(signal)

print(f"Signal: {signal}")
print(f"FFT Magnitudes: {np.abs(fft_res).round(2)}")
print("The spike at index 4 represents the fundamental frequency of the oscillation.")

Parseval’s Theorem

Parseval’s Theorem states that the total energy of a signal in the time domain is equal to the total energy in the frequency domain: $\int |f(x)|^2 dx = \int |\hat{f}(\xi)|^2 d\xi$

Functional Analysis

Functional analysis studies vector spaces endowed with a limit-related structure (like a metric or topology) and the linear operators acting upon them.

1. Normed Linear Spaces

A Normed Linear Space $(V, \|\cdot\|)$ is a vector space $V$ over a field (usually $\mathbb{C}$ or $\mathbb{R}$) with a norm function $\|\cdot\|: V \to [0, \infty)$ satisfying:

1. Definiteness: $\|v\| = 0 \iff v = 0$.
2. Homogeneity: $\|\alpha v\| = |\alpha| \cdot \|v\|$.
3. Triangle Inequality: $\|u + v\| \le \|u\| + \|v\|$.

A Banach Space is a normed linear space that is complete (every Cauchy sequence converges).

2. Hilbert Spaces and Inner Products

A Hilbert Space is a complete inner product space. The inner product $\langle \cdot, \cdot \rangle$ satisfies:

1. Linearity: $\langle \alpha u + \beta v, w \rangle = \alpha \langle u, w \rangle + \beta \langle v, w \rangle$.
2. Conjugate Symmetry: $\langle u, v \rangle = \overline{\langle v, u \rangle}$.
3. Positive Definiteness: $\langle v, v \rangle \ge 0$, with equality only if $v=0$.

The $L^2$ Space is a primary example, where functions are square-integrable: $\langle f, g \rangle = \int_a^b f(x) \overline{g(x)} \, dx, \quad \|f\|_2 = \sqrt{\int_a^b |f(x)|^2 \, dx}$

import numpy as np

# Numerical verification of the Triangle Inequality
# ||f + g|| <= ||f|| + ||g||
x = np.linspace(0, 1, 500)
f = x**2
g = np.exp(x)

def norm_l2(y, dx):
  return np.sqrt(np.trapz(y**2, dx))

dx = x[1] - x[0]
n_f = norm_l2(f, x)
n_g = norm_l2(g, x)
n_sum = norm_l2(f + g, x)

print(f"||f||: {n_f:.4f}")
print(f"||g||: {n_g:.4f}")
print(f"||f+g||: {n_sum:.4f}")
print(f"Is {n_sum:.4f} <= {n_f + n_g:.4f}? {n_sum <= n_f + n_g + 1e-9}")

3. Linear Operators

1. Partial Differential Equations: Proving that a solution exists even if we can’t find it.
2. Signal Processing: Decomposing complex signals into simple frequencies.
3. Quantum Mechanics: The state of a particle is a vector in an infinite-dimensional Hilbert space.

Exercises

The Fundamental Theorem of Calculus

The Fundamental Theorem of Calculus (FTC) is the crown jewel of mathematics. It connects two seemingly unrelated concepts: the derivative (slopes) and the integral (areas). Before this discovery, calculating areas was a tedious process of infinite sums. After the FTC, it became a simple matter of finding an anti-derivative.

Part 1: The Area Function

The first part of the theorem states that if we define an “area function” $F(x)$ that accumulates the area under a function $f(t)$ from a fixed point $a$ to $x$: $F(x) = \int_a^x f(t) \, dt$ Then the derivative of this area function is simply the original function: $F'(x) = f(x)$

The Intuition: The rate at which area is being added at point $x$ is exactly equal to the “height” of the function at that point. If the function is tall, the area grows quickly. If it is zero, the area stops growing.

Part 2: The Shortcut to Integration

The second part provides the formula we use in practice to evaluate definite integrals. If $F$ is any anti-derivative of $f$ (meaning $F' = f$), then: $\int_a^b f(x) \, dx = F(b) - F(a)$

This is revolutionary. To find the area under a curve, you don’t need to draw rectangles. You just need to find the “inverse” of the derivative and plug in the endpoints.

import sympy as sp

# Define the variable and the function
x = sp.Symbol('x')
f = x**2

# Find the anti-derivative (Indefinite Integral)
F = sp.integrate(f, x)
print(f"Function: f(x) = {f}")
print(f"Anti-derivative: F(x) = {F}")

# Evaluate the definite integral from 0 to 3 using FTC Part 2
# Integral = F(3) - F(0)
a, b = 0, 3
result = F.subs(x, b) - F.subs(x, a)

print(f"\nIntegral from {a} to {b} of x^2 dx:")
print(f"F({b}) - F({a}) = {F.subs(x, b)} - {F.subs(x, a)} = {result}")

# Verify with sympy's definite integral tool
assert result == sp.integrate(f, (x, a, b))

Why “Fundamental”?

The FTC is fundamental because it shows that differentiation and integration are inverse operations. It turns a geometry problem (area) into an algebra problem (anti-derivatives). Every time a physicist calculates the energy lost by a falling object or an engineer calculates the total stress on a bridge, they are relying on this 300-year-old bridge between two worlds.

Exercises

Integration: The Art of Accumulation

While derivatives break a function down into its local rates of change, integration builds it back up. Integration is the process of adding up infinitely many tiny pieces to find a whole.

1. The Definite Integral as a Sum

We define the definite integral $\int_a^b f(x) \, dx$ as the signed area under the curve $f(x)$ from $a$ to $b$. Formally, this is reached through a Riemann Sum: we divide the area into $n$ rectangles and take the limit as $n \to \infty$.

Instead of just looking at the formula, let’s watch the approximation get better as we add more rectangles.

import numpy as np
import matplotlib.pyplot as plt

def f(x): return x * np.sin(x) + 2

a, b = 0, 10
n = 10 # Change this to 5, 20, or 100 to see convergence

x = np.linspace(a, b, 1000)
x_rect = np.linspace(a, b, n, endpoint=False)
width = (b - a) / n

plt.plot(x, f(x), 'k', linewidth=2)
plt.bar(x_rect, f(x_rect), width=width, align='edge', alpha=0.3, color='blue', edgecolor='b')
plt.title(f"Riemann Sum Approximation (n={n})")
plt.show()

2. Numerical Integration: Trapezoidal Rule

Computer scientists rarely integrate by hand. They use algorithms like the Trapezoidal Rule, which fits a line (trapezoid) between points instead of a flat rectangle. This usually converges much faster.

import numpy as np

def f(x): return np.exp(-x**2) # The Gaussian bell curve

a, b = 0, 1
n = 5
h = (b - a) / n

x = np.linspace(a, b, n + 1)
y = f(x)

# Trapezoidal rule: h/2 * (y0 + 2y1 + ... + yn)
area = (h/2) * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])

print(f"Estimated area under e^(-x^2) from 0 to 1: {area:.6f}")
print("This is related to the Error Function (erf) used in Statistics.")

3. The Fundamental Theorem of Calculus (FTC)

The FTC states that integration and differentiation are inverse operations. If $F(x)$ is the antiderivative of $f(x)$, then: $\int_a^b f(x) \, dx = F(b) - F(a)$

This turns a hard problem of “infinite summing” into a simple problem of “subtraction.”

4. Improper Integrals

Sometimes we need to integrate over an infinite interval, such as $\int_0^\infty e^{-x} \, dx$. These are vital for calculating total energy or probabilities in bell curves.

5. Summary Check

Laplace Transforms

The Laplace Transform is a powerful tool used in engineering and physics to solve linear differential equations. Its primary “trick” is to transform a function from the Time Domain ($t$) to the Complex Frequency Domain ($s$), where integration and differentiation become simple algebraic multiplication and division.

The Definition

The Laplace transform of a function $f(t)$ is defined as:

$\mathcal{L}\{f(t)\} = F(s) = \int_0^\infty e^{-st} f(t) dt$

Why Transform?

In the time domain, a system might be described by a messy differential equation: $a \frac{d^2y}{dt^2} + b \frac{dy}{dt} + cy = f(t)$

In the Laplace (s) domain, this becomes a simple equation: $(as^2 + bs + c)Y(s) = F(s)$

You can then solve for $Y(s)$ using basic algebra and “Inverse Transform” back to the time domain.

Transfer Functions and Stability

In control theory, the Transfer Function $H(s) = \frac{Y(s)}{X(s)}$ describes how a system responds to an input.

• If the “poles” (the values of $s$ that make the denominator zero) have negative real parts, the system is Stable (vibrations die down).
• If any pole has a positive real part, the system is Unstable (vibrations grow until the system breaks).

import sympy as sp

# Define symbols
t, s = sp.symbols('t s')
# Define a function to transform: sin(t)
f = sp.sin(t)

# Calculate Laplace Transform
F = sp.laplace_transform(f, t, s)

print(f"Time Domain: f(t) = {f}")
print(f"S-Domain (Laplace): F(s) = {F[0]}")
print("\nNotice how the trigonometric function became a rational algebraic fraction.")

Use Cases

1. Circuit Analysis: Capacitors and Inductors become simple algebraic impedances ($1/sC$ and $sL$).
2. Control Systems: Tuning a thermostat or a drone’s flight controller.
3. Signal Processing: Filtering noise out of audio or sensor data.

Exercises

Limits and Continuity

Calculus is built on the concept of change over zero. To make sense of “instantaneous” change, we need a way to talk about what a function does as it gets closer and closer to a point, even if it never actually reaches it. This is the Limit.

1. The Epsilon-Delta Definition

Intuition: A limit $\lim_{x \to c} f(x) = L$ means we can force $f(x)$ to be as close to $L$ as we want (within $\epsilon$), just by making $x$ sufficiently close to $c$ (within $\delta$).

For every $\epsilon > 0$, there exists a $\delta > 0$ such that $0 < |x - c| < \delta \implies |f(x) - L| < \epsilon$.

Let’s vizualize this. If $f(x) = 2x$, and we want to be within $\epsilon = 0.1$ of the output $L=4$ at $x=2$, how close must $x$ be?

import numpy as np
import matplotlib.pyplot as plt

def f(x): return 2 * x

c, L = 2, 4
epsilon = 0.1
delta = epsilon / 2 # For f(x)=2x, delta is exactly epsilon/2

x = np.linspace(c - 2*delta, c + 2*delta, 500)
y = f(x)

plt.plot(x, y, label='f(x)=2x')
plt.axhline(L + epsilon, color='r', linestyle='--', label='L + ε')
plt.axhline(L - epsilon, color='r', linestyle='--', label='L - ε')
plt.axvline(c + delta, color='g', linestyle='--', label='c + δ')
plt.axvline(c - delta, color='g', linestyle='--', label='c - δ')

plt.fill_between([c-delta, c+delta], L-epsilon, L+epsilon, color='yellow', alpha=0.3, label='Safe Zone')
plt.legend()
plt.title("The Epsilon-Delta Game")
plt.show()

2. When Limits Fail

A limit only exists if it matches from both directions. If the “safe zone” can’t be established because the function jumps or oscillates, the limit does not exist.

3. Continuity: The Glue of Calculus

A function is continuous at $c$ if the limit exists and exactly matches the function’s value: $\lim_{x \to c} f(x) = f(c)$.

Let’s test for a “Broken” function: $f(x) = \begin{cases} x^2 & x < 1 \\ 2 & x = 1 \\ x & x > 1 \end{cases}$

import numpy as np

def f(x):
  if x < 1: return x**2
  if x == 1: return 2
  return x

# Check limit from left and right
left_lim = f(0.999999)
right_lim = f(1.000001)
value = f(1.0)

print(f"Limit from left: {left_lim}")
print(f"Limit from right: {right_lim}")
print(f"Actual value at 1: {value}")

is_continuous = np.isclose(left_lim, right_lim) and np.isclose(left_lim, value)
print(f"Is continuous? {is_continuous}")

4. The Intermediate Value Theorem (IVT)

If $f$ is continuous on $[a, b]$, it must hit every value between $f(a)$ and $f(b)$. This is why we can use the Bisection Method to find roots.

5. Summary Check

Mean Value Theorem: The Guarantee of Speed

The Mean Value Theorem (MVT) is often viewed as a “theoretical” result used by mathematicians for proofs. However, it has a very intuitive and practical meaning: if you travel 100 miles in one hour, there must have been at least one instant where your speedometer read exactly 100 mph.

Rolle’s Theorem: The Foundation

Rolle’s Theorem is a special case of the MVT. It states that if a continuous and differentiable function $f(x)$ starts and ends at the same value ($f(a) = f(b)$), there must be at least one point $c$ in between where the derivative is zero ($f'(c) = 0$).

Think of it like this: if you throw a ball up and it comes back to your hand, there was a moment at the very top of its path where its vertical velocity was exactly zero.

The General Mean Value Theorem

The MVT generalizes Rolle’s Theorem to functions that don’t end where they start. It states that for a differentiable function $f(x)$ on $[a, b]$, there exists a point $c \in (a, b)$ such that: $f'(c) = \frac{f(b) - f(a)}{b - a}$

In plain English: The instantaneous rate of change ($f'(c)$) must equal the average rate of change over the interval at some point.

Why do we care?

The MVT is the bridge between local behavior (the derivative at a point) and global behavior (the function values at the endpoints).

• It allows us to bound the error in numerical approximations.
• It proves that if $f'(x) = 0$ everywhere, $f(x)$ must be a constant.
• It is the primary tool used to prove the Fundamental Theorem of Calculus.

import numpy as np
import matplotlib.pyplot as plt

# Interval [a, b]
a, b = 0, 2

# Function f(x) = x^2
def f(x): return x**2
def df(x): return 2*x

# Average rate of change (Secant slope)
avg_rate = (f(b) - f(a)) / (b - a)

# Find 'c' where f'(c) == avg_rate => 2c = 2 => c = 1
c = 1

x = np.linspace(-0.5, 2.5, 100)
plt.plot(x, f(x), label='f(x) = x^2')
plt.plot([a, b], [f(a), f(b)], '--', color='red', label='Average Rate (Secant)')
plt.scatter([c], [f(c)], color='green', zorder=5)

# Tangent at c
tangent_x = np.linspace(0.5, 1.5, 50)
tangent_y = f(c) + df(c)*(tangent_x - c)
plt.plot(tangent_x, tangent_y, color='green', label='Instantaneous Rate (Tangent)')

plt.title("Mean Value Theorem: Parallel Slopes")
plt.legend()
plt.grid(True)
plt.show()

print(f"Average Rate over [{a}, {b}]: {avg_rate}")
print(f"Instantaneous Rate at c={c}: {df(c)}")

Taylor’s Theorem: MVT on Steroids

If the MVT tells us how to approximate a function with a line, Taylor’s Theorem tells us how to approximate it with a polynomial of any degree. $f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \dots + R_n(x)$ The “remainder” $R_n(x)$ is defined using a generalization of the Mean Value Theorem. This allows engineers to know exactly how many terms they need to keep in a computer simulation to ensure the error remains below a certain threshold.

Exercises

Multiple Integration

In 1D calculus, we integrate over a line to find an area. In Multivariable Calculus, we integrate over a 2D region (Double Integral) to find a volume, or over a 3D region (Triple Integral) to find a total mass or charge.

Double Integrals: Volumes under Surfaces

The integral $\iint_R f(x, y) dA$ represents the volume between the $xy$-plane and the surface $z = f(x, y)$ over the region $R$.

Fubini’s Theorem: If the function is nice, you can calculate the double integral by doing two “nested” single integrals in any order: $\int_a^b \int_c^d f(x, y) dy \, dx = \int_c^d \int_a^b f(x, y) dx \, dy$

The Jacobian: Scaling Space

When we change variables (e.g., from Cartesian $(x, y)$ to Polar $(r, \theta)$), the “infinitesimal area” $dA$ changes. We use the Jacobian to account for this stretching.

• Cartesian: $dA = dx \, dy$
• Polar: $dA = r \, dr \, d\theta$

If you forget the $r$ in polar coordinates, your areas and volumes will be wrong!

from scipy import integrate
import numpy as np

# Let's calculate the volume of a hemisphere of radius 1
# Surface: z = sqrt(1 - x^2 - y^2)
# We will use SciPy to integrate over a circular region

def f(y, x):
  if x**2 + y**2 <= 1:
      return np.sqrt(1 - x**2 - y**2)
  return 0

# Integrate over the square [-1, 1] x [-1, 1]
# Note: we filter outside the circle in the function f
volume, error = integrate.dblquad(f, -1, 1, lambda x: -1, lambda x: 1)

print(f"Calculated Volume: {volume:.4f}")
print(f"Theoretical (2/3 * pi): { (2/3) * np.pi:.4f}")

Triple Integrals and Density

If $\rho(x, y, z)$ represents the density of an object at point $(x, y, z)$, the total mass of the object is the triple integral of the density over its volume $V$:

$M = \iiint_V \rho(x, y, z) dV$

Exercises

Multivariable Calculus

Multivariable calculus extends the concepts of single-variable calculus to functions of several variables.

Partial Derivatives and the Gradient

A partial derivative measures the rate of change of a multivariable function with respect to one variable, holding others constant. The Gradient $\nabla f$ is a vector of partial derivatives: $\nabla f = \left( \frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \dots, \frac{\partial f}{\partial x_n} \right)$ It points in the direction of the steepest ascent.

def f(x, y):
  return x**2 + 3*y**2

def gradient_f(x, y):
  # Partial w.r.t x: 2x
  # Partial w.r.t y: 6y
  return (2*x, 6*y)

x0, y0 = 1.0, 2.0
grad = gradient_f(x0, y0)
print(f"Gradient at ({x0}, {y0}): {grad}")

Multiple Integrals

A double integral $\iint_D f(x, y) \, dA$ calculates the volume under the surface $z = f(x, y)$ over a region $D$ in the $xy$-plane.

from scipy.integrate import dblquad

# f(x, y) = x*y
# Over [0, 1] x [0, 2]
area, error = dblquad(lambda y, x: x*y, 0, 1, lambda x: 0, lambda x: 2)

print(f"Double integral of x*y over [0,1]x[0,2]: {area:.4f}")

Vector Fields and Line Integrals

A vector field $\mathbf{F}$ assigns a vector to every point in space. A line integral $\int_C \mathbf{F} \cdot d\mathbf{r}$ measures the work done by a force field along a path $C$.

The Jacobian Matrix

The Jacobian matrix of a vector-valued function $\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m$ is the matrix of all first-order partial derivatives. It represents the best linear approximation of the function near a point. $J_{ij} = \frac{\partial f_i}{\partial x_j}$

Mathematical Optimization

Optimization is the selection of a best element (with regard to some criterion) from some set of available alternatives. In its simplest form, this means finding the values that maximize or minimize a function.

Unconstrained Optimization

For a differentiable function $f(x)$, local extrema occur at points where $f'(x) = 0$. The Second Derivative Test determines if a point is a maximum ($f'' < 0$) or a minimum ($f'' > 0$).

import numpy as np

def f(x):
  return x**4 - 4*x**2 + 3

def df(x):
  return 4*x**3 - 8*x

def ddf(x):
  return 12*x**2 - 8

# Roots of df(x) are 0, sqrt(2), -sqrt(2)
critical_points = [0, np.sqrt(2), -np.sqrt(2)]

for p in critical_points:
  second_der = ddf(p)
  kind = "Min" if second_der > 0 else "Max"
  print(f"Point x={p:.2f} is a {kind} (f''={second_der:.2f})")

Gradient Descent

In higher dimensions, we often use iterative methods like Gradient Descent to find a local minimum: $\mathbf{x}_{n+1} = \mathbf{x}_n - \gamma \nabla f(\mathbf{x}_n)$ where $\gamma$ is the learning rate.

def f(x, y):
  return x**2 + y**2 + 10

def grad_f(x, y):
  return np.array([2*x, 2*y])

point = np.array([5.0, 5.0])
gamma = 0.1

for _ in range(10):
  point = point - gamma * grad_f(*point)
  print(f"Current point: {point}")

Constrained Optimization: Lagrange Multipliers

To optimize $f(x, y)$ subject to $g(x, y) = c$, we look for points where the gradients are parallel: $\nabla f = \lambda \nabla g$

Partial Differential Equations (PDEs)

Ordinary Differential Equations (ODEs) deal with functions of one variable (usually time). Partial Differential Equations (PDEs) deal with functions of multiple variables, such as the temperature $u(x, t)$ at a position $x$ and time $t$.

The Mother of All PDEs: The Laplacian

Most physical PDEs involve the Laplacian ($\nabla^2 u = \Delta u$), which is the divergence of the gradient ($\nabla \cdot \nabla u$). In 1D, this is just $u_{xx}$ (the second spatial derivative). It measures how much the value at a point differs from the average of its neighbors.

Three Classic PDEs

1. The Heat Equation (Diffusion)

$u_t = \alpha \nabla^2 u$ Meaning: The rate of change in temperature ($u_t$) is proportional to how much “different” a point is from its neighbors. Heat flows from hot to cold to “smooth out” the distribution.

2. The Wave Equation

$u_{tt} = c^2 \nabla^2 u$ Meaning: Acceleration ($u_{tt}$) is proportional to the local “curvature” of the medium. This models everything from guitar strings to light waves.

3. Laplace’s Equation (Equilibrium)

$\nabla^2 u = 0$ Meaning: The system has reached a state where every point is the average of its neighbors. This describes electric potentials in vacuum or the steady-state temperature of a plate.

Numerical Solution: The Finite Difference Method

Solving PDEs analytically is hard. Engineers use the Finite Difference Method, which replaces derivatives with differences on a grid.

import numpy as np

# Simulate 1D Heat Diffusion in a rod
L = 1.0       # Length
nx = 50       # Number of points
dx = L / (nx - 1)
alpha = 0.01  # Thermal diffusivity
dt = 0.001    # Time step

# Initial condition: hot spot in the middle
u = np.zeros(nx)
u[int(0.4*nx):int(0.6*nx)] = 100.0 

# Time loop: simplified explicit scheme
for _ in range(100):
  u_new = u.copy()
  for i in range(1, nx-1):
      # u_t = alpha * u_xx
      u_new[i] = u[i] + alpha * dt / dx**2 * (u[i+1] - 2*u[i] + u[i-1])
  u = u_new

print(f"Temperature at middle: {u[nx//2]:.2f}")
print("The heat has started to spread from the center to the edges.")

Exercises

Sequences and Series

Calculus is the study of the infinite. While derivatives and integrals use infinite processes to find slopes and areas, sequences and series study infinite processes directly. How can you add up an infinite number of things and get a finite answer? This question is at the heart of Zeno’s paradoxes and modern data compression.

Limits of Sequences

A sequence is an infinite list of numbers $a_1, a_2, a_3, \dots$. We are interested in its convergence: does the sequence eventually settle down to a specific value $L$? $\lim_{n \to \infty} a_n = L$

In computer science, every iterative algorithm (like Newton’s method or Gradient Descent) is a sequence. We need to prove that these sequences converge to the correct solution.

Infinite Series: Summing to Infinity

A series is the sum of the terms of a sequence: $\sum_{n=1}^\infty a_n$. A series converges if the sequence of its partial sums $S_n = a_1 + \dots + a_n$ converges to a finite limit.

The Geometric Series

The most famous convergent series is the geometric series: $1 + r + r^2 + r^3 + \dots = \frac{1}{1-r} \quad \text{for } |r| < 1$ This formula is used to calculate the “Present Value” of future cash flows in finance and the “Multiplier Effect” in economics.

Taylor Series: Polynomials for Everything

The ultimate application of series is the Taylor Series. It allows us to represent transcendental functions ($e^x, \sin x, \ln x$) as infinite polynomials. $f(x) = \sum_{n=0}^\infty \frac{f^{(n)}(a)}{n!}(x-a)^n$

This is how calculators and computers compute $\sin(1)$. They don’t have a giant lookup table; they use the first few terms of the Taylor series to get as much precision as needed.

import numpy as np
import matplotlib.pyplot as plt

def taylor_exp(x, n):
  """Approximate e^x using first n terms of Taylor series at a=0."""
  approx = 0
  factorial = 1
  for i in range(n):
      if i > 0: factorial *= i
      approx += (x**i) / factorial
  return approx

x_val = 1.0
print(f"Exact e^{x_val}: {np.exp(x_val)}")
for n in [1, 3, 5, 10]:
  print(f"Taylor (n={n}): {taylor_exp(x_val, n)}")

# Visualizing approximation
x = np.linspace(-2, 2, 100)
plt.plot(x, np.exp(x), 'k', label='Exact exp(x)')
plt.plot(x, taylor_exp(x, 3), '--', label='n=3 Approximation')
plt.plot(x, taylor_exp(x, 5), ':', label='n=5 Approximation')
plt.ylim(-1, 8)
plt.legend()
plt.title("Taylor Expansion of e^x")
plt.grid(True)
plt.show()