Numerical Simulation Methods

Introduction

Mechatronic systems combine mechanical motion, electrical circuits, actuators, sensors, controllers, and digital computation. Their behaviour is commonly represented by ordinary differential equations (ODEs), difference equations, algebraic constraints, or combinations of these in a differential-algebraic equation (DAE) model. Although an analytical solution is valuable, realistic systems rarely permit one. Numerical simulation therefore provides the practical bridge between a mathematical model and engineering decisions.

A numerical simulation replaces continuous evolution with a sequence of computed approximations. The quality of the result depends not only on the integration formula, but also on model formulation, step size, event handling, solver selection, error control, and verification. A simulation that produces a smooth graph may still be physically wrong if it violates stability, conservation, actuator limits, or event timing.

This chapter develops the principal numerical methods used in mechatronics and explains how to assess convergence, stability, stiffness, computational cost, and failure modes.

Learning Objectives

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

• formulate a continuous-time mechatronic model for numerical integration;
• distinguish explicit and implicit integration methods;
• apply Euler, Runge–Kutta, and adaptive-step methods;
• explain local and global truncation error, convergence, and stability;
• select an appropriate solver for nonstiff, stiff, hybrid, and constrained systems;
• detect and handle events such as impacts, switching, saturation, and limit crossings;
• diagnose numerical failure and verify simulation results;
• evaluate accuracy against computational cost in an engineering application.

1. From a Physical Model to a Numerical Problem

A common state-space model is

where t is time, x is the state vector, u is the input, p is a parameter vector, f is the state derivative, y is the output, and g is the output function. The initial-value problem requires an initial state x(t₀) = x₀. If x contains position, velocity, current, and temperature, its components may have different units and widely different magnitudes.

A numerical method computes approximations x_n ≈ x(t_n) at times t_n. For a fixed step size, t_n+1 = t_n + h, where h has units of time. Unit consistency is essential: if stiffness is expressed in N/m and mass in kg, the resulting natural frequency is in s⁻¹, and the integration step must be expressed in seconds.

Second-order mechanical equations can be converted to first-order form. For example, a mass–spring–damper system

becomes

with x₁ = q and x₂ = q˙. Here m is mass, c is viscous damping, k is stiffness, and F is force. This transformation is valid when m ≠ 0 and the force law is sufficiently well defined.

A simulation model must also specify initial conditions, parameter values, input signals, output definitions, discontinuities, constraints, and validity limits. Numerical integration cannot correct an incorrectly formulated physical model.

2. One-Step Integration Methods

2.1 Forward Euler method

The forward Euler method evaluates the derivative at the beginning of the interval:

It is explicit, simple, and inexpensive. Its local truncation error is O(h²) and its global error is O(h), provided f is sufficiently smooth. Local error refers to the error introduced over one step assuming the preceding value is exact; global error is the accumulated error over the complete simulation.

For the scalar linear equation x˙ = λx, Euler gives

The factor R(z) = 1 + z, where z = hλ, is the amplification factor. Stability requires the numerical solution not to grow when the exact solution decays. For real negative λ, the condition is

or equivalently 0 ≤ h|λ| ≤ 2. Thus a small physical time constant can impose a very small step even when the output of interest changes slowly.

2.2 Backward Euler method

Backward Euler evaluates the derivative at the unknown end point:

It is implicit because x_n+1 appears on both sides. A nonlinear equation must generally be solved at each step, often using Newton iteration. For x˙ = λx, its amplification factor is

Backward Euler is first-order accurate but A-stable: all modes with Re(λ) < 0 remain numerically stable for any positive h. It is therefore useful for stiff systems, although excessive step sizes can still produce poor transient accuracy and numerical damping.

2.3 Trapezoidal rule

The implicit trapezoidal method averages the beginning and ending derivatives:

It is second-order accurate and A-stable. Compared with backward Euler, it introduces less numerical damping, but it may produce oscillations when applied to very stiff transients or discontinuous inputs. The method is particularly useful when phase accuracy and smooth dynamics are important.

2.4 Runge–Kutta methods

The classical fourth-order Runge–Kutta method, RK4, uses four derivative evaluations:

Here each kᵢ has units of state per unit time, so multiplication by h produces state units. RK4 has global error O(h⁴) for sufficiently smooth problems. It is explicit and often effective for nonstiff mechanical and control models. Its high order does not make it universally superior: a stiff problem can force RK4 to use prohibitively small steps for stability.

3. Adaptive Step-Size Integration

