Calculus Essentials

A limit asks: as an input value approaches some target, what value does the function approach? It does not ask what the function equals at that point — it asks about the trend as we get arbitrarily close.

Formally, we write: $$\lim_{x \to a} f(x) = L$$ This means: for every $\varepsilon > 0$, there exists a $\delta > 0$ such that whenever $0 < |x - a| < \delta$, we have $|f(x) - L| < \varepsilon$. This is the epsilon-delta definition (Weierstrass, 1850s), and it is the rigorous underpinning of everything in calculus.

The epsilon-delta definition sounds mechanical but embodies a profound idea: we can reason about infinite processes by making the error provably small. This is exactly what numerical analysts do when they prove a floating-point algorithm is accurate to within machine epsilon.

Intuition via a CS analogy. Imagine binary search. As you double the number of iterations, your error halves. The limit of the error as iterations go to infinity is 0. You never literally run infinitely many iterations, but the limit tells you the trend. Calculus formalizes this kind of reasoning.

Key limit laws:

• Limits of sums, products, and quotients factor through (when the denominator limit is nonzero)
• The Squeeze Theorem: if $g(x) \le f(x) \le h(x)$ near $a$, and both $g$ and $h$ approach $L$, then $f$ also approaches $L$
• The famous limit $\lim_{x \to 0} \frac{\sin x}{x} = 1$ (proved by squeezing between $\cos x$ and 1)

L'Hopital's Rule (previewing derivatives): if both numerator and denominator approach 0, the limit of their ratio equals the limit of the ratio of their derivatives — a huge computational shortcut.

Why CS cares. Every asymptotic analysis in algorithms (Big-O, $\Theta$, $\Omega$) is implicitly about limits: $f(n) = O(g(n))$ means $\limsup_{n\to\infty} f(n)/g(n) < \infty$. You use limits every time you talk about convergence of an iterative algorithm.

import math

def limit_sin_over_x(epsilon=1e-9):
    """Numerically verify lim_{x->0} sin(x)/x = 1.
    We approach from both sides with decreasing step."""
    results = []
    for exp in range(1, 15):
        x = 10 ** (-exp)
        val = math.sin(x) / x
        results.append((x, val))
        print(f"x={x:.1e}  sin(x)/x = {val:.15f}")
    return results

limit_sin_over_x()

A function $f$ is continuous at a point $a$ if three conditions hold simultaneously:

1. $f(a)$ is defined.
2. $\lim_{x \to a} f(x)$ exists.
3. The limit equals the value: $\lim_{x \to a} f(x) = f(a)$.

Intuitively: you can draw the function without lifting your pencil. Continuity means no sudden jumps, holes, or vertical asymptotes at that point.

Types of discontinuities:

• Removable: the limit exists but the function has a hole (e.g., $f(x) = \frac{\sin x}{x}$ at $x = 0$ before we define $f(0) = 1$)
• Jump: left and right limits both exist but differ (e.g., the floor function at integers)
• Infinite: the function blows up (e.g., $1/x$ at 0)

The Intermediate Value Theorem (IVT): if $f$ is continuous on $[a,b]$ and $f(a) < y < f(b)$, then there exists $c \in (a,b)$ with $f(c) = y$. This sounds obvious but is deep: it depends on the completeness of the real numbers. The IVT is the theoretical engine behind bisection root-finding — you know a root exists in an interval, so you keep halving. Every numerical root finder in a scientific computing library is implicitly applying the IVT.

The Extreme Value Theorem: a continuous function on a closed bounded interval $[a,b]$ attains both its maximum and minimum. This is why gradient descent can find a minimum for a continuous loss function on compact domain.

CS relevance. Activation functions in neural networks must be continuous (and usually differentiable) to allow smooth gradient flow. The step function is discontinuous — this is why it was abandoned for sigmoid and ReLU. Hash functions are deliberately discontinuous (a tiny input change gives a totally different output).

