Computational Mathematics and Scientific Computing

Most equations that describe the real world, from airflow over a wing to the spread of heat in a solar panel, cannot be solved with a neat formula. Computational mathematics is the art of solving them approximately but reliably, using algorithms that a computer can run millions of times per second. In this article we explore the central ideas of the field, and work through calculations that show both the power and the pitfalls of numerical methods.

What Is Scientific Computing?

Scientific computing combines three ingredients:

1. A mathematical model, usually a differential equation, an optimisation problem or a large linear system.
2. A numerical method that replaces exact operations (limits, derivatives, integrals) with finite arithmetic.
3. An implementation on real hardware, from a laptop running Python to a supercomputer with millions of cores.

The core questions are always the same: How accurate is the answer? How fast does the error shrink as we work harder? Will small errors grow and ruin the result? The answers rely on calculus, linear algebra and careful error analysis.

Root Finding with Newton's Method

To solve \( f(x) = 0 \), Newton's method replaces the curve by its tangent line at the current guess and takes the tangent's root as the next guess:

$$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}. $$

When the starting guess is close enough to a simple root, Newton's method converges quadratically: the number of correct digits roughly doubles each step.

Worked Example 1: Computing \( \sqrt{2} \)

Solve \( f(x) = x^2 - 2 = 0 \), so \( f'(x) = 2x \). The iteration simplifies to

$$ x_{n+1} = x_n - \frac{x_n^2 - 2}{2x_n} = \frac{1}{2}\left( x_n + \frac{2}{x_n} \right). $$

Start with \( x_0 = 1 \).

Step 1: \( x_1 = \tfrac{1}{2}(1 + 2) = 1.5 \).

Step 2: \( x_2 = \tfrac{1}{2}\left(1.5 + \tfrac{2}{1.5}\right) = \tfrac{1}{2}(1.5 + 1.333333) = 1.416667 \).

Step 3: \( x_3 = \tfrac{1}{2}\left(1.416667 + \tfrac{2}{1.416667}\right) = \tfrac{1}{2}(1.416667 + 1.411765) = 1.414216 \).

The true value is \( \sqrt{2} = 1.41421356\ldots \) The errors are about \( 0.086 \), \( 0.0025 \) and \( 0.0000021 \): each error is roughly proportional to the square of the previous one, exactly as quadratic convergence predicts. Interestingly, this same formula was used by Babylonian scribes nearly 4,000 years ago, long before calculus.

Solving Differential Equations: Euler's Method

Consider an initial value problem \( y' = f(t, y) \), \( y(t_0) = y_0 \). Euler's method steps forward with step size \( h \) using the slope at the current point:

$$ y_{n+1} = y_n + h\, f(t_n, y_n). $$

Its global error is proportional to \( h \), so it is called a first-order method. Halving the step size roughly halves the error.

Worked Example 2: Radioactive-Style Decay

Solve \( y' = -2y \), \( y(0) = 1 \), up to \( t = 0.3 \) with \( h = 0.1 \). The exact solution is \( y(t) = e^{-2t} \).

Here \( f(t, y) = -2y \), so each step gives \( y_{n+1} = y_n + 0.1(-2y_n) = 0.8\, y_n \).

Step 1: \( y_1 = 0.8 \times 1 = 0.8 \) at \( t = 0.1 \).

Step 2: \( y_2 = 0.8 \times 0.8 = 0.64 \) at \( t = 0.2 \).

Step 3: \( y_3 = 0.8 \times 0.64 = 0.512 \) at \( t = 0.3 \).

Compare with the exact value: \( e^{-0.6} \approx 0.548812 \). The error is \( 0.548812 - 0.512 = 0.036812 \), about 6.7%. Repeating the calculation with \( h = 0.05 \) (six steps with factor 0.9) gives \( 0.9^6 \approx 0.531441 \), with error about 0.0174, roughly half, confirming first-order behaviour. In practice, scientists use higher-order methods such as the classical fourth-order Runge–Kutta method, whose error shrinks like \( h^4 \).

Stability: When Small Steps Matter

For the test equation \( y' = \lambda y \) with \( \lambda < 0 \), Euler's method gives \( y_{n+1} = (1 + h\lambda) y_n \). The numerical solution decays only if \( |1 + h\lambda| < 1 \). With \( \lambda = -2 \), that requires \( h < 1 \).

Try \( h = 1.5 \): the factor becomes \( 1 + 1.5(-2) = -2 \), so the values go \( 1, -2, 4, -8, \ldots \) The true solution decays smoothly to zero, but the computed one explodes and oscillates. This is numerical instability, and it explains why "stiff" problems in chemistry and circuit simulation need implicit methods, such as backward Euler, that remain stable for any step size.

Floating-Point Arithmetic and Rounding Error

Computers store real numbers in binary floating point, following the IEEE 754 standard. A double-precision number carries about 16 significant decimal digits, and many simple decimals cannot be represented exactly. For example, 0.1 in binary is a repeating fraction, so in most programming languages

$$ 0.1 + 0.2 = 0.30000000000000004. $$

The machine epsilon for double precision is \( 2^{-52} \approx 2.2 \times 10^{-16} \), the gap between 1 and the next representable number. Rounding errors are tiny individually, but they can be amplified by:

• Catastrophic cancellation: subtracting two nearly equal numbers destroys significant digits.
• Ill-conditioned problems: where the answer is extremely sensitive to the input.

Conditioning of Linear Systems

Huge linear systems \( A\mathbf{x} = \mathbf{b} \) appear in almost every simulation. The condition number \( \kappa(A) = \|A\|\,\|A^{-1}\| \) measures sensitivity: roughly, you can lose up to \( \log_{10} \kappa(A) \) digits of accuracy. For instance, the matrix

$$ A = \begin{pmatrix} 1 & 1 \\ 1 & 1.0001 \end{pmatrix} $$

has determinant \( 0.0001 \) and a condition number of roughly 40,000, so about four or five of your sixteen digits are at risk. Good algorithms (such as Gaussian elimination with partial pivoting, QR factorisation and iterative methods like conjugate gradients) are designed to avoid making things worse than the problem's inherent conditioning.

Research Frontiers

• Exascale computing: algorithms that scale to more than \( 10^{18} \) operations per second while minimising data movement.
• Mixed precision: using fast low-precision arithmetic on GPUs while correcting to full accuracy.
• Uncertainty quantification: propagating uncertainty in inputs through simulations.
• Scientific machine learning: physics-informed neural networks and learned surrogates that accelerate expensive solvers.
• Randomised numerical linear algebra: using random sampling to approximate huge matrices quickly.

Did You Know?

• In 1996, the Ariane 5 rocket was destroyed shortly after launch partly because a 64-bit floating-point number was converted to a 16-bit integer that overflowed.
• The finite element method, now used to design everything from bridges to bones, grew out of structural engineering work in the 1940s and 1950s.
• Weather forecasts are among the largest routine scientific computations on Earth, solving fluid equations on global grids every few hours.

Key Takeaways

• Scientific computing replaces exact mathematics with algorithms whose errors are carefully controlled.
• Newton's method converges quadratically near a root; Euler's method converges only linearly in \( h \).
• Stability matters: a method can produce nonsense if the step size is too large.
• Floating-point rounding and ill-conditioning limit accuracy, so algorithm design must respect them.

Need step-by-step help with numerical methods, MATLAB or Python assignments, or a computational thesis? Amathify's experts and PhD scholars are available via WhatsApp or the Ask a Question form.