2.6 Numerical Methods: Linear Systems, Interpolation, Curve Fitting, and ODE Solvers

2.6 Numerical Methods: Linear Systems, Interpolation, Curve Fitting, and ODE Solvers

The AMSTHC TOS area "Numerical Methods for Engineers" has five competencies:

1. Recall linear algebra and the physical meaning of derivatives and integrals.
2. Solve systems of linear equations.
3. Use curve fitting and interpolation.
4. Apply numerical differentiation and integration.
5. Evaluate ordinary differential equations and boundary-value problems.

Root finding and Simpson's rules are in Section 2.4. This section covers the rest. Most items can be done on an allowed scientific calculator, many of which solve 2×2 and 3×3 systems directly. You must still know the method behind the answer.

---

Matrix Essentials

A 3×3 determinant is found by cofactor expansion or by the diagonal rule (Sarrus).

---

Solving Simultaneous Linear Equations

Cramer's rule. For $A\mathbf{x} = \mathbf{b}$, $x_i = \det A_i / \det A$, where $A_i$ is $A$ with column $i$ replaced by $\mathbf{b}$. It is practical for 2 or 3 unknowns.

Gaussian elimination. Use row operations to make the matrix upper-triangular, then back-substitute. Partial pivoting, which swaps rows to put the largest coefficient on the diagonal, reduces round-off error.

Example. Solve $2x + y - z = 3$, $x + 3y + 2z = 13$, $3x + y + 3z = 14$.

• Eliminate $x$: R2 − ½R1 gives $2.5y + 2.5z = 11.5$. R3 − 1.5R1 gives $-0.5y + 4.5z = 9.5$.
• Eliminate $y$: R3 + 0.2R2 gives $5z = 11.8$, so $z = 2.36$.
• Back-substitute: $y = (11.5 - 2.5 \times 2.36)/2.5 = 2.24$, and $x = (3 - 2.24 + 2.36)/2 = 1.56$.

Check in equation 3: $3(1.56) + 2.24 + 3(2.36) = 4.68 + 2.24 + 7.08 = 14.00$.

Iterative methods. Jacobi updates every unknown from the previous iteration's values. Gauss-Seidel uses each new value as soon as it is available, so it usually converges faster. Both converge when the matrix is diagonally dominant, meaning each diagonal coefficient exceeds the sum of the other coefficients in its row. Large structural stiffness and pipe-network systems are solved this way.

---

Interpolation

Lagrange form through $(x_0, y_0), \dots, (x_n, y_n)$:

$P(x) = \sum_{i=0}^{n} y_i \prod_{j \ne i} \frac{x - x_j}{x_i - x_j}$

Newton divided differences build the same polynomial incrementally:

$P(x) = f[x_0] + f[x_0,x_1](x - x_0) + f[x_0,x_1,x_2](x - x_0)(x - x_1) + \cdots$

The benefit is that adding a new data point only adds one more term.

Example. A rating curve gives $Q = 2.0$, $5.0$ and $10.0\text{ m}^3/\text{s}$ at stages $1.0$, $2.0$ and $3.0\text{ m}$. Estimate $Q$ at $2.5\text{ m}$.

• Divided differences: $f[x_0,x_1] = 3.0$, $f[x_1,x_2] = 5.0$, $f[x_0,x_1,x_2] = (5.0 - 3.0)/2 = 1.0$.
• $P(2.5) = 2.0 + 3.0(1.5) + 1.0(1.5)(0.5) = 2.0 + 4.5 + 0.75 = 7.25\text{ m}^3/\text{s}$.
• Linear interpolation between the last two points would give $7.5$, so the curvature matters.

---

Curve Fitting by Least Squares

When data contain scatter, fit a curve that minimizes the sum of squared vertical errors rather than passing through every point. The straight-line case, with slope and intercept formulas, correlation and $R^2$, is treated in Section 3.2.

Nonlinear relationships are often linearized:

Fit a line to the transformed data, then convert back.

---

Numerical Differentiation

Example. Water-surface elevations of $10.40$, $10.25$ and $10.06\text{ m}$ are read at $x = 0$, $100$ and $200\text{ m}$. The central-difference slope at $x = 100$ is $(10.06 - 10.40)/200 = -0.0017$.

Reducing $h$ lowers truncation error, but very small $h$ magnifies round-off error.

---

Numerical Solution of Ordinary Differential Equations

For $y' = f(x, y)$ with $y(x_0) = y_0$ and step $h$:

The RK4 slopes are:

$k_1 = f(x_n, y_n), \quad k_2 = f\!\left(x_n + \tfrac{h}{2}, y_n + \tfrac{h}{2}k_1\right), \quad k_3 = f\!\left(x_n + \tfrac{h}{2}, y_n + \tfrac{h}{2}k_2\right), \quad k_4 = f(x_n + h, y_n + hk_3)$

Example. For $y' = x + y$ with $y(0) = 1$ and $h = 0.1$:

• Euler gives $y(0.1) = 1 + 0.1(0 + 1) = 1.100$.
• Heun: the predictor is $1.100$, and $f(0.1, 1.100) = 1.200$. So $y = 1 + 0.05(1 + 1.2) = 1.110$.
• The exact solution $y = 2e^x - x - 1$ gives $1.11034$. Heun is much closer than Euler.

Boundary-value problems. These specify conditions at two ends, like a beam's deflection at both supports. There are two common approaches:

• The shooting method guesses the missing initial slope, integrates, and adjusts the guess until the far-end condition is met.
• The finite-difference method replaces derivatives with difference formulas at interior nodes, producing a system of linear equations.