Scientific Computing with SciPy

Authors

Louis Moresi

Andrew Valentine

NoteSummary

SciPy is a collection of scientific computing tools built on NumPy. It provides robust implementations of common scientific algorithms for interpolation, optimization, integration, signal processing, statistics, and much more.

What is SciPy?

SciPy builds on NumPy to provide a comprehensive ecosystem for scientific computing. It’s organized into submodules, each addressing specific scientific domains:

  • scipy.constants - Physical and mathematical constants
  • scipy.integrate - Integration and differential equations
  • scipy.interpolate - Data interpolation
  • scipy.optimize - Optimization and root finding
  • scipy.signal - Signal processing
  • scipy.spatial - Spatial algorithms and data structures
  • scipy.stats - Statistical functions
  • scipy.linalg - Linear algebra (extends NumPy)
  • scipy.fft - Fast Fourier transforms
  • scipy.ndimage - N-dimensional image processing

Physical Constants

The scipy.constants module provides precise physical constants:

You can search for constants:

Find and print constants relevant to Earth sciences: - Gravitational acceleration (standard) - Stefan-Boltzmann constant - Gas constant (R)

Then calculate Earth’s blackbody temperature using Stefan-Boltzmann law: F = σT⁴, where solar constant F ≈ 1361 W/m²

TipSolution 
from scipy import constants
import numpy as np

g = constants.g
sigma = constants.Stefan_Boltzmann
R = constants.R

print(f"g = {g} m/s²")
print(f"σ = {sigma} W/m²/K⁴")
print(f"R = {R} J/mol/K")

F = 1361 / 4
T_earth = (F / sigma) ** 0.25

print(f"\nEarth's blackbody temperature: {T_earth:.1f} K ({T_earth-273.15:.1f}°C)")

Interpolation

Interpolation fills in values between known data points. This is crucial for working with scattered or sparse data.

2D Interpolation

For gridded data:

Optimization

Find minima, maxima, and zeros of functions:

Finding Roots

Model isostatic balance: A continental crust of density ρc = 2700 kg/m³ and thickness h floats on mantle of density ρm = 3300 kg/m³.

Find the crustal thickness h that results in 1 km of topographic elevation using the isostatic balance equation:

ρc × h = ρm × (h - elevation)

TipSolution 
from scipy import optimize

rho_crust = 2700
rho_mantle = 3300
elevation = 1000

def isostasy(h):
    return rho_crust * h - rho_mantle * (h - elevation)

h_solution = optimize.fsolve(isostasy, 30000)[0]

print(f"Crustal thickness: {h_solution/1000:.2f} km")
print(f"Elevation: {elevation/1000:.2f} km")
print(f"Root depth: {(h_solution-elevation)/1000:.2f} km")

Integration

Numerical integration for definite integrals:

Double Integrals

Statistics

The scipy.stats module provides probability distributions and statistical tests:

Statistical Tests

The Gutenberg-Richter law says earthquake frequency follows: log₁₀(N) = a - b×M

Where N is number of earthquakes ≥ magnitude M, and typically b ≈ 1.

Generate synthetic earthquake data and fit the distribution:

Signal Processing with FFT

Fast Fourier Transform for frequency analysis:

Summary

SciPy provides essential scientific computing tools:

  • Constants - Physical and mathematical values
  • Interpolation - Fill in missing data points
  • Optimization - Find minima/maxima and roots
  • Integration - Numerical integration
  • Statistics - Distributions and hypothesis tests
  • Signal processing - FFT and filtering

SciPy builds on NumPy to provide the algorithms that power scientific Python!

WarningCoding scratch space