When symbolic mathematics reaches its limit, numerical mathematics takes over.
This series goes through SciPy and NumPy, solving what SymPy cannot: equations
with no closed-form solution, integrals with no antiderivative, derivatives from
measured data, curve fits with uncertainty, systems of differential equations
and the Fourier transform, always starting from concrete problems in physics,
chemistry and engineering.
The articles are listed below, in the order they were written.
Most equations, integrals, and differential equations in real problems have no closed-form solution, and SymPy can’t solve them. This 13-article series shows how NumPy and SciPy tackle these problems: roots, polynomials, numerical calculus, optimization, and FFT.
In Kepler’s equation, the unknown appears both inside and outside a sine, and no algebra isolates it. For cases like this, SciPy approximates the root. We will see bisection and brentq, Newton’s method and why it can fail, and nonlinear systems with fsolve and root.
Where is Halley’s Comet on any given date? The answer requires solving Kepler’s equation, which has no closed form. With SciPy’s newton and brentq, we locate the comet, make the law of areas visible, and compute the range of a projectile with air resistance.
From the pH of a weak acid to the volume of a real gas and the friction factor of a pipe: three equations that chemistry and engineering use every day and that only have a useful numerical solution. We solve all three with brentq and numpy.roots and rebuild the Moody chart.
Stored in float64, the coefficients of Wilkinson’s polynomial no longer return the roots 1 to 20: two of them become a complex pair. Before getting there, we look at roots, derivatives, integrals and division with numpy.polynomial, and at the coefficient-order trap.
Experimental data rarely falls on a curve. We will see how NumPy fits polynomials to noisy data, what SymPy cannot do, why extrapolating is dangerous, and why Runge’s phenomenon comes from equally spaced nodes, not from high degree.
Is the reactor temperature rising too fast? The question is a derivative, but all we have are numbers in a table. We will see finite differences, numpy.gradient and scipy.differentiate, and why shrinking the step does not always improve the derivative, because of rounding and noise.
The Gaussian function has no elementary antiderivative, and that is a theorem of Liouville. When the antiderivative does not exist, or all you have is a table of measurements, the way forward is numerical: quad, trapezoid and Simpson, from scipy.integrate, with the care each one demands.
The classic pendulum period formula only holds for small oscillations, and the exact equation has no closed-form solution. With SciPy’s solve_ivp, we measure the real period and solve stiff systems, Lotka-Volterra, falling with drag, and consecutive reactions.
Differentiating and setting the result to zero solves textbook optimization. With several variables, constraints, or a function without a tractable derivative, that calculation breaks down. We will see minimize, SLSQP and bounds from scipy.optimize, local minima, and global search with differential_evolution.
