The text below is copied word for word from the report. The table it discusses is in part (d) above.
The three non-stiff solvers consistently outperform the stiff solvers at all three tolerance levels. At rtol = 10−10, ode113 achieves 7.71 correct digits in just 0.017 s with 2,239 function evaluations. In comparison, ode23s requires 11.2 s and approximately 1.82 million function evaluations, while achieving only 5.23 correct digits. This makes ode23s roughly 650 times slower, despite producing about 2.5 fewer correct digits.
We note that this difference does not necessarily mean that the stiff solvers are less effective. They perform additional computations, including Jacobian evaluations, matrix factorizations, and Newton iterations, that provide little benefit for this non-stiff problem. In the Pleiades system, the step size is primarily determined by the accuracy required during close encounters between the stars, instead of numerical stability restrictions. As a result, solvers designed for non-stiff problems are better suited to this system and achieve greater accuracy with significantly less computational effort.
Reducing the tolerance by six orders of magnitude leads to a substantial increase in the number of steps required by most solvers. The step counts increase by the following factors:
ode23: 98.7×,ode23t: 100.6×,
ode23tb: 99.9×,ode23s: 99.2×,
ode45: 11.5×,ode113: 3.1×,
ode15s: 7.9×.
These results are consistent with the classical relationship steps ∼ tol−1/(q+1), where q is the order of the numerical method. For a second-order method (q = 2), this relationship predicts a 100-fold increase in the step count, which closely matches the results for ode23, ode23t, ode23tb, and ode23s. For ode45, using q = 4 predicts an increase of approximately 15.8 times, compared with the observed factor of 11.5.
Among all the solvers, ode113 shows the smallest increase in step count, at only 3.1 times. Its ability to adjust the method order allows it to maintain efficiency as the tolerance becomes more restrictive. At rtol = 10−10, it is also the fastest solver in the comparison. Meanwhile, ode45 achieves the highest accuracy, with scd = 8.84, while requiring only 0.052 s of computation time.
At rtol = 10−4, none of the solvers achieves two correct digits, with scd values ranging from 0.61 to 1.82. This indicates that the numerical errors affect even the leading digits of the final positions. The three close encounters involving body 1, shown in the figure in part c, can amplify errors accumulated during integration and contribute to the loss of accuracy.
This behavior highlights an important distinction between local and global error. The specified tolerance controls the estimated local error at each integration step, but does not guarantee the same level of accuracy in the final solution. Because the Pleiades system is sensitive to accumulated numerical errors, a tolerance of 10−4 does not necessarily produce four correct digits. Therefore, the accuracy of each solver should be evaluated by comparing its results with a reliable reference solution.
From rtol = 10−4 to 10−10, ode45 improves from scd = 1.14 to 8.84, about 1.3 digits per decade, and ode23 from 1.82 to 7.84. The stiff solvers gain less: ode23t goes from 0.61 to 4.73, about 0.7 digits per decade, and stays two to four digits behind the explicit solvers at equal tolerance.
The mescd values are higher than the scd values in all 21 runs, with differences ranging from 0.59 to 0.98 digits. Because the two measures in Equation 3 use different components and scaling methods, this difference is not constant.
The mescd measure considers the largest absolute error across all 28 components and scales it by a single constant, 6.2251. In every run, the largest absolute error occurs in a velocity component, most commonly x′5 or x′1. These components tend to have relatively large values and change rapidly throughout the simulation.
In contrast, scd considers only the 14 position components and divides each error by its corresponding reference value. The largest relative error occurs in x1, with a reference value of 0.3706, in most runs, and in x5, with a reference value of 0.3426, in the remaining runs. These are the two smallest reference positions in absolute value, making their relative errors more sensitive to a given absolute error.
As a result, the difference between mescd and scd depends on how the numerical errors are distributed between the position and velocity components. Since this distribution changes across solvers and tolerance levels, the gap between the two measures also varies. Despite these differences, both measures give the same ranking of solvers in this experiment. However, their results could differ substantially if a reference position were close to zero at t = 3.
The solvers ode15s, ode23t, and ode23tb each evaluate the Jacobian only once during the entire integration, regardless of the tolerance. They also perform relatively few LU factorizations, ranging from 112 to 351. This suggests that the Newton iterations converge efficiently for this non-stiff problem and allow MATLAB to reuse the same Jacobian without repeatedly recomputing it. For example, at rtol = 10−10, ode23t performs 81,467 linear solves using only 159 LU factorizations, with 81,498 function evaluations, approximately one per step.
In contrast, ode23s requires considerably more computation. Its Rosenbrock method uses the Jacobian directly in its numerical formula, and the results show that a new Jacobian and LU factorization are computed at every step. At rtol = 10−10, the solver performs 58,702 Jacobian evaluations, 58,702 LU factorizations, and 176,106 linear solves, corresponding to three solves per step, one for each stage.
For this 28-dimensional system, approximating each Jacobian numerically requires approximately 28 additional evaluations of the right-hand-side function. This accounts for much of the approximately 1.82 million function evaluations reported for ode23s. In addition, repeatedly performing LU factorizations, which require O(n3) operations for a dense system of size n, adds substantial computational cost. These factors explain why ode23s is the slowest solver in the comparison.
Providing an analytical Jacobian could reduce the number of function evaluations by avoiding the additional evaluations needed for numerical approximation. However, it would not eliminate the cost of repeated LU factorizations. Therefore, ode23s would likely remain less efficient than the other solvers for this non-stiff problem.
The difference between attempted and accepted steps is largest at rtol = 10−4. At this tolerance, ode45 rejects 21 out of 127 attempted steps, while ode15s rejects 59 out of 522. As the tolerance becomes more restrictive, the number of rejected steps decreases, reaching zero for most solvers at rtol = 10−10.
We can explain this behavior through how the solvers adjust their step sizes. With a larger tolerance, the solvers initially attempt larger steps, which may produce excessive errors during close encounters between the stars. These steps are then rejected, forcing the solvers to reduce their step sizes and try again. With tighter tolerances, the solvers already take smaller steps which makes it less likely that the estimated error will exceed the allowed limit.
The initial step sizes of 10−6, 10−9, and 10−12 are much smaller than necessary at t = 0, where the solution is relatively smooth. However, the adaptive solvers adjust their step sizes within the first few steps such that the small initial values have little effect on the overall computation. Although they may introduce a few additional steps, they also help prevent the solver from starting with a step that is too large. Therefore, the initial step size is unlikely to be a major factor in the performance differences observed in the table.
8.1 9. Comparison with the test set (Table II.6.3):
Only the scd values can be directly compared between the two tables, for the reasons explained in Section 7. At rtol = 10−10, the reference test set reports scd values of 5.72 for DASSL, 7.39 for GAMD, 6.94 for MEBDFI, 7.55 for PSIDE-1, 9.20 for RADAU, and 5.20 for VODE. In comparison, our MATLAB results are 8.84 for ode45, 7.84 for ode23, 7.71 for ode113, 6.04 for ode15s, 5.23 for ode23s, 4.73 for ode23t, and 4.68 for ode23tb.
The MATLAB non-stiff solvers perform well compared with those in the reference test set. Among the reported methods, ode45 is surpassed only by RADAU, while ode113 and ode23 achieve higher scd values than both MEBDFI and PSIDE-1. In contrast, the MATLAB stiff solvers generally produce lower accuracy, with several performing near or below the lower end of the reference results. This is consistent with our earlier observation that non-stiff solvers are better suited to the Pleiades problem.
However, the CPU times from the two tables should not be compared directly because the experiments were conducted on different hardware and used different timing procedures. The reference test set reports the median of ten runs, whereas part d requires the minimum time from five runs.