Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Working with Arbitrary Precision and NDSOLVE

• To: mathgroup at smc.vnet.net
• Subject: [mg43245] Working with Arbitrary Precision and NDSOLVE
• From: "Christos Argyropoulos M.D." <charg at med.upatras.gr>
• Date: Sat, 23 Aug 2003 08:08:20 -0400 (EDT)
• Sender: owner-wri-mathgroup at wolfram.com

Hi,

I have the following problem; I want to test the sensitivity of pulse
coupled oscillators to the working precision of numerical integration
routines.

First a little bit of background:
Pulse coupled oscillators (known also as Spiking Neurons or Integrate and
Fire Neurons) are models of real neurons described by a first order ODE of
the form x'(t) = -K x(t) + I(t), x(0)= xrest  if x(t) < threshold; once x
crosses the threshold it is reset back to the resting level  xrest, and
integration restarts. The crossing of a threshold is referred to as a
"spike"; these models code information in the interspike interval and
information reconstruction is based on analysis of the timings of spike
events. I close the background section by mentioning that these models are
used to 1) fit data on biological neural systems 2) used for embedded signal
processing applications 3) are implemented in biomimetic robotic
applications.

The question :
If I set the option WorkingPrecision->n in NDSolve to arbitrary levels (i.e.
8,10,12,14 etc bits) will the internal algorithm carry all operations in the
stated precision AND generate InterpolatingFunction objects of the same
precision ?

Comments on the Code:
I used StoppingCondition to break numerical integration at spike crossing
times, and then restarted the integration using the new initial condition
(this is the best way to model discontinuous dynamics I could come up with).
The timing of spike events is fished out the InterpolatingFunction object
generated at the moment StoppingCondition becomes true.
The code produces a graph of the actual input (inp[t]) with a reconstruction
of the "spiking history" of the model neuron in the form of a SpikeTime,
Interspike Interval for various working precision levels.

The Code:

Needs["Utilities`FilterOptions`"];
Needs["Statistics`NormalDistribution`"];
Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DiscreteDistributions`"];
Needs["LinearAlgebra`MatrixManipulation`"];
Needs["Statistics`DataManipulation`"];
Needs["Statistics`DescriptiveStatistics`"];
Needs["Statistics`LinearRegression`"];
Needs["Statistics`DataSmoothing`"];
Needs["Graphics`Graphics`"];
Needs["NumericalMath`SplineFit`"];

inp[t_]:=2.Exp[-((t-25)/100)^2]+.5 Exp[-((t-250)/25)^2]+
      Abs[Sin[t/20]]+0.05Random[NormalDistribution[0,1]];

workprec={8,10,12,14,16,18,20,22,24,26,28};
graph=ReleaseHold[
      Map[Hold[Module[{},
              results=Module[{time=0.0,timestart=0.0,sol,pnt={},t},
                    While[timestart<500,
                      sol=First[
                          NDSolve[{x'[t]\[Equal]-2.5 x[t]+ inp[t],
                              x[timestart]\[Equal]-2.},x,{t,timestart,500},
                            StoppingTest->x>0,WorkingPrecision\[Rule]#,
                            MaxSteps\[Rule]30000]];AppendTo[pnt,sol];
                      timestart=sol[[1,2,1,1,2]]];pnt];];


            dummy1=Map[t\[GreaterEqual] #[[1]] && t< #[[2]]&,
                Map[Part[#,1,2,1,1]&,results]];
            func=Map[Part[#,1,2]&,results];

            asd[t_]:=Evaluate[
                Evaluate[Inner[Sequence,dummy1,func,Which]][t]];
            dummy2=Map[Part[#,1,2,1,1,2]&,Drop[results,-1]];

            DisplayTogether[
              ListPlot[
                Inner[List,Drop[dummy2,-1],1./(ListConvolve[{1,-1},dummy2]),
                  List],PlotJoined\[Rule]True,PlotRange\[Rule]{All,All},
                PlotStyle\[Rule]{RGBColor[0,0,1]}],

Plot[inp[t],{t,0,500},PlotRange\[Rule]All],PlotRange\[Rule]All,
              Frame\[Rule]True,
              FrameLabel\[Rule]{StyleForm["Time",FontSlant->"Italic",
                    FontFamily->"Courier",FontWeight->"Bold",
                    FontSize\[Rule]14],

                  StyleForm["Analog Intensity",FontSlant->"Italic",
                    FontFamily->"Courier",FontWeight->"Bold",
                    FontSize\[Rule]14],

                  StyleForm["Number of Spikes "<>ToString[Length[dummy2]-1],
                    FontSlant->"Italic",FontFamily->"Courier",
                    FontWeight->"Bold",FontSize\[Rule]14],

                  StyleForm["InterSpike Interval",FontSlant->"Italic",
                    FontFamily->"Courier",FontWeight->"Bold",
                    FontSize\[Rule]14,FontColor\[Rule]RGBColor[0,0,1]]
                  },
              PlotLabel->
                StyleForm["WorkingPrecision "<>ToString[#]<>"\n",
                  FontSlant->"Italic",FontFamily->"Courier",
                  FontWeight->"Bold",FontSize\[Rule]12,
                  FontColor\[Rule]RGBColor[1,0,1]]]]&,workprec]];

Thanks for reading all the way down.

Regards,

Christos Argyropoulos