Study guide

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

PRC Civil Engineer Licensure Examination CELE

Track your progress as you read

Finished sections get a check mark.

1 Chapter 1: Examination Overview, RA 544 & Professional Practice

1.1 CELE Structure, Format, Scoring, and 2026 Subject Order 1.2 Republic Act No. 544: The Civil Engineering Law and Board Regulations 1.3 Civil Engineering Ethics, Contracts, and Professional Practice

2 Chapter 2: Engineering Mathematics, Numerical Methods & Physics

2.1 Advanced Algebra, Progressions & Word Problems 2.2 Trigonometry and Analytical Geometry 2.3 Differential and Integral Calculus 2.4 Differential Equations and Numerical Methods 2.5 Series, Hyperbolic Functions, Multivariable Calculus, and Laplace Transforms 2.6 Numerical Methods: Linear Systems, Interpolation, Curve Fitting, and ODE Solvers 2.7 Physics for Engineers: Kinematics, Newton's Laws, Work-Energy, and Momentum

3 Chapter 3: Engineering Data Analysis & Engineering Economy

3.1 Engineering Data Analysis and Statistical Applications 3.2 Parameter Estimation, Linear Regression, Correlation, and Queuing Models 3.3 Engineering Economy: Time Value of Money and Interest Formulas 3.4 Project Evaluation, Depreciation, and Benefit-Cost Analysis

4 Chapter 4: Surveying Fundamentals, Construction Layout & Route Surveying

4.1 Errors, Corrections, and Pacing/Taping 4.2 Leveling, Earthworks, and Mass Haul Diagrams 4.3 Traverse Computations, Balancing, and Area Determinations 4.4 Route Surveying: Horizontal and Vertical Curves

5 Chapter 5: Highway & Transportation Engineering

5.1 Highway Geometric Design and Sight Distances 5.2 Traffic Engineering, Flow Fundamentals, and Pavement Design 5.3 Highway Policy, Route Location, Drainage, and Earthwork Slopes 5.4 Transportation Planning, Network Analysis, and Intermodal Systems

6 Chapter 6: Construction Materials, Methods, Management & Safety

6.1 Materials for Construction: Concrete, Formwork, Masonry, Steel, and Wood Products 6.2 Construction Equipment, Formwork, and Construction Methods 6.3 Construction Occupational Safety and Health: RA 11058 and DOLE DO 13 6.4 Construction Project Management, CPM/PERT, and Resource Scheduling 6.5 Construction Principles: Design Requirements, Project Organization, Cost Control, and Value Engineering 6.6 Quantity Surveying, Unit Cost Analysis, and Project Cost Estimates 6.7 Construction Trades and Specialties: Plumbing, Plastics, Doors and Windows, Finishes, and Electrical

7 Chapter 7: Fluid Properties, Hydrostatics & Buoyancy

7.1 Fluid Properties, Units, and Pressure Measurement 7.2 Hydrostatic Pressure on Plane and Curved Surfaces 7.3 Buoyancy, Flotation, and Accelerated Fluid Masses

8 Chapter 8: Hydrodynamics, Pipe Flow, Open Channels & Water Supply

8.1 Fluid Dynamics: Continuity and Energy Equations 8.2 Flow Measurement: Orifices, Weirs, and Venturi Meters 8.3 Fluid Flow in Pipes, Friction Losses, and Pipe Networks 8.4 Open Channel Flow, Manning's Equation, and Hydraulic Jump 8.5 Water Supply: Demand, Storage, Pumping, Distribution, and Treatment

9 Chapter 9: Soil Properties, Seepage & Stresses in Soil Mass

9.1 Soil Composition, Weight-Volume Phase Relationships 9.2 Soil Classification Systems (USCS & AASHTO) and Atterberg Limits 9.3 Soil Permeability, Seepage, Flow Nets, and Effective Stress 9.4 Stresses in Soil Mass, Stress Distribution, and Immediate and Secondary Settlement

