Dear all
I am reading the NEURON book and try to understand the numerical methods in NEURON. In my understanding, when NEURON integrate one step, it first setup a matrix GV=I, where G is the Hines matrix representing the conductance, I is the current vector (including axial current, injected current and current of ion channels) and V is the vector of v_new (membrane voltages in next time step). I think V is the vector of v_new because when calculating current, the equation i=g*(ve) is used, and v in this equation is membrane voltage.
However, in impelementation of NEURON, it seems that after we solve the matrix equation GV=I, we get a vector of dv (i.e. v_newv_old) rather than v_new. Could you please tell me what's wrong with my understanding and why we get dv not v_new after solving GV=I in NEURON? Thank you so much.
Numerical methods in NEURON

 Site Admin
 Posts: 5784
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Numerical methods in NEURON
Where did you find that? If you have the book, please read 3.3 Cables and note that spatial discretization of the cable equation produces a set of ODEs, each of which involves dv/dt at a particular location in space.zyc wrote: ↑Sat Jun 09, 2018 2:24 amIn my understanding, when NEURON integrate one step, it first setup a matrix GV=I, where G is the Hines matrix representing the conductance, I is the current vector (including axial current, injected current and current of ion channels) and V is the vector of v_new (membrane voltages in next time step).
Re: Numerical methods in NEURON
Thank you for replying. I watched the source code of NEURON, and I found that when NEURON computes the ion channels (e.g. hh.c), there are two statements:
so I think the matrix represent the conductance and rhs represent current. When current is computed by i=g*(ve), and chapter 4 of NEURON book introduces how to compute v_new from v_old (numerical integration methods), so I guess the solution of the matrix equation are voltages rather than dv.
Code: Select all
NODERHS(_nd) = _rhs;
NODED(_nd) += _g;

 Site Admin
 Posts: 5784
 Joined: Wed May 18, 2005 4:50 pm
 Location: Yale University School of Medicine
 Contact:
Re: Numerical methods in NEURON
OK, we're talking about different things.
The model is described by the set of ODEs
y' = f(y,t)
(bold terms are vectors).
The simplest method for calculating a(n approximate) numerical solution to this is the explicit Euler method
y(t+dt) = y(t) + dt * f(y,t)
where the first term on the right hand side is the value of y at the current time t, and the second term dt * f(y,t) is the "delta y vector" i.e. the amount by which y changes from time t to time t + dt.
NEURON never calculates a "delta y vector." Instead, by default NEURON uses the implicit Euler method, which requires solving
y(t+dt) = y(t) + dt * f(y(t+dt),t)
Gathering all terms at time t+dt onto one side of the equation, and all terms at time t onto the other, we have
y(t+dt) * (I  dt * f(y(t+dt),t)) = y(t)
(the I here is the identity matrix). Multiplying both sides by the inverse of (I  dt * f(y(t+dt),t)) gives us
y(t+dt) = y(t) * (I  dt * f(y(t+dt),t))^1
The model is described by the set of ODEs
y' = f(y,t)
(bold terms are vectors).
The simplest method for calculating a(n approximate) numerical solution to this is the explicit Euler method
y(t+dt) = y(t) + dt * f(y,t)
where the first term on the right hand side is the value of y at the current time t, and the second term dt * f(y,t) is the "delta y vector" i.e. the amount by which y changes from time t to time t + dt.
NEURON never calculates a "delta y vector." Instead, by default NEURON uses the implicit Euler method, which requires solving
y(t+dt) = y(t) + dt * f(y(t+dt),t)
Gathering all terms at time t+dt onto one side of the equation, and all terms at time t onto the other, we have
y(t+dt) * (I  dt * f(y(t+dt),t)) = y(t)
(the I here is the identity matrix). Multiplying both sides by the inverse of (I  dt * f(y(t+dt),t)) gives us
y(t+dt) = y(t) * (I  dt * f(y(t+dt),t))^1
Re: Numerical methods in NEURON
Thank you, I'm more clear now. I see that in "nrn_fixed_step_thread" function, mainly three steps are done:
1) setup_tree_matrix
2) solve the tree matrix
3) update.
So step 1) is to setup the matrix (I  dt * f(y(t+dt),t)), step 2) is to get the inverse of matrix, is that true? If it's true, why in update step, the operation is vec_v(i) += vec_rhs(i), not vec_v = M*vec_v (M is the inverse of matrix)? That's also why I thought NEURON calculates a "delta y vector"
1) setup_tree_matrix
2) solve the tree matrix
3) update.
So step 1) is to setup the matrix (I  dt * f(y(t+dt),t)), step 2) is to get the inverse of matrix, is that true? If it's true, why in update step, the operation is vec_v(i) += vec_rhs(i), not vec_v = M*vec_v (M is the inverse of matrix)? That's also why I thought NEURON calculates a "delta y vector"
Re: Numerical methods in NEURON
In nrn/src/nrnoc/treeset.c there is a lengthy comment (beginning line 81) that discusses how NEURON is doing the calculation.
Re: Numerical methods in NEURON
Thank you so much!