# Numerical IVT: bisection root-finding for f(x) = x^3 - x - 2
# The IVT guarantees a root in [1, 2] because f(1) = -2 < 0 < f(2) = 4

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

def bisect(a, b, tol=1e-10):
    assert f(a) * f(b) < 0, "IVT requires opposite signs at endpoints"
    while (b - a) > tol:
        mid = (a + b) / 2
        if f(mid) == 0:
            return mid
        elif f(a) * f(mid) < 0:
            b = mid
        else:
            a = mid
    return (a + b) / 2

root = bisect(1, 2)
print(f"Root: {root:.10f}")  # ≈ 1.5213797068

The derivative of $f$ at $a$ is defined as the limit of the difference quotient: $$f'(a) = \lim_{h \to 0} \frac{f(a+h) - f(a)}{h}$$

Geometrically, $(f(a+h)-f(a))/h$ is the slope of the secant line through $(a, f(a))$ and $(a+h, f(a+h))$. As $h \to 0$, the secant line pivots until it becomes the tangent line at $a$. The derivative is the slope of that tangent line.

Physically: if $f(t)$ is position at time $t$, then $f'(t)$ is instantaneous velocity — not average velocity over an interval, but the exact rate at one instant.

Differentiability implies continuity (but not vice versa — the absolute value function $|x|$ is continuous everywhere but not differentiable at 0).

Basic differentiation rules (derived from the limit definition):

• Power rule: $(x^n)' = nx^{n-1}$
• Sum rule: $(f+g)' = f' + g'$
• Product rule: $(fg)' = f'g + fg'$
• Quotient rule: $(f/g)' = (f'g - fg') / g^2$
• Chain rule: see next concept

Notation. Leibniz's $\frac{df}{dx}$ (or $\frac{dy}{dx}$) is preferred in applied math and physics because it makes the chain rule look like fraction cancellation: $\frac{dy}{dt} = \frac{dy}{dx} \cdot \frac{dx}{dt}$. Newton's dot notation $\dot{f}$ is common in physics for time derivatives. Lagrange's $f'$ is common in pure math.

Common derivatives to memorize:

That last one — the sigmoid — is something every ML practitioner derives. It is a beautiful special case of the chain rule and quotient rule combined.

Numerical differentiation approximates the derivative using the centered difference: $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$ for small $h$. This has error $O(h^2)$ vs $O(h)$ for the one-sided formula.

# Compare analytical vs numerical derivative of f(x) = x^3 - 3x
# f'(x) = 3x^2 - 3  (by power rule)

def f(x):
    return x**3 - 3*x

def f_prime_analytical(x):
    return 3*x**2 - 3

def f_prime_numerical(x, h=1e-7):
    return (f(x + h) - f(x - h)) / (2 * h)  # centered difference

for x in [-2, 0, 1, 2]:
    analytic = f_prime_analytical(x)
    numeric  = f_prime_numerical(x)
    print(f"x={x:4d}  analytic={analytic:8.4f}  numeric={numeric:8.4f}  err={abs(analytic-numeric):.2e}")

The chain rule tells us how to differentiate a composition of functions. If $y = f(g(x))$, then: $$\frac{dy}{dx} = \frac{dy}{du} \cdot \frac{du}{dx} \quad \text{where } u = g(x)$$

Or equivalently: $(f \circ g)'(x) = f'(g(x)) \cdot g'(x)$.

In Leibniz notation, it looks like fractions canceling — though this is a useful mnemonic, not a rigorous proof. The proof uses the limit definition and a careful squeeze argument.

Intuition. Imagine you're converting temperature: $T_{\text{Kelvin}} = f(\text{Celsius}) = C + 273$, and $C$ changes with altitude $h$ via $C = g(h)$. How fast does Kelvin change with altitude? The chain rule says: multiply the rates. If altitude rises 1 m and that drops Celsius by 0.0065°, then Kelvin drops by exactly the same 0.0065°/m because $dT/dC = 1$.

