From a60cc54e43f3369c3d525ac33981da567ca651e9 Mon Sep 17 00:00:00 2001
From: Sam Yates <halfflat@gmail.com>
Date: Thu, 20 Oct 2016 21:08:00 +0200
Subject: [PATCH] Finalize cable computation document.
---
docs/passive_cable/cable.bib | 24 +++-
docs/passive_cable/cable_computation.tex | 168 +++++++++++++++++++++--
docs/passive_cable/makefile | 14 ++
scripts/PassiveCable.jl | 4 +-
4 files changed, 198 insertions(+), 12 deletions(-)
create mode 100644 docs/passive_cable/makefile
diff --git a/docs/passive_cable/cable.bib b/docs/passive_cable/cable.bib
index 879731a0..1f47e75b 100644
--- a/docs/passive_cable/cable.bib
+++ b/docs/passive_cable/cable.bib
@@ -1,11 +1,33 @@
@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}
+}
diff --git a/docs/passive_cable/cable_computation.tex b/docs/passive_cable/cable_computation.tex
index 3700b823..f01259f2 100644
--- a/docs/passive_cable/cable_computation.tex
+++ b/docs/passive_cable/cable_computation.tex
@@ -1,6 +1,7 @@
\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}
diff --git a/docs/passive_cable/makefile b/docs/passive_cable/makefile
new file mode 100644
index 00000000..8756af0b
--- /dev/null
+++ b/docs/passive_cable/makefile
@@ -0,0 +1,14 @@
+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
diff --git a/scripts/PassiveCable.jl b/scripts/PassiveCable.jl
index a5d714c8..76450490 100644
--- a/scripts/PassiveCable.jl
+++ b/scripts/PassiveCable.jl
@@ -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
--
GitLab