A fixed step is easy to reproduce but inefficient when the system contains both slow and fast intervals. Adaptive methods estimate the local error and change h. A common embedded Runge–Kutta pair computes two approximations of different order, for example orders p and p − 1. Their difference estimates local error:

The norm must account for component scaling. A practical weighted root-mean-square error is

where N is the number of state components, eᵢ is the estimated error in component i, atolᵢ is an absolute tolerance, and rtol is a relative tolerance. Both tolerances are dimensionally consistent with the corresponding state component. A step is accepted when E ≤ 1; otherwise it is rejected and repeated with a smaller step.

A typical proposed step is

for an embedded method whose accepted solution has order p. The safety factor is less than one, and practical solvers impose minimum and maximum step limits. This formula is a control heuristic rather than a guarantee; discontinuities, failed nonlinear iterations, or rapidly changing Jacobians require additional safeguards.

Absolute tolerance is important when a state is near zero; relative tolerance is important when its magnitude is large. In a robot model, position in metres, velocity in m/s, and motor current in amperes should not be given identical unscaled tolerances without considering engineering significance.

4. Stability, Convergence, and Accuracy

4.1 Consistency and convergence

A method is consistent if its discrete equation approaches the differential equation as h → 0. A method is convergent if its numerical solution approaches the exact solution as the step size tends to zero. For well-posed ordinary differential equations, consistency plus suitable stability leads to convergence; this is the practical meaning of the equivalence commonly associated with the Lax principle for linear initial-value problems.

Order describes the rate of error reduction. If the global error is approximately C h^p, reducing h by two should reduce error by approximately 2^p when the asymptotic regime has been reached. This relationship may fail near discontinuities, solver tolerances, round-off limits, or instability.

4.2 Absolute and relative stability

For a linear test equation x˙ = λx, a one-step method has an amplification function R(z). Its absolute stability region is the set of complex values z satisfying

For physical dissipative modes, Re(λ) is negative. A-stability means the complete left half-plane is stable. L-stability additionally requires strong damping as z tends to negative infinity, a useful property for very fast decaying stiff modes. Backward Euler is L-stable; trapezoidal integration is A-stable but not L-stable.

Stability is not the same as accuracy. A stable solution may have unacceptable phase lag, amplitude error, or event timing error. Conversely, an accurate method used with a step outside its stability region can diverge.

4.3 Step-size convergence test

A practical verification procedure is to run the same model with step sizes h, h/2, and h/4. Compare a quantity of interest, such as final position or peak current. If differences decrease at the expected rate, the calculation is likely in a convergent regime. This is evidence, not proof: model errors and unresolved events can remain hidden.

5. Stiffness and Solver Selection

A system is stiff when stable physical modes have widely separated time scales, so an explicit method is restricted by the fastest mode even though the engineer is interested mainly in slower behaviour. In a linear system x˙ = A x, eigenvalues of A provide a first indication. If two stable eigenvalues satisfy |λ_fast| ≫ |λ_slow|, the ratio of time scales is large.

A classic mechatronic example is a motor with electrical and mechanical dynamics. Electrical inductance may produce a time constant of milliseconds, whereas mechanical acceleration may evolve over tenths of a second. An explicit solver must resolve the electrical mode for stability. An implicit solver can step over much of the fast transient while retaining its influence on the slow state.

No table can replace analysis of the model. The selection must consider accuracy, available derivatives or Jacobians, memory, real-time deadlines, reproducibility, and the cost of evaluating the model.

For an implicit method, the nonlinear solve may use Newton's method. With residual

Newton updates are

where J_G is the Jacobian of G. For a smooth model, J_G = I − h ∂f/∂x. A poor initial guess, singular Jacobian, or excessively large step may cause failure. Analytic or automatically generated Jacobians can substantially reduce cost for large systems.

6. Discontinuities and Event Detection

Mechatronic systems are hybrid: continuous states interact with discrete logic. Events include a wheel contacting the ground, a relay switching, a limit switch changing state, a controller saturating, or a safety threshold being crossed.

Represent an event by a guard function

An event occurs when the guard changes sign or reaches zero, subject to a direction condition. If a solver steps from tₙ to tₙ₊₁ and finds opposite signs, it brackets the event. A root-finding method such as bisection or safeguarded secant interpolation locates the event time within the bracket. The solver should then stop at the event, apply the discrete state transition, modify continuous states if required, and restart integration.

For an impact, position may remain continuous while velocity changes according to a restitution law. A simple one-dimensional model is

where v⁻ and v⁺ are velocities immediately before and after impact and e is the coefficient of restitution, with 0 ≤ e ≤ 1 for a passive impact. This equation is a model assumption; it does not apply to every contact mechanism.

