A MATLAB benchmark framework for studying, implementing, and comparing numerical integration methods for differential equations.
This repository started as an ODE integrator benchmark suite, but it has been reorganized into a broader numerical-solver catalog covering:
- classical and modern ODE initial-value integrators,
- stiff, implicit, multistep, Rosenbrock, exponential, IMEX, splitting, and geometric methods,
- nonlinear solvers and linear algebra machinery used inside implicit methods,
- DAE, PDE, SDE, DDE, fractional, BVP, hybrid, multirate, and parallel-in-time benchmark directions,
- MATLAB implementations and future implementation scaffolds organized by method and benchmark family.
The repository is intended for numerical-analysis education, control engineering, robotics simulation, computational physics, chemical kinetics, and scientific-computing benchmark development.
The repository contains two related layers:
-
Runnable MATLAB benchmark framework
Implemented methods and benchmarks can be executed throughrun_all.mandrun_one_benchmark.m. -
Expanded documentation and implementation roadmap
Many additional methods and benchmark problems are documented and organized in category folders. Some of these are planned/scaffolded and are not executed by default until their MATLAB implementations are completed.
By default, the runner excludes planned methods and planned benchmarks. This keeps the repository runnable while the catalog grows.
Open MATLAB in the repository root and run:
run_allRun a single benchmark by registry name:
run_one_benchmark('LinearTestEquation')
run_one_benchmark('LorenzSystem')
run_one_benchmark('RobertsonKinetics')
run_one_benchmark('KeplerTwoBody')Use options explicitly:
opts = default_options();
opts.h = 1e-2;
opts.write_plots = true;
opts.write_tables = true;
opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;
results = run_all(opts);Results are written under:
results/
or, depending on the benchmark runner, under benchmark-specific result folders such as:
results/<BenchmarkName>/
.
├── README.md
├── LICENSE
├── CONTRIBUTING.md
├── CITATION.cff
├── run_all.m
├── run_one_benchmark.m
├── docs/
│ ├── methods/
│ │ ├── 01_ODE_initial_value_integrators/
│ │ ├── 02_DAE_integrators_and_constraint_methods/
│ │ ├── 03_PDE_time_integration_and_spatial_discretization_methods/
│ │ ├── 04_stochastic_differential_equation_methods/
│ │ ├── 05_delay_differential_equation_methods/
│ │ ├── 06_fractional_differential_equation_methods/
│ │ ├── 07_boundary_value_problem_methods/
│ │ └── 08_multirate_and_parallel_in_time_methods/
│ ├── benchmarks/
│ │ ├── 01_ODE_nonstiff_accuracy_and_nonlinear_dynamics/
│ │ ├── 02_ODE_stiff_multiscale_and_chemical_kinetics/
│ │ ├── 03_Hamiltonian_geometric_and_long_time_mechanics/
│ │ ├── 04_DAE_constrained_systems_and_circuits/
│ │ ├── 05_PDE_method_of_lines_and_conservation_laws/
│ │ ├── 06_SDE_stochastic_differential_equations/
│ │ ├── 07_DDE_delay_differential_equations/
│ │ └── 08_Fractional_BVP_hybrid_and_event_driven_systems/
│ ├── metrics/
│ │ ├── README.md
│ │ ├── final_error.md
│ │ ├── max_error.md
│ │ ├── rmse_error.md
│ │ ├── cpu_time.md
│ │ ├── n_steps.md
│ │ ├── nfev.md
│ │ ├── n_rejected.md
│ │ ├── n_jacobian.md
│ │ ├── n_newton.md
│ │ ├── n_linear_solve.md
│ │ ├── max_invariant_error.md
│ │ ├── positivity_violation.md
│ │ ├── constraint_error.md
│ │ └── status.md
│ ├── tutorials/
│ └── references.md
├── src/
│ ├── config/
│ │ ├── default_options.m
│ │ ├── method_registry.m
│ │ └── benchmark_registry.m
│ ├── core/
│ │ ├── run_benchmark.m
│ │ ├── compute_metrics.m
│ │ ├── newton_solve.m
│ │ └── other shared execution utilities
│ ├── methods/
│ │ ├── 01_ODE_initial_value_integrators/
│ │ ├── 02_DAE_integrators_and_constraint_methods/
│ │ ├── 03_PDE_time_integration_and_spatial_discretization_methods/
│ │ ├── 04_stochastic_differential_equation_methods/
│ │ ├── 05_delay_differential_equation_methods/
│ │ ├── 06_fractional_differential_equation_methods/
│ │ ├── 07_boundary_value_problem_methods/
│ │ └── 08_multirate_and_parallel_in_time_methods/
│ ├── benchmarks/
│ │ ├── 01_ODE_nonstiff_accuracy_and_nonlinear_dynamics/
│ │ ├── 02_ODE_stiff_multiscale_and_chemical_kinetics/
│ │ ├── 03_Hamiltonian_geometric_and_long_time_mechanics/
│ │ ├── 04_DAE_constrained_systems_and_circuits/
│ │ ├── 05_PDE_method_of_lines_and_conservation_laws/
│ │ ├── 06_SDE_stochastic_differential_equations/
│ │ ├── 07_DDE_delay_differential_equations/
│ │ └── 08_Fractional_BVP_hybrid_and_event_driven_systems/
│ └── plotting/
├── tests/
│ └── test_smoke.m
└── results/
The docs/ and src/ folders intentionally mirror each other: method documentation categories correspond to method implementation categories, benchmark documentation categories correspond to benchmark implementation categories, and metric documentation explains the quantities reported by the framework.
The methods catalog is organized by numerical method family, not by implementation status. The left column uses human-readable category names, while the right column links directly to the available Markdown documentation files.
| Category | Method / Documentation |
|---|---|
| ODE initial value integrators | Explicit one-step and low-order RK: category README, Forward Euler, Explicit Midpoint, Heun, Ralston, Kutta RK3, Classical RK4 Adaptive, high-order, extrapolation, and stabilized RK: category README, RKF45, Dormand-Prince RK45, Bogacki-Shampine RK23, Cash-Karp RK45, DOP853, RKC, ROCK, SSPRK, Tsitouras 5(4) Implicit one-step, collocation, and DIRK: category README, Backward Euler, Implicit Midpoint, Trapezoidal Rule, DIRK, SDIRK, ESDIRK, Radau IIA 3, Radau IIA 5, Gauss-Legendre IRK2, Gauss-Legendre IRK4 Linear multistep, BDF/NDF, Adams, and predictor-corrector: category README, Adams-Bashforth 2, Adams-Bashforth 4, Adams-Moulton 2, ABM4 predictor-corrector, BDF1, BDF2, BDF3, CVODE Adams, CVODE BDF Rosenbrock-Wanner and linearly implicit stiff methods: category README, Rosenbrock-Euler, ROS2, ROS3, ROS4, RODAS3, RODAS4 Exponential and integrating-factor methods: category README, Exponential Euler, Exponential Midpoint, ETD2, ETDRK4, EPIRK IMEX, additive, and operator-splitting methods: category README, IMEX Euler, IMEX Runge-Kutta, ARK2, ARK3, ARK4, Lie-Trotter splitting Geometric, symplectic, and energy-preserving methods: category README, Symplectic Euler, Verlet, Velocity Verlet, Leapfrog, Yoshida 4, Strang splitting, Discrete Gradient Oscillator Second-order mechanical and structural dynamics methods: category README, Newmark-beta, Central difference, Generalized-alpha, HHT-alpha, Runge-Kutta-Nyström |
| DAE integrators and constraint methods | Category README, Implicit Euler DAE, BDF with mass matrix, DASSL, DASPK, IDA-style BDF, SHAKE, RATTLE, Baumgarte stabilization. |
| PDE time integration and spatial discretization methods | Category README, Method of lines, Finite difference method, Finite volume method, Finite element method, Fourier spectral method, ADI, Crank-Nicolson PDE, Lax-Friedrichs, Lax-Wendroff, WENO. |
| Stochastic differential equation methods | Category README for SDE and stochastic-reaction method documentation. |
| Delay differential equation methods | Category README for method-of-steps, history interpolation, neutral DDE, and state-dependent-delay documentation. |
| Fractional differential equation methods | Category README for fractional ODE/PDE, L1, Grünwald-Letnikov, convolution quadrature, and fast-memory documentation. |
| Boundary value problem methods | Category README for shooting, multiple shooting, collocation, finite-difference BVP, continuation, and pseudo-arclength documentation. |
| Multirate and parallel-in-time methods | Category README for multirate RK, MRI-GARK, waveform relaxation, Parareal, PFASST, and parallel-in-time documentation. |
The ODE folder contains additional subfolders for more precise method families:
01_explicit_one_step_and_low_order_RK
02_adaptive_high_order_extrapolation_and_stabilized_RK
03_implicit_one_step_collocation_and_DIRK
04_linear_multistep_BDF_NDF_Adams_and_predictor_corrector
05_Rosenbrock_Wanner_and_linearly_implicit_stiff_methods
06_exponential_and_integrating_factor_methods
07_IMEX_additive_and_operator_splitting_methods
08_geometric_symplectic_and_energy_preserving_methods
09_second_order_mechanical_and_structural_dynamics_methods
10_adaptivity_events_jacobians_and_mass_matrix_support
11_nonlinear_solvers_for_implicit_methods
12_linear_algebra_backsolves_and_preconditioners
The benchmark catalog is also organized by problem family. The table below uses human-readable category names and links directly to the available benchmark Markdown documentation files.
Registry entries distinguish between runnable implementations and planned/scaffolded entries.
- Implemented means the MATLAB runner can execute the method or benchmark by default.
- Planned means the documentation and file structure exist, but the method or benchmark should not be treated as a validated implementation yet.
- Solver-internal components such as Newton iterations, Jacobian strategies, LU factorizations, Krylov methods, and preconditioners are not always time integrators by themselves. They are included because they are essential for serious stiff, implicit, DAE, and PDE solver benchmarking.
To inspect the current registry from MATLAB:
methods = method_registry();
benchmarks = benchmark_registry();
implementedMethods = methods([methods.implemented]);
implementedBenchmarks = benchmarks([benchmarks.implemented]);
{implementedMethods.name}'
{implementedBenchmarks.name}'run_all.m performs the following steps:
- Loads default options from
default_options.m. - Adds the repository root and subfolders to the MATLAB path.
- Loads method and benchmark registries.
- Filters out planned methods and planned benchmarks by default.
- Runs applicable method/benchmark combinations.
- Writes tables and plots when enabled.
- Optionally writes a global summary table.
Default behavior:
opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;This protects normal runs from failing on future scaffold files.
| Type | Example form | Typical benchmark use |
|---|---|---|
| ODE | y' = f(t,y) |
control, mechanics, circuits, chemical kinetics, chaos |
| Stiff ODE | y' = f(t,y) with separated time scales |
implicit methods, BDF/NDF, Rosenbrock, Radau, SDIRK |
| Hamiltonian ODE | q' = ∂H/∂p, p' = -∂H/∂q |
symplectic methods, energy drift, long-time orbits |
| DAE | F(t,y,y') = 0 |
constrained mechanics, circuits, process systems |
| PDE after discretization | u_t = L(u) or M y' = f(y) |
method of lines, stiffness, conservation laws |
| SDE | dX = a(X)dt + b(X)dW |
stochastic stability, weak/strong convergence |
| DDE | x'(t)=f(t,x(t),x(t-τ)) |
delay systems, history interpolation |
| Fractional DE | Caputo or Grünwald-Letnikov derivatives | memory effects and fractional diffusion |
| BVP | y' = f(x,y), boundary constraints |
shooting, collocation, continuation |
| Hybrid/event-driven system | continuous flow plus jumps/switching | event detection and discontinuity handling |
The benchmark framework is designed to compare both numerical accuracy and solver mechanics. Each metric below has its own documentation file under docs/metrics/.
| Metric | Meaning | Documentation |
|---|---|---|
final_error |
Final-time error against exact or high-accuracy reference solution. | final_error |
max_error |
Maximum trajectory error over the comparison window. | max_error |
rmse_error |
Root-mean-square trajectory error. | rmse_error |
cpu_time |
MATLAB wall-clock execution time measured around one method/benchmark run. | cpu_time |
n_steps |
Accepted step count or output-step count, depending on the solver interface. | n_steps |
nfev |
Number of right-hand-side evaluations. | nfev |
n_rejected |
Rejected adaptive steps. | n_rejected |
n_jacobian |
Jacobian evaluations. | n_jacobian |
n_newton |
Newton or nonlinear iterations for implicit methods. | n_newton |
n_linear_solve |
Linear solves/backsolves used inside implicit or linearly implicit methods. | n_linear_solve |
max_invariant_error |
Maximum drift in mass, energy, angular momentum, or another invariant. | max_invariant_error |
positivity_violation |
Violation of nonnegative state constraints, when applicable. | positivity_violation |
constraint_error |
DAE or geometric constraint residual, when applicable. | constraint_error |
status |
Success, skipped, non-applicable, blow-up, or failure status. | status |
CPU time alone is not enough for solver comparison. For stiff and implicit methods, Jacobian evaluations, Newton iterations, rejected steps, linear solves, and factorization/backsolve counts are often more informative.
Depending on benchmark and method applicability, plotting utilities may generate:
- solution trajectories,
- error versus time,
- final-error bar charts,
- CPU-time bar charts,
- function-evaluation bar charts,
- work-precision diagrams,
- invariant drift plots,
- positivity or constraint-error plots,
- phase portraits,
- state-space projections,
- orbit plots,
- attractor projections.
There is no universal best integrator.
A meaningful comparison is problem-dependent:
- Explicit RK methods are often strong for smooth non-stiff ODEs.
- Adaptive RK methods are useful when automatic local-error control is needed.
- BDF/NDF, Radau, SDIRK, and Rosenbrock-type methods are important for stiff systems.
- Symplectic and variational methods are important for long-time Hamiltonian mechanics.
- Exponential and IMEX methods are important for split stiff/nonstiff systems and PDE semi-discretizations.
- DAE methods require constraint-residual and consistent-initialization diagnostics.
- SDE methods require strong/weak convergence and statistical diagnostics, not only deterministic trajectory error.
- PDE benchmarks require attention to both spatial and temporal discretization error.
To add a method consistently:
- Add or update the Markdown file under the correct folder in
docs/methods/. - Add the MATLAB implementation under the matching folder in
src/methods/. - Register the method in
src/config/method_registry.m. - Set
implemented = trueonly after the method is runnable and tested. - Define method applicability clearly, for example:
all_first_order,stiff_first_order,mechanical_separable,dae,sde,dde,pde_mol.
- Run smoke tests and at least one representative benchmark.
To add a benchmark consistently:
- Add or update the Markdown file under the correct folder in
docs/benchmarks/. - Add the MATLAB benchmark definition under the matching folder in
src/benchmarks/. - Register it in
src/config/benchmark_registry.m. - Define:
- right-hand side or residual,
- default initial condition,
- time interval,
- parameter regimes,
- exact/reference solution strategy,
- invariants or constraints,
- recommended plots,
- method applicability.
- Set
implemented = trueonly when the benchmark can run through the framework.
Run a quick smoke test:
runtests('tests')Run all currently implemented combinations:
run_allRun with custom options:
opts = default_options();
opts.h = 5e-3;
opts.n_output = 2001;
opts.write_plots = true;
opts.write_tables = true;
results = run_one_benchmark('VanDerPolOscillator', opts);Inspect registries:
methods = method_registry();
benchmarks = benchmark_registry();
sum([methods.implemented])
sum([benchmarks.implemented])Each method document should describe:
- category,
- intended MATLAB file,
- mathematical formulation,
- update formula or residual form,
- accuracy/order,
- stability behavior,
- solver requirements,
- strengths and weaknesses,
- best-use cases,
- metrics relevant to the method.
Each benchmark document should describe:
- governing equation,
- initial/boundary conditions,
- exact or reference solution if available,
- stiffness, invariants, constraints, positivity, or qualitative behavior,
- pseudocode for the right-hand side or residual,
- proposed/implemented MATLAB file,
- recommended metrics and plots.
Each metric document should describe:
- what the metric measures,
- the mathematical definition,
- when the metric is meaningful,
- how to interpret good and bad values,
- limitations and common pitfalls.
A simulation video for the repository is available on YouTube:
The code is written in plain MATLAB style and avoids Simulink dependencies. It should run on recent MATLAB releases. Some plotting details may require minor adjustment on very old MATLAB versions.
This project is released under the MIT License. See LICENSE.
If this repository helps your teaching, benchmarking, or research workflow, please cite it using the metadata in CITATION.cff.
The framework includes the following additions based on technical feedback about benchmark design, trajectory-error measurement, timing reliability, and solver-work accounting.
A runnable Swift-Hohenberg equation benchmark has been added. It uses:
- the one-dimensional periodic Swift-Hohenberg equation,
- 128 equispaced spatial points,
- Fourier pseudospectral differentiation,
- a fourth-order stiff spatial operator,
- a deterministic multimode initial condition,
- a high-accuracy stiff reference trajectory.
The spatial discretization produces a 128-dimensional stiff ODE system. This complements the existing Korteweg-de Vries and Kuramoto-Sivashinsky cases by testing scalability, spectral spatial discretization, stiffness handling, and nonlinear pattern-forming dynamics.
Run it with:
results = run_one_benchmark('SwiftHohenbergEquation');Final-time error can be misleading when an inaccurate numerical trajectory eventually approaches the correct equilibrium. The framework therefore reports errors over the full integration interval:
| Metric | Definition and purpose |
|---|---|
final_error |
Error at the final reported time. |
max_error |
Largest trajectory error over the integration interval. |
time_average_error |
Time integral of the pointwise error divided by total integration time. See time_average_error.md. |
rmse_error |
Time-weighted RMS trajectory error computed by trapezoidal integration over the actual output times. See rmse_error.md. |
Using time integration rather than an unweighted average prevents adaptive solvers from being favored or penalized merely because their output times are nonuniform.
A single elapsed-time measurement can vary because of MATLAB JIT compilation, memory allocation, operating-system scheduling, and background applications. Timing now supports:
- an optional untimed warm-up run,
- multiple measured repetitions,
- median elapsed time in
cpu_time, - minimum elapsed time in
cpu_time_min, - timing variability in
cpu_time_std, - the number of repetitions in
timing_repeats.
Default settings are:
opts.timing_warmup = true;
opts.timing_repeats = 3;The median is used as the primary timing statistic because it is less sensitive to occasional operating-system or background-process delays than a single measurement.
The framework can independently count calls to the benchmark right-hand side or mechanical acceleration callback. It uses a closure-based wrapper rather than a global variable.
The counting run is performed separately and is not included in the reported elapsed time, so instrumentation overhead does not distort cpu_time. Results include:
nfev, the measured number of callback evaluations,nfev_source, indicating whether the value came from the instrumented wrapper or solver-provided statistics.
The default setting is:
opts.count_rhs_evaluations = true;See nfev.md and cpu_time.md for interpretation guidance.
Benchmark summaries can now generate bar charts for:
- final-time error,
- time-averaged trajectory error,
- time-weighted RMS trajectory error,
- maximum trajectory error,
- right-hand-side evaluation count,
- median elapsed time with timing-variability error bars.
Warm-up, repeated timing, and the independent counting pass improve measurement quality but require additional solver executions. For a faster exploratory run, use:
opts = default_options();
opts.timing_warmup = false;
opts.timing_repeats = 1;
opts.count_rhs_evaluations = false;
results = run_one_benchmark('SwiftHohenbergEquation', opts);For a more stable timing comparison, increase the number of repetitions:
opts = default_options();
opts.timing_repeats = 7;
results = run_one_benchmark('LinearTestEquation', opts);The feedback-driven additions are covered by:
tests/test_feedback_improvements.m
Run all tests from the repository root with:
runtests('tests')