DEV Community: Abolfazl Mohammadijoo

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

Abolfazl Mohammadijoo — Mon, 15 Jun 2026 01:56:26 +0000

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,
chemical reaction stiffness,
conservation laws,
long-time Hamiltonian behavior,
positivity constraints,
PDE stiffness after spatial discretization,
and failure modes such as blow-up or non-applicability.

This article is based on my MATLAB repository:

GitHub Repository:

mohammadijoo / Differential_Equation_Integrator_Benchmarks_MATLAB

MATLAB benchmark suite for comparing numerical integrators on classical differential-equation problems, including linear decay, harmonic oscillator, Van der Pol, Lorenz, Robertson kinetics, and Kepler two-body dynamics. The repository includes performance plots, accuracy metrics, conservation analysis, and a simulation video available on YouTube.

Differential Equation Integrator Benchmarks in MATLAB

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.

Current status

The repository contains two related layers:

Runnable MATLAB benchmark framework
Implemented methods and benchmarks can be executed through run_all.m and run_one_benchmark.m.
Expanded documentation and implementation roadmap
Many additional methods and benchmark problems are documented and organized in…

View on GitHub

The repository started as a MATLAB ODE integrator benchmark suite, but it has now been reorganized into a broader numerical-solver catalog. It includes runnable MATLAB benchmark code, implemented integrators, benchmark documentation, generated results, and a much larger roadmap for future DAE, PDE, SDE, DDE, fractional, BVP, hybrid, multirate, and parallel-in-time benchmark directions.

1. The Main Mathematical Object: Initial-Value Problems

Most of the currently runnable benchmark framework is based on initial-value problems of the form:

\frac{d y}{d t} = f (t, y), y (t_{0}) = y_{0}

where

y (t) \in R^{n}

is the state vector, 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

and produces numerical states

y_{n} \approx y (t_{n})

A one-step numerical method can be written abstractly as

y_{n + 1} = Φ_{h} (t_{n}, y_{n})

where Phi_h is the numerical flow map generated by the method.

A 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 preserve positivity?
Does it require many function evaluations?
Does it require many Newton or linear-solve iterations?
Does it fail, blow up, or return non-applicable?
Does it preserve qualitative behavior such as oscillations, orbits, phase portraits, or attractors?

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

2. What Changed in the Repository

The original version of the project focused mainly on a smaller set of classical ODE benchmarks and common integrators.

The current repository is broader. It now has two related layers.

Layer 1: Runnable MATLAB Benchmark Framework

This is the part that can be executed directly from MATLAB.

The main entry points are:

run_all

and

run_one_benchmark('LinearTestEquation')
run_one_benchmark('LorenzSystem')
run_one_benchmark('RobertsonKinetics')
run_one_benchmark('KeplerTwoBody')

The default options are obtained from:

opts = default_options();

For example:

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);

By default, planned methods and planned benchmarks are excluded:

opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;

This is important. The repository contains many documented and scaffolded future methods, but the default runner only executes methods and benchmarks marked as implemented.

Layer 2: Expanded Documentation and Roadmap

The repository has also become a structured catalog of numerical methods and benchmark problems.

The documentation and source folders are organized by method and benchmark families. The method catalog includes:

classical ODE initial-value integrators,
stiff, implicit, multistep, Rosenbrock, exponential, IMEX, splitting, and geometric methods,
DAE integrators and constraint methods,
PDE time integration and spatial discretization methods,
stochastic differential equation methods,
delay differential equation methods,
fractional differential equation methods,
boundary-value problem methods,
multirate and parallel-in-time methods.

This means the repository is not only a runnable MATLAB benchmark suite, but also an educational numerical-analysis roadmap.

3. Benchmarking Is Not Just About Final-Time Error

The most direct metric is the final-time error:

e_{final} ∣ y_{N} - y_{ref} (T) ∣_{2}

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

However, final-time error alone can be misleading. A method may accidentally have small error at the final time but large error during the trajectory. Therefore, a trajectory-level maximum error is also useful:

e_{m a x} 0 \leq n \leq N max ∣ y_{n} - y_{ref} (t_{n}) ∣_{2}

A root-mean-square trajectory error gives an average global measure:

e_{RMSE} \frac{1}{N + 1} n = 0 \sum N ∣ y_{n} - y_{ref} (t_{n}) ∣_{22}

For computational cost, the benchmark records or emphasizes quantities such as:

CPU time,
number of accepted steps,
number of right-hand-side evaluations,
number of rejected adaptive steps,
number of Jacobian evaluations,
number of Newton iterations,
number of linear solves,
method status.

For implicit and stiff methods, CPU time alone is not enough. A solver with more steps may still be fast if each step is cheap. A solver with fewer steps may be expensive if each step requires Jacobians, Newton iterations, and matrix factorizations.

For conservative or mechanical systems, invariant drift is also essential. If a system has an invariant I(y), such as energy, mass, or angular momentum, then the invariant error can be measured as:

e_{I} (t_{n}) ∣ I (y_{n}) - I (y_{0}) ∣

and the maximum invariant drift is

e_{I, m a x} 0 \leq n \leq N max ∣ I (y_{n}) - I (y_{0}) ∣

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

4. Stability and the Dahlquist Test Equation

The linear scalar test equation is one of the simplest ways to study basic stability behavior:

