Numerical Integration of Differential Equations in MATLAB: Benchmarking Accuracy, Stability, Stiffness, and Conservation

Numerical integration is one of the central tools of scientific computing. Whenever a physical, biological, chemical, mechanical, electrical, or control system is modeled by differential equations, the next practical question is usually not only “what is the equation?” but also “which numerical method should be trusted to simulate it?”

This question is more subtle than it first appears. A numerical integrator that performs very well on a smooth non-stiff control example may fail completely on a stiff chemical kinetics model. A method that gives excellent short-time accuracy may slowly destroy the energy of a mechanical system. A method that is computationally cheap may require such a small time step that it becomes inefficient. A method that is stable for one equation may be unstable for another.

For this reason, numerical integrators should be studied through benchmark equations that expose different mathematical behaviors: linear stability, oscillation, nonlinear stiffness, chaos, conservation laws, and long-time orbital dynamics.

This article presents a comprehensive discussion of differential equation integrator benchmarking in MATLAB. The discussion is based on the project:

GitHub Repository:
https://github.com/mohammadijoo/DifferentialEquationIntegratorBenchmarks_MATLAB

The repository implements and compares a broad collection of numerical ODE integrators on classical benchmark problems. A companion YouTube simulation video demonstrates the MATLAB implementation, generated plots, MATLAB source codes, and exported CSV metric files.

Video: MATLAB implementation of differential equation integrator benchmarks, including repository code execution, generated plots, and CSV metric files.

---

1. Initial-Value Problems and the Role of Numerical Integrators

The main mathematical object in this benchmark framework is the initial-value problem:

$$\frac{dy}{dt} = f(t,y), \qquad y(t_0)=y_0$$

Here, the state vector belongs to an n-dimensional real space:

$$y(t)\in \mathbb{R}^n$$

The function f(t,y) is the right-hand-side function, and y_0 is the initial condition.

A numerical integrator approximates the exact solution at discrete time points:

$$t_n = t_0 + nh$$

where h is the time step. The numerical approximation is denoted by:

$$y_n \approx y(t_n)$$

The purpose of a benchmark is to compare different update rules of the form:

$$y_{n+1} = \Phi_h(t_n,y_n)$$

where the map Phi_h is the numerical flow map produced by a specific method.

A good benchmark should not only ask whether y_n is close to the exact solution. It should also ask:

• Is the method stable?
• Does it preserve physical invariants?
• Does it handle stiffness?
• Does it remain efficient?
• Does it fail, blow up, or produce NaN values?
• Does it preserve qualitative dynamics such as oscillations, attractors, or orbits?

These questions are especially important in MATLAB, where numerical simulation is frequently used in control engineering, robotics, mechanical dynamics, electrical circuits, chemical kinetics, and computational physics.

---

2. Benchmarking Is Not Just About Error

The most direct error metric is the final-time error:

$$e_{\text{final}} = |y_N - y_{\text{ref}}(T)|$$

where y_ref(T) is either the exact solution or a high-accuracy reference solution at the final time T.

However, this metric alone can be misleading. A method may have small final error but large intermediate error. Therefore, a trajectory-level maximum error is also useful:

$$e_{\max} = \max_{0\le n\le N} |y_n - y_{\text{ref}}(t_n)|$$

A root-mean-square error gives a global average measure:

$$e_{\mathrm{RMSE}} \sqrt{ \frac{1}{N+1} \sum_{n=0}^{N} |y_n-y_{\text{ref}}(t_n)|^2 }$$

For computational efficiency, the benchmark should also record CPU time:

$$T_{\mathrm{CPU}}$$

number of right-hand-side evaluations:

$$N_{\mathrm{fev}}$$

number of accepted steps:

$$N_{\mathrm{steps}}$$

number of rejected adaptive steps:

$$N_{\mathrm{rejected}}$$

number of Newton iterations for implicit methods:

$$N_{\mathrm{Newton}}$$

and number of Jacobian evaluations:

$$N_{\mathrm{Jacobian}}$$

For conservative systems, an additional invariant error is critical. If a system has an invariant such as energy, mass, or angular momentum, then the invariant drift can be measured as:

$$e_I(t_n) = |I(y_n)-I(y_0)|$$

The maximum invariant error is:

$$e_{I,\max} \max_{0\le n\le N} |I(y_n)-I(y_0)|$$