Simply evaluating the event at the end of a large step can produce substantial timing error, penetrate a mechanical stop, or miss a short pulse. After an event, the solver must prevent repeated detection of the same zero by using event direction, hysteresis, a mode flag, or a carefully defined reset rule. Algebraic loops and simultaneous events require an explicit priority or consistency policy.

Saturation deserves special care. If a command is clipped, the clipped value should enter the plant model, while the controller may also require anti-windup logic. Treating saturation only as a plotted output produces an incorrect simulation.

7. Discretization of Control and Physical Models

Digital controllers operate at sampling times. A continuous plant can be approximated by a discrete state model

For a linear continuous plant x˙ = A x + B u with zero-order-held input over a sample period T,

The matrix exponential formulation is exact for the stated linear, constant-input assumptions. If A is nonsingular, B_d = A^(−1)(A_d − I)B; the integral form is safer for singular A.

For nonlinear systems, discretization may use numerical integration over each sample interval. The sample period must be small enough to represent relevant dynamics and switching behaviour, but real-time computation must finish before the next sample. A controller simulation can therefore fail either numerically or operationally: a mathematically accurate algorithm is unusable if it misses its deadline.

8. Worked Example: Euler Stability for a Motor–Load Approximation

Consider a first-order current mode

where τ = 2 ms is the decay time constant. For the homogeneous part, λ = −1/τ = −500 s⁻¹. Forward Euler is stable only if

Therefore

A 1 ms step satisfies the stability condition. A 5 ms step does not: h|λ| = 0.005 × 500 = 2.5, so the amplification factor is 1 − 2.5 = −1.5, whose magnitude exceeds one. The computed current alternates sign and grows, even though the physical current decays.

This result does not mean that 4 ms is an accurate step limit. Stability only prevents unbounded growth. To resolve a 2 ms transient accurately, a step substantially smaller than 2 ms may be needed. An implicit solver could remain stable at 5 ms or more, but it might still fail to reproduce the initial current transient accurately unless tolerances and output requirements permit that approximation.

9. Detailed Applied Case: Servo Axis with a Hard Travel Limit

Background

Consider an illustrative single-axis servo consisting of a motor, gearbox, carriage, encoder, position controller, and a hard mechanical stop. The objective is to predict whether the carriage reaches its target without excessive impact force and whether the controller remains stable.

Problem

The model contains continuous position and velocity, motor current dynamics, actuator saturation, and a discontinuous contact event. A fixed-step simulation with 10 ms sampling appears smooth, but the carriage penetrates the stop by several millimetres. The design team must determine whether this is a physical fault or a numerical artefact.

Analysis

The mechanical position is integrated with velocity. The motor current has a much faster time constant than the carriage. Thus the model is mildly stiff. The stop is represented by the guard g(q) = q − q_max. When g = 0 with positive velocity, contact occurs. If the stop is treated merely as a large spring evaluated at sample instants, the model can generate unrealistic energy and missed contact timing.

The first simulation is repeated with an adaptive nonstiff solver and event location. The event is bracketed and located by interpolation, then the contact law is applied. A second run uses a stiff-capable implicit solver. Tolerances are scaled separately for position, velocity, and current. Step-size refinement is performed around the impact.

Approach and solution

The simulation workflow is:

1. formulate the motor and load equations in consistent SI units;
2. include the actual saturated voltage or torque in the plant equations;
3. define the stop guard and specify event direction;
4. stop at the guard root and apply a contact or restitution model;
5. restart the solver in contact mode;
6. compare explicit and implicit solutions;
7. refine tolerances and maximum step near contact;
8. check energy, travel limits, and controller output against independent calculations.

The coarse fixed-step result is rejected because it is not event-converged. The event-aware solution predicts contact 6.8 ms earlier than the coarse result and reduces computed penetration to the permitted numerical tolerance. The stiff solver takes fewer model evaluations during the fast current transient, while the explicit solver is cheaper per step. Both agree on peak position and impact velocity after refinement, providing stronger evidence of correctness.

Outcome and learning points

The apparent excessive penetration was primarily a timing and event-discretization error, not a confirmed mechanical failure. However, the refined simulation also shows that saturation causes controller windup unless anti-windup is enabled. The case demonstrates that solver choice, event handling, and control logic are inseparable in mechatronic simulation.

10. Computational Cost and Verification

A useful cost measure is the number of right-hand-side evaluations, Jacobian evaluations, nonlinear iterations, rejected steps, and total wall-clock time. An implicit method may perform fewer large steps but require expensive matrix factorizations. An explicit method may be preferable for a small nonstiff model or a hard real-time target.

