Exercises — SciPy — overview of submodules

A verb hides in each task; match verb → drawer.

1. solve a linear system → ==scipy.linalg== (see Linear algebra — solving Ax=b).
2. minimize → ==scipy.optimize== (see Optimization — gradient descent).
3. integrate → ==scipy.integrate== (see Numerical integration (trapezoid & Simpson)).
4. interpolate (curve through points) → ==scipy.interpolate==.
5. distribution / statistics → ==scipy.stats== (see Probability distributions).

It returns the tuple (value, abserr) — the estimate of the integral and an absolute error bound. Two numbers because a numerical answer is worthless without knowing how much to trust it. Always unpack both: val, err = quad(...).

By hand. From the first equation $y=5-2x$. Substitute into the second: $x+3(5-2x)=10\Rightarrow x+15-6x=10\Rightarrow -5x=-5\Rightarrow x=1.$ Then $y=5-2(1)=3$. So $(x,y)=(1,3)$.

from scipy import linalg
import numpy as np
A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([5, 10], dtype=float)
print(linalg.solve(A, b))   # [1. 3.]

Why solve and not inv(A) @ b? solve uses LU decomposition — it factors $A$ once and back-substitutes, which is faster and numerically steadier than forming $A^{-1}$ explicitly.

By hand. ${\int }_0^{\pi }\sin xdx=[-\cos x{]}_0^{\pi }=-\cos \pi +\cos 0=-(-1)+1=2.$

from scipy import integrate
import numpy as np
val, err = integrate.quad(np.sin, 0, np.pi)
print(val, err)   # 2.0  ~2e-14

The estimate is $2.0$ and the error bound is around $10^{-14}$ — essentially exact.

By hand. The minimum of a squared term is where the inside is zero, so $x=3$. The added $+4$ shifts the value up but not the location: $f(3)=0+4=4$. Check with calculus: $\frac{d}{dx}[(x-3{)}^2+4]=2(x-3)=0\Rightarrow x=3$.

from scipy import optimize
res = optimize.minimize(lambda x: (x-3)**2 + 4, x0=0.0)
print(res.x, res.fun)   # [3.0000...] [4.0000...]

A float64 takes 8 bytes. A dense array stores every cell: $10^6\times 10^6=10^{12}\text{cells}\times 8\text{bytes}=8\times 10^{12}\text{bytes}\approx 8\text{terabytes}.$ That will never fit in RAM. Because it is mostly zeros, use `scipy.sparse`, > which stores only the nonzero entries plus their positions. Roughly $5\times 10^6$ nonzeros need on the order of tens of megabytes — a factor of about $10^5$ smaller.

No — agreeing on one tiny problem proves nothing about the general case. scipy.linalg is always compiled against BLAS/LAPACK and offers more routines (lu, qr, schur, expm, …) that numpy.linalg simply lacks. For serious numerical work — decompositions, matrix functions, large systems — prefer scipy.linalg. See Linear algebra — solving Ax=b.

For a curve that bends upward, the straight chord joining the two endpoints lies above the curve (red chord above the black curve in the figure). The trapezoid's flat top therefore covers more area than the region under the curve, so a single strip over-estimates. This is exactly the wiggliness integrate.quad fixes by refining adaptively — it adds more, thinner strips (or higher-order pieces) where the curve bends, driving the error down. See Numerical integration (trapezoid & Simpson).

Step 1 — is it a density? A density must integrate to $1$: ${\int }_0^{1}3x^{2}dx=[x^3{]}_0^1=1.✓$ Step 2 — the median. The running area from $0$ to $m$ is ${\int }_0^{m}3x^{2}dx=m^3.$ Set it to $\frac{1}{2}$: $m^3=\frac{1}{2}\Rightarrow m=(1/2{)}^{1/3}\approx \mathrm{0.7937.}$

from scipy import integrate, optimize
f = lambda x: 3*x**2
area, area_err = integrate.quad(f, 0, 1)
assert area_err < 1e-10          # honour L1.2: read the error bound
print(area, area_err)            # 1.0  ~1e-14  -> valid density
 