This metric is especially important for oscillators, Hamiltonian systems, celestial mechanics, and long-time mechanical simulation.

---

3. Benchmark Equation 1: Linear Decay / Dahlquist Test Equation

The simplest but most important benchmark is the scalar linear test equation:

$$x'(t) = \lambda x(t), \qquad x(0)=x_0$$

Its exact solution is:

$$x(t) = x_0 e^{\lambda t}$$

For negative real parts of lambda, the exact solution decays to zero. A numerical method should reproduce that decay without artificial blow-up.

The key dimensionless quantity is:

$$z = h\lambda$$

For a one-step method applied to the test equation, the update usually takes the form:

$$x_{n+1} = R(z)x_n$$

where R(z) is the stability function of the method. The numerical solution is stable when:

$$|R(z)| \le 1$$

This benchmark is essential because it reveals the absolute stability region of a method.

Forward Euler on the Dahlquist Equation

Forward Euler is:

$$x_{n+1} = x_n + h\lambda x_n$$

Therefore:

$$x_{n+1} = (1+h\lambda)x_n$$

So the stability function is:

$$R(z)=1+z$$

The method is stable only when:

$$|1+z|\le 1$$

For real negative lambda, this implies:

$$-2 \le h\lambda \le 0$$

Thus, even for a simple decaying physical mode, Forward Euler can become unstable if the step size is too large.

Backward Euler on the Dahlquist Equation

Backward Euler is:

$$x_{n+1}=x_n+h\lambda x_{n+1}$$

Solving for the next state gives:

$$x_{n+1} \frac{1}{1-h\lambda}x_n$$

The stability function is:

$$R(z)=\frac{1}{1-z}$$

For negative real parts of z, Backward Euler is stable. This is why Backward Euler is a fundamental stiff integrator: it damps fast stable modes without requiring extremely small time steps.

This benchmark therefore separates conditionally stable explicit methods from stable implicit methods.

---

4. Benchmark Equation 2: Harmonic Oscillator

The harmonic oscillator is:

$$q' = v$$

$$v' = -\omega^2 q$$

Equivalently:

$$q''+\omega^2 q=0$$

The exact solution is periodic:

$$q(t)=A\cos(\omega t)+B\sin(\omega t)$$

The conserved energy is:

$$E(q,v) \frac{1}{2}v^2 + \frac{1}{2}\omega^2 q^2$$

For the exact solution:

$$E(t)=E(0)$$

This benchmark is not only about pointwise accuracy. It is also about phase error and energy drift.

The exact phase-space trajectory satisfies:

$$\frac{v^2}{2E} + \frac{\omega^2 q^2}{2E} 1$$

This is an ellipse in the phase plane. A dissipative method causes the ellipse to spiral inward. An unstable method causes it to spiral outward. A good geometric method should keep the trajectory near the correct energy level.

---

5. Benchmark Equation 3: Van der Pol Oscillator

The Van der Pol oscillator is a nonlinear second-order oscillator. In first-order form:

$$x' = y$$

$$y' = \mu(1-x^2)y - x$$

The parameter mu controls the nonlinear damping.

For small mu, the system behaves like a smooth nonlinear oscillator. For large mu, it becomes a fast-slow system with relaxation oscillations. This introduces stiffness.

The stiffness appears because the solution contains both slow segments and fast transitions. A non-stiff explicit method may need extremely small step sizes to resolve the fast layers.

The Jacobian of the vector field is:

$$J(x,y) \begin{bmatrix} 0 & 1 \ -2\mu xy - 1 & \mu(1-x^2) \end{bmatrix}$$

When mu is large, eigenvalues of this Jacobian can have very different magnitudes. This makes the problem a useful test for stiffness handling.

The Van der Pol benchmark asks:

• Can explicit methods remain stable?
• Do adaptive methods reduce the step size during fast transitions?
• Do stiff methods handle the relaxation regime efficiently?
• How does CPU time scale as mu increases?
• Does the numerical trajectory preserve the limit-cycle structure?

---

6. Benchmark Equation 4: Lorenz System

The Lorenz system is:

$$x' = \sigma(y-x)$$

$$y' = x(\rho-z)-y$$

$$z' = xy-\beta z$$

For the classical parameters:

$$\sigma=10, \qquad \rho=28, \qquad \beta=\frac{8}{3}$$