Verification should address at least four questions:

1. Are the equations implemented correctly? Compare derivatives with hand calculations, limiting cases, dimensional analysis, and conservation laws.
2. Is the numerical solution resolved? Perform step-size and tolerance studies, including event-time convergence.
3. Is the physical model valid? Compare parameters and assumptions with experiments or manufacturer data.
4. Are outputs interpreted correctly? Distinguish solver interpolation from computed states and report units, tolerances, and solver settings.

A simulation should record solver status, rejected steps, minimum step, maximum step, event times, nonlinear convergence, and any warning messages. Reproducibility requires versioned model files, parameter sets, initial conditions, solver settings, and random seeds where stochastic elements exist.

11. Failure Modes and Limitations

Common failure modes include:

• Divergence: the step is outside the stability region, or the model itself is unstable.
• Excessive step rejection: tolerances are too strict, the model is discontinuous, or the Jacobian is poorly conditioned.
• Missed events: the step spans multiple crossings or the guard function is not monitored.
• Chattering: noise or numerical error repeatedly changes a switching condition; hysteresis or a minimum dwell time may be needed.
• Constraint drift: integrating constrained coordinates independently violates geometric or algebraic constraints.
• Algebraic-loop failure: instantaneous dependencies cannot be resolved with the selected formulation.
• Artificial damping or oscillation: the integration method changes amplitude or phase.
• Underflow, overflow, or loss of significance: state scaling and nondimensionalization may be required.
• False confidence: a visually smooth plot conceals incorrect units, parameters, saturation, or event timing.

DAE models require consistent initial conditions satisfying both differential and algebraic equations. Index reduction, coordinate selection, or a DAE solver may be necessary. Contact and friction models can be nonsmooth; a smooth numerical result is not evidence that a nonsmooth physical law has been represented correctly.

12. Best Practices and Contemporary Tools

Use a solver with adaptive error control for exploratory offline work, and use fixed-step simulation only when its step is justified by stability, accuracy, and real-time requirements. Scale states and tolerances by physical importance. Limit the maximum step near known events and rapidly changing inputs. Supply Jacobians or sparsity information when available. Separate model equations from solver configuration so that the model can be tested independently.

Automatic differentiation and symbolic tools can generate consistent derivatives and Jacobians, reducing implementation errors. Hardware-in-the-loop systems increasingly use real-time targets, where deterministic execution and overruns are critical. Co-simulation standards allow separate mechanical, electrical, and control models to exchange variables, but communication step size, interpolation, and power exchange can introduce additional numerical error. AI-based surrogate models may reduce computational cost, but they require validation over the operating envelope and must not silently replace safety-critical physical constraints.

Ethically and professionally, engineers should not report only the solver output that supports a design decision. Numerical uncertainty, failed runs, parameter uncertainty, and model limitations should be documented. In safety-related systems, simulation evidence must be complemented by testing, traceability, independent review, and applicable standards and regulatory requirements.

Chapter Summary

Numerical simulation converts continuous and hybrid mechatronic models into computable approximations. Euler, trapezoidal, Runge–Kutta, and implicit methods differ in order, stability, cost, and suitability for stiff systems. Consistency, convergence, and stability must be assessed separately from physical-model validity. Adaptive step control balances accuracy and computation, while scaled absolute and relative tolerances are essential for multivariable systems.

Stiff systems often require implicit methods; nonstiff systems commonly benefit from explicit Runge–Kutta methods. Events such as impacts, switching, saturation, and limit crossings must be located and handled explicitly. Solver selection should reflect dynamics, accuracy, real-time constraints, event structure, and available derivatives. Step refinement, tolerance studies, dimensional checks, conservation checks, and independent validation are central to trustworthy simulation. Numerical results are engineering evidence, not automatically engineering truth.

Key Terms

Numerical integration: Approximation of a differential equation using discrete computational steps.

State: Minimum set of variables needed to describe system evolution.

Explicit method: Method whose next state is computed directly from known quantities.

Implicit method: Method requiring solution of an algebraic or nonlinear equation for the next state.

Local truncation error: Error introduced during one step assuming exact prior data.

Global error: Accumulated error over a simulation interval.

Consistency: Property that the discrete method approaches the governing equation as step size tends to zero.

Convergence: Property that the numerical solution approaches the exact solution as step size tends to zero.

Stability region: Set of step-size and system-eigenvalue combinations for which numerical modes remain bounded.

A-stability: Stability for all continuous modes with negative real part.

Stiffness: Presence of time scales that strongly restrict explicit step size despite stable physical dynamics.

