7. Ordinary Differential and Differential-Algebraic Equations

Introduction and Context

Many mechanical systems are governed by laws that relate a quantity to its rate of change. Newton’s second law, for example, converts forces into accelerations, while energy balances, heat-transfer relations, and constitutive equations generate evolution models. Once these models are discretized in space or expressed through lumped parameters, they commonly take the form of ordinary differential equations (ODEs). When algebraic constraints must simultaneously be satisfied, the model becomes a differential-algebraic equation (DAE).

The numerical treatment of these equations is central to computational mechanics. A method that is excellent for a slowly varying oscillator may fail for a stiff thermal system, a constrained multibody mechanism, or a highly dissipative control model. Numerical solution therefore requires more than applying a formula: the analyst must identify the problem type, select an appropriate discretization, control error, assess stability, and verify that physical constraints are respected.

This chapter develops methods for initial-value problems (IVPs), boundary-value problems (BVPs), and DAEs. Explicit and implicit one-step methods, multistep methods, shooting, finite differences, stiffness, stability, adaptivity, and constraint handling are considered together so that method selection can be connected to engineering judgment.

Learning Objectives

After studying this chapter, the reader should be able to:

1. Formulate mechanical models as ODEs, BVPs, or DAEs.

2. Derive and apply explicit and implicit time-integration methods.

3. Explain consistency, zero-stability, convergence, and absolute stability.

4. Distinguish stiff and non-stiff problems and select suitable solvers.

5. Use multistep, shooting, and finite-difference methods for engineering problems.

6. Apply adaptive step-size control and interpret numerical error.

7. Identify DAE index, consistent initial conditions, and constraint drift.

8. Evaluate numerical results using residuals, conservation laws, and physical reasoning.

From Mechanical Laws to Computational Problems

A first-order ODE system is written as

\[

\dot{\mathbf y}=\mathbf f(t,\mathbf y), \qquad \mathbf y(t_0)=\mathbf y_0.

\]

Higher-order equations are usually converted to this form by defining additional state variables. For a forced single-degree-of-freedom oscillator,

\[

m\ddot q+c\dot q+kq=F(t),

\]

let \(y_1=q\) and \(y_2=\dot q\). Then

\[

\dot y_1=y_2,\qquad \dot y_2=\frac{F(t)-cy_2-ky_1}{m}.

\]

The initial conditions specify the state at one time. This is an IVP.

A BVP specifies conditions at different locations or times. Steady heat conduction in a rod, for example, may satisfy

\[

\frac{d}{dx}\left(k(T)\frac{dT}{dx}\right)+s(x,T)=0,

\]

with temperature prescribed at one end and heat flux prescribed at the other. A second-order equation needs two independent conditions, but they need not be imposed at the same point.

A DAE has the semi-explicit form

\[

\dot{\mathbf y}=\mathbf f(t,\mathbf y,\mathbf z),\qquad

\mathbf 0=\mathbf g(t,\mathbf y,\mathbf z),

\]

where \(\mathbf y\) contains differential variables and \(\mathbf z\) contains algebraic variables. In a constrained pendulum, positions may be subject to a geometric relation while tension acts as an unknown multiplier enforcing that relation.

The distinction is computationally important. An ODE solver advances every state through a differential equation. A DAE solver must advance differential variables while solving algebraic equations and preserving constraints.

One-Step Time Integration

For a uniform step \(h\), the forward Euler method is

\[

\mathbf y_{n+1}=\mathbf y_n+h\mathbf f(t_n,\mathbf y_n).

\]

It is explicit: the new state is obtained without solving equations. Its local truncation error is \(O(h^2)\), and its global error is \(O(h)\). Although simple, it can be unsuitable for oscillatory or stiff systems because its stability region is limited.

Backward Euler evaluates the derivative at the new state:

\[

\mathbf y_{n+1}=\mathbf y_n+h\mathbf f(t_{n+1},\mathbf y_{n+1}).

\]

This requires a nonlinear solve. For a linear system \(\dot y=\lambda y\), it gives the amplification factor

