How to choose a solver

How to choose a solver#

Jaxley provides four different solvers for the differential equations: Backward Euler (bwd_euler, the default), Forward Euler (fwd_euler), Exponential Euler (exp_euler) and the Crank-Nicolson solver (crank_nicolson). The solver can be changed as follows:

jx.integrate(..., solver="exp_euler")

Below, we will describe the efficiency, accuracy, and numerical robustness of these solvers. We note that these solvers only impact the way in which the voltage update is performed; they do not impact the ion channels. Because of this, any efficiency difference of the solvers can become insignicant in simulations with many different ion channels.

⚠️ IMPORTANT!
We try to be general here, but different morphologies and different hardware might lead to different conclusions. If you find that our advice is different from your experience, we would appreciate if you would let us know by creating an issue on GitHub.

Single-compartment model#

For single compartment neurons, all solvers should be relatively stable, give similar results, and have similar computational speed. We note that exponential Euler is not yet optimized for network simulations, and will be very slow for them.

Multi-compartment models#

For multi-compartment models, the solver can have a big impact on accuracy and speed.

Forward Euler#

fwd_euler can be very efficient, but it will typically be unstable and explode. Avoid fwd_euler for multicompartment simulations.

Backward Euler and Crank Nicolson#

bwd_euler and crank_nicolson are numerically stable. On CPU, they are typically the most efficient solvers. On GPU, they can also be efficient, but you might get speed-ups with exponential Euler (see below).

Exponential Euler#

exp_euler is also numerically stable. On CPU, it will typically be signicantly slower than bwd_euler or crank_nicolson. On GPU, you might be able to get speed-ups with exp_euler, even with default settings (in particular for small morphologies of <500 compartments and small batchsizes):

v = jx.integrate(..., solver="exp_euler")

If you run many simulations with the same axial_resistivity, length, radius, and capacitance, the exp_euler can give very large speedups on GPU (up to five times faster than bwd_euler or crank_nicolson). To achieve this, precompute the transition matrix once (before all of your simulations)

cell.customize_solver_exp_euler(
    exp_euler_transition=cell.build_exp_euler_transition_matrix(delta_t)
)

and then run the integration loop (which will automatically detect that the transition matrix has been precomputed and will not compute it again):

v1 = jx.integrate(cell, solver="exp_euler")
cell.set("gLeak", 1e-4)
v2 = jx.integrate(cell, solver="exp_euler")

We note that exponential Euler is not yet optimized for network simulations, and can be very slow for them.

Navigation