10 Chapter 10: Consolidation, Soil Strength, Earth Pressures & Foundations

10.1 Soil Compaction and Soil Compressibility/Consolidation 10.2 Shear Strength of Soils and Mohr-Coulomb Failure Criteria 10.3 Lateral Earth Pressures (Rankine & Coulomb) and Retaining Walls 10.4 Shallow and Deep Foundation Bearing Capacity and Slope Stability

11 Chapter 11: Engineering Mechanics & Statics

11.1 Force Systems, Resultants, and Equilibrium of Rigid Bodies 11.2 Centroids, Moments of Inertia, and Friction 11.3 Cables, Trusses, and Internal Forces in Members

12 Chapter 12: Strength of Materials

12.1 Axial Stress, Strain, and Deformation 12.2 Torsion, Flanged Bolt Couplings, and Thin-Walled Pressure Vessels 12.3 Beam Stresses: Flexural and Transverse Shear 12.4 Beam Deflections, Combined Stresses, and Mohr's Circle

13 Chapter 13: Structural Analysis & NSCP Load Criteria

13.1 Indeterminate Structures: Approximate and Classical Methods 13.2 Moment Distribution Method and Slope Deflection 13.3 Influence Lines for Determinate and Indeterminate Structures 13.4 NSCP Minimum Design Loads: Dead, Live, Wind, and Earthquake

14 Chapter 14: Reinforced Concrete Design (NSCP 2015 / USD)

14.1 Ultimate Strength Design (USD) Principles and Singly Reinforced Beams 14.2 Doubly Reinforced Beams and T-Beams 14.3 Shear and Diagonal Tension, Stirrup Design, and Development Length 14.4 Reinforced Concrete Columns and Isolated Footing Design

15 Chapter 15: Steel, Timber & Prestressed Design and Materials Testing

15.1 Structural Steel Tension Members and Connections 15.2 Structural Steel Compression Members and Beam Design 15.3 Timber Design Principles, Working Stress, and Sawn Lumber Beams 15.4 Prestressed Concrete Fundamentals and Prestress Losses 15.5 Construction Materials Testing and Concrete Mix Proportioning

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

Key Takeaways Explain This Topic AI

Gaussian elimination reduces a linear system to upper-triangular form and then solves it by back substitution.

Gauss-Seidel iteration converges reliably when the coefficient matrix is diagonally dominant.

A Lagrange or Newton polynomial through n + 1 data points has degree at most n and reproduces every point exactly.

The central difference [f(x+h) − f(x−h)]/(2h) has error of order h², better than the order-h forward difference.

The classical fourth-order Runge-Kutta method has global error of order h⁴, while Euler's method is only first order.

Last updated: October 2026

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

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

Recall linear algebra and the physical meaning of derivatives and integrals.

Solve systems of linear equations.

Use curve fitting and interpolation.

Apply numerical differentiation and integration.

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

Operation Rule

Product A B AB A B Defined when columns of A A A = rows of B B B ; ( A B ) i j = ∑ k a i k b k j (AB)_{ij} = \sum_k a_{ik} b_{kj} ( A B ) ij ​ = ∑ k ​ a ik ​ b kj ​

Determinant (2×2) det ⁡ [ a b c d ] = a d − b c \det\begin{bmatrix} a & b \\ c & d \end{bmatrix} = ad - bc det [ a c ​ b d ​ ] = a d − b c

Inverse (2×2) A − 1 = 1 a d − b c [ d − b − c a ] A^{-1} = \dfrac{1}{ad - bc}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix} A − 1 = a d − b c 1 ​ [ d − c ​ − b a ​ ]

  • Singular matrix det ⁡ A = 0 \det A = 0 det A = 0 : no unique solution
  • Transpose ( A B ) T = B T A T (AB)^T = B^T A^T ( A B ) T = B T A T

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

Solving Simultaneous Linear Equations

