Skip to content
Snippets Groups Projects
Commit a60cc54e authored by Sam Yates's avatar Sam Yates
Browse files

Finalize cable computation document.

parent d3f81690
No related branches found
No related tags found
No related merge requests found
@article{bhalla1992,
title = {Rallpacks: a set of benchmarks for neuronal simulators},
author = {Upinder S. Bhalla and David H. Bilitch and James M. Bower},
journal = {Trends in Neurosciences},
volume = {15},
number = {11},
pages = {453--458},
year = {1992},
month = {11},
author = {Upinder S. Bhalla and David H. Bilitch and James M. Bower},
doi = {10.1016/0166-2236(92)90009-W}
}
@article{churchill1937,
title = {The inversion of the Laplace transformation by a direct expansion in series and its application to boundary-value problems},
author = {R. V. Churchill},
journal = {Mathematische Zeitschrift},
volume = {42},
number = {1},
pages = {567--579},
year = {1937},
doi = {10.1007/BF01160095}
}
@article{borwein2009,
title = {Uniform bounds for the complementary incomplete gamma function},
author = {Jonathan M Borwein and O-Yeat Chan},
journal = {Mathematical Inequalities and Applications},
volume = {12},
number = {1},
pages = {115--121},
year = {2009},
doi = {10.7153/mia-12-10}
}
\documentclass[parskip=half]{scrartcl}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[citestyle=authoryear-icomp,bibstyle=authoryear]{biblatex}
\usepackage{booktabs}
\usepackage[style=british]{csquotes}
......@@ -17,6 +18,8 @@
\addbibresource{cable.bib}
\title{Analytic solution to passive cable model}
\author{Sam Yates}
\date{October 20, 2016}
\begin{document}
......@@ -24,7 +27,7 @@
\section{Background}
The Rallpack suite \cite{bhalla1992} is a set of three tests for cable-based neuronal
The Rallpack suite \parencite{bhalla1992} is a set of three tests for cable-based neuronal
simulators, for validation and benchmarking. Validation data is provided with the
distribution of the suite, which at the time of writing can now be found within the
\textsc{Genesis} simulator simulator distribution.\footnote{see \url{http://genesis-sim.org}.}
......@@ -86,7 +89,7 @@ With the model boundary conditions,
\begin{align}
v(x, 0) &= E, \\
\left.\frac{\partial v}{\partial x}\right\vert_{x=0} & = -Ir, \\
\left.\frac{\partial v}{\partial x}\right\vert_{x=L} & = 0,
\left.\frac{\partial v}{\partial x}\right\vert_{x=L} & = 0,
\end{align}
\end{subequations}
where $I$ is the injected current and $L$ is the cable length.
......@@ -116,7 +119,11 @@ by
\section{Solution to the normalized cable equation}
Let $G(x, s)$ be the Laplace transform of $g(x, t)$ with respect to $t$. From
Let $G(x, s)$ be the Laplace transform of $g(x, t)$ with respect to $t$.
If the inverse transform of $G$ exists, it must agree with $g$ almost everywhere
for $t>0$, or for all $t>0$ given the continuity of $g$.
From
\eqref{eq:normcable} and \eqref{eq:normcableinitial},
\begin{equation}
\label{eq:lap}
......@@ -155,19 +162,162 @@ Consequently,
\end{aligned}
\end{equation}
$G(x, s)$ has poles when $s=0$, $m=0$, or when $mL = k\pi i$ for some non-zero
integer $k$.
\subsection*{Inversion}
Sufficient conditions for the inverse transform of $G(x,s)$ to exist
and be representable in series form are as follows
\parencite[][Theorem 4]{churchill1937}):
\begin{enumerate}
\item
$G(x, s)$ is analytic in some right half-plane $H$,
\item
the singularities of $G(x, s)$ are all poles, and
\item for some $k>1$, $|s^k G(x, s)|$ is bounded in $H$ and
on the circles $|s|=\rho_i$, where $\rho_i$ is an unbounded monotonically
increasing sequence.
\end{enumerate}
If in addition all the poles $\sigma_j$ of $G$ are simple, then
the series expansion of the inverse transform is given by
\begin{equation}
g(x,t) = \sum_j e^{\sigma_j t}\, \Res(G;\sigma_j),\quad t>0.
\end{equation}
$G(x,s)$ has poles when $s=0$, $m=0$, or $\sinh mL=0$.
In the following, let
\begin{gather*}
a_k = k\pi /L\\
m_k = a_k i\\
s_k = m_k^2 -1 = -k^2\pi^2/L^2-1.
\end{gather*}
for $k\in\mathbb{Z}$.
The poles are then $0$, $-1$, and $s_k$ for $k\geq 1$.
There are no branch points
arising from $m=\sqrt{1+s}$, as letting $m=-\sqrt{1+s}$ leaves $G$ unchanged.
For $|s+1|>\epsilon$,
\begin{equation}
\begin{aligned}
|s^{3/2}G(x,s)|^2
& \leq (1+\epsilon)^{-1} \left| \frac{\cosh m(L-x)}{\sinh mL} \right|^2
\\
& \leq (1+\epsilon)^{-1} (1+|\coth mL|)^2
\label{eq:gbounds}
\end{aligned}
\end{equation}
which is bounded in the half-plane $\Re(s) \geq -1+\epsilon$.
As $\coth u+iv$ is periodic in $v$, \eqref{eq:gbounds} is also bounded outside
the regions $\{s: m\in D_k\}$, where $D_k$ are non-overlapping $\delta$-neighborhoods of $m_k$
for some fixed small $\delta>0$. These are contained within disks of radius
$2k\pi\delta/L+\delta^2$ about $s_k$ for $k\geq 1$, which then are separated by
circles $|s|=\rho_k$ about the origin for some $\rho_k$ in $(|s_k|, |s_k+1|)$.
\textbf{Pole at $s=0$}. Recalling $m=\sqrt{1+s}$, $m$ and $\sinh mL$ are non-zero in
a neighbourhood of $s=0$, and so
\textbf{Pole at $s=0$}.\\
Recalling $m=\sqrt{1+s}$, $m$ and $\sinh mL$ are non-zero in\\
a neighbourhood of $s=0$, and so the pole is simple and
\begin{equation}
\begin{aligned}
\Res(G, 0) & = - \frac{1}{m}\cdot\left.\frac{\cosh m(L-x)}{\sinh mL}\right|_{s=0}\\
\Res(G; 0) & = - \frac{1}{m}\cdot\left.\frac{\cosh m(L-x)}{\sinh mL}\right|_{s=0}\\
& = - \frac{\cosh (L-x)}{\sinh L}.
\end{aligned}
\end{equation}
%
\textbf{Pole at $s=-1$}.\\
Let $G(x,s)=f(x,s)/h(s)$, where
\begin{equation}
\begin{aligned}
f(x, s) &= -\frac{1}{s}\cosh m(L-x)\\
h(s) &= m\sinh mL.
\end{aligned}
\label{eq:hf}
\end{equation}
Noting that $dm/ds = \frac{1}{2}m^{-1}$,
\begin{equation}
\begin{aligned}
h'(s) &= \frac{1}{2}m^{-1}\sinh mL + \frac{1}{2}L\cosh mL \\
&= \frac{1}{2}L + \frac{1}{2}L + O(m^2) \quad(m\to 0) \\
&= L + O(s+1) \quad(s\to -1).
\label{eq:hprime}
\end{aligned}
\end{equation}
The pole is therefore simple, and
\begin{equation}
\Res(G; -1) = f(x, -1)/h'(-1) = -\frac{1}{L}.
\end{equation}
\textbf{Pole at $s=s_k$, $k \geq 1$}.\\
Taking $h$ and $f$ as above \eqref{eq:hf},
$m_k$ is non-zero for $k\geq 1$ and
\[
h'(s_k) = \frac{1}{2}L\cosh m_kL.
\]
Consequently the pole is simple and
\begin{equation}
\begin{aligned}
\Res(G; s_k)
& = f(x, s_k)/h'(s_k)\\
& = -\frac{2}{s_k L}\frac{\cosh m_k(L-x)}{\cosh m_kL} \\
& = -\frac{2}{s_k L}\frac{\cosh m_kL\cosh m_kx-\sinh m_kL\sinh m_kL)}{\cosh m_kL} \\
& = -\frac{2}{s_k L}\cosh m_k x,
\end{aligned}
\end{equation}
as $\sinh m_k=0$.
In terms of $a_k$,
\begin{equation}
\Res(G; s_k) = \frac{2}{L}\cdot\frac{1}{1+a_k^2}\cdot\cos a_kx.
\end{equation}
The series exapnsion for $g(x, t)$ therefore is
\begin{equation}
g(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{
1+2\sum_{k=1}^\infty \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}.
\label{eq:theg}
\end{equation}
\section{Computation of series}
As $t$ approaches zero, the series \eqref{eq:theg} converges increasingly slowly.
For computational purposes, an estimate on the series residual allows the determination
of stopping criteria for a given tolerance.
Let $g_n$ be the partial sum
\begin{equation}
g_n(x, t) = -\frac{\cosh(L-x)}{\sinh L} + \frac{1}{L}e^{-t}\left\{
1+2\sum_{k=1}^n \frac{e^{-ta_k^2}}{1+a_k^2}\cos a_k x\right\}.
\end{equation}
so that $g(x, t) =\lim_{n\to\infty} g_n(x,t)$. Let $\bar{g}_n = |g-g_n|$ be the
residual. The $a_k$ form an increasing sequence, so
\begin{equation}
\begin{aligned}
\bar{g}_n(x,t)
& \leq \frac{2}{L}e^{-t}\sum_{n+1}^\infty\frac{e^{-ta_k^2}}{1+a_k^2}\\
& \leq \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{1+u^2}\,du\\
& < \frac{2}{L}e^{-t}\int_{a_n}^\infty \frac{e^{-tu^2}}{u^2}\,du.
\end{aligned}
\label{eq:gbar}
\end{equation}
Substituting $u=v/\sqrt{t}$ gives the identity
\begin{equation}
\int_{a}^\infty \frac{e^{-tu^2}}{u^2}\,du
= \frac{1}{2}\sqrt{t}\int_{a^2t}^\infty e^{-v} v^{\frac{-3}{2}}\,dv
= \frac{1}{2}\sqrt{t}\,\Gamma(-\frac{1}{2},a^2 t),
\end{equation}
where $\Gamma(\alpha, z)$ is the upper incomplete gamma function.
For real $\alpha<1$ and $z>0$, \textcite[][Theorem 2.3]{borwein2009} give the upper bound
\[
\Gamma(\alpha, z) \leq z^{\alpha-1} e^{-z}.
\]
Substituting into \eqref{eq:gbar} gives
\begin{equation}
\begin{aligned}
\bar{g}_n(x,t)
& < \frac{1}{L}e^{-t}\sqrt{t}\,\Gamma(-\frac{1}{2},a_n^2 t) \\
& \leq \frac{1}{L}e^{-t}\sqrt{t}\,(a_n^2 t)^{-\frac{3}{2}}e^{-a_n^2t} \\
& = \frac{e^{-t(1+a_n^2)}}{L t a_n^3}.
\end{aligned}
\end{equation}
\printbibliography
\end{document}
SOURCES := cable_computation.tex
LATEXMK := latexmk -e '$$clean_ext=q/bbl run.xml/' -xelatex -use-make -halt-on-error
all: cable_computation.pdf
cable_computation.pdf: cable_computation.tex cable.bib
$(LATEXMK) $<
clean:
for s in $(SOURCES); do $(LATEXMK) -c "$$s"; done
realclean:
for s in $(SOURCES); do $(LATEXMK) -C "$$s"; done
......@@ -86,7 +86,7 @@ end
# computed potential at (x,t) on cable.
function cable(x, t, L, lambda, tau, r, V, I; tol=1e-8)
scale = -I*r*lambda;
scale = I*r*lambda;
if scale == 0
return V
else
......@@ -135,7 +135,7 @@ function rallpack1(x, t; tol=1e-8)
lambda = sqrt(d/4 * RM/RA)
tau = CM*RM
return cable(x, t, L, lambda, tau, r, EM, I, tol=tol)
return cable(x, t, L, lambda, tau, r, EM, -I, tol=tol)
end
end # module PassiveCable
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment