2.4 Differential Equations and Numerical Methods

2.4 Differential Equations and Numerical Methods

Differential equations and numerical methods bridge pure mathematics and applied engineering mechanics in civil engineering. While first-order and second-order ordinary differential equations (ODEs) model hydraulic seepage, thermal dissipation in curing mass concrete, and structural dynamics, practical field problems often resist closed-form analytical solutions. Numerical methods—including root-finding iterations (Bisection, Newton-Raphson) and numerical integration (Trapezoidal, Simpson's 1/3 and 3/8 rules)—provide accurate numerical approximations required on the board examination.

---

First-Order Ordinary Differential Equations (ODEs)

A first-order differential equation involves independent variable $x$, dependent variable $y$, and the first derivative $\frac{dy}{dx}$.

1. Variable Separable ODEs

The equation can be manipulated algebraically into isolated single-variable differentials: $f(x)\,dx + g(y)\,dy = 0 \implies \int f(x)\,dx + \int g(y)\,dy = C$

2. Homogeneous Differential Equations

An equation in the differential form $M(x, y)\,dx + N(x, y)\,dy = 0$ is homogeneous of degree $n$ if $M(kx, ky) = k^n M(x, y)$ and $N(kx, ky) = k^n N(x, y)$.

• Standard Substitution: Let $y = vx \implies dy = v\,dx + x\,dv$ (or $x = vy \implies dx = v\,dy + y\,dv$).
• Substituting transforms the equation into a separable ODE in terms of $v$ and $x$.

3. Exact Differential Equations

The equation $M(x, y)\,dx + N(x, y)\,dy = 0$ is exact if and only if: $\frac{\partial M}{\partial y} = \frac{\partial N}{\partial x}$ The general solution is $F(x, y) = C$, found by: $F(x, y) = \int M(x, y)\,\partial x + g(y), \quad \text{where } g'(y) = N(x, y) - \frac{\partial}{\partial y}\left[\int M(x, y)\,\partial x\right]$

4. First-Order Linear Differential Equations

The standard canonical form is: $\frac{dy}{dx} + P(x)y = Q(x)$ Where $P(x)$ and $Q(x)$ are continuous functions of $x$ alone.

• Integrating Factor ($I(x)$): $I(x) = e^{\int P(x)\,dx}$
• General Closed-Form Solution: $y \cdot I(x) = \int Q(x) I(x)\,dx + C$

---

Engineering Applications of First-Order Models

1. Newton's Law of Cooling (Thermal Curing of Concrete)

Mass concrete structures (such as gravity dams and bridge piers) generate substantial internal hydration heat that must dissipate to the ambient environment. Newton's Law states that the rate of temperature change of a body is proportional to the difference between its temperature $T$ and the surrounding ambient temperature $T_m$: $\frac{dT}{dt} = -k(T - T_m)$ Separating variables and integrating yields: $T(t) = T_m + (T_0 - T_m)e^{-kt}$ Where $T_0$ is the initial temperature at $t = 0$, and $k$ is the cooling rate constant ($\text{time}^{-1}$).

2. Exponential Growth and Decay

Governed by $\frac{dN}{dt} = kN \implies N(t) = N_0 e^{kt}$. If $k < 0$, it models radioisotope decay or contaminant breakdown; half-life is $t_{1/2} = \frac{\ln 2}{|k|}$.

---

Second-Order Linear Homogeneous ODEs with Constant Coefficients

Second-order equations describe dynamic mechanical vibrations, structural frame oscillations, and beam deflections: $a \frac{d^2 y}{dx^2} + b \frac{dy}{dx} + c y = 0 \quad (a \ne 0)$

The Characteristic (Auxiliary) Equation

Assuming trial solution $y = e^{rx}$ yields the characteristic quadratic equation: $a r^2 + b r + c = 0$

Roots are obtained via $r = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$. The physical and mathematical nature of the response is classified by the discriminant $\Delta = b^2 - 4ac$:

---

Numerical Root-Finding Algorithms

When non-linear transcendental equations $f(x) = 0$ (such as pipe friction factor formulas or open channel critical depths) cannot be solved analytically, numerical iterative algorithms are employed.

1. The Bisection Method

Based on the Intermediate Value Theorem: If a continuous function $f(x)$ satisfies $f(a) \cdot f(b) < 0$, at least one real root exists within the bracket $[a, b]$.

• Iteration Formula: Midpoint $c = \frac{a + b}{2}$. If $f(a)f(c) < 0$, root lies in $[a, c]$; otherwise root lies in $[c, b]$.
• Convergence Rate: Linear (error halves each iteration, convergence factor $= 0.5$).
• Number of Iterations for Tolerance $\epsilon$: $n \ge \frac{\ln\left(\frac{b - a}{\epsilon}\right)}{\ln 2}$
• Strengths & Weaknesses: Absolutely guaranteed convergence, but slow compared to gradient methods.

2. The Newton-Raphson Method

Uses the local tangent line approximation at trial point $x_n$ to project to the next root estimate: $x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

• Convergence Rate: Quadratic near a simple root ($|e_{n+1}| \approx M |e_n|^2$), doubling the number of significant correct decimal places with each iteration.
• Common Failure Modes:
  1. Zero Derivative ($f'(x_n) = 0$): Horizontal tangent causes division by zero.
  2. Oscillatory Cycles: Diverges or enters an infinite periodic loop around points of inflection.
  3. Poor Starting Value: Converges to an unintended distant root if $x_0$ is poorly chosen.
  4. Multiple Roots: If a root has multiplicity $m > 1$, convergence slows from quadratic to linear.

---

Numerical Integration Methods

Numerical quadrature approximates the definite integral $I = \int_a^b f(x)\,dx$ using $n$ equal panels of step size $h = \frac{b - a}{n}$ across discrete ordinates $y_i = f(x_i)$.

1. The Trapezoidal Rule

Approximates the region under each subinterval with a linear trapezoid: $\int_a^b f(x)\,dx \approx \frac{h}{2}\left[f(x_0) + 2\sum_{i=1}^{n-1} f(x_i) + f(x_n)\right]$

• Polynomial Exactness: Degree $\le 1$ (exact for linear profiles).
• Global Truncation Error: $E_T = -\frac{(b - a)h^2}{12} f''(\xi) = O(h^2)$

2. Simpson's 1/3 Rule

Approximates consecutive pairs of panels with quadratic parabolas.

• Strict Constraint: Requires an even number of intervals $n$ (an odd number of ordinates $n + 1$). $\int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(x_0) + 4\sum_{\text{odd } i} f(x_i) + 2\sum_{\text{even } i} f(x_i) + f(x_n)\right]$
• Polynomial Exactness: Degree $\le 3$ (exact for cubics as well as quadratics due to error cancellation!).
• Global Truncation Error: $E_S = -\frac{(b - a)h^4}{180} f^{(4)}(\xi) = O(h^4)$

3. Simpson's 3/8 Rule

Approximates sets of three panels with cubic polynomials.

• Strict Constraint: Requires the number of intervals $n$ to be a multiple of 3. $\int_a^b f(x)\,dx \approx \frac{3h}{8}\left[f(x_0) + 3f(x_1) + 3f(x_2) + 2f(x_3) + 3f(x_4) + 3f(x_5) + 2f(x_6) + \dots + f(x_n)\right]$
• Polynomial Exactness: Degree $\le 3$.
• Global Truncation Error: $E_{3/8} = -\frac{(b - a)h^4}{80} f^{(4)}(\xi) = O(h^4)$

---

Step-by-Step Worked Problems

Worked Example: Thermal Dissipation in Curing Mass Concrete

Problem: A freshly cast mass concrete bridge pier foundation is placed at an initial core temperature of $T_0 = 65^\circ\text{C}$. The ambient surrounding air temperature is maintained at $T_m = 25^\circ\text{C}$. After $4\text{ hours}$ of cooling, the core temperature drops to $55^\circ\text{C}$. Assuming Newton's Law of Cooling governs heat transfer:

1. Determine the cooling rate constant $k$ (in $\text{hr}^{-1}$).
2. What will be the concrete core temperature after a total elapsed time of $12\text{ hours}$?

Solution:

1. Apply Newton's Law of Cooling formula: $T(t) = T_m + (T_0 - T_m)e^{-kt}$ $T(t) = 25 + (65 - 25)e^{-kt} = 25 + 40e^{-kt}$
2. Use the condition at $t = 4\text{ hours}$ ($T = 55^\circ\text{C}$): $55 = 25 + 40e^{-4k}$ $30 = 40e^{-4k} \implies e^{-4k} = \frac{30}{40} = 0.75$ $-4k = \ln(0.75) \approx -0.28768$ $k = \frac{0.28768}{4} \approx 0.07192\text{ hr}^{-1}$
3. Compute the core temperature at $t = 12\text{ hours}$: $T(12) = 25 + 40e^{-0.07192(12)} = 25 + 40e^{-0.86304}$ $e^{-0.86304} = (e^{-4k})^3 = (0.75)^3 = 0.421875$ $T(12) = 25 + 40(0.421875) = 25 + 16.875 = 41.88^\circ\text{C}$
4. The core temperature after $12\text{ hours}$ is $41.88^\circ\text{C}$.

---

CELE Board Exam Traps & Strategic Checklists