Extended chain rule (for deep compositions): $$\frac{d}{dx} f_n(f_{n-1}(\cdots f_1(x) \cdots)) = f_n'(\cdots) \cdot f_{n-1}'(\cdots) \cdots f_1'(x)$$ This is literally backpropagation in a neural network with $n$ layers. Each layer contributes a multiplicative factor to the gradient.

Why gradients vanish or explode. If each $|f_i'| < 1$, that product of $n$ terms shrinks exponentially — vanishing gradient. If each $|f_i'| > 1$, it grows exponentially — exploding gradient. The entire research area of residual networks (ResNets), gradient clipping, and careful weight initialization exists to manage this chain-rule phenomenon.

Automatic Differentiation (AD) is the algorithmic application of the chain rule. Modern ML frameworks (PyTorch, JAX, TensorFlow) maintain a computational graph of operations and apply the chain rule symbolically via either forward-mode AD or reverse-mode AD (which is exactly backpropagation). This gives exact (to floating-point precision) gradients, not numerical approximations.

# Manual chain rule: d/dx[ sin(x^2 + 1) ]
# Let u = x^2 + 1, y = sin(u)
# dy/dx = cos(u) * 2x = cos(x^2 + 1) * 2x

import math

def chain_rule_example(x):
    u = x**2 + 1
    # dy/du = cos(u), du/dx = 2x
    return math.cos(u) * 2 * x

# Verify numerically
def f(x): return math.sin(x**2 + 1)
def numeric_deriv(x, h=1e-7): return (f(x+h) - f(x-h)) / (2*h)

for x in [0.5, 1.0, 2.0]:
    analytic = chain_rule_example(x)
    numeric  = numeric_deriv(x)
    print(f"x={x}  chain-rule={analytic:.8f}  numerical={numeric:.8f}")

# --- PyTorch shows how AD does this automatically ---
# import torch
# x = torch.tensor(1.0, requires_grad=True)
# y = torch.sin(x**2 + 1)
# y.backward()  # applies chain rule via reverse-mode AD
# print(x.grad)  # = cos(x^2+1) * 2x at x=1

The derivative of a derivative is the second derivative, written $f''(x)$ or $\frac{d^2f}{dx^2}$. It measures how the rate of change is itself changing — i.e., concavity and acceleration.

Critical points and optimization. If $f'(a) = 0$, then $a$ is a critical point — the tangent is horizontal. By the First Derivative Test: if $f'$ changes from positive to negative at $a$, then $a$ is a local maximum; if from negative to positive, a local minimum.

Second Derivative Test: if $f'(a) = 0$ and $f''(a) > 0$, then $a$ is a local minimum (concave up — like a bowl). If $f''(a) < 0$, local maximum (concave down — like a hill). If $f''(a) = 0$, the test is inconclusive.

Gradient Descent. In machine learning with a loss function $L(\theta)$ over parameters $\theta$: $$\theta_{t+1} = \theta_t - \eta \cdot L'(\theta_t)$$ We move opposite to the gradient because we want to descend. $\eta$ is the learning rate. The second derivative (or the Hessian in multiple dimensions) tells us the curvature — too large a step on a highly curved region will overshoot. Adaptive optimizers like Adam estimate the second moment of gradients to adjust step size per-parameter.

Taylor's Theorem — the crown jewel. Any smooth function can be approximated near $a$ by a polynomial: $$f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \cdots + \frac{f^{(n)}(a)}{n!}(x-a)^n + R_n$$ The remainder $R_n$ is bounded by $(x-a)^{n+1}$ times some constant. This is how computers evaluate $\sin$, $\cos$, $e^x$ — via truncated Taylor/Maclaurin series. It is also why Newton's method for root-finding converges quadratically: it fits a linear (tangent line) approximation at each step, using $f$ and $f'$.

# Newton's Method: finds roots of f using f and f'
# Quadratic convergence: error roughly squares each step