Adaptive step size: Variable integration interval selected using an error estimate.

Tolerance: User-specified bound controlling numerical error.

Event function: Guard function used to identify a discrete transition or discontinuity.

Hybrid system: System combining continuous dynamics and discrete modes.

Saturation: Limitation that clips an actuator or signal to an allowable range.

DAE: Differential-algebraic equation containing differential and algebraic relationships.

Jacobian: Matrix of partial derivatives used in nonlinear solution and sensitivity calculations.

Numerical damping: Artificial reduction of amplitude introduced by a discretization method.

Co-simulation: Coupled simulation in which separate subsystem models exchange variables.

Review Questions

Short Answer Questions

1. What is the difference between local truncation error and global error?
2. Why can an explicit solver be inefficient for a stiff motor model?
3. Define an event function and explain how a solver locates an event.
4. What is the purpose of using both absolute and relative tolerances?
5. Why does a stable numerical solution not necessarily imply an accurate solution?

Long Answer Questions

1. Derive the forward Euler update for a first-order state-space model and discuss its accuracy and stability limitations.
2. Compare forward Euler, backward Euler, trapezoidal integration, and RK4 in terms of implicitness, order, stability, and computational cost.
3. Explain stiffness using eigenvalues and describe why implicit methods are advantageous for stiff systems.
4. Describe a complete verification procedure for a numerical simulation of a mechatronic system.
5. Explain how impacts, saturation, and switching should be represented in a hybrid simulation.

Analytical / Application Questions

1. For x˙ = −250x, determine the largest stable forward-Euler step and explain whether that step should automatically be considered accurate.
2. A robot joint model contains a 1 ms electrical time constant and a 0.5 s mechanical time constant. Recommend a solver strategy for offline simulation and for hard real-time simulation, justifying both choices.
3. A fixed-step simulation of a servo stop shows penetration and oscillatory contact force. Propose a diagnostic and correction plan involving event detection, step refinement, contact modelling, and convergence testing.

Multiple Choice Questions with Answers

Question 1

For the forward Euler method applied to x˙ = λx, which amplification factor is correct?

A. R(z) = 1 + z, where z = hλ

B. R(z) = 1/(1 − z)

C. R(z) = e^z

D. R(z) = 1 − z²

Correct Answer: A. R(z) = 1 + z, where z = hλ

Question 2

What is the global order of the classical RK4 method for a sufficiently smooth nonstiff problem?

A. First order

B. Second order

C. Third order

D. Fourth order

Correct Answer: D. Fourth order

Question 3

Which characteristic most strongly indicates stiffness?

A. All state variables have identical units

B. Stable modes have widely separated time scales

C. The input is constant

D. The model has only one state

Correct Answer: B. Stable modes have widely separated time scales

Question 4

A numerical step for an adaptive solver is normally accepted when the scaled error measure is:

A. Exactly zero

B. Greater than one

C. Less than or equal to one

D. Equal to the maximum step size

Correct Answer: C. Less than or equal to one

Question 5

For forward Euler applied to x˙ = −500x, what is the largest stable positive step size?

A. 0.5 ms

B. 2 ms

C. 4 ms

D. 5 ms

Correct Answer: C. 4 ms

Question 6

What is the main purpose of event root finding?

A. To increase the model order

B. To locate the time at which a guard function reaches zero

C. To remove all discontinuities from the physical model

D. To guarantee that the plant is stable

Correct Answer: B. To locate the time at which a guard function reaches zero

Question 7

Which method is generally most appropriate as a first choice for a stiff continuous model?

A. Large-step forward Euler

B. Implicit BDF or Radau-type integration

C. A plotting routine without error control

D. A purely algebraic solver

Correct Answer: B. Implicit BDF or Radau-type integration

Question 8

Why are state-scaled tolerances useful in a robot model?

A. They eliminate the need for initial conditions

B. They make every state dimensionless automatically

C. They account for different magnitudes and engineering units of states

D. They guarantee exact event timing

Correct Answer: C. They account for different magnitudes and engineering units of states

Question 9

At a passive one-dimensional impact, the restitution relation in the chapter is:

A. v^+ = e v^−

B. v^+ = −e v^−

C. v^+ = v^−/e

D. v^+ = −v^−/e

Correct Answer: B. v^+ = −e v^−

Question 10

Which observation provides evidence of numerical convergence?

A. The plot uses many colours

B. The result changes substantially when the step is halved

C. Results approach a limiting value under systematic step refinement

D. The solver completes without displaying a warning

Correct Answer: C. Results approach a limiting value under systematic step refinement