ode45gpu: GPU-Accelerated ODE Solver

• ode45gpu is a specialized GPU implementation of MATLAB’s ode45 using the Dormand–Prince 4(5) scheme for parallel integration of independent ODE systems.
• It restructures the solver with unrolled scalar operations to comply with MATLAB’s arrayfun limitations, achieving significant speed-ups over conventional methods.
• Benchmark results show up to 22× faster performance and near-ideal scaling on GPUs like Tesla K80, P100, and V100 for large-scale atom-ion trajectory simulations.

Searching arXiv for recent and foundational papers on ode45gpu and GPU ODE integration. ode45gpu is an optimized and specialized implementation of the Runge–Kutta algorithm used in MATLAB’s ODE solver ode45, introduced in connection with atomiongpu.m as a GPU-accelerated MATLAB workflow for atom-ion dynamics. In the available description, its defining use case is the massively parallel simulation of trajectories of a trapped ion and an atom starting far away, with one trajectory executed per GPU thread through MATLAB’s gpuArray/arrayfun machinery. The implementation retains adaptive embedded Runge–Kutta step control, but restructures the solver so that each thread integrates its entire trajectory independently and returns only final states and selected observables rather than a dense trajectory history (Chowdhury et al., 15 Sep 2025).

1. Definition and problem setting

In the cited usage, ode45gpu is not presented as a built-in MATLAB solver. It is described instead as “our optimized and specialized implementation of the Runge-Kutta algorithm used in MATLAB’s ODE solver ode45,” embedded in atomiongpu.m and designed to parallelize molecular-dynamics-style trajectory calculations for atom-ion collisions (Chowdhury et al., 15 Sep 2025).

The reported benchmark problem is a 12-D atom-ion ODE with

$V(r)= -C_4/r^4 + C_8/r^8$

in a Paul trap. The 12 components correspond to “3 coordinates and 3 velocities for each of two particles,” and the solver is used to obtain final states and observables such as “the probability of complex formation, the distribution of observables such as the scattering angle and complex lifetime, and plots of specific trajectories” (Chowdhury et al., 15 Sep 2025).

A common misconception is to treat ode45gpu as a general synonym for any GPU-resident ode45-style integrator. The published description is narrower: it is a particular MATLAB implementation specialized for a specific class of trajectory simulations, even though its structure is representative of a broader literature on GPU integration of many independent ODE systems (Chowdhury et al., 15 Sep 2025).

2. Embedded Runge–Kutta formulation

ode45gpu implements “exactly the same Dormand–Prince embedded Runge–Kutta 4(5) scheme that underlies MATLAB’s ode45.” For

$\dot y = f(t,y), \qquad y(t_0)=y_0,$

the method advances one step with seven stages,

$$k_1 = f(t_n,y_n), \quad \ldots, \quad k_7 = f\!\left(t_n+c_7 h,\; y_n + h\sum_{j=1}^6 a_{7j}k_j\right),$$

and then forms both a fifth-order and a fourth-order update,

$$y_{n+1}^{(5)} = y_n + h\sum_{i=1}^7 b_i k_i, \qquad y_{n+1}^{(4)} = y_n + h\sum_{i=1}^7 \hat b_i k_i.$$

Their difference provides the embedded local error estimate

$e_{n+1} = y_{n+1}^{(5)} - y_{n+1}^{(4)}.$

The paper explicitly lists the Dormand–Prince Butcher tableau, including

$$b_i^{(5)} = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]$$

and

$$b_i^{(4)} = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40].$$

A scalar or weighted norm of $e_{n+1}$ is compared against the user’s absolute and relative tolerances, atol and rtol; accepted steps satisfy $\|e_{n+1}\|\le 1$, and rejected steps are retried with a smaller $h$ (Chowdhury et al., 15 Sep 2025).

The step-size controller is given in the classical form

$\dot y = f(t,y), \qquad y(t_0)=y_0,$0

Within the broader GPU ODE literature, this places ode45gpu alongside other explicit embedded Runge–Kutta strategies, although adjacent packages also emphasize alternative pairs such as Runge–Kutta–Cash–Karp or, for moderate stiffness, Runge–Kutta–Chebyshev (Chowdhury et al., 15 Sep 2025).

