CyRK - Change Log

2026

v0.20.X

v0.20.0 (2026-10-10)

New Features

  • Tsit5: the 5(4) pair of Tsitouras (2011). It costs the same six diffeq calls per step as RK45 and on the benchmark problems needed 7% to 49% fewer calls than RK45 to reach errors of 1e-6 and below (it is not always cheaper: on a Kepler orbit with eccentricity 0.9 it needed 3% to 6% more). Its fourth-order interpolant is built from the step’s own stages, so dense output, t_eval, and events cost no extra calls.

  • Vern7 and Vern8: Verner’s (2010) “most efficient” 7(6) and 8(7) pairs, at 10 and 13 diffeq calls per step. Vern8 needed 7% to 29% fewer calls than DOP853 for relative errors of 1e-6 to 1e-10 on the benchmark problems. Vern7 needed up to about 40% more than DOP853 on most of them and up to 19% fewer on the Pleiades problem. Their interpolants match the order of the method (7 and 8) and need 5 and 7 extra stages, which are only evaluated on steps that build an interpolant.

    • Vern7 uses the interpolant coefficients as Verner corrected them in 2024.

    • Vern8’s error stops improving near a relative error of 1e-12 on long integrations, where rounding in its large weights accumulates. DOP853 and Vern7 keep improving to rtol=1e-13.

  • All three are available from cysolve_ivp (ODEMethod.TSIT5, ODEMethod.VERN7, ODEMethod.VERN8, appended to the enum as 9, 10, and 11), pysolve_ivp, and nbsolve2_ivp (method="Tsit5", "Vern7", or "Vern8"), and support everything the other explicit methods do. The legacy nbsolve_ivp does not provide them.

  • The tableaus, error weights, extra stages, and interpolants were checked against the Runge-Kutta order conditions in high-precision arithmetic (every condition holds to better than 1e-33 for the source coefficients, and to about 3e-32 for the values written to “rk.cpp” with about 36 significant digits). The Tsit5 tableau comes from OrdinaryDiffEq.jl (MIT license, noted in “LICENSE.md”) and its interpolant weights were derived from it. The Verner coefficients are from Verner’s published files.

  • The interpolants are stored as their values at Chebyshev-Lobatto nodes and evaluated with the barycentric formula (“dense.cpp”). Evaluated from its monomial coefficients, which reach 3e5, Vern8’s interpolant would lose about two digits at tight tolerances (7e-11 instead of 1e-11 on a Kepler orbit at rtol=1e-13). With the barycentric form the dense output is as accurate as the steps.

  • Updates and fixes GitHub actions.

Performance

  • “Performance/performance.py” now tracks Tsit5, Vern7, and Vern8 (“cyrk_performance-Tsit5.csv”, “-Vern7.csv”, and “-Vern8.csv”).

Fixes

  • Dense output of a backward integration used the interpolant of the neighboring step (“cysolution.cpp”, CySolverResult::call). Between the step times its values were wrong at the level of the step size (errors from 1e-8 to 3e-5 where forward dense output gave 1e-10), and at t_start and anywhere in the first step it used the placeholder interpolant, returning NaN (0 for LSODA). Every method was affected. Each time is now evaluated with the interpolant of the step that contains it, and backward dense output is as accurate as forward.

  • pysolve_ivp: an event with terminal = False or terminal = 0 ended the integration at its first occurrence, because the value was used as a count of allowed occurrences. As in SciPy and the “Events” documentation, False and 0 now mean a non-terminal event (max_allowed still applies), True and a positive int n end the integration at the n-th occurrence, and a negative value raises a ValueError.

  • LSODA failed with “Required step size is less than spacing between numbers” at the first step no larger than the spacing between numbers at t, a check added in v0.19.2. ODEPACK (and SciPy’s LSODA) take such steps and continue, and they can be temporary: far from t = 0 repeated error test failures can cut the step to the spacing before it grows back, and after a switch to the Adams formulas a stale Lipschitz estimate from the BDF steps can hold the step at the spacing until the next stiffness test, 20 steps later. LSODA now fails only after more than LSODA_MAX_TINY_STEPS (100, “lsoda.hpp”) such steps in a row, which still ends the creep that non-finite rates cause. A step that left t unchanged counts as one too: ODEPACK tests the step size before its error test cuts it further, so its own count misses some of them.

  • A failed step lost its error status when the problem had events: the event check that followed the step replaced the status with NO_ERROR. BDF, LSODA, and Radau then ended with success = False, status 0, and the message “No errors were encountered”, and RK23, RK45, and DOP853 stepped on past the failure until the step limit. Events are no longer checked after a failed step, so the step’s own status (“cysolver.cpp”, CySolverBase::take_step) is reported.

  • The dense LU solve behind BDF, LSODA, and Radau (“c_lu.cpp”, c_dense_lu_solve and c_complex_dense_lu_solve) gave wrong results whenever the factorization pivoted in more than one column. The factorization swaps whole rows, as LAPACK’s dgetf2 does, but the solve applied each interchange during the forward substitution, as LINPACK’s dgesl does for a factorization that leaves the earlier columns of L in place. The solve now applies every interchange first, as LAPACK’s dgetrs does. A problem whose Jacobian has entries far larger than its diagonal (a stiff variable strongly coupled to slow ones) got wrong Newton corrections, so BDF and Radau took many times the steps they needed, or failed with step-size errors, and Radau’s error estimates were wrong.

  • LSODA with lband and uband and no analytic Jacobian could fail or slow down when the banded factorization pivoted. The banded LU factorization (“c_lu.cpp”, c_banded_lu_factor) did not clear the fill-in rows of each column it reached, as LAPACK’s dgbtf2 does, and LSODA’s finite-difference Jacobian only rewrites the band, so the previous factorization’s fill-in entered U through the row interchanges. On a stiff lower-bidiagonal chain with lband=1, uband=0, LSODA failed after 310 steps where the dense Jacobian takes 492. It now takes the same 492.

  • The first step size search (“cysolver.cpp”, CySolverBase::p_calc_first_step_size, used by every method but LSODA) called the differential equation beyond the end of a span shorter than its probe step (at t = 0.014 for a span of 1e-9). The probe now stays inside the span and the first step is capped at the span and at max_step, as in SciPy’s select_initial_step. The search also restores the initial state when it is done instead of leaving the solver at the probe point.

  • Radau grew its step by 9 (the safety factor times MAX_FACTOR) instead of SciPy’s 10 after a step with a zero error estimate, for example when nothing changes. It took 18 steps where SciPy takes 17 to cover 1e10 from a constant state.

  • pysolve_ivp raised an IndexError for an rtol or atol given as a numpy scalar (such as np.float64(1e-6)) and an AttributeError for one given as an int or a list. Any real scalar or array-like is now accepted.

  • BDF called the differential equation once more than it needed to each time it refreshed an analytic Jacobian (jac_ptr) after a failed Newton iteration. Only the finite-difference estimate needs that derivative.

  • LSODA failed at its first step (“The Newton iteration of an implicit method failed to converge”) on a stiff problem started on its slow manifold, such as a fast variable at its quasi-static value. Its first step size comes from f(t0), which then has no fast component, and the Adams functional iteration cannot converge at it. ODEPACK (and SciPy’s LSODA) cuts the step by 4 after each failure and gives up after 10, a factor of about 3e5 in all, which fell short from a stiffness of about 1e6 up (“c_lsoda.cpp”, stoda_handle_corrector_failure). The first step now gets 10 more retries, which bound the step by the Lipschitz constant that the failed iteration measured, so a stiffness of up to at least 1e20 now starts. Problems whose first step ODEPACK starts take the same steps as before.

  • LSODA could hold its step at the Adams stability bound without end: 5.3e-7 for a fast variable of stiffness 1e6 at rtol 1e-9, and 7.5e-12 near the singular center of a planet’s structure integration (the gravity equation’s 2/r). There the corrector converges at roundoff, so it forms no new Lipschitz estimate, and the step selection holds h through its order increase branch, which skips ODEPACK’s stability test and so never sets irflag. The stiffness test then kept the Adams formulas. The hold now sets irflag when the step is at the stability bound, and LSODA switches to BDF. ODEPACK and SciPy’s LSODA hold the step in both cases.

