Before selecting a numerical method, the mathematical nature of the problem must be defined by how the constraints are applied. For an IVP, all auxiliary conditions are specified at a single point, allowing the numerical solution to "march" forward step-by-step. For a BVP, auxiliary conditions are specified at multiple distinct points (e.g., $y(a) = y_a$ and $y(b) = y_b$). Since you do not have all information at one starting point, the solution must simultaneously satisfy all boundary constraints.

• Local Truncation Error ($E_a$): The pure mathematical error introduced during a single isolated step of the method, assuming the value at the start of that specific step was exactly correct. For Euler's method, the local truncation error is proportional to the step size squared ($O(h^2)$).
• Global Truncation Error: The total accumulated error at the end of the entire simulation. Because the method must take approximately $1/h$ individual steps to traverse a total distance, the accumulated global truncation error is typically one full order lower than the local error. For Euler's method, it degrades to $O(h)$.

Heun's method (Improved Euler) is a classic predictor-corrector method. It first uses standard explicit Euler's method to loosely estimate (predict) a preliminary value $y_{i+1}^0$ at the end of the step. It then evaluates the slope at this predicted future point and averages it with the initial slope to correct the final estimate.

The general RK2 formula computes a weighted average of two intermediate slopes. Different algebraic choices of these weighting parameters yield specific named methods:

• Heun's Method (or Improved Euler): Uses the beginning and end slopes.
• Midpoint Method: Uses the slope at the middle of the interval.
• Ralston's Method: Evaluates the second slope at $x_i + \frac{3}{4}h$, yielding a theoretically minimized local truncation error bound for a second-order method.

A highly efficient adaptive algorithm that calculates two separate RK estimates—one of order 4 and one of order 5—using the exact same intermediate function evaluations (slopes). The absolute difference between the two computed estimates provides a continuous, real-time measure of the local truncation error:

• If the calculated error strictly exceeds a specified tolerance, the step size $h$ is decreased, and the step is completely repeated.
• If the calculated error is much smaller than the tolerance, $h$ is mathematically increased for the next step, drastically speeding up the simulation.

This dual-estimate strategy forms the computational core for standard ODE solver functions like MATLAB's ode45 and Python's scipy.integrate.solve_ivp.

• Adams-Bashforth: Explicit predictor formulas that use previous points ($y_{i}, y_{i-1}, y_{i-2}$) to fit an extrapolating polynomial to the derivative.
• Adams-Moulton: Implicit corrector formulas that are systematically paired with Adams-Bashforth predictors to create highly accurate predictor-corrector pairs. They use the predicted future point to form an interpolating polynomial.

For a linear or nonlinear $n$-th order ODE, define $n$ new state variables. For example, for the damped oscillator $y'' + 3y' + 2y = 0$:

1. Let $y_1 = y$ (Position)
2. Let $y_2 = y'$ (Velocity) This yields a coupled system of two first-order ODEs:

• $y_1' = y_2$
• $y_2' = -3y_2 - 2y_1$

Methods like Euler's or RK4 are then applied simultaneously to the entire system state vector $\{Y\} = [y_1, y_2]^T$ at each discrete time step.

When analyzing 2D autonomous systems (where time $t$ does not appear explicitly, $y_1' = f(y_1, y_2)$ and $y_2' = g(y_1, y_2)$), we can plot $y_2$ against $y_1$ directly, ignoring time. This creates a phase plane. Plotting the numerical solutions as trajectories in this plane reveals the fundamental stability characteristics of the system, such as convergence to a stable equilibrium point (a node or focus) or a continuous periodic orbit (a limit cycle). Phase plane analysis is critical for understanding nonlinear dynamics without requiring exact analytical time-domain solutions.

• Stability Region: The set of complex numbers $h\lambda$ in the complex plane for which the numerical method produces a bounded, non-diverging solution (i.e., $|y_{i+1}/y_i| \le 1$).
• A-stability: A numerical method is A-stable if its stability region encompasses the entire left half of the complex plane ($Re(h\lambda) \le 0$). This mathematically guarantees the numerical solution decays whenever the exact analytical solution decays, regardless of how large the step size $h$ is chosen.

Explicit methods (like standard RK4) have finite, bounded stability regions, meaning $h$ must be strictly constrained to stay within the region. Implicit methods (like Backward Euler or Trapezoidal rule) are often unconditionally A-stable.

Instead of approximating the derivatives directly like FDM, weighted residual methods construct an approximate continuous global solution $u(x) \approx \sum c_i \phi_i(x)$ using predefined trial functions $\phi_i(x)$.

• Collocation Method: Forces the differential equation to be satisfied exactly at specific isolated points (collocation points), yielding a system of linear equations for the unknown coefficients $c_i$.
• Galerkin Method: Forces the error (residual) of the differential equation to be mathematically orthogonal to the trial functions integrated over the entire domain. This powerful method forms the rigorous mathematical foundation for the Finite Element Method (FEM).

An iterative numerical approach used to selectively find the dominant (largest magnitude) eigenvalue and its corresponding eigenvector. It is highly efficient for massive engineering systems where computing all eigenvalues is computationally prohibitive.

Numerical solutions for ODEs are central to analyzing civil engineering systems driven by differential equations:

• Structural Dynamics (IVPs): Calculating the time-history response (displacement, velocity, acceleration) of mass-spring-damper models representing buildings subjected to earthquake ground motions or wind gusts.
• Beam Deflection (BVPs): Solving the Euler-Bernoulli fourth-order differential equation to find the continuous deflection curve of a beam subjected to complex distributed loads, constrained by fixed or pinned boundary conditions.
• Environmental Engineering: Modeling the concentration decay of pollutants in a river over time using coupled first-order reaction equations.