the system exhibits chaotic behavior.

The Lorenz system is a warning against relying only on long-time pointwise error. Because the system is chaotic, nearby trajectories separate exponentially:

$$|\delta y(t)| \approx |\delta y(0)|e^{\lambda_L t}$$

where lambda_L is a positive Lyapunov exponent.

This means that even a very accurate numerical method will eventually deviate from the reference trajectory. Therefore, for chaotic systems, the benchmark should emphasize:

• short-time trajectory accuracy,
• qualitative attractor geometry,
• sensitivity to time step,
• stability of the computed trajectory,
• preservation of the butterfly-shaped attractor.

The Lorenz system is not a conservation-law benchmark. It is a qualitative-dynamics benchmark. It tests whether the method gives a credible numerical representation of chaotic flow over a meaningful time window.

---

7. Benchmark Equation 5: Robertson Stiff Chemical Kinetics

The Robertson problem is a classic stiff chemical kinetics benchmark. A standard form is:

$$y_1' = -0.04y_1 + 10^4 y_2y_3$$

$$y_2' = 0.04y_1 - 10^4 y_2y_3 - 3\times 10^7 y_2^2$$

$$y_3' = 3\times 10^7 y_2^2$$

The usual initial condition is:

$$y_1(0)=1, \qquad y_2(0)=0, \qquad y_3(0)=0$$

The system conserves total mass:

$$y_1(t)+y_2(t)+y_3(t)=1$$

Therefore, an important invariant error is:

$$e_m(t_n) |y_{1,n}+y_{2,n}+y_{3,n}-1|$$

The Robertson problem is difficult because the reaction rates differ by many orders of magnitude. The term below creates a very fast transient:

$$3\times 10^7 y_2^2$$

This is exactly the type of system where explicit methods become inefficient or unstable. Stiff solvers, implicit methods, BDF methods, and Rosenbrock-type methods are much more appropriate.

The benchmark should also check positivity:

$$y_i(t)\ge 0, \qquad i=1,2,3$$

A numerical method that produces negative chemical concentrations may be mathematically inaccurate and physically meaningless.

---

8. Benchmark Equation 6: Kepler Two-Body Problem

The Kepler two-body problem describes motion under an inverse-square central force:

$$r' = v$$

$$v' = -\mu \frac{r}{|r|^3}$$

Here:

$$r= \begin{bmatrix} x\ y \end{bmatrix}, \qquad v= \begin{bmatrix} v_x\ v_y \end{bmatrix}$$

The total energy is:

$$E(r,v) \frac{1}{2}|v|^2 \frac{\mu}{|r|}$$

The angular momentum in two dimensions is:

$$L = x v_y - y v_x$$

For the exact Kepler problem:

$$E(t)=E(0)$$

$$L(t)=L(0)$$

This benchmark is a strong test of long-time geometric behavior. A method may have good local accuracy but slowly distort the orbit. The orbit may precess, expand, shrink, or collapse depending on the numerical method.

For orbital mechanics and robotics simulation, this distinction is important. Long-time qualitative fidelity is often more important than small short-time pointwise error.

Symplectic methods, splitting methods, and energy-aware methods are especially meaningful for this benchmark.

---

9. Numerical Integrator Families

The benchmark framework compares many families of methods. Each family represents a different philosophy of numerical integration.

---

10. Explicit One-Step Methods

Explicit methods compute the next numerical state directly from already known values.

Forward Euler

Forward Euler is the simplest method:

$$y_{n+1}=y_n+h f(t_n,y_n)$$

It is first order:

$$e_{\text{global}} = O(h)$$

Its advantages are simplicity and low computational cost. Its disadvantages are poor stability and low accuracy.

Forward Euler is useful as a baseline, not usually as a final method for serious simulation.

---

Explicit Midpoint Method

The explicit midpoint method uses a half-step slope:

$$k_1=f(t_n,y_n)$$

$$k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right)$$

$$y_{n+1}=y_n+h k_2$$

It is second order:

$$e_{\text{global}}=O(h^2)$$

Compared with Forward Euler, it improves accuracy by sampling the vector field inside the interval.

---

Heun’s Method

Heun’s method, also called the improved Euler method, averages the slope at the beginning and predicted end of the step:

$$k_1=f(t_n,y_n)$$

$$k_2=f(t_n+h,y_n+hk_1)$$