3. GPU mapping in MATLAB

The GPU implementation is driven by restrictions and opportunities in MATLAB’s arrayfun model. The published description states that ordinary MATLAB ode45 “builds and manipulates small vectors and employs many temporary allocations,” which is acceptable on a CPU but incompatible with efficient GPU execution under arrayfun, because arrayfun “forbids dynamic array creation inside the kernel.” The implementation therefore “completely unroll[s] the 7-stage RK in scalar form.” Instead of manipulating a 12-component state as a MATLAB vector, it uses twelve separate local variables, and the stage derivatives become scalar temporaries rather than dynamically allocated arrays (Chowdhury et al., 15 Sep 2025).

A single call of the form arrayfun(@ode45gpu, C4Mat, C8Mat, …, y0Mat, rtolMat, atolMat) dispatches one GPU thread per matrix element. Inputs are stored as gpuArray matrices of size ntheta × nphi, and each thread executes the entire Dormand–Prince integration loop for one trajectory. “Inputs are read from global memory; all temporaries (y₁…y₁₂, k-stages, local step-size) go into registers or thread-local memory; results (final y₁…y₁₂, observables) are written back to global memory exactly once.” The same function contains the error estimate, the adaptive controller, and a hard cap on internal stages: “Break out if the number of internal stages ever exceeds a hard cap (one million)—again fully in-thread” (Chowdhury et al., 15 Sep 2025).

The execution model contains “zero inter-thread communication or reduction.” Step-size adaptation, error estimates, bounce counts, lifetime tracking, and related quantities are all local to the thread, so “there is no need for shared memory or explicit synchronization.” MATLAB and the driver map the threads “into blocks/warps with its own heuristics (typically 128–256 threads per block, subject to occupancy), but that is under the hood” (Chowdhury et al., 15 Sep 2025).

4. Interface and workflow in atomiongpu.m

The front end to ode45gpu is atomiongpu.m, which is described as a helper script exposing a pure-MATLAB parameter block. The user edits physical parameters, simulation parameters, and output or parallelization options; “no MEX or CUDA code is written by the user” (Chowdhury et al., 15 Sep 2025).

Category	Parameters stated in the description
Physical parameters	mion, matom, Tion, Tatom, collisiontype, potential, n,m, Cn,Cm, De,Re, ax,ay,az,qx,qy,qz, OmegaRF
Simulation parameters	tmax, maxsteps, r0, ntrajectories, ntheta, nphi, positions, rtol, atol
Parallelization and output	processor, ncores, longest-lived or custom trajectory flags, saveworkspace, savecsvs, onlyonecsv

Internally, the workflow is described in five steps. First, the script “constructs gpuArray matrices for each y₀(i) and constant parameter.” Second, it calls [tend, yend1…yend12, other] = arrayfun(@ode45gpu, …) so that each GPU thread integrates one trajectory. Third, it gathers results back to the host. Fourth, it computes observables including “scattering angles, bounce counts, lifetimes, kinetic-energy changes.” Fifth, it “generates CSVs, heat-maps over (θ,φ), and plots the longest-lived or custom trajectories” (Chowdhury et al., 15 Sep 2025).

The per-thread function signature is also explicitly characterized. Inputs include C4, C8, the trap parameters ax, ay, az, qx, qy, qz, OmegaRF, the time interval t0, tf, the initial-state components y01…y012, and the error-control tolerances rtol, atol. Outputs include tend, the final state components yend1…yend12, and other, described as “a small vector of on-the-fly observables (e.g. complex lifetime, bounce count, number of substeps, final ion KE)” (Chowdhury et al., 15 Sep 2025).

5. Performance characteristics