def newton(f, f_prime, x0, tol=1e-12, max_iter=50):
    x = x0
    for i in range(max_iter):
        fx = f(x)
        if abs(fx) < tol:
            print(f"Converged in {i} iterations")
            return x
        x = x - fx / f_prime(x)  # x_{n+1} = x_n - f/f'
    return x

import math
root = newton(
    f       = lambda x: math.cos(x) - x,    # find x = cos(x)
    f_prime = lambda x: -math.sin(x) - 1,
    x0      = 0.5
)
print(f"Fixed point of cos: {root:.15f}")  # Dottie number ≈ 0.739085

The definite integral $\int_a^b f(x),dx$ is the signed area between the curve $y = f(x)$ and the $x$-axis, from $a$ to $b$. The definition is constructive:

1. Divide $[a,b]$ into $n$ subintervals of width $\Delta x = (b-a)/n$.
2. Pick a sample point $x_i^*$ in each subinterval.
3. Form the Riemann sum: $\sum_{i=1}^n f(x_i^*) \Delta x$.
4. The integral is the limit of this sum as $n \to \infty$ (and $\Delta x \to 0$).

This limit exists (and is the same regardless of where we pick $x_i^*$) whenever $f$ is Riemann integrable — which includes all continuous functions and many discontinuous ones.

Interpretations beyond area:

• If $v(t)$ is velocity, then $\int_a^b v(t),dt$ is net displacement.
• If $p(x)$ is a probability density function, then $\int_a^b p(x),dx$ is the probability the random variable falls in $[a,b]$.
• If $r(t)$ is a data rate (bytes/second), then $\int_0^T r(t),dt$ is total bytes transferred.

Numerical integration. On a computer, you approximate integrals using quadrature rules:

• Left/right Riemann sums: error $O(h)$
• Midpoint rule: error $O(h^2)$
• Trapezoid rule: $\int_a^b f(x),dx \approx \frac{h}{2}[f(a) + 2f(a+h) + \cdots + f(b)]$, error $O(h^2)$
• Simpson's rule: fits parabolas, error $O(h^4)$
• Gaussian quadrature: chooses sample points and weights optimally, spectacularly accurate for smooth functions

The antiderivative (indefinite integral). $\int f(x),dx = F(x) + C$ means $F'(x) = f(x)$. Finding antiderivatives is in general harder than differentiation — there is no product rule analogue (there are integration by parts and substitution instead), and many elementary functions have no closed-form antiderivative (e.g., $e^{-x^2}$, which defines the error function $\text{erf}$).

# Trapezoid rule integration: estimate int_0^pi sin(x) dx = 2
import math

def trapezoid(f, a, b, n):
    h = (b - a) / n
    total = f(a) + f(b)
    for i in range(1, n):
        total += 2 * f(a + i * h)
    return total * h / 2

for n in [10, 100, 1000, 10000]:
    approx = trapezoid(math.sin, 0, math.pi, n)
    error  = abs(approx - 2.0)
    print(f"n={n:6d}  approx={approx:.10f}  error={error:.2e}")
# Error decreases as O(h^2) = O(1/n^2): 10x more points -> 100x less error

The Fundamental Theorem of Calculus (FTC) is one of the greatest theorems in all of mathematics. It reveals that differentiation and integration are inverse operations — two seemingly different geometric ideas (slope vs. area) are in fact two sides of the same coin.

FTC Part 1 (The derivative of an integral). Define the accumulation function $A(x) = \int_a^x f(t),dt$. Then $A'(x) = f(x)$.

In words: the rate at which area accumulates is exactly the height of the function. If you pile up a little more area from $x$ to $x + \Delta x$, it's approximately $f(x) \cdot \Delta x$, and dividing by $\Delta x$ and taking the limit gives exactly $f(x)$.