Documentation

  • “Implicit Methods” page: how accurate a user Jacobian has to be, and how one tight relative tolerance tightens the BDF and Radau Newton iterations for every variable.

  • “Implicit Methods” page: where CyRK’s implicit methods deliberately differ from SciPy’s.

  • New “Integration Methods” page: the order, stages, error estimate, interpolant, and memory of every method, guidance on which methods suit which problems, and work-precision measurements (diffeq calls and wall time to reach a given error) on the Lorenz, Arenstorf, Pleiades, Van der Pol, and linear oscillator problems. The script behind the measurements is “Benchmarks/explicit_work_precision.py”.

  • README, “C++ API”, “Numba”, and “Dense Output and t_eval” pages list the new methods and their memory use and interpolant costs.

v0.19.X

v0.19.4 (2026-10-08)

Fixes

  • An event whose root finder failed stopped Python interpreter (or, from C++, threw std::out_of_range out of the solver). CyrkErrorMessages had no message for ARGUMENT_ERROR, OPTIMIZE_SIGN_ERROR, OPTIMIZE_CONVERGENCE_ERROR, or EVENT_SETUP_FAILED, and CySolverResult::update_status looked messages up with .at. These codes now have messages, and the integration ends unsuccessfully with the root finder’s code as its status.

    • New c_get_error_message(code) (“c_common.hpp”, cimport from CyRK.cy.common) looks a message up without throwing, falling back to a generic “Unknown error code” message. CyRK’s own lookups use it; C++ and Cython callers should prefer it over CyrkErrorMessages.at.

  • CyRK.get_error_message(code) raised a KeyError for an integer code and never returned a message. It now returns the code’s message.

v0.19.3 (2026-10-06)

Fixes

  • Reading t or y of a pysolve_ivp result whose Python diffeq raised an error would crash Python with an access violation. WrapCySolverResult.t and y now read the result’s current storage on every access, return empty arrays ((0,) and (num_dy, 0)) when it holds nothing.

v0.19.2 (2026-10-04)

Fixes

  • pysolve_ivp: an exception raised by the Python diffeq or by an event function is now raised to the caller. The C++ solvers call these functions through a noexcept hook, so before this fix Python printed “Exception ignored in: ‘CyRK.cy.pysolver.PySolver.diffeq’” and the integration went on with the previous step’s derivative, reporting success with a wrong solution (an event that raised was counted as having fired). PySolver now stores the exception, stops the solver (status PYTHON_FUNCTION_ERROR, the step in progress finishes on NaN rates without saving data), and raises it once the solver returns; the same holds for a KeyboardInterrupt, which now stops a long integration, and for the dense-output call and call_vectorize of a solution with extra outputs, which evaluate the diffeq again. A diffeq that returns fewer than num_y + num_extra values raises a ValueError instead of being read past its end.

    • CySolverBase::set_external_error(error_code) is the C++ hook behind this: it stops the integration from outside the solver, and the status it sets is not overwritten by the step in progress.

  • LSODA accepted steps taken with NaN rates and reported success with NaN states. Its weighted norms (vmnorm, fnorm, bnorm) used fmax, which drops a NaN argument, so a NaN correction measured zero and passed both the corrector’s convergence test and the local error test. The norms now propagate NaN, so such a step is rejected as with every other method.

  • LSODA no longer creeps forward without end once its step falls to the spacing between numbers at t (which NaN or inf rates force, and which inf rates already did before the fix above, until the step limit set by the available memory). ODEPACK only counts such a step (nhnil) and goes on. CyRK’s LSODA now fails with STEP_SIZE_ERROR_SPACING, as the other methods do.

Tests

  • Tests/D_PySolver_Tests/test_f_pysolve_failures.py: exceptions from the diffeq (every method, at the initial state, with pass_dy_as_arg and args, during a dense call with extra outputs, a KeyboardInterrupt, a reused solver afterwards) and from events reach the caller; a short return raises; NaN and inf rates end every method’s integration unsuccessfully, and LSODA’s saved states stay finite.

v0.19.1 (2026-09-15)

Fixes

  • Reading a solution’s dense output is now safe from several threads at once when the problem has extra outputs. The interpolant does not carry extra outputs, so each read (CySolverResult::call and call_vectorize) evaluates the diffeq once more at the interpolated state, and it did so by borrowing the retained solver’s current-state arrays (offload_to_temp, write y_now and t_now, run the diffeq, read dy_now, load_back_from_temp). Reads stay thread safe only if the diffeq and its pre-eval function also only read their arguments. Solutions of Python diffeqs (pysolve_ivp) still go through the solver’s arrays and must not be read from several threads.

Tests

  • CyRK.cy.parallel_test.run_dense_output_threads_test(num_threads, num_repeats) solves a Lorenz problem with extra outputs once and reads its dense output from several C++ threads at once, counting reads that differ from a serial read. test_dense_output_read_from_several_threads requires zero.

v0.19.0 (2026-09-08)

New

  • c_cysolve_batch (C++, “CyRK/cy/c_parallel.hpp”; cimport from CyRK.cy.parallel): solves a vector of c_CySolveJobs, each with its own CySolverResult, across std::thread workers. The “Parallelization” documentation page provides more information.

    • CyRK.cy.parallel_test.run_parallel_benchmark(num_threads, num_runs, t_end) times a batch of Lotka-Volterra solves so the threading crossover can be measured on any machine. The “Parallelization” page reports the measurements from one machine and the rules of thumb drawn from them (threading pays once time per job times number of jobs reaches a few milliseconds; scaling to many threads needs jobs of roughly 100 us or longer).

Fixes

  • Fixed the macOS wheels on PyPI, which could not be imported unless an OpenMP run time happened to exist at the path used on the build machine (Library not loaded: @rpath/libomp.dylib). The wheels were built outside of cibuildwheel, so nothing bundled the library they linked to. Every earlier macOS wheel back to at least v0.10.1 was affected.

Build

  • Dropped the OpenMP compile and link flags. CyRK never used it. The solvers are nogil and thread safe so it can be parallized externally without the OpenMP support. The wheels no longer link to (or have to bundle) libomp, libgomp, or vcomp140, and a source build on macOS no longer needs brew’s llvm and libomp.

    • CyRK.cy.prange_test was replaced by CyRK.cy.parallel_test, which drives a batch of solves through the new C++ thread helper described under “New”.

  • Merged _build_cyrk.py into setup.py. The cmdclass hook was a leftover from the 2022 move to pyproject.toml; setup.py has declared the extensions since v0.2.4 (the only way setuptools can mark a wheel as platform specific), so the hook was redundant and, through py-modules, was installing _build_cyrk as a top-level module in users’ site-packages.

  • The macOS arm64 wheels are now built by cibuildwheel alongside the Linux and Windows wheels. Every wheel is installed into a fresh environment and CyRK’s built-in solver checks are run before the upload job can start. Free-threaded CPython wheels are skipped because numba does not ship them currently.

  • The wheel workflow only uploads to PyPI for a published release, or when its manual upload input is set, so it can be run against a branch first.

Tests

  • The macOS test workflow now installs CyRK with clang through setup-python, same as linux and windows workflows do (and the same toolchain the wheels are built with). The conda environment and the brew llvm/libomp steps have been removed.

Conda-Forge

  • The recipe no longer lists vcomp14, llvm-openmp, or libgomp in its host or run requirements.

Documentation

  • README: the macOS install section no longer asks for Homebrew’s llvm and libomp; each platform’s default compiler works.

  • “Parallelization” page: explains that it is the module containing the prange loop, not CyRK, that has to be compiled with OpenMP, and lists the flags for MSVC, GCC, and Apple clang.

  • Demos/jupyter_cyhack.py adds OpenMP flags to %%cython cells only where the compiler supports them (Homebrew libomp on macOS; MSVC and GCC out of the box) and honors CYRK_DEMO_NO_OPENMP; without OpenMP the cells still compile with serial prange loops.

v0.18.X

v0.18.1 (2026-08-18)

Dependencies

  • Bumped NumPy upper version pin from <2.4 to <2.6 for both runtime and build.

v0.18.0 (2026-08-17)

New Implicit Integration Methods:

  • Radau: an implicit Runge-Kutta method of the Radau IIA family of order 5. It is a C++ port of SciPy’s scipy/integrate/_ivp/radau.py.

    • Like BDF, it is available from cysolve_ivp, pysolve_ivp, and nbsolve2_ivp, and supports dense output, t_eval, events, extra output, solution reuse, backward integration, and per-variable tolerance arrays.

    • Unlike BDF and LSODA this is a single-step method, so it carries no solution history and changes its step size without having to rescale anything. That makes it the strongest of the three when a stiff problem changes character sharply: on the new Robertson benchmark it reaches the end in 165 steps where BDF needs 215 and LSODA needs 250.

    • Added complex dense LU factorization to “c_lu.hpp(cpp)” (the equivalents of LAPACK’s zgetrf and zgetrs) for Radau’s complex iteration matrix. The real and complex routines now share one implementation rather than duplicating the algorithm.

  • LSODA: an Adams / BDF method that monitors the problem as it integrates and switches between the non-stiff Adams formulas and the stiff BDF formulas on its own.

    • It is built on a lightly modified copy of the C translation of ODEPACK’s LSODA that ships with SciPy, driven one step at a time so that CyRK keeps control of the solution storage, the events, and the t_eval interpolation. See the new Third-Party Code section of “LICENSE.md” for the notices, and please cite ODEPACK if you use this method.

    • Like BDF, it is available from cysolve_ivp, pysolve_ivp, and nbsolve2_ivp, and supports dense output, t_eval, events, extra output, solution reuse, backward integration, and per-variable tolerance arrays.

    • Added three options that only LSODA understands. Passing them to any other method now raises an AttributeError rather than being ignored.

      • min_step - the smallest step size LSODA is allowed to take.

      • lband and uband: the bandwidth of the Jacobian. Declaring a bandwidth on a large system whose variables only couple to their neighbors is a substantial speed up.

    • Two bugs in the upstream C translation of LSODA were corrected in CyRK’s copy (both are noted at the top of “c_lsoda.hpp”): the user’s minimum step size was stored in a local variable rather than in the common block, so the option was silently discarded, and bnorm was handed the banded Jacobian without the row offset that it was written to.

  • BDF: an implicit multi-step method based on the backward differentiation formulas, with the order varying automatically between 1 and 5. It is a C++ port of SciPy’s scipy/integrate/_ivp/bdf.py and includes the accuracy enhancement from the modified (NDF) formulas. This method is ideal for stiff ODEs (an example shows a decrease from 46,451 to 348 steps when compared to RK45).

    • It is available everywhere the existing methods are (cysolve_ivp, pysolve_ivp, and nbsolve2_ivp) and supports everything they do: dense output, t_eval, events, extra output, solution reuse, backward integration, and per-variable tolerance arrays.

    • Each step solves its algebraic system with a simplified Newton iteration. The iteration matrix is factorized with CyRK’s own LU routines and reused across steps until the step size, the order, or a convergence failure forces a refactorization.

    • BDF holds the Jacobian and its factorization as dense matrices, so its memory grows with the square of the number of dependent variables. That footprint is checked against max_ram_MB during setup, which reports a memory allocation error rather than starting a solve that could not finish in reasonable time.

C++ Backend

  • Prepared the solver base classes for implicit integration methods.

    • Moved the tolerances (rtols, atols) and the step size limits (max_step_size, first_step_size) from RKConfig up to ProblemConfig, since every adaptive method needs them. RKConfig is retained as an alias so existing code keeps working, and its long constructor and update_properties overload now live on ProblemConfig.

    • Moved the tolerance parsing (p_setup_error_control) and the first step size estimator (p_calc_first_step_size) from RKSolver to CySolverBase for the same reason.

    • Added CySolverBase::p_finalize_setup, called at the end of setup once the first step size is known. Multi-step methods need it to build their solution history.

    • CySolverDense now refreshes its interpolation order and its step size from the solver every time its state is set, instead of only when it is constructed. Two new virtuals, set_Q_order_max and get_dense_base_y_ptr, cover the differences between an interpolant that works forward from the start of a step and one that works backward from the end of it.

    • A non-positive max_step is now rejected during setup instead of producing an invalid step size clamp.

  • Added an optional analytic Jacobian, JacobianFuncType. cysolve_ivp and baseline_cysolve_ivp take a new jac_ptr argument at the end of their argument lists; it is unused until an implicit method is selected. When it is null, CySolverBase::p_estimate_jacobian builds the Jacobian with adaptive finite differences, following SciPy’s num_jac.

  • Added “c_lu.hpp(cpp)”: dense and banded LU factorization (the equivalents of LAPACK’s dgetrf, dgetrs, dgbtrf, and dgbtrs). These are implemented in CyRK so that it still does not need to link against an external BLAS or LAPACK library.

  • Added the JACOBIAN_IS_SINGULAR, NEWTON_CONVERGENCE_ERROR, and LSODA_INTERNAL_ERROR error codes. They are all added together so that the error code values stay stable.

Fixes

  • Fixed get_method_str in “numba_solver.py” reporting the wrong method name; its integers did not match the ODEMethod enum.

  • Fixed CySolverResult::solve reporting success for a solve that never ran.

  • Fixed the wheel installing a top level Tests package into the user’s site-packages.

  • Fixed every exclude in “MANIFEST.in” silently matching nothing.

  • Fixed a %d format specifier being handed a size_t in “events_test.pyx”.

  • Removed three cached pointers in “rk.cpp” that were never read (l_K_ptr in p_compute_stages and p_estimate_error, l_K_ptr_index_ptr in set_Q_array).

  • CyRK’s x86-64 builds target AVX2 deliberately, but a CPU without it used to fault with an illegal instruction somewhere inside a compiled extension. Importing CyRK now checks for AVX2 first and raises an ImportError that says what happened and how to get a build that runs. The check is fail-open: if the CPU’s capabilities cannot be determined the import proceeds.

  • Added CYRK_NO_AVX2. Setting it builds without AVX2, and leaving it set at run time turns the check above off, so a CPU without AVX2 has a supported route to a working install.

  • Fixed -mavx2 and -mfma being passed to the linker rather than the compiler on Linux, where they had no effect on code generation at all.

  • Silenced MSVC’s C4551 and C4018 on Windows builds.

Documentation

  • Added an “Agreement with SciPy” section to the implicit methods page, with two new figures. The first is a work-precision diagram showing that CyRK and SciPy buy accuracy at the same rate for every implicit method and problem tested.

  • Added “Benchmarks/scipy_implicit_comparison.py”, which regenerates those figures and statistics.

  • Added a new “Implicit Methods” documentation page covering how to choose between BDF and LSODA, the Jacobian, the cost of these methods, and the LSODA-only options.

  • Updated “C++_API.md” with the new cysolve_ivp argument list and the new “c_lu”, “bdf”, “lsoda”, and “c_lsoda” modules.

  • Documented the new error codes in “Status_and_Error_Codes.md”.

  • Added a Third-Party Code section to “LICENSE.md” with the notices for the vendored LSODA code and for the algorithms that are ported from SciPy, along with SciPy’s license. Added the ODEPACK citations to the README.