The published performance data compare several implementations for 2,500 trajectories on one Haswell core under MATLAB R2023a, and then examine GPU scaling on Tesla K80, P100, and V100 devices. The reported CPU timings are: ode45 with a naive force function f(t,y) at “~10 hours,” ode45 with an optimized MATLAB vector-free force g(t,y) at “~3.5 hours (3× speed-up),” and ode45gpu (single-threaded) with the unrolled RKF4(5)” at “~10 minutes—22× speed-up over standard ode45.” Multi-core CPU scaling withparfor` is reported to follow “almost ideally as ∼ ncores^-0.935” (Chowdhury et al., 15 Sep 2025).

Configuration	Reported result
ode45 with naive force	~10 hours for 2 500 trajectories
ode45 with optimized vector-free force	~3.5 hours (3× speed-up)
ode45gpu with unrolled RKF4(5)	~10 minutes—22× speed-up over standard ode45

For GPU throughput, the tests use “n identical copies of the chosen chaotic trajectory” with $\dot y = f(t,y), \qquad y(t_0)=y_0,$1 from $\dot y = f(t,y), \qquad y(t_0)=y_0,$2 to $\dot y = f(t,y), \qquad y(t_0)=y_0,$3. The runtime is reported to be “essentially flat” below $\dot y = f(t,y), \qquad y(t_0)=y_0,$4, because “all GPU SMs [are] saturated.” Above that threshold, runtime grows linearly, with asymptotic costs of “≃0.00008 s/trajectory” on K80, “≃0.00003 s/trajectory” on P100, and “≃0.000025 s/trajectory” on V100. The paper further reports that “uniformly sampled trajectories (different initial angles) give the same asymptotic slope, with only modest scatter for small n” (Chowdhury et al., 15 Sep 2025).

At larger scales, the stated example is “∼10 million trajectories on 8 GPUs in under 15 h (∼23 000 traj/s per GPU).” A plausible implication is that the implementation is primarily throughput-oriented: once the device is saturated, performance is governed less by the control overhead of adaptive stepping than by the aggregate cost of many independent right-hand-side evaluations and per-thread step-control loops (Chowdhury et al., 15 Sep 2025).

6. Relation to broader GPU ODE integration research

ode45gpu belongs to a broader family of GPU ODE solvers for “large number of independent ordinary differential equation systems.” Hegedüs described a “general purpose, modular program package” capable of using professional graphics cards, with numerical schemes given as “the explicit and adaptive Runge--Kutta--Cash--Karp algorithm and the explicit fourth order Runge--Kutta method with fixed time step.” That package deliberately does not store intermediate trajectory points; instead, “with pre-declared device functions, the required special features or properties of a solution can be easily extracted and stored each into a dedicated variable,” and event handling is incorporated to detect special points and to support “non-smooth dynamics---e.g. impact dynamics” (Hegedűs, 2018).

The implementation strategy in the surrounding literature is consistent with the model used by ode45gpu: one thread integrates one independent ODE system. Niemeyer and Sung’s account of “GPU-Based Parallel Integration of Large Numbers of Independent ODE Systems” states explicitly that “one GPU thread (CUDA or OpenCL work-item) integrates one independent ODE system,” emphasizes structure-of-arrays storage

$\dot y = f(t,y), \qquad y(t_0)=y_0,$5

and treats divergence mitigation by grouping systems with similar stiffness or initial conditions. In that literature, nonstiff explicit RKCK on GPU outperforms multithreaded CPU implementations once $\dot y = f(t,y), \qquad y(t_0)=y_0,$6, while moderately stiff RKC can outperform both explicit CPU and implicit CPU codes for large ensembles (Niemeyer et al., 2016).

More general GPU time-integration infrastructures extend the same ideas beyond specialized MATLAB kernels. Balos et al. describe GPU-enabled SUNDIALS data structures and show that in ARKode the “integrator control (step-size logic, tableau storage, error test) lives entirely on the host,” while vector operations and user callbacks operate on GPU data; for explicit ARKODE, Dormand–Prince 4(5) is the default ERK method (Balos et al., 2020). Singh et al. combine Boost.Odeint with OpenFPM through a distributed algebra layer that makes a GPU/MPI Dormand–Prince 4(5) stepper available in a concise template-expression language, with the state never leaving the GPU except for initial setup or final output (Singh et al., 2023). This suggests that ode45gpu occupies the MATLAB-specialized end of a wider design space whose common theme is the efficient parallel integration of many independent or distributed ODE problems on GPU hardware.