6.2 Interpolation, Numerical Integration & Numerical Solution of ODEs

6.2 Interpolation, Numerical Integration & Numerical Solution of ODEs

The second half of the PRC Enhanced TOS topic Numerical Methods & Analysis deals with data that arrives as a table rather than a formula, and with differential equations that have no closed-form solution. Electrical engineers meet both daily: conductor ampacity tables that need interpolation, measured load curves that need integration into energy, and transient-stability swing equations that need step-by-step solution.

---

1. Interpolation

Linear Interpolation

The workhorse for reading between two rows of a PEC or manufacturer's table:

$y = y_0 + \frac{(x - x_0)(y_1 - y_0)}{x_1 - x_0}$

Lagrange Interpolating Polynomial

For $n+1$ data points, the unique polynomial of degree $\le n$ through them is

$P_n(x) = \sum_{i=0}^{n} y_i L_i(x), \qquad L_i(x) = \prod_{\substack{j=0 \\ j \ne i}}^{n} \frac{x - x_j}{x_i - x_j}$

Newton's Divided-Difference Form

$P_n(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$

with divided differences built recursively as

$f[x_i, x_{i+1}] = \frac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}$

Examinable fact: Lagrange and Newton forms give the same polynomial through the same points — the interpolating polynomial of degree $\le n$ through $n+1$ distinct points is unique. Newton's form is preferred in practice because adding a new data point only appends one term instead of rebuilding everything.

Interpolation vs. extrapolation: estimating inside the data range is interpolation and is comparatively safe; going outside it is extrapolation, and high-degree polynomials diverge violently there (Runge's phenomenon). Spline interpolation, which fits low-degree polynomials piecewise with matched slopes and curvatures, avoids this oscillation.

---

2. Finite Differences & Numerical Differentiation

The central difference is second-order accurate at the same cost as the one-sided forms, so it is the default. Be aware that numerical differentiation amplifies measurement noise (the subtraction of nearly equal numbers is divided by a small $h$), which is why differentiating a noisy current signal is far less reliable than integrating it.

---

3. Numerical Integration

With $h = (b-a)/n$:

Trapezoidal Rule

$\int_a^b f(x)\,dx \approx \frac{h}{2}\left[f_0 + 2(f_1 + f_2 + \cdots + f_{n-1}) + f_n\right], \qquad E = O(h^{2})$

Simpson's 1/3 Rule — $n$ must be even

$\int_a^b f(x)\,dx \approx \frac{h}{3}\left[f_0 + 4(f_1 + f_3 + \cdots) + 2(f_2 + f_4 + \cdots) + f_n\right], \qquad E = O(h^{4})$

The coefficient pattern is $1, 4, 2, 4, 2, \dots, 4, 1$. Although it fits parabolas, Simpson's 1/3 rule integrates cubics exactly — a favourite examination fact.

Simpson's 3/8 Rule — $n$ must be a multiple of three

$\int_a^b f(x)\,dx \approx \frac{3h}{8}\left[f_0 + 3f_1 + 3f_2 + 2f_3 + 3f_4 + \cdots + f_n\right]$

---

4. Numerical Solution of Ordinary Differential Equations

For $\dfrac{dy}{dx} = f(x,y)$ with $y(x_0) = y_0$:

Euler's Method

$y_{k+1} = y_k + h\,f(x_k, y_k), \qquad \text{global error } O(h)$

Simple, but only first-order accurate and conditionally stable — halving $h$ only halves the error while doubling the work.

Heun's Method (Improved Euler, RK2)

$y_{k+1} = y_k + \frac{h}{2}\left[f(x_k, y_k) + f(x_{k+1}, y_k + h f(x_k,y_k))\right], \qquad O(h^{2})$

Fourth-Order Runge-Kutta (RK4)

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

$y_{k+1} = y_k + \frac{1}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right), \qquad \text{global error } O(h^{4})$

RK4 costs four function evaluations per step but delivers fourth-order accuracy, which is why it is the standard integrator for the transient-stability swing equation and for machine dynamic models.

---

Solved Board Exam Examples

Example 1: Simpson's 1/3 Rule on a Load Curve

A feeder demand is recorded every 2 hours over an 8-hour shift: 4.0, 5.2, 6.8, 6.0 and 4.4 MW. Estimate the energy delivered.

Solution. Here $h = 2$ h and $n = 4$ intervals — even, so Simpson's 1/3 rule applies:

$E \approx \frac{2}{3}\left[4.0 + 4(5.2) + 2(6.8) + 4(6.0) + 4.4\right]$

$= \frac{2}{3}\left[4.0 + 20.8 + 13.6 + 24.0 + 4.4\right] = \frac{2}{3}(66.8) = \boxed{44.53 \text{ MWh}}$

The trapezoidal rule on the same data gives 44.40 MWh — close here, but with error of order $h^2$ rather than $h^4$.

Example 2: One RK4 Step on an RL Transient

For $\dfrac{di}{dt} = 10 - 2i$ with $i(0) = 0$, take one RK4 step of $h = 0.1$ s.

Solution.

$k_1 = 0.1(10 - 0) = 1.0000$ $k_2 = 0.1\left[10 - 2(0 + 0.5)\right] = 0.1(9.0) = 0.9000$ $k_3 = 0.1\left[10 - 2(0 + 0.45)\right] = 0.1(9.1) = 0.9100$ $k_4 = 0.1\left[10 - 2(0 + 0.91)\right] = 0.1(8.18) = 0.8180$

$i(0.1) = 0 + \frac{1}{6}\left[1.0000 + 2(0.9000) + 2(0.9100) + 0.8180\right] = \frac{5.4380}{6} = \boxed{0.90633 \text{ A}}$

The exact solution is $i(t) = 5(1 - e^{-2t})$, giving $i(0.1) = 0.906346$ A — RK4 is correct to five decimal places in a single step, while plain Euler would have returned 1.0000 A, an error of 10%.

Example 3: Lagrange Interpolation of Conductor Ampacity

Ampacity is tabulated as 115 A at 30 °C, 103 A at 40 °C and 89 A at 50 °C. Estimate the ampacity at 35 °C using a second-degree Lagrange polynomial.

Solution. With $x_0 = 30, x_1 = 40, x_2 = 50$ and $x = 35$:

$L_0 = \frac{(35-40)(35-50)}{(30-40)(30-50)} = \frac{75}{200} = 0.375$ $L_1 = \frac{(35-30)(35-50)}{(40-30)(40-50)} = \frac{-75}{-100} = 0.750$ $L_2 = \frac{(35-30)(35-40)}{(50-30)(50-40)} = \frac{-25}{200} = -0.125$

$P(35) = 0.375(115) + 0.750(103) - 0.125(89) = 43.125 + 77.25 - 11.125 = \boxed{109.25 \text{ A}}$

Simple linear interpolation between the first two rows would give 109.0 A; the quadratic form captures the slight curvature of the derating curve. Note that the coefficients sum to exactly 1.000, a useful arithmetic check.