$$y_{n+1}=y_n+\frac{h}{2}(k_1+k_2)$$

This method is also second order. It is a predictor-corrector style explicit Runge-Kutta method.

---

Ralston’s Method

Ralston’s method chooses coefficients to reduce the local truncation error:

$$k_1=f(t_n,y_n)$$

$$k_2=f\left(t_n+\frac{2h}{3},y_n+\frac{2h}{3}k_1\right)$$

$$y_{n+1}=y_n+h\left(\frac{1}{4}k_1+\frac{3}{4}k_2\right)$$

It is another second-order method, but with optimized weighting.

---

11. Classical Runge-Kutta Methods

Runge-Kutta methods improve accuracy by combining several internal stages.

Classical RK4

The classical fourth-order Runge-Kutta method is:

$$k_1=f(t_n,y_n)$$

$$k_2=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_1\right)$$

$$k_3=f\left(t_n+\frac{h}{2},y_n+\frac{h}{2}k_2\right)$$

$$k_4=f(t_n+h,y_n+hk_3)$$

$$y_{n+1} y_n + \frac{h}{6} (k_1+2k_2+2k_3+k_4)$$

RK4 has global error:

$$e_{\text{global}}=O(h^4)$$

It is one of the most important fixed-step methods in engineering simulation. It performs very well on smooth non-stiff problems, but it is not a stiff solver and it is not symplectic.

This distinction matters. RK4 may be excellent for short-time non-stiff control examples but less appropriate for very stiff chemical kinetics or long-time Hamiltonian mechanics.

---

12. Adaptive Embedded Runge-Kutta Methods

Adaptive methods estimate their own local error and adjust the step size.

An embedded Runge-Kutta method computes two approximations of different orders:

$$y_{n+1}^{(p)}$$

and:

$$y_{n+1}^{(p-1)}$$

The error estimate is:

$$\epsilon_n \left| y_{n+1}^{(p)} y_{n+1}^{(p-1)} \right|$$

The method accepts the step if:

$$\epsilon_n \le \mathrm{tol}$$

A common step-size update rule is:

$$h_{\text{new}} s h \left( \frac{\mathrm{tol}}{\epsilon_n} \right)^{1/(p+1)}$$

where the safety factor satisfies:

$$0<s<1$$

RKF45 and Dormand-Prince RK45

RKF45 and Dormand-Prince RK45 are adaptive embedded Runge-Kutta methods. They are especially useful when the solution has regions of different difficulty.

For example:

• In the Van der Pol oscillator, small steps are needed near fast transitions.
• In smooth regions, larger steps are acceptable.
• In non-stiff engineering simulations, adaptive RK methods often give excellent error-per-cost performance.

Dormand-Prince RK45 is closely related to the style of algorithm used in many ode45-type solvers.

However, adaptivity alone does not solve stiffness. A non-stiff adaptive explicit method may reduce the step size repeatedly and still become inefficient on very stiff problems.

---

13. Implicit One-Step Methods

Implicit methods define the next numerical state through an equation involving the unknown future state.

Backward Euler

Backward Euler is:

$$y_{n+1} y_n + h f(t_{n+1},y_{n+1})$$

This requires solving the nonlinear equation:

$$G(y_{n+1}) y_{n+1} y_n h f(t_{n+1},y_{n+1}) 0$$

Newton’s method can be used:

$$y_{n+1}^{(m+1)} y_{n+1}^{(m)} \left[ I hJ_f(t_{n+1},y_{n+1}^{(m)}) \right]^{-1} G(y_{n+1}^{(m)})$$

where:

$$J_f = \frac{\partial f}{\partial y}$$

Backward Euler is first order but strongly stable for stiff decay. It is dissipative, which is good for stiff chemical transients but bad for conservative mechanics if energy preservation is desired.

---

Implicit Midpoint Method

The implicit midpoint method is:

$$y_{n+1} y_n + h f\left( t_n+\frac{h}{2}, \frac{y_n+y_{n+1}}{2} \right)$$

It is second order. For Hamiltonian systems, it has favorable geometric properties and often behaves much better than purely dissipative methods.

---

Trapezoidal Rule

The trapezoidal rule is:

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

It is second order and implicit. It often gives good accuracy but may show oscillatory behavior for some stiff problems because it is not as strongly damping as Backward Euler.

---