Tests

  • Added “Tests/H_Implicit_Tests” covering accuracy against analytic solutions, stiff performance, dense output, t_eval, events, extra output, backward integration, tolerance arrays, solution reuse, step size limits, and the LSODA banded Jacobian.

  • Added the new methods to the existing accuracy test suites.

  • Made the trimmed RK23 tests deterministic.

  • Extended the RK23 trimming to “test_a_cysolve_ivp.py”.

  • Replaced diffeq_stiff in “test_a_pysolve_ivp.py” with the Robertson problem, which is actually stiff, and used it to check that an explicit method exhausts a step budget that an implicit method finishes comfortably inside.

Performance

  • Added a stiff problem to the performance stack: “Performance/robertson.py” holds the Robertson chemical kinetics problem, the standard benchmark for stiff solvers. None of the existing benchmark problems are stiff, so none of them showed what the implicit methods are for. On this one the explicit methods are held to tiny steps by stability: over the larger time span RK45 needs about 34,600 steps where BDF needs 109 and LSODA needs 140. It is also available to cysolve_ivp through cytester as differential equation number 11.

  • Added the new integration methods to “Performance/performance.py”.

    • The legacy nbsolve_ivp does not implement the implicit methods, so its column is left blank for them.

    • The two problems with 10,000 dependent variables are skipped for the implicit methods, since a dense Jacobian of that size would need a 10,000 by 10,000 factorization on every step.

  • Noted in “Performance.md” that a stiff problem’s step count is set by stability rather than by the error tolerances, which is the case where an implicit method wins despite its higher per step cost.

Benchmarks

  • Made the integration method selectable in “Benchmarks/CyRK - SciPy Comparison.ipynb” through a new integration_method setting. It was hard coded to RK45; it now accepts any of the six methods and drops nbsolve_ivp from the comparison when an implicit one is chosen, since that solver only implements the explicit Runge-Kutta methods. Figures for a method other than RK45 are saved with the method in their file name so they do not overwrite the existing ones.

  • Added the Robertson problem to “Benchmarks/diffeq_builder.py” as 'robertson' so that the notebook has a stiff problem to compare the implicit methods on. Its middle species sits around 1e-5 for the whole integration, so the plot scales it up by 1e4 to keep it visible.

  • Added event functions for the Robertson problem to cytester, which previously raised NotImplementedError when asked for events on differential equation number 11.

  • Fixed the notebook plotting the nbsolve_ivp result a second time in place of the nbsolve2_ivp result.

v0.17.X

v0.17.2 (2026-08-15)

Fixes

  • Fixed reused solutions silently keeping the previous run’s problem shape.

    • CySolverResult::setup only rebuilt the configuration’s derived properties (num_y, num_dy, capture_extra, …) the first time it ran. Callers are allowed to set the source properties directly, which pysolve_ivp does, so on a reuse those derived properties went stale and the second run silently solved the wrong problem. They are now rebuilt on every setup.

    • The most visible symptom was that changing num_extra on a reused pysolve_ivp solution dropped the extra output entirely. That now works.

    • cysolve_ivp was never affected; it builds its configuration through update_properties, which has always re-initialized.

  • pysolve_ivp now raises an AttributeError when a reused solution is given a y0 with a different number of dependent variables, instead of silently returning garbage. The solver’s dependent variable storage is shared with numpy arrays that were built for the previous run, so changing that size is not supported; build a new solver instead. cysolve_ivp shares no storage with Python objects and keeps its ability to be reused across problems of different sizes.

Documentation

  • Documented in “CySolverResult_Reuses.md” what can and can not change between reuses.

Tests

  • Added “Tests/D_PySolver_Tests/test_e_pysolve_reuse_shape.py” covering reuse with a changed num_y, num_extra, and t_eval, for both pysolve_ivp and cysolve_ivp.

v0.17.1 (2026-02-23)

Package

  • Adds new “numba_c_api_.cpp” to manifest.

v0.17.0 (2026-02-23)

Fixes

  • Fixed issue where events may not be checked in pysolve_ivp.

  • Fixed issue where dense output calls were not working when performing backwards integration.

C++ Backend

  • Changed how dense outputs are constructed so they retain a pointer to the solution instance not the solver.

    • This will allow dense to survive even if the solver is released. The exception is if capture_extra is set to true In that case we need to call the diffeq from dense::call so we must keep the reference (all dense are destroyed if the solver is released).

  • Exposed the force_retain_solver parameter in cysolve_ivp.

    • If set to True (the default) then the solver will be retained at the end of integration, even if it is not necessarily needed. This can avoid issues where the solver is accessed by the user after it has been released.

    • If set to False (and capture_extra = False or capture_dense = False) then the solver’s memory will be released at the end of integration. This can reduce memory overhead particularly if your problem has a large number of dependent y variables.

  • Unrolled most of the RK calculations to improve performance.

nbsolve_ivp Changes

  • Implemented a new nbsolve2_ivp which is a numba-safe wrapper to the CySolver C++ backend. This will greatly improve future development by putting all functionality into C++ and just having a light numba wrapper around it. The new nbsolve2_ivp can have cysolve_ivp functionality (like dense outputs) that the current nbsolve_ivp lacks.

    • The current nbsolve_ivp will be replaced by this new nbsolve2_ivp in a future release (treating it as experimental for now).

    • Added deprecation warnings and a new warnings=True bool to the current nbsolve_ivp.

    • The result of integration using nbsolve2_ivp is different than CyRK’s other solvers and there are some important considerations to keep in mind. Please review the new “Numba” documentation.

