Numerical imprecision when number of segments is changed

[briandepasquale]

I'm running a pretty simple model and found what I think must be numerical errors in the integrator. Just wondering if the developers have any suggestions about how to avoid it.

branch = jx.Branch(jx.Compartment(), nseg=nseg)
cell = jx.Cell(branch, parents=jnp.asarray([-1]))
cell.delete_recordings()
cell.branch(0).loc(0.0).record("v")
voltages = jx.integrate(cell, delta_t=0.025, t_max=10.0)
plt.plot(voltages[0], c="b")
plt.ylim(-80, 0)

Running the following code gives a result of -70 for the entire simulation, as it should.

Increasing nseg to 3 causes the voltage to drift to 0. Using fwd_euler prevents this for nseg=3. (Although my understanding is that bwd_euler should generally be preferred over fwd_euler.) Decreasing delta_t and using bwd_euler or CN improves the situation, but there is still some error. Generally increasing nseg (with the default solver and delta_t) reduces the error (e.g. solution is perfect for nseg=4, better but not great for nseg=5,6,7 and back to no error for nseg=13 or 20.

I'm somewhat familiar with numerical methods for solving PDEs, so I'm guessing there is some spatial-temporal tradeoff going on. What is somewhat surprising is how nseg=2 or 4 seem to be outliers here.

Any general guidelines for choosing a default nseg value (assuming a user uses default delta_t and solver) would be appreciated! (Seems like the answer is use nseg=2, but any other guidance would be welcome!)

[michaeldeistler]

Hi there,

thanks for creating this! The default solver for jaxley

voltages = jx.integrate(..., voltage_solver="jaxley.stone")

only supports nseg that are powers of 2, i.e. nseg=1,2,4,8.

If you want other nseg, then please use

voltages = jx.integrate(..., voltage_solver="jaxley.thomas")

I should have done this before but forgot, so will do the following now: 1) change the default to jaxley.thomas and 2) raise a NotImplementedError for nseg=3,5,6,7,9,...

Thanks a lot bringing this up!
Michael

[michaeldeistler]

Hi again,

I just made a new release of tridiax (which jaxley uses for its solvers) to fix this. If you do

pip install tridiax --upgrade

then things should work with any solver.

Sorry for the bug, and thanks a lot for reporting!

Let us know if anything is unclear!
Michael