14. Linear Multistep Methods

Linear multistep methods reuse previous solution values and previous function evaluations.

A general form is:

$$\sum_{j=0}^{k}\alpha_j y_{n+j} h \sum_{j=0}^{k}\beta_j f_{n+j}$$

If the final beta coefficient is zero, the method is explicit. If it is nonzero, the method is implicit.

---

Adams-Bashforth Methods

The two-step Adams-Bashforth method is:

$$y_{n+1} y_n + h \left( \frac{3}{2}f_n \frac{1}{2}f_{n-1} \right)$$

The four-step Adams-Bashforth method is:

$$y_{n+1} y_n + \frac{h}{24} \left( 55f_n 59f_{n-1} + 37f_{n-2} 9f_{n-3} \right)$$

These methods can be efficient because they reuse previous slopes. However, they require startup values and can have stability limitations.

---

Adams-Moulton and Predictor-Corrector Methods

A second-order Adams-Moulton method can be written as:

$$y_{n+1} y_n + \frac{h}{2} (f_{n+1}+f_n)$$

Because f_{n+1} depends on y_{n+1}, this is implicit.

A predictor-corrector method first predicts the next state using an explicit formula, then corrects it using an implicit formula.

For example, the predictor can be:

$$y_{n+1}^{P} y_n + h \left( \frac{3}{2}f_n \frac{1}{2}f_{n-1} \right)$$

The corrector can be:

$$y_{n+1}^{C} y_n + \frac{h}{2} \left[ f(t_{n+1},y_{n+1}^{P}) + f_n \right]$$

ABM4 predictor-corrector methods extend this idea to higher order.

---

BDF2

The second-order backward differentiation formula is:

$$\frac{3y_{n+1}-4y_n+y_{n-1}}{2h} f(t_{n+1},y_{n+1})$$

Equivalently:

$$3y_{n+1} 4y_n + y_{n-1} 2h f(t_{n+1},y_{n+1}) 0$$

BDF methods are important in stiff integration. They are widely used in chemical kinetics, process systems, and differential-algebraic equation solvers.

---

15. Rosenbrock and Linearly Implicit Methods

Fully implicit methods may require nonlinear Newton iterations. Rosenbrock methods reduce this cost by solving linear systems involving the Jacobian.

A simple Rosenbrock-Euler step has the form:

$$\left(I-h\gamma J_n\right)k_1 f(t_n,y_n)$$

$$y_{n+1}=y_n+h k_1$$

where:

$$J_n=\frac{\partial f}{\partial y}(t_n,y_n)$$

Instead of solving a nonlinear equation, the method solves a linear system. This makes Rosenbrock methods attractive for stiff systems.

For stiff equations such as Robertson kinetics, linearly implicit methods can offer a practical compromise between stability and cost.

---

16. Exponential Euler

For systems where the local linear part is important, exponential integrators use the matrix exponential.

Suppose:

$$y'=Ay+g(t,y)$$

The exponential Euler method can be written as:

$$y_{n+1} e^{hA}y_n + h\varphi_1(hA)g(t_n,y_n)$$

where:

$$\varphi_1(z)=\frac{e^z-1}{z}$$

Exponential methods can be effective when the linear dynamics dominate the behavior, especially when the matrix A contains stiff linear modes.

---

17. Symplectic and Geometric Methods

Mechanical systems often have Hamiltonian structure. For a separable Hamiltonian:

$$H(q,p)=T(p)+V(q)$$

the equations are:

$$q' = \frac{\partial H}{\partial p}$$

$$p' = -\frac{\partial H}{\partial q}$$

A symplectic method preserves the geometric structure of phase space. It may not conserve energy exactly at every step, but it often keeps energy bounded over long time intervals.

This is why symplectic methods are valuable for oscillators, orbital mechanics, molecular dynamics, and robotics.

---

Symplectic Euler

One form of symplectic Euler is:

$$p_{n+1} p_n h\nabla V(q_n)$$

$$q_{n+1} q_n + h\nabla T(p_{n+1})$$

For quadratic kinetic energy:

$$T(p)=\frac{1}{2}p^T M^{-1}p$$

the position update becomes:

$$q_{n+1} q_n + hM^{-1}p_{n+1}$$

This method is simple but has much better long-time behavior than Forward Euler on Hamiltonian systems.

---

Verlet Method

For a second-order system:

$$q''=a(q)$$

the Verlet method is:

$$q_{n+1} 2q_n q_{n-1} + h^2 a(q_n)$$

It is second order and symplectic for appropriate mechanical systems.

---

Velocity Verlet

Velocity Verlet updates position and velocity at the same time grid:

$$q_{n+1} q_n + h v_n + \frac{h^2}{2}a(q_n)$$

$$v_{n+1} v_n + \frac{h}{2} \left[ a(q_n)+a(q_{n+1}) \right]$$

This method is especially useful for particle dynamics, pendulum-like systems, and Kepler motion.

---

Leapfrog

Leapfrog staggers position and velocity:

$$v_{n+\frac{1}{2}} v_{n-\frac{1}{2}} + h a(q_n)$$

$$q_{n+1} q_n + h v_{n+\frac{1}{2}}$$

The name comes from the fact that position and velocity updates “leap over” each other.

---

Yoshida Fourth-Order Method

Higher-order symplectic methods can be built by composing lower-order symplectic maps.

If S_h is a second-order symplectic step, a fourth-order Yoshida composition has the form:

$$S_h^{(4)} S_{w_1h} S_{w_0h} S_{w_1h}$$

with coefficients:

$$w_1 \frac{1}{2-2^{1/3}}$$

$$w_0 -\frac{2^{1/3}}{2-2^{1/3}}$$

This improves accuracy while preserving the symplectic structure.

---

18. Splitting Methods

Splitting methods decompose the vector field into parts that are easier to integrate.

Suppose:

$$y' = (A+B)y$$

If the flows of A and B can be computed separately, Strang splitting uses:

$$\Phi_h \approx \Phi_{A,h/2} \circ \Phi_{B,h} \circ \Phi_{A,h/2}$$

This gives second-order accuracy:

$$\Phi_h \Phi_{A,h/2} \Phi_{B,h} \Phi_{A,h/2} + O(h^3)$$

Splitting methods are useful in Hamiltonian dynamics, reaction-diffusion systems, and problems where different physical effects can be integrated separately.

---

19. Newmark-Beta Method

The Newmark-beta method is widely used in structural dynamics. For:

$$M q'' + C q' + Kq = F(t)$$

the method updates displacement and velocity using parameters beta and gamma:

$$q_{n+1} q_n + h v_n + h^2 \left[ \left(\frac{1}{2}-\beta\right)a_n + \beta a_{n+1} \right]$$

$$v_{n+1} v_n + h \left[ (1-\gamma)a_n + \gamma a_{n+1} \right]$$

A common average-acceleration choice is:

$$\gamma=\frac{1}{2}, \qquad \beta=\frac{1}{4}$$

This method is important in mechanical engineering, vibration analysis, and finite element structural simulation.

---

20. Energy-Oriented and Discrete Gradient Methods

Energy-preserving methods are designed to control invariant drift.

For a system with energy:

$$E(y)$$

a discrete gradient method uses a discrete gradient satisfying:

$$E(y_{n+1})-E(y_n) \bar{\nabla}E(y_n,y_{n+1})^T (y_{n+1}-y_n)$$

If the update is designed so that:

$$y_{n+1}-y_n h S \bar{\nabla}E(y_n,y_{n+1})$$

where the matrix S is skew-symmetric:

$$S^T=-S$$

then:

$$E(y_{n+1})-E(y_n) h \bar{\nabla}E^T S \bar{\nabla}E 0$$

The equality follows because for any vector a:

$$a^T S a = 0$$

when S is skew-symmetric.

This is the key idea behind energy-preserving integration. It is especially useful when invariant preservation is more important than exact pointwise tracking.

---

21. Work-Precision Analysis

A work-precision diagram compares computational cost against numerical error.

The horizontal axis often represents CPU time or number of function evaluations:

$$\text{work} T_{\mathrm{CPU}} \quad \text{or} \quad N_{\mathrm{fev}}$$

The vertical axis represents error:

$$\text{precision} e_{\text{final}}, \quad e_{\max}, \quad \text{or} \quad e_{\mathrm{RMSE}}$$

A method is attractive if it lies toward the lower-left region of the plot:

$$\text{low error} + \text{low cost}$$

However, interpretation depends on the problem.

