Check/Output (discretizsation) matrix

Anything that doesn't fit elsewhere.
Post Reply
danarwed
Posts: 6
Joined: Fri Jan 08, 2010 8:05 am

Check/Output (discretizsation) matrix

Post by danarwed »

Hello dear forum users,
is there any possibility to check or print the matrix which is assembled when solving a model with implicit time
integration methods (CN, BE)? I know the Matrix class, which has a print method, but I am not sure whether this datatype
is used for matrices evolving in implicit time stepping schemes.

Why do I want to check the matrix?
I tried to solve a simple model:
Cylindrical cell with 1 passive compartment, no reaction term, potential initially set to
some constant value, alpha-synapse firing after t=1ms, model time t=1ms.
To solve the model I used nseg=16384, dt=0.1ms and the Backward Euler method:

Code: Select all

secondorder = 0
create soma
nsegments=16384

soma {
    nseg = nsegments
    diam = 4
    L = 158.113885
    Ra = 200.
    /*
    // passive current
    insert pas
    g_pas = 0.0005      // [S/cm^2], conductive
    e_pas = 0              // [mv], equilibrium potential
    */
    cm = 2.                 // [uf/cm^2]
}

objectvar syn
soma syn    = new AlphaSynapse(0.0)
syn.tau        = 1.0 
syn.onset    =  1.0
syn.e          =  50
syn.gmax    = 10

dt    = 0.0005
tstop = 10
finitialize(0.00362067)

proc integrate() {
    print soma.nseg
    while(t < tstop) {
        fadvance()
        printf("%g %.14g\n", t, soma.v(1.0))
    }
}
This is the output I get:

Code: Select all

    1   
16384
0 0.00362067
0.1 0.00362067
0.2 0.00362067
0.3 0.00362067
0.4 0.00362067
0.5 0.00362067
0.6 0.00362067
0.7 0.00362067
0.8 0.00362067
0.9 0.00362067
1 0.00362067
1.1 3.8569371857645
1.2 10.221250121198
1.3 16.991046904552
1.4 23.102472056881
1.5 28.258439789191
1.6 32.486941084025
...up to 10 ms
Astonishing for me is that the initial value from t=0.0ms up to t=1ms, when the synapse is activated,
the potential is not changed. If the system is solved numerically, I should at least see some
round-off errors induced by numerics -- or isn't the system solved at all, that is,
is Neuron able to detect the cell to be inactive? As far as I understand the code,
'#define ELIMINATE_T_ROUNDOFF 0' says that round-off errors are not eliminated.
That is why I wanted to check the system matrix for its entries, to understand the numerics.

Hopefully my questions makes sense as I do not now how Neuron solves the system internally.
Many thanks for your help!
Dan
ted
Site Admin
Posts: 6394
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Check/Output (discretizsation) matrix

Post by ted »

Astonishing for me is that the initial value from t=0.0ms up to t=1ms, when the synapse is activated,
the potential is not changed.
With NEURON, when a simulation violates one's expectations, it usually means that either (1) there is a mismatch between the conceptual model (which is in the mind of the modeler) and the computational model (the model in the computer, which is created by executing the modeler's source code), or (2) the modeler's expectations were incorrect.

The model consists of a ladder network of perfect capacitors that are connected at their ungrounded ends by axial resistors. At t = 0, all capacitors are initialized to the same potential. Before the synaptic conductance is activated (i.e. before t = 1 ms), there is no longitudinal potential gradient, and no pathway for charge to escape to ground. Consequently the potentials of all ungrounded nodes remain unaltered until after t = 1 ms.

Which is exactly what the simulation did.

You can verify this for yourself by writing out the ODEs for such a ladder network and solving them analytically. The task might be made easier and insight might more readily be gained if you start with a single compartment, then two compartments, then three, and finally apply inductive reasoning to leap to "as many compartments as you like."
danarwed
Posts: 6
Joined: Fri Jan 08, 2010 8:05 am

Re: Check/Output (discretizsation) matrix

Post by danarwed »

Dear ted,

thank you for your fast reply. Of course, I agree with your statement

Code: Select all

The model consists of a ladder network of perfect capacitors that are connected at their ungrounded ends by axial resistors. At t = 0, all capacitors are initialized to the same potential. Before the synaptic conductance is activated (i.e. before t = 1 ms), there is no longitudinal potential gradient, and no pathway for charge to escape to ground. Consequently the potentials of all ungrounded nodes remain unaltered until after t = 1 ms.

Which is exactly what the simulation did.
since it is clear that the potential should not change if the nodes are ungrounded, no longitudinal gradient exists, no sink/source is activated and the Helmholtz reaction term is switched off. But still, an equation system has to
be solved (from the Mathematician's point of view), which might introduce numerical cancellation or round-off errors, which I expected to detect with an output of 16 digits for the potential. As the computational simulation
gives the analytical value with 100% conformity, I just wondered if there is any system solved or if Neuron is clever enough to detect that 'knothing has to be done'.

Anyway, putting stress on the numerics of my computations it would be useful to be able to have a closer look to any involved matrices -- especially for solving more complcated models.
Is there any way to output the matrix to file or stdout?

Thank you in advance for a short answer.
Best Regards,
Dan
danarwed
Posts: 6
Joined: Fri Jan 08, 2010 8:05 am

Re: Check/Output (discretizsation) matrix

Post by danarwed »

Hello Ted,

I had a look in the source code and found 'nnr_print_matrix' and similar functions (spFileMatrix).
After compiling Neuron with 'nrn_print_matrix' in the solve functions I can have a look at the questioned
matrix.

Anyway thanks for your hints above.
Best regards,
Dan
hines
Site Admin
Posts: 1711
Joined: Wed May 18, 2005 3:32 pm

Re: Check/Output (discretizsation) matrix

Post by hines »

There is no roundoff error in this case because each axial current is identically 0 for the first ms since a compartment v and an adjacent compartment v
are identical. So the equation being solved is dv/dt = 0. There will be round-off errors after t=1ms. You can see this if you replace the AlphaSynapse with
an IClamp (and finitialize(0)). Then sum of cm*v*compartment_area will undoubtedly not be exactly equal to IClamp.amp*IClamp.dur
eg wth nsegment = 1000
c = 0
for (x) c += v(x)*cm(x)*area(x)*(1e-5)
print c, c-1
gives 1 6.3504757e-14
and for 100 and 10000 segments respectively
1 0
1 4.1304737e-12
Post Reply