Introduction to SciPy, Constants & Precision
Hey, it’s Anand here. You already know NumPy gives you fast arrays and basic vector math. SciPy takes it a step further: it is a production-grade scientific toolbox built directly on top of NumPy, leveraging compiled C and Fortran libraries (LAPACK, BLAS, QUADPACK, MINPACK, ODEPACK) under the hood.
In numerical computing, getting an answer is easy; knowing whether the answer is correct is what separates engineering from guesswork.
SciPy assumed version: 1.14+ / 1.18+. In modern SciPy, always import submodules explicitly (e.g. from scipy import constants, optimize, stats). Never rely on importing the top-level scipy alone.
1. What Even Is SciPy?
SciPy (pronounced “Sigh Pie”) provides optimized routines across science and engineering:
| Submodule | Purpose & Key Algorithms |
|---|---|
scipy.constants | Physical constants (CODATA 2018/2022) & metric conversion factors |
scipy.linalg | Dense linear algebra (LU, Cholesky, QR, SVD, Schur, matrix exponentials) |
scipy.sparse | Sparse matrix formats (CSR, CSC, DIA) and iterative solvers (cg, gmres, splu) |
scipy.optimize | Root-finding (brentq, root), unconstrained & constrained minimization (SLSQP), curve fitting |
scipy.integrate | Numerical quadrature (quad, simpson) and initial/boundary value ODE solvers (solve_ivp, solve_bvp) |
scipy.interpolate | 1D/2D splines (CubicSpline, PchipInterpolator, RegularGridInterpolator) |
scipy.stats | Probability distributions, hypothesis testing (Welch t-test, Mann-Whitney), bootstrap resampling |
scipy.signal | Filter design (Butterworth, Elliptic, FIR), Welch PSD, spectrograms, peak detection |
scipy.fft | Fast Fourier Transforms (fft, rfft, frequencies, zero-padding) |
scipy.spatial | KD-Trees (cKDTree), Convex Hulls, Voronoi diagrams, 3D Rotations |
scipy.special | Bessel, Gamma, error functions, elliptic integrals, log-sum-exp |
2. Numbers on a Computer: Tolerances & Cancellation
Computers store real numbers using IEEE 754 double precision floats: 1 sign bit, 11 exponent bits, and 52 fraction bits (53 significant bits including the hidden bit).
Two fundamental consequences follow:
- Most decimal fractions (like
0.1) cannot be stored exactly in binary. - The gap between 1.0 and the next representable float is machine epsilon ().
Never compare floating-point numbers with ==. Always check closeness with both relative and absolute tolerance:
|a - b| <= atol + rtol * |b|
Inspecting Machine Epsilon & Catastrophic Cancellation
When subtracting two nearly equal numbers, their leading bits cancel, amplifying tiny rounding errors into massive relative errors. For example:
f(x) = (1 - cos(x)) / x^2
As x -> 0, the analytical limit is 0.5. But computing 1 - cos(x) naively for x = 1e-8 cancels out to zero! Using the trigonometric identity 2 * sin^2(x/2) / x^2 eliminates subtraction entirely.
Try running the precision check below:
Example
Output:
3. Physical Constants & Unit Conversion (scipy.constants)
scipy.constants stores the internationally adopted CODATA physical constants. Since the 2019 SI redefinition, fundamental constants like the speed of light, Planck’s constant, and Boltzmann’s constant are exact integers or exact decimal values by definition.
Working with CODATA Values & Conversions
Physical Constants
Access standard constants directly:
- Speed of light:
constants.speed_of_light() - Planck’s constant:
constants.Planck() - Gravitational constant:
constants.gravitational_constant() - Gas constant:
constants.R()
Test constant lookup and conversion interactively:
Example
Output:
4. Try It Yourself
Problem: Stable Log & Exponential Evaluation
When evaluating exp(x) - 1 or log(1 + x) for small values (e.g. x = 1e-12), naive computation loses precision. Test np.expm1(x) and np.log1p(x) to preserve 16 digits of precision.
Show solution
import numpy as np
x = 1e-12
naive_exp = (np.exp(x) - 1.0) / x
stable_exp = np.expm1(x) / x
print(f"Naive (exp(x)-1)/x: {naive_exp:.12f}")
print(f"Stable expm1(x)/x : {stable_exp:.12f}")