FTC Part 2 (Evaluation theorem). If $F$ is any antiderivative of $f$ (meaning $F' = f$), then: $$\int_a^b f(x),dx = F(b) - F(a)$$ This is the computational engine: to evaluate a definite integral exactly, find the antiderivative and evaluate it at the two endpoints. This reduces integration from an infinite limiting process (Riemann sums) to algebraic evaluation.

Example. $\int_0^1 x^2 ,dx$. The antiderivative of $x^2$ is $\frac{x^3}{3}$. So the integral is $\frac{1^3}{3} - \frac{0^3}{3} = \frac{1}{3}$.

Why this is miraculous. These are very different operations. Integration counts area via limits of sums; differentiation computes local slope via a limit of a ratio. That they are exact inverses is not at all obvious geometrically, and took millennia of mathematical development to nail down rigorously (the rigorous proof had to wait for Cauchy and Riemann in the 1800s). The theorem also provides the proof that antiderivatives are unique up to a constant.

CS implication. The FTC justifies analytical formulas for integrals, which appear everywhere: the normalization constant of a Gaussian distribution, the cumulative distribution function, the expected value of a continuous random variable, the energy in a signal (Parseval's theorem via Fourier analysis), and more. Numerical integration is needed when there's no closed-form antiderivative.

# FTC Part 1: accumulation function and its derivative
import math

def accumulate(f, a, x, n=10000):
    """Numerically compute A(x) = integral_a^x f(t) dt"""
    if x == a: return 0.0
    h = (x - a) / n
    total = sum(f(a + (i + 0.5) * h) for i in range(n))
    return total * h  # midpoint Riemann sum

f = math.sin  # f(x) = sin(x), antiderivative is -cos(x)
a = 0.0

# FTC Part 2: exact answer via antiderivative F(x) = -cos(x)
def F(x): return -math.cos(x)

for b in [0.5, 1.0, 1.5, math.pi]:
    numerical  = accumulate(f, a, b)
    analytical = F(b) - F(a)  # FTC Part 2
    print(f"int_0^{b:.4f} sin(x)dx: numerical={numerical:.8f}  FTC={analytical:.8f}")

The FTC reduces integration to finding antiderivatives. But how do you find antiderivatives? There is no single algorithm (unlike differentiation), but there are powerful techniques.

U-Substitution is the chain rule in reverse. If the integral has the form $\int f(g(x)) g'(x),dx$, substitute $u = g(x)$, $du = g'(x),dx$: $$\int f(g(x)) g'(x),dx = \int f(u),du$$

Example: $\int 2x \cos(x^2),dx$. Let $u = x^2$, $du = 2x,dx$. Becomes $\int \cos(u),du = \sin(u) + C = \sin(x^2) + C$.

Integration by Parts is the product rule in reverse: $$\int u,dv = uv - \int v,du$$ Choose $u$ and $dv$ such that $\int v,du$ is easier. The mnemonic LIATE helps choose $u$ (Logarithm, Inverse trig, Algebraic, Trigonometric, Exponential — pick the first applicable type).

Example: $\int x e^x,dx$. Let $u = x$ (algebraic, so it will simplify: $du = dx$), $dv = e^x dx$ ($v = e^x$). Then $\int xe^x,dx = xe^x - \int e^x,dx = xe^x - e^x + C = e^x(x-1) + C$.

Important integrals to know for CS/ML:

• $\int e^{-x^2},dx$ has no elementary closed form — it defines $\text{erf}(x)$, crucial for the Gaussian distribution
• $\int_0^\infty x^{n-1} e^{-x},dx = \Gamma(n)$ — the Gamma function, the continuous extension of factorial
• $\int_0^1 x^{a-1}(1-x)^{b-1},dx = B(a,b)$ — the Beta function, the normalization of the Beta distribution (conjugate prior for Bernoulli trials)

Symbolic integration in software (SymPy, Mathematica, WolframAlpha) automates these techniques. Knowing them helps you verify, recognize patterns, and understand what the computer is doing.

# SymPy: symbolic integration demonstrating substitution and IBP
from sympy import *

x = symbols('x')

# U-substitution: int 2x*cos(x^2) dx
expr1 = 2*x * cos(x**2)
print("int 2x*cos(x^2) dx =", integrate(expr1, x))

# Integration by parts: int x*exp(x) dx
expr2 = x * exp(x)
print("int x*exp(x) dx   =", integrate(expr2, x))

# Definite integral: Gaussian normalization int_{-inf}^{inf} exp(-x^2) dx
expr3 = exp(-x**2)
result = integrate(expr3, (x, -oo, oo))
print("int_{-inf}^{inf} exp(-x^2) dx =", result)  # sqrt(pi)

This concept ties everything together in the context most CS students will use calculus most intensively: training neural networks.

Loss function and gradient. Given a neural network $\hat{y} = f(x; \theta)$ parameterized by weights $\theta$, and a loss $L(\hat{y}, y)$ measuring error, we want to minimize $L$ over $\theta$. The gradient $\nabla_\theta L$ is the vector of all partial derivatives — the generalization of $f'$ to many variables. It points in the direction of steepest increase; we step in the opposite direction.

Forward pass and reverse-mode AD. In a neural network with layers $z_1 = f_1(x)$, $z_2 = f_2(z_1)$, ..., $\hat{y} = f_n(z_{n-1})$, computing $\nabla_\theta L$ requires the chain rule across all layers: $$\frac{\partial L}{\partial \theta_k} = \frac{\partial L}{\partial \hat{y}} \cdot \frac{\partial \hat{y}}{\partial z_{n-1}} \cdots \frac{\partial z_{k+1}}{\partial z_k} \cdot \frac{\partial z_k}{\partial \theta_k}$$ Backpropagation is the dynamic-programming algorithm that computes all these partial derivatives in one backward pass by reusing intermediate computations. It has $O(n)$ complexity instead of $O(n^2)$ from naive re-computation.

Deriving the sigmoid derivative (a must-know exercise): $$\sigma(x) = \frac{1}{1+e^{-x}}$$ $$\sigma'(x) = \frac{e^{-x}}{(1+e^{-x})^2} = \sigma(x)(1 - \sigma(x))$$ This elegant form means you can compute the derivative from the output alone — no need to store the input during forward pass.

ReLU's derivative is the Heaviside step function: $\text{ReLU}'(x) = 1$ if $x > 0$, $0$ if $x < 0$ (undefined at 0, handled by convention). The clean derivative is why ReLU almost entirely replaced sigmoid for hidden-layer activations.

Integrals in ML. Expected values, KL divergence, entropy — all integrals (or sums for discrete distributions). The reparameterization trick in VAEs converts a non-differentiable sampling step into a differentiable one by pushing randomness into a parameter-free noise variable, preserving the gradient flow.

import math

def sigmoid(x):
    return 1 / (1 + math.exp(-x))

def sigmoid_deriv(x):
    s = sigmoid(x)
    return s * (1 - s)  # chain rule + quotient rule combined

def relu(x):
    return max(0.0, x)

def relu_deriv(x):
    return 1.0 if x > 0 else 0.0

# Mini backprop through 2-layer net: y = sigmoid(w2 * relu(w1 * x))
# dL/dw1 via chain rule (L = (y - y_true)^2 / 2)
def forward_backward(x, w1, w2, y_true):
    # Forward
    z1 = w1 * x
    a1 = relu(z1)
    z2 = w2 * a1
    y_hat = sigmoid(z2)
    loss = 0.5 * (y_hat - y_true) ** 2
    # Backward (chain rule)
    dL_dy  = y_hat - y_true
    dL_dz2 = dL_dy * sigmoid_deriv(z2)
    dL_dw2 = dL_dz2 * a1
    dL_da1 = dL_dz2 * w2
    dL_dz1 = dL_da1 * relu_deriv(z1)
    dL_dw1 = dL_dz1 * x
    return loss, dL_dw1, dL_dw2

loss, g1, g2 = forward_backward(x=1.0, w1=0.5, w2=0.8, y_true=1.0)
print(f"Loss={loss:.4f}  dL/dw1={g1:.6f}  dL/dw2={g2:.6f}")