NMODL SympySolver - cnexp
This notebook describes the implementation of the cnexp
part of the SympySolverVisitor
, which solves the systems of ODEs defined in DERIVATIVE
blocks when these ODEs are linear and independent.
For a higher level overview of the approach to solving ODEs in NMODL, please see the nmodl-odes-overview notebook.
For a more general tutorial on using the NMODL python interface, please see the tutorial notebook.
Implementation
The SympySolverVisitor
for solver method cnexp
does the following:
Get list of all global scope variables from the Symbol Table, as well as any local variables in DERIVATIVE block
For each differential equation in DERIVATIVE block:
Parse equation into SymPy, giving it the list of variables
This gives us a differential equation of the form:
where the function depends on , as well as possibly other variables reprensented by which we assume do not depend on or
Solve equation analytically using sympy.dsolve to give a solution of the form:
where is some function that depends on the value of at time t, the timestep , and the other variables ().
Replace ODE with analytic solution as C code using sympy.printing.ccode
If we fail to find a solution then do nothing - so currently NMODL reverts to using the legacy CNEXP solver routine (same as mod2c or nocmodl)
Pade approximant
There is an option use_pade_approx
which if enabled does the following extra step:
Given the analytic solution :
Expand the solution in a Taylor series in
dt
, extract the coefficientsConstruct the (1,1) Pade approximant to the solution using these Taylor coefficients
Return this approximate solution (correct to second order in ) as C code
(Replacing the exponential with a Pade aproximant here was suggested in sec 5.2 of (https://www.eccomas2016.org/proceedings/pdf/7366.pdf) - since the overall numerical integration scheme in NEURON is only correct to first or second order in , it is valid to expand the analytic solution here to the same order and so avoid evaluating the exponential function)
Implementation Tests
The unit tests may be helpful to understand what these functions are doing - SympySolverVisitor
tests are located in test/visitor/sympy_solver.cpp, and tests involving cnexp
have the tag “[cnexp]
”
Examples
[1]:
%%capture
! pip install nmodl
[2]:
import nmodl.dsl as nmodl
def run_sympy_solver(mod_string, pade=False):
# parse NMDOL file (supplied as a string) into AST
driver = nmodl.NmodlDriver()
AST = driver.parse_string(mod_string)
# run SymtabVisitor to generate Symbol Table
nmodl.symtab.SymtabVisitor().visit_program(AST)
# constant folding, inlining & local variable renaming passes
nmodl.visitor.ConstantFolderVisitor().visit_program(AST)
nmodl.visitor.InlineVisitor().visit_program(AST)
nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)
# run SympySolver visitor
nmodl.visitor.SympySolverVisitor(use_pade_approx=pade).visit_program(AST)
# return contents of new DERIVATIVE block as list of strings
return nmodl.to_nmodl(
nmodl.visitor.AstLookupVisitor().lookup(
AST, nmodl.ast.AstNodeType.DERIVATIVE_BLOCK
)[0]
).splitlines()[1:-1]
Ex. 1
Single constant ODE
[3]:
ex1 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = 4
}
"""
print("exact solution:\t", run_sympy_solver(ex1, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex1, pade=True)[0])
exact solution: m = 4.0*dt+m
pade approx: m = 4.0*dt+m
Ex. 2
Single linear ODE
[4]:
ex2 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = a*m
}
"""
print("exact solution:\t", run_sympy_solver(ex2, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex2, pade=True)[0])
exact solution: m = m*exp(a*dt)
pade approx: m = m*(-a*dt-2.0)/(a*dt-2.0)
Ex. 3
Single linear ODE
[5]:
ex3 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = (minf-m)/mtau
}
"""
print("exact solution:\t", run_sympy_solver(ex3, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex3, pade=True)[0])
exact solution: m = minf-(-m+minf)*exp(-dt/mtau)
pade approx: m = (-dt*m+2.0*dt*minf+2.0*m*mtau)/(dt+2.0*mtau)
Ex. 4
Single linear ODE that can be simplified
[6]:
ex4 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = (minf-m)/mtau - m/mtau - 2*minf/mtau + 3*m/mtau
}
"""
print("exact solution:\t", run_sympy_solver(ex4, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex4, pade=True)[0])
exact solution: m = minf+(m-minf)*exp(dt/mtau)
pade approx: m = (-dt*m+2.0*dt*minf-2.0*m*mtau)/(dt-2.0*mtau)
Ex. 5
Single nonlinear ODE with analytic solution
[7]:
ex5 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = m^3
}
"""
print("exact solution:\t", run_sympy_solver(ex5, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex5, pade=True)[0])
exact solution: m = sqrt(-pow(m, 2)/(2.0*dt*pow(m, 2)-1.0))
pade approx: m = m*(dt*pow(m, 2)-2.0)/(3.0*dt*pow(m, 2)-2.0)
Ex. 6
Single nonlinear ODE (more complicated and not handled yet): unsupported, so does not modify DERIVATIVE block, leaves it to a later visitor pass to deal with
[8]:
ex6 = """
BREAKPOINT {
SOLVE states METHOD cnexp
}
DERIVATIVE states {
m' = exp(m)^2
}
"""
print("exact solution:\t", run_sympy_solver(ex6, pade=False)[0])
print("pade approx:\t", run_sympy_solver(ex6, pade=True)[0])
exact solution: m = 0.5*log(-exp(2*m)/(2*dt*exp(2*m)-1))
pade approx: m = 0.5*(dt*(log(exp(2*m))-2.0)*exp(2*m)-log(exp(2*m)))/(dt*exp(2*m)-1.0)