Benchmarks

  • Added in a repeater inside cysolve_ivp tester so we can more accurately benchmark the integrators performance without also tracking the python wrappers overhead.

    • This led to a huge increase in performance on the benchmarks. cysolve_ivp is now consistently 100 to 500x faster than scipy solve_ivp`.

  • Added a benchmark for large number of dependent y variables (y=10_000) to cytester and the SciPy Comparison notebook.

Demos

  • Fixed demos so they properly use cross-platform headers (required ipython hack which is included in the new “Demos/jupyter_cyhack.py”).

  • Added demo to show cysolve prange.

Documentation

  • Added resets/restarts to the CySolver examples notebook.

  • Adjusts change log section numbering.

  • Added brief note in documentation about parallelizing cysolve.

  • Added document describing memory usage of the C++ structures.

  • Added document describing performance tips and tricks.

  • Added a new Numba.md documentation file describing nbsolve_ivp and nbsolve2_ivp.

Dependencies

  • Bumps max numpy pins to 2.4.x

  • Adds support for Python 3.14.x

Tests

  • Added tests for nbsolve2_ivp.

  • Turned back on backwards integration tests that were being skipped.

2025

v0.16.X

v0.16.1 (2025-12-10)

Dependencies

  • CyRK now supports numpy 2.4.x

Repository Changes

  • Changed license to Apache 2.0

  • Added citation.cff file

  • Added codemeta.json file

  • Add GitHub Workflow to automatically update version numbers

Tests

  • Retired MacOS 13 tests.

v0.16.0 (2025-09-24)

CyRK now supports Events!

  • New C++ Event class which wraps Cython functions and is integrated into CySolver.

  • pysolve_ivp and cysolve_ivp have new arguments for the user to provide event functions.

  • Added tests and documentation for these new features.

Other New Features

  • Added new CyRK.optimize package with will house various numerical/scientific optimizations tools.

    • To this package, added a Cython wrapped brentq method based off of SciPy’s implementation. (Note this is not properly exposed to python (don’t need to for CyRK but might be a nice thing to add in the future)

Changes

  • Moved Cython wrappers for “common.cpp” to its own file “common.pyx/pxd”; refactored previous “common.cpp/hpp” to “c_common.cpp/hpp”

  • Improved performance of pysolve_ivp.

Fixes

  • Fixed some incorrect types in cysolver_api wrapper.

v0.15.X

v0.15.1 (2025-08-15)

Fixes

  • Fixed issue where cysolver would not properly reset in-between reuse runs unless setup was explicitly called.

v0.15.0 (2025-08-13)

Changes

  • Moved away from using std::optional in integrator configs. It is just not currently supported by cython enough. Causing more headaches than its worth.

    • This led to a small performance improvement for cysolve_ivp for problems with a low number of integration steps.

Tests

  • Added tests to check on the reuse functionality introduced in v0.14.0.

Fixes

  • Fixed some initialization of C++ object properties to be more efficient.

  • Other compile warnings were resolved.

v0.14.X

v0.14.1 (2025-08-13)

  • Minor documentation fix.

v0.14.0 (2025-08-13)

New

  • Reuse mechanism

    • The major change with this version is a “reuse” functionality that allows users to reuse previous ODE result structures to avoid memory allocations.

    • Added demos and documentation to showcase how this new reuse functionality works.

  • Added a demo cell in “1 - Getting Started” notebook to showcase the various functionality of the WrapCySolverResult object.

  • Added ReadTheDocs documentation website at: https://cyrk.readthedocs.io/en/latest.

Changes

  • C++ Backend:

    • Moved to ErrorCodes enum for tracking error codes.

    • Changed how CySolverResult works:

      • Transitioned away from using shared pointers.

      • Instantiating classes requires less arguments. Moving more of the initializations to setups so that these setup methods can be called later on to change parameters without reallocating the entire object.

    • Added more checks to look out for memory allocation problems.

    • Added more checks and error catches to avoid other problems.

    • Optimized solution storage, solver, and dense output to better fit in typical cache sizes.

    • Vectors are now pre-allocated with guesstimated sizes at objection creation to reduce the need for subsequent reallocations.

  • CySolver:

    • WarpCySolverResult.call now reshapes the output array.

    • Other arrays are now passed to cysolve_ivp as vectors. These changes remove the need to pass the size of the vectors to cysolve_ivp.

      • Additional arguments are now passed to cysolve_ivp as vector[char] now.

      • Former y0_ptr is now y0_vec and expects a vector[double].

      • Former t_eval_ptr is now t_eval_vec and expects a vector[double].

      • Former r_tols_ptr is now rtols_vec and expects a vector[double].

      • Former a_tols_ptr is now atols_vec and expects a vector[double].

    • Improved how additional arguments are displayed when WrapCySolverResult.print_diagnostics is called.

  • Other:

    • Updated benchmark plot to show the new reset results.

Fixes

  • CySolver Backend:

    • Fixed potential issue where provided t-eval pointer could stop pointing to correct block of memory.

    • Fixed issue where extra outputs were not getting saved on the first time step.

    • Fixed potential problem where t_eval indices were being converted to ints. For very large arrays there could be a wrap around bug.

    • Moved to one vector for dependent y arrays storage and one vector for dy arrays. Pointers added to respective sections of each vector.

    • Most vectors are now resized and values are changed via assignment rather than being cleared and appended to. This allows for resets to retain the original memory and avoid memory allocations unless absolutely necessary.

  • PySolver Backend:

    • Fixed potential issue where the python module may not have gotten decremented.

    • Removed the PySolver methods from the C++ codebase. Don’t see any reason they can’t live in pure python land.

v0.13.X

v0.13.5 (2025-04-09)

Fixes

  • Missing files in manifest.

v0.13.4 (2025-04-09)

Fixes

  • Fixed *.pyx and *.pxd files not being included in distributions.

v0.13.3 (2025-04-09)

Fixes

  • Fixed issue in manifest.

v0.13.2 (2025-04-09)

Changes

  • Modified GitHub workflows so that tests are only run on commits when asked.

Fixes

  • Fixed cythonized files that were being included with sdist, causing problems with conda-forge.

v0.13.1 (2025-04-08)

Changes

  • Added more debug info during installation.

Fixes

  • Fixed const errors at compile time.

v0.13.0 (2025-04-03)

New & Changes:

  • Added support for Python 3.13 and Numpy 2.x

Fixes

  • Fixes memory leak in CySolver (this was mostly effecting pure cython calls; python calls may not have had an issue. Also only when dense output was on).

    • The fix also led to a fairly big performance boost to pysolve.

  • Fixed error in CySolver tester that could cause crashes or incorrect results (only effected some tests and benchmarks).

v0.12.X

v0.12.2 (2025-03-27)

  • Re-release to trigger conda build that included dependency fixes.

2024

v0.12.X (cont’d)

v0.12.1 (2024-12-02)

Fixes

  • Added fix so that arrays returned by WrapCySolverResult class are not destroyed when class is. Instead they are managed by python’s garbage collector once they are held by a python variable. This addresses Github issues #80 and #78

Tests

  • Added test to check that arrays stay alive after underlying WrapCySolverResult class is destroyed.

v0.12.0 (2024-12-02)

New & Changes:

  • MAX_STEP can now be cimported from the top-level of CyRK: from CyRK cimport MAX_STEP

  • Changed uses of void* to char* for both diffeq additional args and pre-eval outputs. The signature of these functions has changed, please review documentation for correct usage.

  • Added new diagnostic info display tool to cysolve_ivp and pysolve_ivp output that you can access with <result>.print_diagnostics().

Fixes

  • Fixed issue with cysolve_ivp (pysolve_ivp did not have this bug) where additional args are passed to diffeq and dense output is on and extra output is captured.

    • Calling the dense output when extra output is on requires additional calls to the diffeq. However, after integration there is no gurantee that the args pointer is pointing to the same memory or that the values in that memory have not changed.

    • Best case, this could cause unexpected results as new values are used for additional args; worst case it could cause access violations if the diffeq tries to access released memory.

    • Now the CySolverBase makes a copy of the additional argument structure so it always retains those values as long as it is alive. This requires an additional memory allocation at the start of integration. And it requires the user to provide the size of the additional argument structure to cysolve_ivp.

Tests

  • Fixed tests where additional args were not being used.

  • Fixed issue with diffeq test 5.

    • This fixes GitHub Issue #67

  • Added tests to solve large number of diffeqs simultaneously to try to catch issues related to GitHub issue #78.

Documentation:

  • Updated the “Advanced CySolver.md” documentation that was out of date.

v0.11.X

v0.11.5 (2024-11-27)

New

  • Added a steps_taken tracking variable to the C++ class CySolverResult and the Cython wrapped WrapCySolverResult so that users can see how many steps were taken during integration even when using t_eval.

v0.11.4 (2024-11-25)

C++ Changes

  • Moved away from stack to heap allocated vectors for arrays that depend on the size of y0. Prior stack allocation required a hard upper limit which was a limitation for some use cases (See GitHub Issue #74).

v0.11.3 (2024-11-18)

New

  • Added helper function to allow users to include CyRK headers during installation. CyRK.get_include() this returns a list which must be appended to your include list in “setup.py”.

  • Some minor performance improvements when using dense output or t_eval.

  • Added new cysolve_ivp C++ function, baseline_cysolve_ivp_noreturn, and cython wrapped version, cysolve_ivp_noreturn, that require the storage class as input (then modify it) rather than provide it as output.

C++ Changes

  • Major rework of the C++ Solver backend

    • The CySolverResult storage structure is now the owner of: the solution data, the dense outputs, and the solver itself (before the solver was stand alone). This ensures that the solver is alive for the dense outputs which, depending on settings, may need to make calls to the solver even after the integration is complete.

    • Removed many pointers in favor of direct access of variables. This was due to some variables moving and hanging pointers causing crashes.

    • Drastically reduced the size of CySolverDense and moved its state data to a heap allocated array. This has led to large performance gains when dense_output=True for integrations requiring a large number of steps.

  • num_y and num_extra are now size_t type instead of unsigned ints.

  • integration method is not int instead of unsigned int.

Fixes

  • Fixed issue where if t_eval was less than t_end it could cause an access violation.

  • Fixed issue where dense output was sometimes being created twice when t_eval was provided.

  • Addressed several compile warnings.

  • Fixed issue where when t_eval is provided and dense output = true, cysolver may call a dense solution that has moved (hanging pointer).

v0.11.2 (2024-11-12)

New

  • Using extra output (via num_extra) with cysolve_ivp and pysolve_ivp now works when dense_output is set to true. CySolverSolution will now make additional calls to the differential equation to determine correct values for extra outputs which are provided alongside the interpolated y values.

    • Added relevant tests.

Other

  • Refactored some misspellings in the cysolver c++ backend.

Fix

  • Fixed missing np.import_array in cysolver backend.

v0.11.1 (2024-11-11)

Fixes

  • MacOS was not correctly installing openmp version of the cython object files. Fixed this but…

Issues

  • Cython prange was failing (for the array interps) only on MacOS on some versions of Python. Couldn’t figure it out so removing prange for now.

v0.11.0 (2024-11-09)

New

  • WrapCySolverResult result class now provides user access to attribute num_y.

Removed

  • Removed previous cyrk_ode and older version of the CySolver class-based solver.

    • The functionality of cyrk_ode is now handled by the new (as of v0.10.0) pysolve_ivp function.

    • The functionality of CySolver is partly handled by the new (as of v0.10.0) cysolve_ivp function.

      • Note that the new cysolve_ivp is a functional approach. A class based approach like the older CySolver class supported is no longer available but could be easy to implement. If there is interest please create a Github issue for it.

Refactors

  • Refactored the new cysolver and pysolver files to remove “New”. This will break imports based on previous versions.

  • Broke up cysolver and pysolver into different files to isolate each other’s code.

Other

  • Changed the default ordering for diffeq function inputs to follow the scheme dydt(dy, t, y); previously it was dydt(t, y, dy). This affects the cy2nb and nb2cy helper functions.

  • Updated performance module to use new methods over old.

Demos

  • Fixed typo in the type of the mixed-type args container.

  • Updated to work with new refactoring.

Tests

  • Updated tests to use pysolver where cyrk_ode was used.

  • Changed tolerances and other inputs to try to make some tests faster.

Dependencies

  • Tested that CyRK works with numpy v2.X; but a lot of other packages don’t right now. So setting it as upper limit.

  • Tested that CyRK can not work with Python 3.13 yet due to numba dependence. See issue

v0.10.X

v0.10.2 (2024-11-05)

New

  • Added new interpolate_from_solution_list function to help interpolate between multiple solutions across a domain.

Bugs

  • Fixed issue where CyRK.test_cysolver() used incorrect kind and order of arguments.

  • Fixed MacOS compile issues when using OpenMP (for both x86-64 and arm64 macos).

  • Fixed issue where MacOS was failing certain tests.

  • Building new wheels to fix Github issue #62.

Tests

  • Added tests to check all built in testers.

v0.10.1 (2024-07-25)

C++ Backend

  • Changed optional args from double pointer to void pointers to allow for arbitrary objects to be passed in.

    • Added description of this feature to “Documentation/Advanced CySolver.md” documentation and “Demos/Advanced CySolver Examples.ipynb” jupyter notebook.

  • Allow users to specify a “Pre-Eval” function that can be passed to the differential equation. This function should take in time, y, and args and update an output pointer which can then be used by the diffeq to solve for dydt.

cysolve_ivp

  • Change call signature to accept new pre_eval_func function.

  • Added more differential equations to tests.

  • Added tests to check new void arg feature.

  • Added tests to check new pre-eval function feature.

MacOS

  • Going back to GCC for C and C++ compile instead of clang (ran into inconsistent test failures with clang).

v0.10.0 (2024-07-17)

C++ Backend

  • This version of CyRK introduces a major rework of the backend integrator which is now written in pure C++.

  • CySolver is now a Cython wrapper to this C++ integrator which can be accessed via Python.

    • Access this function by using from CyRK cimport cysolve_ivp (this must be done within Cython).

    • The plan is to replace CyRK’s older CySolver with this function.

  • There is now a new PySolver version of this wrapper that allows a user to pass a python function to the C++ backend.

    • Access this function by using from CyRK import pysolve_ivp.

    • This is designed as a drop-in-place replacement for SciPy’s solve_ivp.

    • The plan is to replace CyRK’s older cyrk_ode with this function.

Implemented Dense Output and Improved t_eval for new C++ backend

  • Both pysolve_ivp and cysolve_ivp now utilize a much more accurate interpolator when t_eval is provided.

  • Users can now also request the interpolators be saved with the data, enabling Dense Output functional calls.

  • This closes #45.

    • Note that these improvement was not made for nbsolve_ivp, cyrk_ode, or CySolver methods. See below to learn about these methods’ deprecation.

  • Added tests, documentation, and demos describing these features.

Deprecating Older CyRK Methods

  • The new C++ backend is more flexible, faster, and allows for easy additions of new features. It is common across the cython-only, python, and njit-safe numba solvers. Therefore making a change to it propagates to all three solvers - allowing for easier maintenance and new features. For these reasons, the older cyrk_ode, CySolver, and nbrk_ode are now marked as deprecated. No new features will be implemented for those functions and they will be removed in the next major release of CyRK.

  • Deprecated cyrk_ode

  • Deprecated CySolver

  • Warnings will be issued if these functions are used in this release. To suppress these warnings set raise_warnings to False in the respective function calls.

CySolver

  • Changed error message to use a stack-allocated char-array and associated pointer.

  • Added new argument to constructor raise_warnings (default: True) to allow users to suppress warnings.

cyrk_ode

  • Added new argument to constructor raise_warnings (default: True) to allow users to suppress warnings.

WrapCySolverResult

  • cysolve_ivp and pysolve_ivp now return a class structure that stores the result of the integration along with some meta data. The accessible attributes are:

    • cysolve_ivps CySolverResult:

      • success (bool): Flag if the integration was successful.

      • message (str): Message to give a hint on what happened during integration.

      • error_code (int): Additional error/status code that hints on what happened during integration.

      • size (int): Length of time domain.

      • y (float[:, :]): 2D Float Array of y solutions (and any extra output).

      • t (float[:]): 1D Float Array of time domain at which y is defined.

nbsolve_ivp

  • The older nbrk_ode has been refactored to nbsolve_ivp to match the signature of the new cython-based functions (and scipy’s solve_ivp).

  • The output of nbsolve_ivp is now a named tuple that functions similar to the WrapCySolverResult

Memory Management

  • Changed exit code when memory can not be allocated.

  • Changed some heap allocated arrays in CySolver to be stack allocated

    • This change limits the total number of y-dependent variables and extra output that is tracked to 50. This size is easy to change. If your use case requires a larger size then open an issue and an alternative can be discussed.

  • Converted the underlying storage arrays for CySolver to LinkedLists arrays.

Bug Fixes

  • Fixed issue where the Cython-based solvers might use the incorrect memory freeing function.

Other Changes

  • Moved from GCC to Clang on MacOS builds. There was a new problem that appeared with GCC’s linker and could not find a working solution. The original move away from clang was done to support openMP multiprocessing. CyRK does not currently use that so the switch back should be okay.

Known Issues

  • There is an occasional bug with backwards integration on pysolve_ivp and cysolve_ivp. See Github Issue #56.

v0.9.X

v0.9.0 (2024-05-22)

Major Changes

  • Shifted from using the Python-based PyMem_Alloc, PyMem_Free to c-based malloc, free.

  • CySolver _solve method is now gil-free.

    • This has led to a 35%–230% speed boost at low values of steps (faster start up).

Other Changes

  • CI will now build x64-86 and arm64 wheels for MacOS (change suggested by @dihm in #49).

    • Did have to use this nomkl workaround which may cause problems. TBD.

v0.8.X

v0.8.8 (2024-04-30)

New Features

  • Added utils.free_mem function to free memory so that future changes to the memory allocation system will call the proper free function that works with the utils.allocs.

Changes

  • Changed all instances of PyMem_Free to new free_mem from the utils.

v0.8.7 (2024-04-28)

  • Updated manifest and rebuilt wheels.

v0.8.6 (2024-02-13)

Major Changes

  • Added support for Python 3.12.

  • Converted CySolver’s rk_step method into a pure-c implementation to allow for further optimizations.

  • Changed all files to compile with c rather than c++.

    • Had to change cpp_bools to bints to make this change.

Bug Fixes

  • Fixed issue where CyRK was not installing on MacOS due to issue with LLVM and OpenMP.

    • Have opted to go to gcc for macOS install due to issues with OpenMP and clang.

  • Fixed incorrect type for rk method in CySolver (should eliminate some compile warnings).

  • Fixed issue in benchmark where incorrect results were being displayed for CySolver.

2023

v0.8.X (cont’d)

v0.8.5 (2023-10-27)

Other Changes

  • Improved t-eval setter and resetter in CySolver

Bug Fixes

  • Fixed bug where incorrect memory was being accessed whenever CySolver integration failed (leading to seg fault).

v0.8.4 (2023-10-18)

Performance

  • Removed some try/finally blocks, they were largely not needed in CySolver as allocated memory is released on the class destruction.

  • Cleaned up unused variables in cython-based solvers.

Other Changes

  • Added “force_fail” parameter to CySolver to force the integrator to fail to test if memory is released properly.

  • CySolver class pointers now initialize to NULL at the start of init.

  • CySolver now owns all data that is heap-allocated (via class attributes). This allows better management of data in the event of crashes or integration failures.

Bug Fixes

  • Max number of steps was being performed before extra_output was parsed in CySolver which could lead to incorrect max num steps.

  • Fixed incorrect type for CySolver.user_provided_max_num_steps.

v0.8.3 (2023-10-05)

New Features

  • Added tests to check if memory access violations can occur when CySolver is resolved many times.

  • Added new cdef methods to CySolver for more efficient changes of parameters:

    • CySolver.change_y0_pointer - Changes the y0 pointer without having to pass memoryviews.

    • CySolver.change_t_eval_pointer - Changes the t_eval pointer without having to pass memoryviews.

    • Added way for user to limit RAM usage for cython-bases solvers. This also changed how max number of steps was calculated.

Performance

  • Changing RK variables back to stack-allocated c-arrays rather than malloc arrays.

  • Improved how CySolver and cyrk_ode expected size is calculated and how much it grows with each concat.

  • Files now compile across multiple threads during installation.

Other Changes

  • Moved some common constants for both CySolver and cyrk_ode out of their files and into cy.common.

  • Added more meaningful memory error messages to cython files.

    • Memory allocations (or reallocations) are now performed by helper functions in CyRK.utils.

  • Better future-proofed package structure (mainifests, gitignores, etc.).

  • Converted most Py_ssize_t to size_t.

  • Cleaned up a lot of unused variables and imports.

Bug Fixes

  • Fixed potential memory leaks in cython-based solvers when exceptions are raised.

    • The new safe guards (likely the try/finally blocks) did cause a somewhat sizable hit to performance.

v0.8.2 (2023-09-25)

New Features

  • Added a helper flag to control if CySolver.reset_state is called at the end of initialization.

v0.8.1

New Features

  • Added more interp functions that take pointers as arguments.

Changes

  • Converted interp functions to use each other where possible, rather than having separate definitions.

  • Cleaned up .pxd file formatting.

Performance

  • Moved some backend functionality for CyRK.interp to pure c file for performance gains.

Bug Fixes

  • Fixed issue with “cy/common.pyx” not having the correct cython flags during compilation.

v0.8.0

New Features

  • Added new interp functions that work with c pointers. These can only be cimported.

  • Added new “CyRK.cy.common.pyx” file for functions that are used by both cyrk_ode and CySolver.

    • Moved interpolation functionality into CyRK.cy.common. Restructured cyrk_ode and CySolver to use this new function for interpolations.

Changes

  • Refactored many CySolver internal attributes to reflect to change from memoryviews to pointers. The most important ones for the user are:

    • CySolver.y_new_view -> CySolver.y_ptr

    • CySolver.dy_new_view -> CySolver.dy_ptr

    • CySolver.t_new -> CySolver.t_now

    • CySolver.arg_array_view -> CySolver.args_ptr

  • Changed RK constants back to c arrays initialized with PyMem_Malloc. The memory for these arrays are setup in the cython-based solvers. Afterwards, there are helper functions in CyRK.rk to populate the arrays with correct values.

  • Moved to a more generalized scheme for compiling cython files. See “cython_extensions.json”, “_build_cyrk.py”, and “setup.py” for details.

Performance

  • Transitioned many arrays from numpy to c arrays allocated with PyMem_Malloc, etc. These changes led to a significant performance boost for cython-based solvers.

  • Copied some performance lessons that were learned from the cython-based solvers to the numba-based nbrk_ode.

v0.7.X

v0.7.1

Changes

  • Changed cyrk_ode to match the format used by CySolver for its rk_step.

Performance

  • Minor calculation taken out of tight inner loops in cython-based solvers.

Bug Fixes

  • Added back noexcepts to dabs functions used by cyrk_ode that were mistakenly removed in final dev commit of v0.7.0.

  • Fixed issue where cython-based solvers could overshoot t_span[1].

  • Fixed issue where interp functions would give wrong result when requested x was between x_array[0] and x_array[1].

  • Fixed issue where interp functions would give wrong results if requested x was negative and x_array was positive (or vice versa).

  • The use of carrays for RK constants led to floating point rounding differences that could impact results when step sizes are small.

    • Converted RK constants to numpy arrays which seem to handle the floats much better.

    • Also changed the interaction with these variables to be done solely through constant memoryviews. This may provide a bit of a performance boost.

v0.7.0

Major Changes

  • Added support for Cython 3.0.0

    • Added noexcept to pure cython functions to avoid a potential python error check.

New Features

  • Added the ability to pass arrayed versions of rtol and atol to both the numba and cython-based solvers (cyrk_ode and CySolver).

    • For both solvers, you can pass the optional argument “rtols” and/or “atols”. These must be C-contiguous numpy arrays with float64 dtypes. They must have the same size as y0.

    • Added tests to check functionality for all solvers.

    • This resolves Issue #1.

  • Added new optional argument to all solvers max_num_steps which allows the user to control how many steps the solver is allowed to take.

    • If exceeded the integration with fail (softly).

    • Defaults to 95% of sys.maxsize (depends on system architecture).

  • New CySolver.update_constants method allows for significant speed boosts for certain differential equations.

    • See test diffeqs, which have been updated to use this feature, for examples.

Other Changes

  • Improved documentation for most functions and classes.

  • To make more logical sense with the wording, CySolver.size_growths now gives one less than the solver’s growths attribute.

  • Cleaned up status codes and created new status code description document under “Documentation/Status and Error Codes.md”

  • Fixed compile warning related to NPY_NO_DEPRECATED_API.

  • Converted RK variable lengths to Py_ssize_t types.

  • Changed default tolerances to match scipy: rtol=1.0e-3, atol=1.0e-6.

Performance

  • Various minor performance gains for cython-based solvers.

  • Moved key loops in CySolver into self-contained method so that gil can be released.

  • New CySolver.update_constants method allows for significant speed boosts for certain differential equations.

Bug Fixes:

  • Fixed potential seg fault when accessing CySolver’s arg_array_view.

  • Fixed potential issue where CySolver’s first step size may not be reset when variables that affect it are.

  • Fixed missed declaration in cyrk_ode.

  • Fixed bug where the state reset flag was not being passed from CySolver.solve wrapper method.

v0.6.X

v0.6.2

New Features

  • Added auto_solve key word to CySolver class. This flag defaults to True. If True, then the solver will automatically call self.solve() after initialization.

  • Added new parameter change functions to CySolver so that certain parameters can be changed after the class is initialized for a performance boost.

    • Look for the “self.change_” methods in cysolver.pyx/pxd. There is a main change method, CySolver.change_parameters which allows you to change multiple parameters at once.

Bug Fixes:

  • Fixed issue where CySolver could give incorrect results if the solve() method was called multiple times on the same instance.

  • Removed extraneous code from CySolver.__init__.

  • Changed several cython integer variables to all use Py_ssize_t types. Corrected type conversions.

v0.6.1

New Features

  • Added top level parameters (like MAX_STEP) used in CySolver to cysolver.pxd so they can be cimported.

  • Added new argument to array.interp and array.interp_complex: provided_j the user can provide a j index, allowing the functions to skip the binary search.

  • Added new function interpj to array module that outputs the interpolation result as well as the j index that was found.

Bug Fixes

  • Fixed issue with array tests not actually checking values.

Other Changes

  • Reordered tests since numba tests take the longest.

  • Added initializedcheck=False to the array module compile arguments.

v0.6.0

New Features

  • CyRK now works with python 3.11.

  • Created the CySolver class which is more efficient than the cyrk_ode function.

  • New functions in CyRK.cy.cysolvertest to help test and check performance of CySolver.

Performance

  • Removed python lists from cyrk_ode leading to an increase in performance of 15–20%.

Bug Fixes

  • Fixed compile error with cyrk_ode “complex types are unordered”.

    • This was not a problem before so likely something has changed in newer cython versions.

  • Fixed missing declarations for variables in cyrk_ode.

  • Fixed potential problems during installation where paths may be incorrect depending on OS.

v0.5.X

v0.5.3

Performance

  • Removed dynamic optional arguments from cyrk_ode. Now it checks if those arguments are set to None.

Other Changes

  • Changed cyrk_ode arguments to const to avoid memoryview buffer problems. (Change made by David Meyer)

v0.5.2

  • Fixed issues with MacOS wheel build during CI.

v0.5.1

  • Fixed issues with MacOS wheel build during CI.

v0.5.0

New Features

  • cyrk_ode now supports both float and complex-typed y and dydt functions.

Performance

  • Converted various ints to shorts, chars, or Py_ssize_t. Py_ssize_t is recommended by Cython for loop integers to better support 64-bit architecture.

  • Added custom interpolation functions which, depending on the size of the array, can be up to 10x faster than numpys.

  • Removed unnecessarily variables from cyrk_ode.

  • Had to turn off fastmath for nbrk_ode. See issue 24. This negatively impacted the numba integrator’s performance by around 5%.

Other Changes

  • Refactored, cleaned up, and added comments and docstrings to cyrk_ode.

  • Changed both nbrk_ode and cyrk_ode tests to use pytest parameterization.

  • Changed the accuracy test for both nbrk_ode and cyrk_ode to check against a known function.

  • Added openmp dependence during compile time to allow for the use of prange.

  • Moved cyrk_ode’s Runge-Kutta constants to a separate module CyRK.rk.

Bug Fixes

  • Fixed issue (for nbrk_ode and cyrk_ode) where incorrect step size could be used due to bad minimum step check (see issue 20).

v0.4.X

v0.4.0

New Features

  • Added the ability to save intermediate (non-dependent y) results during integration for nbrk and cyrk ode solver.

    • See Documentation/Extra Output.md for more information.

Performance

  • Minor performance improvements to cyrk_ode (switch to c++ compiler and some related functionality)

  • The new feature that saves intermediate results during integration had a minor impact on performance (even when not using the feature). However, it is within most tests margin of error.

v0.3.X

v0.3.0

Bug Fixes

  • nbrk_ode fixes

    • Improved the storage of results during integration, greatly reducing memory usage. This provides a massive increase in performance when dealing with large time spans that previously required the processor to search outside its cache during integration.

Performance Improvements

  • Various improvements to nbrk_ode make it about 200% faster on small time-spans and over 30x+ faster on large timespans.

  • Improvements to cyrk_ode provided a modest (~5%) performance increase.

Other Changes

  • Helper functions now have an additional optional kwarg cache_njit which is set to False but can be toggled to enable njit caching.

  • Fixed issue in function timing calculation used in the benchmark plot.

v0.2.X

v0.2.4

Bug Fixes
Other Changes

  • Updated CI workflows to utilize cibuildwheel for building binary wheels.

2022

v0.2.X (cont’d)

v0.2.3

Bug Fixes

  • cyrk_ode fixes

    • Bug in doubling up on the time step in the final inter-step diffeq calculation.

v0.2.2

Other Changes

  • Added a performance tracking package to measure CyRK’s performance over time and versions.

v0.2.1

New Features

  • Added helper functions from CyRK import nb2cy, cy2nb which convert differential equation argument signatures between the formats required for cyrk and nbrk ode solvers.

Other Changes

  • Added back some commented out tests that were left over from the bug fixed in v0.2.0.

  • Added tests to check that performs both cyrk and nbrk integrations on larger time domains.

  • Removed the specific test that looked at the underlying issue fixed in v0.2.0 (this is still checked by other tests).

v0.2.0

Bug Fixes

  • Fixed issues with the metadata provided by pyproject.toml.

  • cyrk_ode fixes

    • Fixed bug that was causing the ubuntu slowdown and likely other errors

    • Added a cabs absolute value function to ensure that complex numbers are being properly handled when abs() is called.

  • nbrk_ode fixes

    • Fixed warning during numba integration about contiguous arrays.

    • Fixed issue where variable was referenced before assignment when using nbrk’s DOP853

Performance Improvements

  • cyrk_ode improvements

    • Integrator now selects an output message based on the status code rather than building a string during the integration loop.

    • Switched the loop order for the final list to ndarray conversion. Before the time domain was being redundantly built y_size times.

Other Changes

  • pyproject.toml provides more constrained package list rather than an open search.

  • Added back ubuntu tests and publishes to GitHub workflows.

v0.1.X

v0.1.3

  • Fixed a broken example in the readme documentation.

  • Added better quick tests for both the numba and cython versions.

  • Added SciPy dependency that was missing (required for numba integrator).

  • Increased the lower limit on numpy package version (to fix issue 7) and removed the upper limit version.

  • Removed python 3.7 as a supported version.

  • Updated graphic in readme.

  • Converted over to using pyproject.toml instead of setup.py

    • removed version.py from project folder.

v0.1.2

  • Made the calling argument for the numba solver more consistent with the cython one by letting first_step==0 be equivalent to == None

v0.1.1

  • Corrected issues with installation

  • Improved GitHub workflows

Initial Beta Version 0.1.0 Launched!