y^{'} = λ y, y (0) = y_{0}

The exact solution is

y (t) = y_{0} e^{λ t}

If lambda < 0, the exact solution decays. A numerical method should reproduce that decay without artificial growth.

For a one-step method applied to the linear test equation, the numerical update can often be written as

y_{n + 1} = R (z) y_{n}, z = hλ

where R(z) is the stability function of the method.

The method is stable for a given z if

∣ R (z) ∣ \leq 1

This benchmark is useful because it exposes the absolute-stability region of a method. For example, Forward Euler is simple, but its stability region is limited. Backward Euler is much more robust for decay problems and is often useful for stiff systems, but it can strongly damp dynamics.

5. The Original Six Core Benchmarks

The early version of the article focused on six major ODE benchmarks. These are still important and remain part of the project.

5.1 Linear Test Equation

Purpose:

stability behavior,
decay accuracy,
basic convergence,
step-size sensitivity.

Equation:

y^{'} = λ y

This benchmark is small, but it is not trivial. It is the standard first test for understanding whether a method behaves correctly on exponential growth or decay.

5.2 Harmonic Oscillator

The harmonic oscillator is a basic conservative mechanical system:

q^{'} = v, v^{'} = - ω^{2} q

Its energy is

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

For the exact solution,

E (t) = E (0)

This benchmark tests:

phase accuracy,
amplitude error,
energy drift,
long-time oscillatory behavior,
suitability of symplectic or energy-preserving methods.

A method can have small local error but still slowly increase or decrease the oscillator energy. That is why geometric diagnostics are important.

5.3 Van der Pol Oscillator

The Van der Pol oscillator is a nonlinear oscillator:

x^{'} = v

v^{'} = μ (1 - x^{2}) v - x

For small mu, the problem behaves like a moderately nonlinear oscillator. For large mu, it develops fast and slow time scales and becomes stiff.

This benchmark tests:

nonlinear dynamics,
limit-cycle behavior,
stiffness transition,
robustness under changing time scales.

It is useful because a method that works well for mu = 1 may struggle when mu becomes large.

5.4 Lorenz System

The Lorenz system is a classical chaotic system:

x ˙ = σ (y - x)

y ˙ = x (ρ - z) - y

z ˙ = x y - β z

This benchmark tests:

nonlinear instability,
sensitive dependence on initial conditions,
qualitative trajectory behavior,
attractor visualization,
reliability over finite time windows.

For chaotic systems, a numerical trajectory should not be expected to match a reference solution forever. Instead, the benchmark should focus on reasonable finite-time accuracy and qualitative behavior.

5.5 Robertson Stiff Chemical Kinetics

The Robertson problem is a classical stiff chemical kinetics benchmark:

y_{1'} = - 0.04 y_{1} + 1 0^{4} y_{2} y_{3}

y_{2'} = 0.04 y_{1} - 1 0^{4} y_{2} y_{3} - 3 \cdot 1 0^{7} y_{22}

y_{3'} = 3 \cdot 1 0^{7} y_{22}

It also has a mass-conservation property:

y_{1} + y_{2} + y_{3} = 1

and the species concentrations should remain nonnegative:

y_{i} (t) \geq 0

This benchmark tests:

stiffness,
positivity preservation,
conservation of total mass,
robustness of implicit and linearly implicit methods,
cost of Jacobians, Newton iterations, and linear solves.

Explicit methods can require extremely small time steps on this problem. Stiff solvers are usually more appropriate.

5.6 Kepler Two-Body Problem

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

r^{'} = v

v^{'} = - μ \frac{r}{∣ r ∣ ^{3}}

where

r = (x, y), v = (v_{x}, v_{y})

The total energy is

E (r, v) \frac{1}{2} ∣ v ∣^{2} \frac{μ}{∣ 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 tests:

long-time orbital behavior,
energy drift,
angular momentum drift,
phase error,
orbit distortion,
suitability of symplectic methods.

A method may have good short-time accuracy but slowly distort the orbit. In orbital mechanics and robotics simulation, long-time qualitative fidelity is often more important than only pointwise short-time error.

6. Expanded Implemented Benchmark Families

The current repository goes beyond the original six examples. It now includes a wider set of implemented benchmarks grouped by mathematical behavior.

6.1 Non-Stiff Accuracy and Nonlinear Dynamics

This group includes smooth and nonlinear ODE examples such as:

Brusselator ODE,
Duffing oscillator,
forced damped pendulum,
logistic growth,
Lotka-Volterra predator-prey model,
rigid-body Euler equations,
SIR epidemic model.

These benchmarks are useful for testing:

nonlinear behavior,
oscillations,
phase-space geometry,
population dynamics,
sensitivity to step size,
qualitative trajectory preservation.

For example, the logistic equation has the form:

y^{'} = ry (1 - \frac{y}{K})

This benchmark is useful because it has simple nonlinear saturation and a known qualitative behavior: solutions should approach the carrying capacity K.

The Lotka-Volterra predator-prey model has the form:

x^{'} = αx - β x y

y^{'} = δ x y - γ y

It tests whether a method can represent coupled nonlinear oscillations without artificial damping or growth.

6.2 Stiff Multiscale and Chemical Kinetics Benchmarks

The stiff benchmark family now includes examples such as:

Akzo Nobel chemical kinetics,
Allen-Cahn method-of-lines stiff ODE,
HIRES problem,
Oregonator / Belousov-Zhabotinsky model,
pollution model ODE,
Prothero-Robinson problem,
ring modulator circuit ODE,
stiff Brusselator regime.

These examples are important because stiffness is not only a theoretical issue. It appears in chemical kinetics, reaction-diffusion systems, circuit models, and multiscale physical systems.

A stiff system usually contains widely separated time scales. Informally, a system may have both fast decaying modes and slow solution features. The numerical method must remain stable on the fast modes while accurately tracking the slow dynamics.

A linearized stiff system may have eigenvalues satisfying

∣ Re (λ_{m a x}) ∣ ≫ ∣ Re (λ_{m i n}) ∣

where the fastest stable mode imposes a severe step-size restriction on explicit methods.

This is why stiff benchmarks should record more than error and CPU time. They should also examine:

rejected steps,
Jacobian evaluations,
Newton iterations,
linear solves,
failure status,
positivity violation,
conservation error.

6.3 Hamiltonian, Geometric, and Long-Time Mechanics Benchmarks

The repository also includes mechanical and Hamiltonian-style benchmarks such as:

charged particle in magnetic field,
double pendulum,
FPUT lattice,
Henon-Heiles Hamiltonian,
N-body gravitational problem,
planar restricted three-body problem,
Toda lattice.

These benchmarks are designed to expose long-time structure preservation.

A Hamiltonian system has the canonical form:

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

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

For the exact continuous system, the Hamiltonian is conserved:

H (q (t), p (t)) = H (q (0), p (0))

However, a general-purpose integrator does not necessarily preserve this structure. A non-symplectic method may introduce artificial damping or artificial energy growth over long simulations.

This is why symplectic methods, variational integrators, discrete-gradient methods, and energy-momentum ideas are important for this family of benchmarks.

6.4 PDE Method-of-Lines and Conservation-Law Benchmarks

The expanded repository also includes PDE-derived problems after spatial discretization, including:

advection equation,
Cahn-Hilliard equation,
Fisher-KPP equation,
Gray-Scott reaction-diffusion,
heat equation,
inviscid Burgers equation,
Korteweg-de Vries equation,
Kuramoto-Sivashinsky equation,
viscous Burgers equation,
wave equation.

A PDE can often be converted into a large ODE system using the method of lines. A general semi-discrete form is:

M y ˙ = F (t, y)

or, when the mass matrix is the identity,

y ˙ = F (t, y)

For diffusion-dominated problems, spatial discretization can generate stiffness. For hyperbolic conservation laws, numerical dissipation, dispersion, shock handling, and stability restrictions become important.

For example, the heat equation

u_{t} = α u_{xx}

after finite-difference spatial discretization becomes a system of ODEs. The eigenvalues of the discrete Laplacian can impose severe time-step restrictions on explicit methods.

The advection equation

u_{t} + c u_{x} = 0

tests phase accuracy, numerical diffusion, and wave propagation behavior.

These PDE-derived benchmarks make the project more useful for computational physics, numerical PDEs, and scientific machine simulation workflows.

7. Planned Benchmark Directions

The repository also contains planned or scaffolded benchmark directions. These should not be treated as fully validated runnable benchmark cases until their implementation status is changed to implemented.

The planned directions include:

DAE constrained systems and circuits,
stochastic differential equations,
delay differential equations,
fractional differential equations,
boundary-value problems,
hybrid and event-driven systems.

This planned layer is useful because it shows how the project can grow beyond ordinary ODE integration.

For a DAE, the problem may have the form:

F (t, y, y ˙) = 0

For a delay differential equation:

x^{'} (t) = f (t, x (t), x (t - τ))

For an SDE:

d X = a (X, t), d t + b (X, t), d W_{t}

For a fractional differential equation with a Caputo derivative:

^{C} D_{t α} y (t) = f (t, y (t)), 0 < α < 1

Each of these problem classes requires different diagnostics. For example, SDEs require statistical metrics, DDEs require history interpolation diagnostics, DAEs require constraint residuals, and fractional equations require memory-aware numerical methods.

8. Implemented Integrator Families

The current runnable method registry includes several important families of numerical methods.

8.1 Explicit One-Step and Low-Order Runge-Kutta Methods

Implemented examples include:

Forward Euler,
Explicit midpoint,
Heun method,
Ralston method,
Kutta RK3,
classical RK4.

Forward Euler is the simplest method:

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

It is first order:

e_{global} = O (h)

It is useful as a baseline, but it is usually not competitive for serious simulation because of its low order and limited stability.

The explicit midpoint method uses an internal half-step:

k_{1} = f (t_{n}, y_{n})

k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1})

y_{n + 1} = y_{n} + h k_{2}

Classical RK4 uses four stages:

k_{1} = f (t_{n}, y_{n})

k_{2} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1})

k_{3} = f (t_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{2})

k_{4} = f (t_{n} + h, y_{n} + h k_{3})

y_{n + 1} y_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

RK4 is a strong baseline for smooth non-stiff ODEs, but it is not a stiff solver and not a structure-preserving geometric method.

8.2 Adaptive Runge-Kutta Methods

Implemented adaptive RK examples include:

Runge-Kutta-Fehlberg RKF45,
Dormand-Prince RK45.

Adaptive methods estimate local error by comparing two approximations of different orders. A typical embedded pair produces

y_{n + 1 [p]} and y_{n + 1 [p + 1]}

and estimates the local error as

err_{n + 1} ∣ ∣ y_{n + 1 [p + 1]} - y_{n + 1 [p]} ∣ ∣

The step size is then adjusted. A common form is

h_{new} s h (\frac{tol}{err _{n + 1}})^{1/ (p + 1)}

where s is a safety factor.

Adaptive RK methods are useful for non-stiff problems with varying solution scales. However, adaptivity does not automatically solve stiffness. A stiff problem may still require implicit, BDF, Rosenbrock, or other stiffly stable methods.

8.3 Implicit One-Step Methods

Implemented implicit one-step examples include:

Backward Euler,
implicit midpoint,
implicit trapezoidal rule.

Backward Euler is

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

This is an implicit equation for y_{n+1}. It usually requires a nonlinear solve.

Define the residual:

G (Y) = Y - y_{n} - h f (t_{n + 1}, Y)

The method solves

G (Y) = 0

typically by Newton's method.

Backward Euler is first order and strongly damping, but it is very useful for stiff decay problems.

The implicit trapezoidal rule is

y_{n + 1} y_{n} + \frac{h}{2} [f (t_{n}, y_{n}) + f (t_{n + 1}, y_{n + 1})]

It is second order and often less dissipative than Backward Euler.

8.4 Linear Multistep and Predictor-Corrector Methods

Implemented multistep examples include:

Adams-Bashforth 2,
Adams-Bashforth 4,
Adams-Moulton 2,
Adams-Bashforth-Moulton 4 predictor-corrector,
BDF2.

A linear multistep method uses multiple previous solution values:

j = 0 \sum s α_{j} y_{n + j} h j = 0 \sum s β_{j} f (t_{n + j}, y_{n + j})

Adams-Bashforth methods are explicit. BDF methods are implicit and are important in stiff ODE integration.

BDF2 has the form

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

Because BDF2 is implicit, it generally requires nonlinear solves, Jacobians, and linear algebra. Therefore, benchmarking BDF-style methods should include more than error and CPU time.

8.5 Rosenbrock and Linearly Implicit Methods

Implemented linearly implicit examples include:

Rosenbrock-Euler,
ROS2.

Rosenbrock methods avoid fully nonlinear solves by linearizing around the current state. A simplified Rosenbrock-Euler form is:

(I - hγ J_{n}) 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})

These methods are useful for stiff systems because they use Jacobian information without requiring a full nonlinear Newton iteration at every step.

The benchmark should record:

Jacobian evaluations,
linear solves,
right-hand-side evaluations,
stability behavior,
and failure status.

8.6 Exponential Methods

The implemented catalog also includes exponential Euler.

Exponential methods are important for systems with a stiff linear part:

y^{'} = A y + g (t, y)

The exact solution formula involves the matrix exponential:

y (t + h) e^{h A} y (t) + \int_{0 h} e^{(h - τ) A} g (t + τ, y (t + τ)), d τ

Exponential Euler approximates the nonlinear part while treating the linear part through exponential propagation.

These methods are especially relevant for semi-discretized PDEs, where the linear stiff part may dominate the stability restriction.

8.7 Geometric, Symplectic, and Energy-Aware Methods

Implemented geometric methods include:

Velocity Verlet,
Stormer-Verlet,
Yoshida fourth-order symplectic composition,
discrete-gradient oscillator method.

For a separable Hamiltonian system

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

the equations are

q^{'} = \nabla_{p} T (p)

p^{'} = - \nabla_{q} V (q)

Velocity Verlet has the form

p_{n + \frac{1}{2}} p_{n} - \frac{h}{2} \nabla V (q_{n})

q_{n + 1} q_{n} + h \nabla T (p_{n + \frac{1}{2}})

p_{n + 1} p_{n + \frac{1}{2}} - \frac{h}{2} \nabla V (q_{n + 1})

Symplectic methods are not designed to minimize final-time error only. Their strength is long-time qualitative behavior. They often keep the numerical trajectory on a more physically meaningful path over long simulations.

9. Work-Precision Behavior

A useful benchmark should compare error against cost.

A work-precision diagram usually plots an error metric against a work metric such as:

CPU time,
number of function evaluations,
number of accepted steps,
number of Jacobian evaluations,
number of linear solves.

A method is not better simply because it has lower error at one step size. It should be compared across different tolerances or step sizes.

A typical work-precision study asks:

How fast does error decrease as work increases?
Which methods fail before reaching the target accuracy?
Which methods become unstable for larger step sizes?
Which methods are accurate but too expensive?
Which methods are cheap but unreliable?

For example, RK4 may be efficient for smooth non-stiff ODEs, but a stiff chemical kinetics problem may require a stiff method even if each step is more expensive.

10. Applicability Matters

One important improvement in the repository is the explicit idea of method applicability.

Not every method should run on every benchmark.

For example:

Velocity Verlet is meaningful for separable mechanical systems, but not for arbitrary chemical kinetics.
Stochastic methods are meaningful for SDEs, not deterministic ODEs.
DAE methods are meaningful for constrained problems, not ordinary explicit ODE-only tests.
PDE spatial discretization methods belong to PDE-derived benchmarks.
BVP methods should not be evaluated like initial-value solvers.

Therefore, the benchmark framework should distinguish between:

success,
failure,
skipped,
not applicable,
planned but not implemented.

This is better than forcing every method to run on every problem. A solver comparison should respect mathematical structure.

11. Repository Structure

The repository now mirrors documentation and source-code organization.

A simplified view is:

.
├── README.md
├── LICENSE
├── CONTRIBUTING.md
├── CITATION.cff
├── run_all.m
├── run_one_benchmark.m
├── docs/
│   ├── methods/
│   ├── benchmarks/
│   ├── tutorials/
│   └── references.md
├── src/
│   ├── config/
│   │   ├── default_options.m
│   │   ├── method_registry.m
│   │   └── benchmark_registry.m
│   ├── core/
│   │   ├── run_benchmark.m
│   │   ├── compute_metrics.m
│   │   └── newton_solve.m
│   ├── methods/
│   ├── benchmarks/
│   └── plotting/
├── tests/
└── results/

The registries are important because they define which methods and benchmarks are implemented, planned, applicable, or skipped.

To inspect implemented methods and benchmarks inside MATLAB:

methods = method_registry();
benchmarks = benchmark_registry();

implementedMethods = methods([methods.implemented]);
implementedBenchmarks = benchmarks([benchmarks.implemented]);

{implementedMethods.name}'
{implementedBenchmarks.name}'

12. Recommended MATLAB Workflow

A typical workflow is:

run_all

To run a single benchmark:

run_one_benchmark('VanDerPolOscillator')

To customize options:

opts = default_options();

opts.h = 5e-3;
opts.n_output = 2001;
opts.write_plots = true;
opts.write_tables = true;
opts.write_global_summary = true;

results = run_one_benchmark('LorenzSystem', opts);

To keep the run restricted to implemented entries:

opts.include_planned_methods = false;
opts.include_planned_benchmarks = false;

Results are written under:

results/

or under benchmark-specific folders such as:

results/<BenchmarkName>/

Typical outputs include:

CSV summary tables,
final-error plots,
CPU-time plots,
function-evaluation plots,
work-precision diagrams,
invariant drift plots,
phase portraits,
orbit plots,
attractor projections,
benchmark-specific solution visualizations.

13. Why the Expanded Benchmark Set Is Useful

The expanded repository is useful because numerical integrators have different strengths.

Explicit Runge-Kutta methods are often strong for smooth non-stiff ODEs.

Adaptive RK methods are useful when local-error control is needed.

Implicit methods are important for stiff decay and stiff chemical kinetics.

BDF and multistep methods are important for stiff and long-time integration.

Rosenbrock methods are useful when Jacobian-based linearized stiffness handling is needed.

Exponential methods are useful when the problem has a strong linear stiff part.

Symplectic and geometric methods are important for long-time Hamiltonian mechanics.

PDE method-of-lines benchmarks expose stiffness, dispersion, diffusion, and conservation issues caused by spatial discretization.

DAE, SDE, DDE, fractional, BVP, and hybrid systems require additional diagnostics beyond ordinary ODE error metrics.

Therefore, the main lesson is:

There is no universal best integrator. There is only the best integrator for a given problem structure, accuracy requirement, stability constraint, and computational budget.

14. Practical Interpretation of Benchmark Results

When interpreting the benchmark outputs, I would avoid ranking methods using a single metric.

A better interpretation is problem-dependent.

For the linear test equation, stability behavior and decay accuracy matter.

For the harmonic oscillator, phase error and energy drift matter.

For Van der Pol, stiffness transition and nonlinear oscillation matter.

For Lorenz, finite-time accuracy and qualitative chaotic behavior matter.

For Robertson kinetics, stiffness, positivity, mass conservation, and solver robustness matter.

For Kepler and other Hamiltonian systems, long-time orbital geometry, energy drift, and angular momentum drift matter.

For PDE method-of-lines benchmarks, both spatial and temporal discretization behavior matter.

For future DAE benchmarks, constraint residuals and consistent initialization will matter.

For future SDE benchmarks, strong and weak convergence should be measured statistically, not only through a single deterministic trajectory.

This is why benchmark results should be interpreted as a numerical profile, not as a single winner list.

15. 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 a geometric energy-preserving 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.

The updated MATLAB repository reflects this idea by separating runnable implementations from a larger numerical-method roadmap. It includes implemented methods, implemented benchmarks, planned/scaffolded future directions, plotting utilities, result tables, and documentation categories for a broad set of differential-equation problems.

The most important conclusion is:

A good benchmark suite should measure accuracy, cost, stability, stiffness behavior, invariant drift, method applicability, and failure modes. Numerical integration is problem-dependent, and solver choice should be guided by the mathematical structure of the model.

Building a Four-Bar Linkage Mechanism Simulator in Haskell

Abolfazl Mohammadijoo — Tue, 09 Jun 2026 21:01:42 +0000

Most developers know Haskell as a language for functional programming, type safety, compilers, parsers, and beautiful mathematical abstractions.

But can Haskell also be used to build an interactive engineering simulator?

That was the motivation behind my project:

Four-Bar Mechanism Haskell Simulator

Repository:

https://github.com/mohammadijoo/Four-Bar-Mechanism-Haskell

This project is a browser-backed desktop-style GUI application written in Haskell. It visualizes, classifies, and animates a planar four-bar linkage mechanism, which is one of the most classical mechanisms in mechanical engineering, kinematics, and machine design.

The GUI is built with Threepenny-GUI, so the interface runs in a local browser window, while the mathematical model and mechanism logic remain written in Haskell.

For me, the interesting part was not only drawing a moving linkage. It was about connecting mechanism design theory, computational geometry, and functional programming in one small educational simulator.

What is a four-bar linkage?

A four-bar linkage is a closed-loop mechanical system made from four rigid links connected by four revolute joints.

In this project, the four links are:

Symbol	Name	Description
`g`	Ground link	Fixed distance between pivots A and B
`a`	Input link	Rotating link from A to moving pivot C
`b`	Output link	Link from fixed pivot B to moving pivot D
`f`	Floating link / coupler	Link connecting moving pivots C and D

The fixed pivots are placed at:

A = (0, 0), B = (g, 0)

The input link rotates by angle $α$ . Therefore, point C can be computed directly as:

C = (a cos α,; a sin α)

Point D is more interesting.

It must satisfy two geometric distance constraints:

∣ D - C ∣ = f

∣ D - B ∣ = b

So the simulator solves the position of point D using a circle-intersection method.

One circle is centered at C with radius f.
The other circle is centered at B with radius b.

Where those two circles intersect, the mechanism can close.

That is the basic geometric heart of the simulator.

Why four-bar linkages matter

Four-bar mechanisms appear in many real mechanical systems.

Some examples include:

Application	Role of the four-bar idea
Engines	Converts rotating and oscillating motion
Suspension systems	Guides constrained motion
Robotic grippers	Creates controlled opening and closing paths
Pumps	Converts rotary motion into reciprocating motion
Walking mechanisms	Generates repeated leg-like motion
Windshield wipers	Converts motor rotation into sweeping motion
Machine tools	Produces repeatable constrained motion

A four-bar linkage is simple enough to understand visually, but rich enough to teach real mechanism design concepts.

That is why it appears so often in mechanical engineering and machine design courses.

Mobility of the mechanism

A planar mechanism can be analyzed using the Kutzbach-Grübler mobility equation.

M = 3 (n - 1) - 2 j_{1} - j_{2}

where:

Symbol	Meaning
`M`	Mobility / degrees of freedom
`n`	Number of links
`j1`	Number of lower pairs, such as revolute joints
`j2`	Number of higher pairs

For an ideal planar four-bar linkage:

n = 4, j_{1} = 4, j_{2} = 0

Substituting these values gives:

M = 3 (4 - 1) - 2 (4) - 0

M = 9 - 8 = 1

So a four-bar linkage has exactly one degree of freedom.

This means that once the input angle is chosen, the rest of the linkage configuration is determined by the loop-closure constraints, except for the selected assembly branch.

That is a nice result: one input motion drives the whole mechanism.

The implemented geometric model

The simulator computes the mechanism position from the link lengths and the input angle.

First, the distance between B and C is calculated:

d = ∣ B - C ∣

For the mechanism to close, the two circles must intersect. The assembly condition is:

∣ f - b ∣ \leq d \leq f + b

If this condition is not satisfied, there is no real point D, and the mechanism cannot be assembled at that input angle.

When the condition is satisfied, there are usually two possible solutions for point D. These represent the two assembly branches of the mechanism.

In the GUI, the user can switch between these two branches using the Flip branch button or the F key.

This is one of the parts I like most about the project: the animation is not just a pre-drawn visual effect. The simulator actually solves the linkage geometry.

Grashof condition

One of the most important ideas in four-bar linkage design is the Grashof condition.

Let the four link lengths be sorted as:

s \leq p \leq q \leq l

where s is the shortest link and l is the longest link.

The Grashof index is:

G = s + l - p - q

The classification is:

Condition	Meaning
`G < 0`	Grashof mechanism
`G = 0`	Change-point mechanism
`G > 0`	Non-Grashof mechanism

This matters because a Grashof mechanism usually has at least one link that can make a full rotation relative to the ground.

In practical terms, this helps answer an important design question:

Can the input link rotate continuously, or will it only rock back and forth?

That is exactly the kind of question a mechanism designer needs to answer early.

Validity index

The simulator also checks whether the mechanism is physically assemblable.

The validity index is:

V = l - s - p - q

A valid four-bar mechanism must satisfy:

V \leq 0

If:

V > 0

then the longest link is longer than the sum of the other three links. In that case, the loop cannot close.

So the simulator reports not only the motion type, but also whether the chosen link lengths describe a valid mechanism.

This is useful because a random set of four link lengths does not always produce a real four-bar linkage.

Full linkage state classification

The simulator computes three excess values:

T_{1} = g + f - b - a

T_{2} = b + g - f - a

T_{3} = f + b - g - a

The signs of these values help classify the detailed linkage state.

The GUI then displays the simplified motion type:

Classification	Meaning
`crank-crank`	Both relevant links can rotate fully
`crank-rocker`	Input rotates fully, output rocks
`rocker-crank`	Input rocks, output rotates fully
`rocker-rocker`	Both links oscillate

This turns the project into more than a moving diagram. It becomes a small mechanism classification tool.

Why use Haskell?

Haskell is probably not the first language most people would choose for a mechanical simulator.

Python, MATLAB, C++, JavaScript, or Julia are more common choices for engineering applications.

But that is exactly why this project was interesting.

Haskell is very good at separating pure mathematical logic from side-effecting application code.

In this project, the mechanism model is separated into a pure Haskell module:

src/FourBar.hs

This module contains the core mechanism logic, including:

Linkage
Pose
excesses
grashofIndex
validityIndex
classifyMotion
basicMotionState
solvePose
nearestValidAlpha

The GUI and event loop are handled separately in:

app/Main.hs

This separation is important.

The mechanism solver can later be reused in tests, command-line tools, notebooks, optimization programs, or even a web backend.

So Haskell helps keep the engineering model clean, testable, and reusable.

That is one of the main reasons I enjoyed building this project.

Threepenny-GUI and the browser-backed interface

The project uses Threepenny-GUI.

Threepenny-GUI allows a Haskell program to create an interface that is rendered in a web browser. The application starts a local GUI server, and the user opens:

http://127.0.0.1:8023

The browser is the visual layer, but the application logic remains in Haskell.

The GUI includes:

Feature	Description
Link sliders	Change `g`, `a`, `b`, and `f` interactively
SVG drawing	Renders the mechanism as scalable browser graphics
Play / Pause	Starts and stops the animation
Flip branch	Switches between the two assembly solutions
Reset	Restores the default mechanism
Angle display	Shows live alpha and beta values
Classification panel	Shows Grashof, validity, and crank/rocker state
Keyboard controls	Space, `F`, and `R` shortcuts

I like this architecture because it combines two worlds:

Haskell \to mechanism logic

Browser + SVG \to visual interface

Threepenny-GUI acts as the bridge between them.

SVG visualization

The mechanism is rendered using SVG.

SVG is a good fit for this kind of project because it is scalable, clean, and works naturally inside a browser interface.

The simulator draws:

fixed pivots
moving pivots
links
labels
angle arcs
the current linkage configuration

Because the GUI is browser-rendered, the interface can use HTML, CSS, fonts, and SVG while still keeping the engineering logic in Haskell.

The result is a simple but effective architecture:

Pure Haskell solver + Threepenny-GUI + SVG rendering

Repository structure

The repository is organized like this:

Four-Bar-Mechanism-Haskell/
├── app/
│   └── Main.hs
├── src/
│   └── FourBar.hs
├── resources/
│   └── fonts/
│       └── Inter-Regular.ttf
├── fourbar-linkage-threepenny.cabal
├── README.md
├── LICENSE
└── .gitignore

The most important files are:

src/FourBar.hs
app/Main.hs

FourBar.hs contains the pure mechanism mathematics.

Main.hs builds the GUI, handles events, updates the state, and renders the current mechanism as SVG.

This structure makes the project easier to understand because the engineering model is not mixed directly into the interface code.

How to run the project

First, install Haskell using GHCup.

Then clone the repository:

git clone https://github.com/mohammadijoo/Four-Bar-Mechanism-Haskell.git
cd Four-Bar-Mechanism-Haskell

Update Cabal:

cabal update

Run the simulator:

cabal run fourbar-linkage-threepenny

If the browser does not open automatically, open:

http://127.0.0.1:8023

Controls

Control	Action
Drag `g` slider	Change ground link length
Drag `a` slider	Change input link length
Drag `b` slider	Change output link length
Drag `f` slider	Change floating/coupler link length
Play / Pause	Start or stop animation
Flip branch	Switch assembly branch
Reset	Restore default mechanism
Space	Toggle Play / Pause
`F`	Flip branch
`R`	Reset

The sliders make the simulator useful for experimentation.

You can change link lengths and immediately see how the mechanism behavior changes.

That is helpful for learning because mechanism theory becomes visual.

What this project demonstrates

This project demonstrates several ideas at the same time:

Area	What the project shows
Mechanical engineering	Four-bar linkage geometry and classification
Kinematics	Loop-closure constraints and assembly branches
Computational geometry	Circle intersection for solving point D
Haskell programming	Pure modeling separated from GUI code
GUI development	Browser-backed interface with Threepenny-GUI
Visualization	SVG-based scalable mechanism drawing

For me, the most useful lesson is that Haskell can be used for more than abstract examples.

It can also be used to build small, practical, interactive engineering tools.

Possible future improvements

Some future improvements could make the simulator even more useful:

Add coupler point tracing
Plot the coupler curve
Export motion data as CSV
Add velocity and acceleration analysis
Add transmission angle plots
Add mechanical advantage visualization
Add automated tests for mechanism classification
Package the app as a standalone desktop executable

The most interesting improvement would probably be coupler curve tracing.

Many real design problems ask:

What link lengths should I choose so that a point on the coupler follows a desired path?

Adding that feature would move the project from simulation toward mechanism synthesis.

Conclusion

The four-bar linkage is simple, but it contains many important engineering ideas.

With only four links and four revolute joints, it can produce rich motion patterns and useful mechanical transformations.

In this project, I used Haskell to model the mechanism, Threepenny-GUI to build the interface, and SVG to render the animation in a browser-backed desktop-style application.

The result is a small educational simulator that connects theory, geometry, and functional programming.

Repository:

https://github.com/mohammadijoo/Four-Bar-Mechanism-Haskell

Thanks for reading.

Building a Double Pendulum Chaos Simulator in Prolog

Abolfazl Mohammadijoo — Mon, 08 Jun 2026 13:27:36 +0000

Most developers know Prolog as a language for symbolic reasoning, expert systems, search, theorem proving, and artificial intelligence. But can Prolog also be used to build a numerical physics simulator with a desktop GUI?

That was the motivation behind my project:

Prolog Double Pendulum Chaos Simulator

Repository:

https://github.com/mohammadijoo/Prolog-Double-Pendulum-Chaos-Simulator

The project is a desktop GUI application written in SWI-Prolog using the XPCE graphical toolkit. It simulates the classical double pendulum, one of the most famous examples of deterministic chaos.

The application lets the user define physical parameters such as initial angles, angular velocities, link lengths, masses, gravity, damping, numerical time step, and perturbation size. Then it animates the double pendulum and draws the trajectory of the outer endpoint.

The interesting part is that the motion is completely deterministic, but for many initial conditions it becomes visually irregular, highly sensitive, and practically unpredictable.

That is the essence of chaos.

What is a double pendulum?

A simple pendulum has one mass attached to one link. For small angles, it behaves almost like a linear oscillator. Its motion is smooth, periodic, and relatively easy to predict.

A double pendulum adds only one more link and one more mass:

Link 1 rotates around a fixed pivot.
Link 2 rotates around the end of link 1.
Each link has its own angle and angular velocity.
Gravity acts on both masses.

That small structural change is enough to make the dynamics much richer.

The state of the system can be written as:

x = [θ_{1}, θ_{2}, ω_{1}, ω_{2}]^{T}

where:

Symbol	Meaning
`theta1` / $θ_{1}$	angle of the first link
`theta2` / $θ_{2}$	angle of the second link
`omega1` / $ω_{1}$	angular velocity of the first link
`omega2` / $ω_{2}$	angular velocity of the second link

In first-order state-space form:

θ ˙_{1} = ω_{1}

θ ˙_{2} = ω_{2}

ω ˙_{1} = f_{1} (θ_{1}, θ_{2}, ω_{1}, ω_{2}, m_{1}, m_{2}, L_{1}, L_{2}, g)

ω ˙_{2} = f_{2} (θ_{1}, θ_{2}, ω_{1}, ω_{2}, m_{1}, m_{2}, L_{1}, L_{2}, g)

This means the angular accelerations depend on the current angles, angular velocities, masses, lengths, and gravity.

Why is the double pendulum chaotic?

Chaos does not mean randomness.

A chaotic system can be fully deterministic. The equations are fixed. There is no random force inside the mathematical model. But tiny differences in the starting condition can grow rapidly and produce very different future behavior.

This is called sensitive dependence on initial conditions, often described as the butterfly effect.

In this project, I added a small perturbation option so that the user can compare two nearly identical initial states:

x_{A} (0) = [θ_{1}, θ_{2}, ω_{1}, ω_{2}]^{T}

x_{B} (0) = [θ_{1} + δ, θ_{2}, ω_{1}, ω_{2}]^{T}

Here, $δ$ may be extremely small, for example:

δ = 0.0 5^{\circ}

At the beginning, both motions may look almost identical. After some time, the trajectories can separate dramatically.

That is not because the simulator is broken. It is exactly what the nonlinear dynamics predict.

The double pendulum becomes chaotic because of several factors.

1. Nonlinear trigonometric coupling

The equations contain terms like:

sin (θ_{1} - θ_{2})

cos (θ_{1} - θ_{2})

sin (θ_{1} - 2 θ_{2})

These terms couple the two links in a nonlinear way.

2. Velocity-squared terms

The accelerations contain terms involving:

ω_{12}, ω_{22}

These terms represent inertial and centripetal effects. When the pendulum moves fast, these nonlinear terms become very important.

3. Energy transfer between links

Energy can move back and forth between the two links. Sometimes the second link swings gently. Sometimes it flips rapidly. This exchange creates complex motion.

4. Bounded but irregular motion

The pendulum is physically bounded by its link lengths, but the endpoint path can become tangled and aperiodic. This is a typical sign of deterministic chaos: the motion is bounded, but it does not settle into a simple repeating pattern.

The implemented dynamic model

For angles measured from the downward vertical direction, the project uses the classical double-pendulum equations.

Define:

Δ = θ_{1} - θ_{2}

D = 2 m_{1} + m_{2} - m_{2} cos (2 θ_{1} - 2 θ_{2})

Then the angular acceleration of the first link is:

ω ˙_{1} = \frac{- g ( 2 m _{1} + m _{2} ) sin ( θ _{1} ) - m _{2} g sin ( θ _{1} - 2 θ _{2} ) - 2 m _{2} sin ( Δ ) ( ω _{22} L _{2} + ω _{12} L _{1} cos ( Δ ))}{L _{1} D} - c ω_{1}

The angular acceleration of the second link is:

ω ˙_{2} = \frac{2 sin ( Δ ) ( ω _{12} L _{1} ( m _{1} + m _{2} ) + g ( m _{1} + m _{2} ) cos ( θ _{1} ) + ω _{22} L _{2} m _{2} cos ( Δ ) )}{L _{2} D} - c ω_{2}

where:

Symbol	Meaning
$m_{1}, m_{2}$	masses
$L_{1}, L_{2}$	link lengths
$g$	gravitational acceleration
$c$	optional damping coefficient

When $c = 0$ , the system is ideal and conservative except for numerical integration error.

When $c > 0$ , a small damping term is applied.

Converting angles to Cartesian positions

To draw the pendulum, the simulator converts angular states into Cartesian coordinates.

The first bob position is:

x_{1} = L_{1} sin (θ_{1})

y_{1} = L_{1} cos (θ_{1})

The second bob position is:

x_{2} = x_{1} + L_{2} sin (θ_{2})

y_{2} = y_{1} + L_{2} cos (θ_{2})

The simulator draws the trajectory of:

(x_{2}, y_{2})

which is the endpoint of the outer link.

That endpoint trajectory is the most visually interesting part of the simulation. It often forms irregular loops, folds, and tangled paths.

Numerical integration with RK4

The system is described by nonlinear differential equations. To animate it, the simulator must repeatedly advance the state forward in time.

For this project I used the classical fourth-order Runge-Kutta method, usually called RK4.

For a state equation:

x ˙ = f (x, t)

RK4 computes:

k_{1} = f (x_{n}, t_{n})

k_{2} = f (x_{n} + \frac{h}{2} k_{1}, t_{n} + \frac{h}{2})

k_{3} = f (x_{n} + \frac{h}{2} k_{2}, t_{n} + \frac{h}{2})

k_{4} = f (x_{n} + h k_{3}, t_{n} + h)

and then:

x_{n + 1} = x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})

In the GUI, the user can choose the time step dt.

A smaller value such as 0.01 usually gives smoother and more accurate behavior.

A larger value such as 0.03 or 0.04 may run faster but can introduce larger numerical error.

This is especially important in chaotic systems, because numerical error can also grow over time.

Why use Prolog for this?

Prolog is usually not the first language people choose for physics simulation. Python, MATLAB, C++, Julia, or JavaScript are more common choices.

That is exactly why this project was interesting.

Prolog is a logic programming language. Instead of writing programs mainly as sequences of commands, you describe facts, rules, relations, and goals.

A Prolog predicate can describe a relationship such as:

deg_to_rad(Degrees, Radians).

or:

rk4_step(State, Params, Dt, NextState).

In this project, Prolog is used for:

simulation state management
parsing user input
numerical derivative computation
RK4 integration
GUI control logic
drawing commands
event handling
timer-based animation
logging and reset behavior

This makes the project a useful experiment in applying a symbolic/logical language to a numerical and graphical problem.

The main lesson is that Prolog is not limited to toy examples or pure symbolic AI. With SWI-Prolog and XPCE, it can also build interactive desktop applications.

SWI-Prolog and XPCE

The project uses SWI-Prolog, a mature open-source Prolog implementation.

For the GUI, it uses XPCE, the graphical toolkit available with SWI-Prolog.

XPCE provides GUI elements such as:

frames
dialogs
buttons
text fields
drawing canvases
graphical objects
timers
event messages

The simulator uses XPCE to create:

the main desktop application window
a MATLAB-like trajectory view
a CAD-like mechanism view
control fields for initial conditions
Play / Pause and Reset buttons
real-time drawing of links, bobs, pivots, grids, and trajectories

GUI concept

The GUI is divided into two main visual areas.

Left pane: trajectory view

The left side shows the endpoint trajectory. It is similar to a scientific plotting or MATLAB-style visualization.

It includes:

coordinate grid
axes
double pendulum shape
endpoint trail
energy drift information

This view emphasizes the mathematical behavior of the system.

Right pane: mechanism view

The right side shows the physical pendulum in a more mechanical or CAD-like style.

It includes:

pivot
first link
second link
first bob
second bob
current mechanism configuration

This view emphasizes the physical structure of the system.

Bottom control pane

The bottom pane contains the input fields and controls:

initial angle 1
initial angle 2
angular velocity 1
angular velocity 2
link lengths
masses
gravity
damping
perturbation
time step
maximum trail points
Play / Pause
Reset
clock
motion log

Repository structure

The repository is organized as follows:

prolog-double-pendulum-chaos-simulator/
│
├── run.pl
├── README.md
├── LICENSE
├── .gitignore
├── CONTRIBUTING.md
│
├── src/
│   ├── double_pendulum_gui.pl
│   └── double_pendulum_physics.pl
│
├── docs/
│   ├── EQUATIONS.md
│   └── DEVELOPMENT_NOTES.md
│
├── examples/
│   └── initial_conditions.txt
│
└── screenshots/
    └── .gitkeep

The two most important source files are:

src/double_pendulum_physics.pl
src/double_pendulum_gui.pl

Physics module

The physics module contains the mathematical model and numerical integration code.

Important predicates include:

deg_to_rad(+Degrees, -Radians)

Converts degrees to radians.

rad_to_deg(+Radians, -Degrees)

Converts radians to degrees.

derivatives(+State, +Params, -Derivative)

Computes the state derivative of the double pendulum.

Example state representation:

state(Theta1, Theta2, Omega1, Omega2)

Example parameter representation:

params(M1, M2, L1, L2, G, Damping)

Derivative representation:

deriv(DTheta1, DTheta2, DOmega1, DOmega2)

The RK4 integrator is exposed through:

rk4_step(+State, +Params, +Dt, -NextState)

The Cartesian position calculation is handled by:

bob_positions(+State, +Params, -X1, -Y1, -X2, -Y2)

The module also estimates total mechanical energy:

total_energy(+State, +Params, -Energy)

This helps observe numerical drift during simulation.

GUI module

The GUI module contains the XPCE interface and animation loop.

Important predicates include:

main.

which builds and opens the desktop GUI.

build_gui.

which creates the main frame, drawing canvas, text fields, buttons, labels, and layout.

toggle_simulation.

which handles Play / Pause behavior.

start_new_simulation.

which reads the input fields, builds the initial state, creates the perturbation state, creates the timer, and starts the animation.

simulation_tick.

which is called repeatedly by the XPCE timer.

Each tick performs the main animation tasks:

advance both pendulums using RK4
update the simulation clock
append new endpoint trajectory points
redraw the GUI
update logs
stop automatically after the maximum simulation time

The drawing helper predicates wrap XPCE drawing objects such as lines, circles, boxes, and text.

How to run the project

First install SWI-Prolog with XPCE support.

Clone the repository:

git clone https://github.com/mohammadijoo/Prolog-Double-Pendulum-Chaos-Simulator.git
cd Prolog-Double-Pendulum-Chaos-Simulator

Run:

swipl -q -f none run.pl

You can also start SWI-Prolog manually:

swipl

Then load and run:

?- [run].
?- main.

To verify that XPCE is available:

?- use_module(library(pce)).
true.

Suggested initial conditions

Here are some useful test cases.

Strong chaotic behavior

theta1 deg      120
theta2 deg      -10
omega1 rad/s    0.0
omega2 rad/s    0.0
L1 m            1.0
L2 m            1.0
m1 kg           1.0
m2 kg           1.0
g m/s^2         9.81
damping         0.000
perturb deg     0.05
dt s            0.02
trail points    1500

Smaller perturbation

theta1 deg      120
theta2 deg      -10
omega1 rad/s    0.0
omega2 rad/s    0.0
perturb deg     0.005
dt s            0.01

High-energy motion

theta1 deg      170
theta2 deg      80
omega1 rad/s    0.0
omega2 rad/s    0.0
perturb deg     0.05
dt s            0.01

Damped motion

theta1 deg      130
theta2 deg      -30
omega1 rad/s    0.0
omega2 rad/s    0.0
damping         0.02
perturb deg     0.05
dt s            0.02

With damping, the system eventually loses energy, but the early transient motion can still show strong sensitivity.

What this project demonstrates

This project is not just a pendulum animation.

It demonstrates several important ideas:

Deterministic chaos
The system follows exact equations, but long-term prediction becomes practically difficult.
Sensitivity to initial conditions
Very small perturbations can eventually produce very different trajectories.
Nonlinear dynamics
Trigonometric coupling and velocity-squared terms make the system behave very differently from a linear oscillator.
Numerical simulation
RK4 integration makes it possible to animate the differential equations step by step.
Prolog beyond symbolic AI
Prolog can also be used for GUI applications, numerical computation, and engineering-style educational software.
Visualization as a learning tool
Seeing the endpoint trajectory makes the concept of chaos much easier to understand than formulas alone.

Screenshots

Possible future improvements

Some possible directions for future versions:

export trajectory data to CSV
add phase-space plots
add Lyapunov exponent estimation
compare Euler, RK2, RK4, and symplectic integration
add real-time parameter sliders
add multi-pendulum ensemble visualization
add GIF/video export
add web version using Prolog backend and JavaScript frontend
add unit tests for the physics predicates
improve energy conservation diagnostics

A particularly interesting next step would be estimating the largest Lyapunov exponent. That would turn the visual idea of “nearby states separating” into a measurable quantity.

Conceptually, if two nearby trajectories separate approximately as:

d (t) \approx d (0) e^{λ t}

then $λ$ is a Lyapunov exponent.

If:

λ > 0

the system shows exponential divergence of nearby trajectories, which is a signature of chaos.

Conclusion

The double pendulum is a beautiful example of how simple mechanical components can produce complex nonlinear behavior.

Two links, two masses, gravity, and deterministic equations are enough to create motion that appears almost random.

Using Prolog for this project made it even more interesting. It shows that Prolog can be used not only for symbolic reasoning and AI-style logic programming, but also for building interactive engineering simulations with a desktop GUI.

Repository:

https://github.com/mohammadijoo/Prolog-Double-Pendulum-Chaos-Simulator

If you are interested in Prolog, numerical simulation, nonlinear dynamics, or educational engineering software, feel free to check the repository, try the simulator, and suggest improvements.

Thanks for reading.