def running_area(m):
    val, err = integrate.quad(f, 0, m)   # unpack BOTH, per L1.2
    return val - 0.5                     # root here => median
m = optimize.brentq(running_area, 0, 1)
print(m)                          # 0.7937005...

Why brentq? We need a root of $g(m)={\int }_0^{m}f-\frac{1}{2}$ on a known bracket $[0,1]$ where $g$ changes sign ($g(0)=-\frac{1}{2}<0$, $g(1)=+\frac{1}{2}>0$). brentq reliably brackets and squeezes onto that sign change. Note we still unpack val, err inside running_area — the advice from L1.2 does not lapse just because quad sits in a helper. This is the optimize verb (root-finding) feeding on the integrate verb — two drawers in one pipeline. See Probability distributions.

The points lie exactly on a line: from $(0,1)$ we read $c=1$; the slope between $(0,1)$ and $(2,5)$ is $\frac{5-1}{2-0}=2$, so $m=2$. Every point satisfies $y=2x+1$.

from scipy import optimize
import numpy as np
x = np.array([0, 1, 2], dtype=float)
y = np.array([1, 3, 5], dtype=float)
popt, _ = optimize.curve_fit(lambda x, m, c: m*x + c, x, y)
print(popt)   # [2. 1.]  -> m=2, c=1

Which drawer is this? curve_fit lives in ==scipy.optimize==, not scipy.linalg — it wraps MINPACK's Levenberg–Marquardt least-squares solver, and least-squares fitting is an optimization verb ("minimize the squared residuals"). scipy.linalg also can do least squares (via linalg.lstsq, a decomposition), but the moment you phrase the task as "fit these parameters," the natural drawer is optimize. See Optimization — gradient descent.

Take $f(x)=x^4-8x^2$. Its derivative $f^{\prime }(x)=4x^3-16x=4x(x^2-4)$ is zero at $x=-2,0,+2$. The two minima are at $x=\pm 2$ (each giving $f=-16$); $x=0$ is a local maximum. A downhill search starting at x0 = -1 slides to $x=-2$; starting at x0 = +1 slides to $x=+2$. Same function, different basin, different reported minimizer.

from scipy import optimize
f = lambda x: x**4 - 8*x**2
print(optimize.minimize(f, x0=-1.0).x)   # ~[-2.]
print(optimize.minimize(f, x0= 1.0).x)   # ~[ 2.]

Lesson: local optimizers find the valley you start in. For global minima you must try multiple starts or a global method. See Optimization — gradient descent.

The red dot at $x=-1$ sits on the left downslope, so the search rolls left into the valley at $x=-2$. The red dot at $x=+1$ sits on the right downslope and rolls right into $x=+2$. The hump at $x=0$ is a local maximum — its slope is zero but the curve bends down on both sides, so any nudge escapes it. A gradient search will never rest there because it always steps in the downhill direction and $x=0$ has downhill on either side.

import numpy as np
from scipy import stats, fft
N = 1024
samples = stats.norm.rvs(loc=0, scale=1, size=N, random_state=0)  # SAMPLE verb
print(round(samples.mean(), 1))   # ~0.0 (law of large numbers)
print(round(samples.std(), 1))    # ~1.0
spectrum = fft.fft(samples)       # TRANSFORM verb
print(spectrum.shape)             # (1024,)  -> same length, now complex
print(spectrum.dtype)             # complex128

(a) With $N=1024$ draws from $\mathcal{N}(0,1)$, the sample mean is very close to $0$ and the sample standard deviation close to $1$ — that is what "sampling from" a distribution means: rvs = random variates (see Probability distributions). (b) fft.fft returns exactly one complex number per input point, so the output has shape $(1024,)$ — same length as the input, but now living in the frequency domain instead of the time/sample domain (see Fourier transform (FFT)). Why different drawers? scipy.stats.rvs answers "give me random numbers shaped like this distribution" (the sample verb); scipy.fft.fft answers "which frequencies is this signal built from?" (the transform verb). Two distinct verbs, two distinct drawers, even though both operate on the very same NumPy array.

Match the verb → pick the drawer → predict by hand → run → read the error / check the bracket / question the local optimum. That loop is the entire skill.