Cramer's rule. For A x = b A\mathbf{x} = \mathbf{b} A x = b , x i = det ⁡ A i / det ⁡ A x_i = \det A_i / \det A x i ​ = det A i ​ / det A , where A i A_i A i ​ is A A A with column i i i replaced by b \mathbf{b} 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 2 x + y − z = 3 2x + y - z = 3 2 x + y − z = 3 , x + 3 y + 2 z = 13 x + 3y + 2z = 13 x + 3 y + 2 z = 13 , 3 x + y + 3 z = 14 3x + y + 3z = 14 3 x + y + 3 z = 14 .

Eliminate x x x : R2 − ½R1 gives 2.5 y + 2.5 z = 11.5 2.5y + 2.5z = 11.5 2.5 y + 2.5 z = 11.5 . R3 − 1.5R1 gives − 0.5 y + 4.5 z = 9.5 -0.5y + 4.5z = 9.5 − 0.5 y + 4.5 z = 9.5 .

Eliminate y y y : R3 + 0.2R2 gives 5 z = 11.8 5z = 11.8 5 z = 11.8 , so z = 2.36 z = 2.36 z = 2.36 .

Back-substitute: y = ( 11.5 − 2.5 × 2.36 ) / 2.5 = 2.24 y = (11.5 - 2.5 \times 2.36)/2.5 = 2.24 y = ( 11.5 − 2.5 × 2.36 ) /2.5 = 2.24 , and x = ( 3 − 2.24 + 2.36 ) / 2 = 1.56 x = (3 - 2.24 + 2.36)/2 = 1.56 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 3(1.56) + 2.24 + 3(2.36) = 4.68 + 2.24 + 7.08 = 14.00 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 ) , … , ( x n , y n ) (x_0, y_0), \dots, (x_n, y_n) ( x 0 ​ , y 0 ​ ) , … , ( x n ​ , y n ​ ) :

P ( x ) = ∑ i = 0 n y i ∏ j ≠ i x − x j x i − x j P(x) = \sum_{i=0}^{n} y_i \prod_{j \ne i} \frac{x - x_j}{x_i - x_j} P ( x ) = ∑ i = 0 n ​ y i ​ ∏ j  = i ​ x i ​ − x j ​ x − 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 ) + ⋯ 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 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 ​ ) + ⋯

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

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

Divided differences: f [ x 0 , x 1 ] = 3.0 f[x_0,x_1] = 3.0 f [ x 0 ​ , x 1 ​ ] = 3.0 , f [ x 1 , x 2 ] = 5.0 f[x_1,x_2] = 5.0 f [ x 1 ​ , x 2 ​ ] = 5.0 , f [ x 0 , x 1 , x 2 ] = ( 5.0 − 3.0 ) / 2 = 1.0 f[x_0,x_1,x_2] = (5.0 - 3.0)/2 = 1.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 m 3 / s 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} P ( 2.5 ) = 2.0 + 3.0 ( 1.5 ) + 1.0 ( 1.5 ) ( 0.5 ) = 2.0 + 4.5 + 0.75 = 7.25 m 3 / s .

Linear interpolation between the last two points would give 7.5 7.5 7.5 , so the curvature matters.

Caution

