6.1 Numerical Root-Finding & Solution of Linear Systems

6.1 Numerical Root-Finding & Solution of Linear Systems

Numerical Methods & Analysis is topic I of the Mathematics subject in the PRC Enhanced Table of Specifications (PRBEE Resolution No. 40, s. 2024), weighted 3.75% of the whole examination and 15 of the 100 Mathematics items. It entered the REE syllabus with CHED Memorandum Order No. 88, s. 2017, which lists Numerical Methods and Analysis among the required BSEE professional courses. Every modern power-system tool — load flow, short-circuit study, transient stability, harmonic analysis — is a numerical method underneath, so this topic is applied as well as examinable.

---

1. Error, Precision & Termination

Absolute error: $E_a = |x_{\text{true}} - x_{\text{approx}}|$

Relative error: $E_r = \left|\dfrac{x_{\text{true}} - x_{\text{approx}}}{x_{\text{true}}}\right| \times 100%$

Approximate relative error between successive iterates, which is what an algorithm can actually compute:

$\varepsilon_a = \left|\frac{x_{k+1} - x_{k}}{x_{k+1}}\right| \times 100\%$

A result is said to be correct to $n$ significant figures when $\varepsilon_a < (0.5 \times 10^{2-n})%$.

Two error sources are always in play: truncation error (from cutting an infinite process short) and round-off error (from finite machine precision). Refining a step size reduces truncation error but eventually amplifies round-off, so there is an optimum step, not a monotonic improvement.

---

2. Bracketing Methods

Bisection Method

Requires a sign change: $f(a)f(b) < 0$. Each pass replaces the bracket with its better half at $c = (a+b)/2$.

$\text{Error after } n \text{ iterations} \le \frac{b-a}{2^{n}} \qquad \Longrightarrow \qquad n \ge \frac{\ln\!\left(\frac{b-a}{\varepsilon}\right)}{\ln 2} = \log_{2}\!\frac{b-a}{\varepsilon}$

Bisection is guaranteed to converge but only linearly. It is the safety net when faster methods diverge.

False Position (Regula Falsi)

Same bracketing guarantee, but the new estimate is the straight-line intercept rather than the midpoint:

$c = b - \frac{f(b)(a-b)}{f(a)-f(b)}$

Usually faster than bisection, but one endpoint can stagnate for strongly curved functions.

---

3. Open Methods

Newton-Raphson Method

$x_{k+1} = x_{k} - \frac{f(x_{k})}{f'(x_{k})}$

Quadratic convergence near a simple root — the number of correct digits roughly doubles each iteration. Its failure modes are examinable:

• $f'(x_k) \to 0$ sends the tangent nearly horizontal and the next iterate to infinity;
• a poor initial guess can converge to a different root or cycle indefinitely;
• convergence degrades to linear at a repeated (multiple) root.

Secant Method

Replaces the analytic derivative with a backward difference:

$x_{k+1} = x_{k} - \frac{f(x_{k})\,(x_{k} - x_{k-1})}{f(x_{k}) - f(x_{k-1})}$

Needs two starting values and no derivative; converges at order $\approx 1.618$ (the golden ratio) — slower than Newton-Raphson but cheaper per step.

Comparison

---

4. Systems of Linear Equations

Direct Methods

Gaussian elimination reduces $[A|b]$ to upper-triangular form by row operations, then back-substitutes. Partial pivoting — swapping in the row with the largest available pivot magnitude — is what keeps the process numerically stable and is required whenever a pivot is small.

LU decomposition factors $A = LU$ once, after which any number of right-hand sides is solved by one forward and one back substitution. This is why fault studies that reuse the same network matrix factor it a single time.

Cramer's rule, $x_i = \det(A_i)/\det(A)$, is examinable for $2\times2$ and $3\times3$ systems but is computationally hopeless beyond that.

Iterative Methods

Jacobi: all unknowns are updated from the previous sweep's values.

$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}\right)$

Gauss-Seidel: each new value is used immediately within the same sweep, typically halving the iteration count.

$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)}\right)$

Both converge if $A$ is strictly diagonally dominant:

$|a_{ii}| > \sum_{j \neq i} |a_{ij}| \quad \text{for every row } i$

A power network's bus admittance matrix is naturally diagonally dominant and extremely sparse, which is exactly why Gauss-Seidel was the first practical load-flow algorithm and why Newton-Raphson load flow later displaced it on convergence speed.

---

Solved Board Exam Examples

Example 1: Newton-Raphson on a Load Equation

Find the root of $f(x) = x^{3} - 2x - 5$ near $x_0 = 2$, performing two Newton-Raphson iterations.

Solution. $f'(x) = 3x^{2} - 2$.

Iteration 1: $f(2) = 8 - 4 - 5 = -1$, $f'(2) = 12 - 2 = 10$.

$x_1 = 2 - \frac{-1}{10} = 2.100$

Iteration 2: $f(2.1) = 9.261 - 4.2 - 5 = 0.061$, $f'(2.1) = 13.23 - 2 = 11.23$.

$x_2 = 2.1 - \frac{0.061}{11.23} = \boxed{2.09457}$

The true root is 2.094551, so two iterations already deliver six-figure accuracy — the signature of quadratic convergence.

Example 2: Bisection Iteration Count

How many bisection iterations are needed to locate a root in $[1, 2]$ to within $10^{-4}$?

Solution.

$n \ge \log_{2}\frac{2-1}{10^{-4}} = \log_{2}(10\,000) = \frac{\ln 10\,000}{\ln 2} = \frac{9.2103}{0.6931} = 13.29$

Rounding up, $\boxed{14 \text{ iterations}}$ are required. Newton-Raphson would reach the same tolerance in about four.

Example 3: One Gauss-Seidel Sweep

Solve the system below with one Gauss-Seidel sweep starting from $x = y = 0$:

$4x + y = 9, \qquad x + 3y = 10$

Solution. Rearranged, $x = (9 - y)/4$ and $y = (10 - x)/3$.

$x^{(1)} = \frac{9 - 0}{4} = 2.25, \qquad y^{(1)} = \frac{10 - 2.25}{3} = 2.5833$

Note that $y^{(1)}$ used the new $x^{(1)}$, not the old zero — that immediate substitution is precisely what distinguishes Gauss-Seidel from Jacobi (Jacobi would have given $y^{(1)} = 10/3 = 3.3333$). The exact solution is $x = 1.5454$, $y = 2.8182$; each row satisfies $|a_{ii}| > \sum_{j \ne i}|a_{ij}|$, so convergence is assured.