For a smooth non-stiff problem, RK4 or adaptive RK methods may dominate. For a stiff chemical kinetics problem, an implicit or Rosenbrock method may dominate. For a long-time orbit problem, a symplectic method may have larger pointwise error but much better invariant behavior.

This is why work-precision diagrams should be interpreted together with invariant plots, phase portraits, failure status, and qualitative solution plots.

---

22. Why No Integrator Is Universally Best

The central lesson of integrator benchmarking is:

$$\text{best method} \ne \text{same method for every equation}$$

A method should be selected according to the mathematical structure of the problem.

For a smooth non-stiff problem, RK4, RKF45, or Dormand-Prince RK45 may be excellent choices.

For stiff decay or stiff chemical kinetics, Backward Euler, BDF2, Rosenbrock, or stiff MATLAB solvers are more appropriate.

For conservative mechanical systems, symplectic, splitting, or energy-oriented methods are often better for long-time simulation.

For chaotic systems, short-time accuracy and qualitative attractor behavior are more meaningful than long-time pointwise agreement.

For adaptive simulations, step-size control is essential, but it should not be confused with stiffness handling.

---

23. MATLAB as a Benchmarking Environment

MATLAB is a strong environment for this type of benchmark because it provides:

• matrix and vector operations,
• fast plotting,
• readable implementation of algorithms,
• built-in timing tools such as tic and toc,
• easy CSV export,
• numerical linear algebra routines,
• well-known ODE solvers for comparison.

A benchmark repository can therefore produce both visual and quantitative outputs.

Typical outputs include:

• solution trajectory plots,
• error-versus-time plots,
• final error bar charts,
• CPU-time bar charts,
• function-evaluation bar charts,
• work-precision diagrams,
• invariant-error plots,
• phase portraits,
• CSV metric files.

The CSV files are especially important because they allow the benchmark to be inspected beyond visual plots. They can be opened in MATLAB, Python, Excel, R, or any data-analysis tool.

---

24. Educational Value

This benchmark framework is useful for several audiences.

For numerical analysis students, it demonstrates order, stability, stiffness, and convergence.

For control engineering students, it shows how simulation quality affects closed-loop and open-loop dynamic analysis.

For robotics students, it explains why mechanical simulation should consider energy and geometric structure.

For computational physics students, it connects numerical methods to conservation laws and phase-space behavior.

For MATLAB users, it provides readable implementations of many important ODE integrators.

For researchers, it gives a reproducible framework for comparing methods across different differential equation types.

---

25. Practical Interpretation of the Six Benchmarks

The six benchmark equations can be interpreted as a structured test suite.

The Dahlquist test equation asks whether the method has acceptable stability behavior.

The harmonic oscillator asks whether the method controls phase error and energy drift.

The Van der Pol oscillator asks whether the method can handle nonlinear oscillation and stiffness transition.

The Lorenz system asks whether the method can represent chaotic dynamics credibly over a useful time window.

The Robertson problem asks whether the method can survive severe stiffness while preserving positivity and mass.

The Kepler problem asks whether the method can preserve orbital geometry, energy, and angular momentum over long times.

Together, these problems form a balanced benchmark suite. They are not redundant. Each reveals a different failure mode.

---

26. Conclusion

Numerical integration is not only a matter of implementing formulas. It is a matter of matching the numerical method to the mathematical structure of the differential equation.

The same method can be excellent on one benchmark and poor on another. RK4 may be an excellent general-purpose baseline for smooth non-stiff systems, but it is not a stiff solver and not an energy-preserving geometric method. Backward Euler may handle stiff decay robustly, but it can overdamp conservative systems. Symplectic methods may preserve long-time mechanical structure, but they are not designed for arbitrary stiff chemical kinetics. Adaptive RK methods are powerful, but adaptivity does not automatically solve stiffness.

A good benchmark must therefore measure accuracy, cost, stability, failure behavior, and invariant drift. It must include equations with different mathematical character: linear decay, oscillation, nonlinear stiffness, chaos, chemical kinetics, and orbital mechanics.

The MATLAB benchmark project DifferentialEquationIntegratorBenchmarks_MATLAB provides a practical environment for studying these issues. It generates plots and CSV metric files, and the accompanying YouTube simulation video shows how the MATLAB codes run and produce visual comparisons and numerical performance outputs.

The main lesson is:

$$\boxed{ \text{There is no universal best integrator; there is only the best integrator for a given problem structure.} }$$