High-degree polynomials through many equally spaced points can oscillate wildly near the ends (Runge's phenomenon). Use piecewise or spline interpolation for long data sets.

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 R^2 R 2 , is treated in Section 3.2.

Nonlinear relationships are often linearized :

Model Transform Linear form

y = a e b x y = a e^{bx} y = a e b x take ln ⁡ y \ln y ln y ln ⁡ y = ln ⁡ a + b x \ln y = \ln a + bx ln y = ln a + b x

y = a x b y = a x^b y = a x b take log ⁡ \log lo g of both sides log ⁡ y = log ⁡ a + b log ⁡ x \log y = \log a + b\log x lo g y = lo g a + b lo g x

y = a + b / x y = a + b/x y = a + b / x let u = 1 / x u = 1/x u = 1/ x y = a + b u y = a + bu y = a + b u

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

Numerical Differentiation

Formula Expression Error order

Forward difference f ( x + h ) − f ( x ) h \dfrac{f(x+h) - f(x)}{h} h f ( x + h ) − f ( x ) ​ O ( h ) O(h) O ( h )

Backward difference f ( x ) − f ( x − h ) h \dfrac{f(x) - f(x-h)}{h} h f ( x ) − f ( x − h ) ​ O ( h ) O(h) O ( h )

Central difference f ( x + h ) − f ( x − h ) 2 h \dfrac{f(x+h) - f(x-h)}{2h} 2 h f ( x + h ) − f ( x − h ) ​ O ( h 2 ) O(h^2) O ( h 2 )

Second derivative (central) f ( x + h ) − 2 f ( x ) + f ( x − h ) h 2 \dfrac{f(x+h) - 2f(x) + f(x-h)}{h^2} h 2 f ( x + h ) − 2 f ( x ) + f ( x − h ) ​ O ( h 2 ) O(h^2) O ( h 2 )

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

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

Numerical Solution of Ordinary Differential Equations

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

Method Update Global error

Euler y n + 1 = y n + h f ( x n , y n ) y_{n+1} = y_n + h f(x_n, y_n) y n + 1 ​ = y n ​ + h f ( x n ​ , y n ​ ) O ( h ) O(h) O ( h )

Heun (improved Euler) predictor y ∗ = y n + h f ( x n , y n ) y^* = y_n + hf(x_n,y_n) y ∗ = y n ​ + h f ( x n ​ , y n ​ ) ; y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y ∗ ) ] y_{n+1} = y_n + \frac{h}{2}[f(x_n,y_n) + f(x_{n+1},y^*)] y n + 1 ​ = y n ​ + 2 h ​ [ f ( x n ​ , y n ​ ) + f ( x n + 1 ​ , y ∗ )] O ( h 2 ) O(h^2) O ( h 2 )

Runge-Kutta 4 y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) y n + 1 ​ = y n ​ + 6 h ​ ( k 1 ​ + 2 k 2 ​ + 2 k 3 ​ + k 4 ​ ) O ( h 4 ) O(h^4) O ( h 4 )

The RK4 slopes are:

k 1 = f ( x n , y n ) , k 2 = f ⁣ ( x n + h 2 , y n + h 2 k 1 ) , k 3 = f ⁣ ( x n + h 2 , y n + h 2 k 2 ) , k 4 = f ( x n + h , y n + h k 3 ) 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) k 1 ​ = f ( x n ​ , y n ​ ) , k 2 ​ = f ( x n ​ + 2 h ​ , y n ​ + 2 h ​ k 1 ​ ) , k 3 ​ = f ( x n ​ + 2 h ​ , y n ​ + 2 h ​ k 2 ​ ) , k 4 ​ = f ( x n ​ + h , y n ​ + h k 3 ​ )

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

Euler gives y ( 0.1 ) = 1 + 0.1 ( 0 + 1 ) = 1.100 y(0.1) = 1 + 0.1(0 + 1) = 1.100 y ( 0.1 ) = 1 + 0.1 ( 0 + 1 ) = 1.100 .

Heun: the predictor is 1.100 1.100 1.100 , and f ( 0.1 , 1.100 ) = 1.200 f(0.1, 1.100) = 1.200 f ( 0.1 , 1.100 ) = 1.200 . So y = 1 + 0.05 ( 1 + 1.2 ) = 1.110 y = 1 + 0.05(1 + 1.2) = 1.110 y = 1 + 0.05 ( 1 + 1.2 ) = 1.110 .

The exact solution y = 2 e x − x − 1 y = 2e^x - x - 1 y = 2 e x − x − 1 gives 1.11034 1.11034 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.

Loading diagram... Picking a Numerical Method