\[

R(z)=\frac{1}{1-z},\qquad z=h\lambda.

\]

Backward Euler is first-order accurate but strongly stable for decaying modes. It is therefore valuable for stiff systems, even though it may introduce noticeable numerical damping.

The trapezoidal rule is

\[

\mathbf y_{n+1}=\mathbf y_n+\frac{h}{2}\left[\mathbf f(t_n,\mathbf y_n)+\mathbf f(t_{n+1},\mathbf y_{n+1})\right].

\]

It is second-order accurate and A-stable for linear test problems. It produces less artificial damping than backward Euler, but rapidly varying stiff components can lead to numerical oscillations when the step is too large.

Runge–Kutta methods obtain several derivative estimates within each step. The classical fourth-order method uses four stages and has global error \(O(h^4)\). Explicit Runge–Kutta methods are attractive for non-stiff dynamics because they require only function evaluations. Embedded pairs, such as a fourth/fifth-order pair, provide two estimates of different order. Their difference estimates local error and supports automatic step-size adjustment.

Stability, Consistency, and Stiffness

A numerical method is consistent if its local truncation error tends to zero as \(h\to0\). For multistep methods, zero-stability is also required: perturbations in the starting values must remain bounded as the number of steps increases. Consistency plus zero-stability implies convergence under standard assumptions.

Absolute stability is examined using the test equation

\[

\dot y=\lambda y,

\]

whose exact solution decays when \(\operatorname{Re}(\lambda)<0\). A numerical method is stable for a particular \(z=h\lambda\) when its amplification factor satisfies \(|R(z)|\le1\). The shape of the stability region matters more than order when the equations contain rapidly decaying modes.

A problem is stiff when stable physical solutions contain components with widely separated time scales, forcing an explicit method to use very small steps for stability rather than accuracy. Consider

\[

\dot y=-1000y+\sin t.

\]

The transient decays on a time scale of approximately \(10^{-3}\), whereas the forcing varies on a time scale of order one. Explicit Euler requires approximately \(|1-1000h|<1\), so \(h<0.002\), even after the transient is no longer important. An implicit method can use much larger steps while retaining stability.

Stiffness occurs in thermal networks, chemical or curing models, contact regularization, feedback control, fluid–structure coupling, and discretized diffusion equations. It is not defined solely by a large coefficient; it is a relation between system time scales, desired simulation interval, and method stability.

Nonlinear Solves and Jacobians

Implicit methods require the solution of a nonlinear residual equation. For backward Euler, define

\[

\mathbf F(\mathbf y_{n+1})=\mathbf y_{n+1}-\mathbf y_n-h\mathbf f(t_{n+1},\mathbf y_{n+1}).

\]

Newton’s method gives

\[

\left[\mathbf I-h\frac{\partial\mathbf f}{\partial\mathbf y}\right]\Delta\mathbf y=-\mathbf F,

\qquad \mathbf y\leftarrow\mathbf y+\Delta\mathbf y.

\]

A good predictor, usually an explicit Euler estimate, improves convergence. The Jacobian may be analytic, approximated by finite differences, or evaluated using automatic differentiation. In large systems, sparse linear algebra and matrix-free Krylov methods can greatly reduce computational cost.

A nonlinear solve should not be confused with time-integration accuracy. A loose nonlinear tolerance can produce an error larger than the requested time-discretization error; an unnecessarily tight tolerance wastes computation. Solver tolerances should be related to absolute and relative solution tolerances.

Multistep Methods

Multistep methods reuse previously computed values. A two-step Adams–Bashforth formula is

\[

y_{n+1}=y_n+h\left(\frac32f_n-\frac12f_{n-1}\right),

\]

and is explicit. The two-step Adams–Moulton formula is

\[

y_{n+1}=y_n+h\left(\frac{5}{12}f_{n+1}+\frac{2}{3}f_n-\frac{1}{12}f_{n-1}\right),

\]

and is implicit. Predictor–corrector schemes commonly use Adams–Bashforth for prediction and Adams–Moulton for correction.

