Evolutionary Optimisation of Runge–Kutta Methods for Oscillatory Problems

We propose a new strategy for constructing Runge–Kutta (RK) methods using evolutionary computation techniques, with the goal of directly minimising global error rather than relying on traditional local properties. This approach is general and applicable to a wide range of differential equations. To highlight its effectiveness, we apply it to two benchmark problems with oscillatory behaviour: the (2+1)-dimensional nonlinear Schrödinger equation and the N -Body problem (the latter over a long interval), which are central in quantum physics and astronomy, respectively. The method optimises four free coefficients of a sixth-order, eight-stage parametric RK scheme using a novel objective function that compares global error against a benchmark method over a range of step lengths. It overcomes challenges such as local minima in the free coefficient search space and the absence of derivative information of the objective function. Notably, the optimisation relaxes standard RK node bounds ( c i ∈ [ 0 , 1 ] ), leading to improved local stability, lower truncation error, and superior global accuracy. The results also reveal structural patterns in coefficient values when targeting high eccentricity and non-sinusoidal problems, offering insight for future RK method design.