---
title: 'ode45gpu: GPU-Accelerated ODE Solver'
url: https://www.emergentmind.com/topics/ode45gpu
type: topic
---

# ode45gpu: GPU-Accelerated ODE Solver

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 [2509.12381].

## 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 [2509.12381].

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” [2509.12381].

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 [2509.12381].

## 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$ [2509.12381].

The step-size controller is given in the classical form
$$
h_{\rm new} = 0.8\,h\left(\frac{1}{\|e_{n+1}\|}\right)^{1/5}.
$$
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 [2509.12381].

## 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 [2509.12381].

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” [2509.12381].

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” [2509.12381].

## 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” [2509.12381].

| 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” [2509.12381].

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)” [2509.12381].

## 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 with `parfor` is reported to follow “almost ideally as ∼ ncores^-0.935” [2509.12381].

| 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 $n$ from $1$ to $10^6$. The runtime is reported to be “essentially flat” below $n \simeq 10^3–10^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” [2509.12381].

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 [2509.12381].

## 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” [1810.03931].

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
$$
Y_{\rm global}[i + N_{\rm ode}\cdot j] = y_i^{(j)},
$$
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 $N_{\rm ODE}\sim 10^3–10^4$, while moderately stiff RKC can outperform both explicit CPU and implicit CPU codes for large ensembles [1611.02274].

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 [2011.12984]. 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 [2309.05331]. 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.

Source: https://www.emergentmind.com/topics/ode45gpu