Backward differentiation formulas (BDFs) are implicit multistep methods designed particularly for stiff problems. Their advantages include good stability for dissipative systems and efficient use of function evaluations after the history has been established. Their limitations include the need for starting values, sensitivity to rapidly changing step sizes, and reduced suitability for strongly oscillatory problems.

Adaptive Integration

A fixed step is rarely optimal over a complete simulation. An adaptive method estimates the local error and chooses a new step. If the estimated error is \(e_n\), a typical update is

\[

h_{new}=s h\left(\frac{tol}{e_n}\right)^{1/(p+1)},

\]

where \(p\) is the order of the accepted method and \(s<1\) is a safety factor. In practice, error is scaled componentwise:

\[

e_i^{scaled}=\frac{|e_i|}{atol_i+rtol_i|y_i|}.

\]

This avoids treating a small displacement and a large temperature with the same absolute criterion. A step is accepted when a norm of the scaled error is below unity. Adaptive algorithms must also limit excessive step growth, reject repeated failures, and detect events such as impacts or valve switching.

Boundary-Value Problems: Shooting and Finite Differences

In the shooting method, an unknown initial condition is guessed, the IVP is integrated, and the terminal boundary error is computed. For

\[

y''=f(x,y,y'),\qquad y(a)=\alpha,\quad y(b)=\beta,

\]

set \(y(a)=\alpha\) and guess \(y'(a)=s\). Integration produces \(y(b;s)\). The scalar residual is

\[

\Phi(s)=y(b;s)-\beta.

\]

A root-finding method adjusts \(s\) until \(\Phi(s)=0\). Secant or Newton iteration may be used. Sensitivity equations can provide the derivative \(d\Phi/ds\).

Shooting is conceptually simple and effective when the IVP is stable. It may fail when small errors in the initial slope grow strongly, as in unstable or stiff boundary-value problems. Multiple shooting divides the domain into subintervals, integrates separately, and imposes continuity conditions. This reduces sensitivity and is often more robust.

Finite differences convert derivatives at grid points into algebraic equations. For a linear problem

\[

y''+p(x)y'+q(x)y=r(x),

\]

central differences give

\[

y''(x_i)\approx\frac{y_{i+1}-2y_i+y_{i-1}}{h^2},

\qquad

y'(x_i)\approx\frac{y_{i+1}-y_{i-1}}{2h}.

\]

Substitution yields a sparse algebraic system for the nodal values. Central differences are generally second-order accurate on a uniform grid. Boundary conditions must be discretized consistently; a derivative condition may use a one-sided formula or a ghost point. Nonlinear BVPs lead to nonlinear algebraic systems, usually solved by Newton iteration.

Finite differences are often preferable to shooting for stiff BVPs, complicated boundary conditions, and solutions with interior layers. Their accuracy depends on mesh quality, and adaptive mesh refinement may be necessary where gradients are large.

Differential-Algebraic Equations

A DAE can be represented as

\[

\mathbf F(t,\mathbf y,\dot{\mathbf y})=\mathbf0.

\]

The index measures, informally, how many differentiations are needed to obtain an explicit ODE form. Index-one semi-explicit systems can often be solved by differentiating or locally eliminating algebraic variables. Higher-index systems are more sensitive to perturbations and require careful treatment.

For a constrained mechanical system with coordinates \(\mathbf q\), mass matrix \(\mathbf M\), constraint \(\boldsymbol\phi(\mathbf q,t)=0\), and Lagrange multipliers \(\boldsymbol\lambda\), a common formulation is

\[

\mathbf M\ddot{\mathbf q}=\mathbf f(\mathbf q,\dot{\mathbf q},t)+\mathbf G^T\boldsymbol\lambda,

\qquad \boldsymbol\phi(\mathbf q,t)=0,

\]

where \(\mathbf G=\partial\boldsymbol\phi/\partial\mathbf q\). The multipliers represent constraint forces. Initial positions and velocities must satisfy both position-level and velocity-level constraints; otherwise the numerical problem begins inconsistently.

Naive integration of the differential equations can cause constraint drift, so that \(\boldsymbol\phi\ne0\) after several steps. Remedies include solving constraints simultaneously at each step, projection of the state back onto the constraint manifold, stabilization terms, and coordinate reduction. Projection improves geometric satisfaction but can alter energy or momentum unless designed carefully.

For index-one DAEs, implicit Euler or BDF methods are often effective. The nonlinear system includes both differential update equations and algebraic constraints. Consistent initialization is a separate computational task: algebraic variables and, where necessary, differential derivatives must be determined so that all equations and differentiated constraints hold.

Applied Case Study: Constrained Four-Bar Mechanism

Background

Consider an illustrative planar four-bar mechanism used to transmit rotary motion. The crank angle is prescribed as a function of time, while the coupler and rocker coordinates are unknown. Link lengths impose two loop-closure equations.

Situation and Problem

If the mechanism is simulated by integrating unconstrained Cartesian coordinates, round-off and truncation errors gradually violate loop closure. Near toggle positions, the constraint Jacobian becomes poorly conditioned, and an explicit solver produces large, nonphysical excursions unless an extremely small step is used.

Analysis

Let \(\mathbf q\) contain the Cartesian coordinates of selected joints and let \(\boldsymbol\phi(\mathbf q,t)=0\) represent the two loop-closure relations. The equations of motion are augmented with Lagrange multipliers. The resulting DAE is more faithful than replacing the links by stiff springs, because spring regularization introduces artificial high frequencies and can create severe stiffness.

The crank motion and initial configuration are first checked for consistency. At each implicit time step, Newton iteration solves the dynamic residual together with \(\boldsymbol\phi=0\). The Jacobian is assembled in block form, including the constraint Jacobian.

Approach and Solution

An implicit BDF integrator is selected because the mechanism becomes numerically stiff near the toggle configuration. Relative and absolute tolerances are scaled separately for positions, velocities, and multipliers. The code monitors loop-closure residuals, Newton convergence, input work, and reaction forces. An event detector identifies the prescribed crank limits and prevents steps from crossing discontinuities unnoticed.

Outcome

The computed positions remain on the constraint manifold to within the prescribed tolerance. Reaction forces are obtained from the multipliers, rather than inferred indirectly from spring forces. A comparison with a reduced-coordinate model is used as an independent verification for configurations where coordinate reduction is nonsingular.

Learning Points

The case demonstrates that model formulation determines numerical behavior. Constraint equations are not optional post-processing conditions; they are part of the governing problem. Implicit integration, consistent initialization, residual monitoring, and independent verification are essential near singular or nearly singular configurations.

Practical Applications and Professional Judgment

ODE solvers support dynamic analysis of suspension systems, rotating machinery, robotics, thermal management, batteries, process equipment, and control systems. BVP methods are used for beam deflection, lubrication films, boundary layers, heat transfer, and optimal trajectories. DAEs occur in multibody dynamics, electrical–mechanical systems, incompressible flow discretizations, and network models.

A professional workflow begins with nondimensionalization and model inspection. Estimate characteristic time scales, identify discontinuities, determine whether Jacobians are sparse, and decide whether conservation or constraint preservation is critical. Select the simplest method that satisfies accuracy and stability requirements, then verify the result through mesh or tolerance refinement.

Common mistakes include using explicit Euler for a stiff system, imposing BVP conditions as if they were IVP data, ignoring inconsistent DAE initialization, accepting a converged nonlinear solve with a large residual, and reporting many digits unsupported by discretization error. Solver success flags do not guarantee a physically correct solution.

Numerical models should also be traceable. Units, tolerances, solver versions, initial conditions, and event rules should be recorded. In safety-critical engineering, independent implementations, benchmark problems, sensitivity studies, and conservative validation procedures are appropriate. Automatic differentiation, sparse direct solvers, parallel integration, and surrogate models increasingly support large simulations, but automation does not replace physical verification. Machine-learning surrogates must preserve relevant constraints and should not be trusted outside their validated training domain.

Current Developments

Modern simulation software increasingly combines adaptive implicit integration, automatic Jacobians, event handling, and sparse linear algebra. Structure-preserving methods aim to conserve energy, momentum, or geometric properties over long simulations. For Hamiltonian mechanical systems, symplectic methods can control long-term energy error better than conventional methods, even when their pointwise order is lower.

DAE formulations are also being integrated with optimization, parameter estimation, and digital twins. In these applications, derivatives of outputs with respect to design variables are important; sensitivity equations and adjoint methods can be more efficient than repeated finite-difference simulations. The central trend is not merely faster integration, but more reliable integration of models containing constraints, multiple scales, uncertainty, and coupled physical domains.

Best-Practice Checklist

• Classify the model as IVP, BVP, ODE, or DAE before selecting a solver.
• Estimate time scales and inspect stiffness.
• Use tolerances that reflect the physical scales of each variable.
• Supply accurate Jacobians or exploit automatic differentiation when possible.
• Treat nonlinear and linear solver tolerances as part of the error budget.
• Check residuals, constraints, conservation laws, and event locations.
• Refine time steps or spatial meshes to establish convergence.
• Use independent formulations or benchmark solutions for verification.
• Record computational settings so results are reproducible.

Chapter Summary

Ordinary differential equations describe evolution through differential laws, whereas differential-algebraic equations combine evolution with instantaneous constraints. Initial-value problems are advanced in time using explicit, implicit, Runge–Kutta, or multistep methods. Explicit methods are efficient for non-stiff systems; implicit methods are usually more reliable for stiff dynamics but require nonlinear and linear solves. Stability, rather than formal order alone, determines whether a method is suitable.

Boundary-value problems can be solved by shooting, multiple shooting, or finite differences. Shooting converts the problem into an IVP plus a root-finding problem, while finite differences produce a global algebraic system and are often more robust for unstable or stiff cases. Adaptive methods regulate local error by changing the step size. DAE simulation requires consistent initialization, appropriate index treatment, and active monitoring of constraint residuals. Across all problem classes, verification through refinement, residuals, physical laws, and independent comparisons is essential.

Key Terms

ODE: An equation involving derivatives of unknown functions with respect to one independent variable.

IVP: A problem in which all required conditions are specified at one initial point.

BVP: A problem whose conditions are imposed at two or more points.

DAE: A coupled system of differential and algebraic equations.

Stiffness: A numerical difficulty caused by widely separated time scales and restrictive stability requirements.

Consistency: The property that local discretization error vanishes as the step tends to zero.

Zero-stability: Bounded propagation of perturbations in a multistep method.

Convergence: Tendency of the numerical solution toward the exact solution as discretization is refined.

Absolute stability: Stability of a numerical method for a specified value of \(h\lambda\).

Amplification factor: The ratio relating successive numerical solution values for a test equation.

Implicit method: A method in which the new solution appears inside the update equation.

Explicit method: A method in which the new solution is obtained without solving an implicit equation.

Shooting method: A BVP technique that adjusts unknown initial data until a terminal condition is met.

Multiple shooting: A shooting method using several subintervals and continuity conditions.

Constraint drift: Accumulation of numerical error that violates algebraic or geometric constraints.

DAE index: An indication of the number of differentiations required to obtain an explicit ODE structure.

Consistent initialization: Initial data satisfying differential, algebraic, and necessary differentiated constraints.

Jacobian: The matrix of partial derivatives used to linearize nonlinear residual equations.

Adaptive integration: Time integration in which the step size changes according to an error estimate.

Residual: The imbalance remaining when an approximate solution is substituted into governing equations.

Review Questions

Short Answer Questions

1. What distinguishes an IVP from a BVP?

2. Why can an explicit method be inefficient for a stiff problem?

3. Define consistency and zero-stability.

4. What is constraint drift in a DAE simulation?

5. Why are Jacobians useful in implicit time integration?

Long Answer Questions

1. Derive the backward Euler method for a first-order ODE and discuss its stability for the test equation.

2. Compare shooting and finite-difference methods for solving a second-order BVP.

3. Explain stiffness and compare the suitability of explicit Runge–Kutta and implicit BDF methods.

4. Describe adaptive step-size control using an embedded error estimate.

5. Explain the formulation, initialization, and numerical treatment of a constrained mechanical DAE.

Analytical / Application Questions

1. For \(\dot y=-500y+\cos t\), estimate the explicit Euler stability restriction and explain why an implicit method may be preferable over a long simulation.

2. A shooting computation converges for one initial-slope guess but diverges for nearby guesses. Diagnose the likely numerical cause and propose two remedies.

3. A multibody simulation shows increasing loop-closure error while Newton iterations converge at every step. Explain why this can occur and design a verification procedure to distinguish integration error from model inconsistency.

22. Multiple Choice Questions with Answers

Question 1

Which formulation most directly represents a constrained mechanical system with algebraic position relations?

A. An explicit Euler formula only

B. A DAE with differential equations and constraint equations

C. A scalar interpolation polynomial

D. A finite-difference equation without boundary conditions

Correct Answer: B. A DAE with differential equations and constraint equations

Question 2

For the test equation y' = λy with negative real λ, explicit Euler is stable when which condition holds for real λ?

A. 0 < hλ < 2

B. -2 < hλ < 0

C. hλ > 2

D. hλ = 1

Correct Answer: B. -2 < hλ < 0

Question 3

Why is backward Euler often selected for stiff systems?

A. It never requires a nonlinear solve

B. It has a large stability region for decaying modes

C. It is always fourth-order accurate

D. It exactly conserves mechanical energy

Correct Answer: B. It has a large stability region for decaying modes

Question 4

In a shooting method, what quantity is normally driven to zero by adjusting the guessed initial condition?

A. The local truncation error at every mesh point

B. The terminal boundary residual

C. The Jacobian determinant

D. The time-step ratio

Correct Answer: B. The terminal boundary residual

Question 5

Which feature commonly indicates stiffness?

A. Only one slowly varying time scale

B. Widely separated decay time scales restricting explicit step size

C. Absence of differential terms

D. A boundary condition imposed at one endpoint

Correct Answer: B. Widely separated decay time scales restricting explicit step size

Question 6

What is the principal purpose of an embedded Runge–Kutta pair?

A. To eliminate all round-off error

B. To estimate local error and adapt the step size

C. To convert every BVP into a DAE

D. To enforce algebraic constraints exactly

Correct Answer: B. To estimate local error and adapt the step size

Question 7

A DAE simulation has converged Newton iterations but increasing constraint violation. Which issue is most directly suggested?

A. Constraint drift during time integration

B. Excessive spatial resolution

C. A guaranteed exact solution

D. Zero numerical damping

Correct Answer: A. Constraint drift during time integration

Question 8

Finite-difference discretization of a second derivative on a uniform grid is commonly based on which expression?

A. (yᵢ₊₁ − yᵢ)/h

B. (yᵢ₊₁ − 2yᵢ + yᵢ₋₁)/h²

C. (yᵢ₊₁ + yᵢ₋₁)/2h

D. h²(yᵢ₊₁ − yᵢ₋₁)

Correct Answer: B. (yᵢ₊₁ − 2yᵢ + yᵢ₋₁)/h²

Question 9

What must be true of DAE initial data for a physically meaningful simulation?

A. They must satisfy the algebraic and necessary differentiated constraints

B. They must maximize the initial residual

C. They must be generated by explicit Euler

D. They must make all multipliers zero

Correct Answer: A. They must satisfy the algebraic and necessary differentiated constraints

Question 10

Which practice provides the strongest evidence that a computed engineering solution is reliable?

A. Reporting more decimal places

B. Checking refinement, residuals, constraints, and physical laws

C. Using the largest possible time step

D. Ignoring solver tolerance information

Correct Answer: B. Checking refinement, residuals, constraints, and physical laws