New analysis for passive cable equation (#842)

* Replace derivation of cylindrical cable equation formula with S–L approach.
* Relate analysis of cylindrical form to Rallpack 1 in a separate section.
* Add analysis of tapered cable case.
* Unify discussion of cable with constant electrical properties but varying geometry.
* Add section describing how one can well-estimate the inter-cell flux given point samples or CV means for the membrane voltage.

.gitignore

@@ -44,12 +44,10 @@
 *.fls
 *.bbl
 
-
 # cmake
 CMakeFiles
 CMakeCache.txt
 cmake_install.cmake
-Makefile
 
 # mechanisms generated from .mod files
 mechanisms/multicore/*.hpp

doc/math/passive_cable/makefile → doc/math/cable_equation/Makefile

-SOURCES := cable_computation.tex
+SOURCES := cable_equation.tex
 
-LATEXMK := latexmk -e '$$clean_ext=q/bbl run.xml/' -xelatex -use-make -halt-on-error
+LATEXMK := latexmk -e '$$clean_ext=q/bbl xdv thm run.xml/' -xelatex -use-make -halt-on-error
 
-all: cable_computation.pdf
+all: cable_equation.pdf
 
-cable_computation.pdf: cable_computation.tex cable.bib
+cable_equation.pdf: cable_equation.tex # cable.bib
 	$(LATEXMK) $<
 
 clean:

doc/math/passive_cable/cable.bib → doc/math/cable_equation/cable.bib

@@ -31,3 +31,40 @@
     year = {2009},
     doi = {10.7153/mia-12-10}
 }
+
+@article{sorolla2013,
+    title = {Algorithm to calculate a large number of roots of the cross-product of {B}essel functions},
+    author = {Edén Sorolla and Juan R. Mosig and Michael Mattes},
+    journal = {{IEEE} Transactions on Antennas and Propagation},
+    volume = {61},
+    number = {4},
+    month = {4},
+    year = {2013},
+    doi = {10.1109/TAP.2012.2231929}
+}
+
+@misc{nistdlmf,
+    title = {{NIST} Digital Library of Mathematical Functions (Release 1.0.23)},
+    editor = {F. W. J. Olver and A. B. {Olde Daalhuis} and D. W. Lozier and B. I. Schneider and R. F. Boisvert and C. W. Clark and B. R. Miller and B. V. Saunders},
+    url = {http://dlmf.nist.gov},
+    urldate = {2019-06-15}
+}
+
+@article{lommel1879,
+    title = {Zur Theorie der Bessel'schen Functionen},
+    author = {E. Lommel},
+    journal = {Mathematische Annalen},
+    volume = {14},
+    issue = {4},
+    year = {1879},
+    pages = {510--536}
+}
+
+@book{reid1980,
+    title = {{S}turmian theory for ordinary differential equations},
+    author = {William Thomas Reid},
+    year = {1980},
+    publisher = {Springer-Verlag},
+    location = {New York},
+    doi = {10.1007/978-1-4612-6110-0}
+}

doc/math/cable_equation/cable_equation.tex

+\documentclass[parskip=half]{scrartcl}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage[citestyle=authoryear-icomp,bibstyle=authoryear]{biblatex}
+\usepackage{booktabs}
+\usepackage[style=british]{csquotes}
+\usepackage{fontspec}
+\usepackage{mathtools}
+\usepackage{ntheorem}
+\usepackage{siunitx}
+
+\usepackage{tikz}
+\usetikzlibrary{angles,calc,intersections,quotes}
+
+\usepackage[p,osf]{cochineal}
+\usepackage[cochineal,vvarbb]{newtxmath}
+\usepackage[cal=boondoxupr,calscaled=1.0]{mathalfa}
+
+\newfontfamily\dispfamily{Cabin}[Scale=0.95]
+\addtokomafont{disposition}{\dispfamily}
+
+\MakeAutoQuote{«}{»}
+\DeclareMathOperator{\Res}{Res}
+
+\addbibresource{cable.bib}
+
+\newcommand{\Int}[2]{\int_{#1}^{#2}\!}
+\newcommand{\D}{\mathop{}\!d}
+
+\newcommand{\cZ}{\mathcal{Z}}
+\newcommand{\cC}{\mathcal{C}}
+
+\theoremstyle{nonumberplain}
+\newtheorem{lemma*}{Lemma}
+\newenvironment{proof}{\textbf{Proof.}}{\hfill$\square$}
+
+\title{Analytic solutions to the cable equation}
+
+\author{Sam Yates}
+%\date{July 29, 2019}
+
+\sisetup{mode=text}
+
+\begin{document}
+
+% Work-around for siunitx vs newtxmath bug
+\ExplSyntaxOn
+    \cs_undefine:N \c__siunitx_minus_tl
+    \tl_const:Nn \c__siunitx_minus_tl { - }
+\ExplSyntaxOff
+
+\maketitle
+
+\section{The cable equation}
+
+The cable equation describes the evolution of the potential
+on a long, thin, conducting cable in a conducting medium, separated
+by a leaky dielectric. It assumes that the behaviour
+can be modelled entirely as a one dimensional problem,
+in terms of the linear conductivity of the cable $\sigma$,
+and the conductance $g$ and capacitance $c$ per unit length of the
+dielectric. The potential $v(x, t)$ then satisfies
+\begin{equation}
+    (\sigma v')' = c \dot v + g v.
+\end{equation}
+where $f'$ denotes the derivative of $f$ with respect to the first
+variable, and $\dot{f}$ the derivative with respect to the second.
+The quantities $\sigma$, $g$, and $c$ are functions of position,
+and will depend on the electrical and geometrical properties of
+the system.
+
+For a cable of radius $r(x)$ with constant bulk resistivity
+$R_L$, areal capacitance $C_M$ and areal surface resistivity
+$R_M$, these parameters are given by
+\begin{align}
+    \sigma(x) &= \pi r(x)^2 / R_L, \\
+    g(x) &= 2 \pi r(x) \sqrt{1 + r'(x)^2} / R_M, \\
+    c(x) &= 2 \pi r(x) \sqrt{1 + r'(x)^2} C_M.
+\end{align}
+Letting $\theta(x) = \arctan r'(x)$, the cable equation becomes
+\begin{equation}
+    \label{eq:constelec}
+    \frac{R_M}{2 R_L}\left(
+        2\sin\theta\cdot v' + r \cos\theta\cdot v''
+    \right) =
+    R_M C_M \dot v + v.
+\end{equation}
+
+\section{Analytic solutions for $v(x, t)$}
+
+Two important special cases for the cable equation are the cylinder,
+with $r(x)$ constant and $\theta(x)=0$, and the conical frustum,
+with $\theta(x)$ constant and $r(x)=x\tan\theta$.
+
+Given constant electrical properties, as in \eqref{eq:constelec},
+a zero voltage at $t=0$, and Neumann boundary conditions, what
+is $v(x, t)$? The Rallpack 1 model \autocite{bhalla1992} is
+an example of such a passive cable model, and is discussed in
+Appendix \ref{ap:rallpack}.
+
+\subsection{The cylinder}
+
+Suppose the cable lies on the interval $[0, B]$, with initial
+voltage zero, zero current
+at $x=0$, and fixed injected current $I$ at $x=B$. The boundary
+conditions are then $v(x,0) = 0$, $v'(0, t) = 0$, and $v'(B, t) = -E$,
+where $-E = I\cdot R_L/\pi r^2$.
+
+A change of variables gives
+\begin{equation}
+    v(x, t) = -\lambda E\cdot u(x/\lambda, t/\tau),
+\end{equation}
+where
+\begin{equation}
+    \lambda = \sqrt{\frac{R_M r}{2 R_L}}, \quad \tau = R_M C_M, \quad b = B/\lambda,
+\end{equation}
+and $u$ satisfies
+\begin{equation}
+    \begin{aligned}
+        u''(x, t) &= \dot u(x, t) + u(x, t),\\
+        u(x, 0) &= 0,\\
+        u'(0, t) &= 0,\\
+        u'(b, t) &= 1.
+    \end{aligned}
+    \label{eq:ucyl}
+\end{equation}
+
+The corresponding time-invariant problem has solution $\hat u(x)$ where
+\begin{equation}
+    \hat u'' = \hat u,\quad \hat u'(0) = 0,\quad \hat u'(b) = 1,
+\end{equation}
+giving
+\begin{equation}
+    \hat u(x) = \frac{\cosh x}{\sinh b}.
+\end{equation}
+
+Consider separable transient solutions $u(x,t)-\hat u(x)$, with zero derivative
+at $x=0$ and $x=b$, of the form $\psi(x)\eta(t)$. From \eqref{eq:ucyl},
+$\eta(t)=e^{-\lambda t}$ for some $\lambda$ and $\psi$ satisfies
+\begin{gather}
+    \label{eq:efncyl}
+    \psi'' + (\lambda-1)\psi = 0,\\
+    \label{eq:efncylbc}
+    \quad \psi'(0)=\psi'(b)=0.
+\end{gather}
+This is a regular Sturm-Liouville problem
+\autocite[see e.g.][Theorem~2.1, p.~146]{reid1980} with discrete non-negative eigenvalues
+$\lambda_k-1$ and corresponding orthogonal eigenfunctions $\psi_k$,
+\begin{equation}
+    \langle \psi_k, \psi_{k'} \rangle =
+    \Int{0}{b} \psi_k(x)\psi_{k'}(x) \D x =0\quad\text{if $k\neq k'$}.
+\end{equation}
+
+The least eigenvalue corresponds to $\lambda_0 = 1$ with eigenfunction
+$\psi_0(x)=1$.
+
+For $\lambda>1$, \eqref{eq:efncyl} gives
+$\psi(x) = c_1\cos\sqrt{\lambda-1}x + c_2\sin\sqrt{\lambda-1}x$. From $\psi'(0)=0$,
+we have $c_2=0$. From $\psi'(b)=1$, it follows $c_1=0$ or $\sqrt{\lambda-1}b = k\pi$
+for some positive integer $k$. Consequently we have eigenvalues and eigenfunctions,
+\begin{gather}
+    \label{eq:lkcyl}
+    \lambda_k = 1+\left(\frac{k\pi}{b}\right)^{\mathrlap{2}},\\
+    \psi_k(x) = \cos \frac{k\pi x}{b}.
+\end{gather}
+for integers $k\geq 0$, including the constant $\lambda=1$ case.
+
+The time-dependent solution $u(x, t)$ then is
+\begin{equation}
+    u(x,t) = \hat u(x) - \sum_{k=0}^{\infty} a_k e^{-\lambda_k t}\psi_k(x),
+\end{equation}
+with $a_k$ determined by the initial conditions,
+\begin{equation}
+    0 = \langle u(x,0), \psi_k(x) \rangle
+    = \langle \hat u(x), \psi_k(x) \rangle
+    - a_k\langle \psi_k(x), \psi_k(x)\rangle,
+\end{equation}
+with
+\begin{align}
+    \langle \hat u(x), \psi_k(x) \rangle
+    = \frac{1}{\sinh b}\Int{0}{b}\cosh x\cos\frac{k\pi}{b}x \D x
+    = (-1)^k\lambda_k^{-1},
+    \intertext{and}
+    \langle \psi_k(x), \psi_k(x) \rangle
+    = \Int{0}{b}\cos^2 \frac{k\pi}{b}x \D x
+    =
+    \begin{cases}
+        b & \text{if $k=0$,}\\
+        \frac{b}{2} & \text{otherwise.}\\
+    \end{cases}
+\end{align}
+Consequently,
+\begin{equation}
+    \label{eq:cylseries}
+    u(x,t) = \frac{\cosh x}{\sinh b} - \frac{1}{b}e^{-t} -
+    \frac{2}{b}\sum_{k=1}^{\infty} (-1)^k \lambda_k^{-1} e^{-\lambda_k t} \cos\frac{k\pi x}{b},
+\end{equation}
+with $\lambda_k$ given in \eqref{eq:lkcyl}. A discussion on the residuals in this series can
+be found in Appendix~\ref{ap:cylcomp}.
+
+Note that the time derivative $\dot u(x,t)$ can be expressed in terms of Jacobi theta functions:
+\begin{equation}
+    \begin{aligned}
+        \dot u(x,t) &= \frac{1}{b}e^{-t}\left(1+2\sum_{k=1}^{\infty} (-1)^ke^{-\frac{\pi^2 k^2 t}{b}}\cos\frac{k\pi x}{b}\right)\\
+        &= \frac{1}{b}e^{-t}\theta_4\left(\frac{\pi x}{2 b}\middle| \frac{i\pi t}{b^2}\right).
+    \end{aligned}
+\end{equation}
+
+\subsection{The tapered cable}
+
+Consider the cable as a conical frustum on the interval $[A, B]$ with radius
+$r(x) = x \tan\theta$ and boundary conditions $v'(A, t) = 0$, $v'(B, t) = -E$,
+and $v(x, 0) = 0$ (see Figure \ref{fig:tapered}). A fixed injecting current $I$
+at $x=B$ corresponds to $-E=I\cdot R_L/\pi r(B)^2$. Reparameterizing gives
+\begin{equation}
+    v(x, t) = -\lambda E\cdot u(x/\lambda, t/\tau),
+\end{equation}
+where
+\begin{equation}
+    \lambda = \frac{R_M\sin\theta}{2 R_L}, \quad \tau = R_M C_M, \quad a = A/\lambda, \quad b = B/\lambda,
+\end{equation}
+and $u$ satisfies
+\begin{equation}
+    \label{eq:conu}
+    \begin{aligned}
+        x u''(x, t) + 2 u'(x, t) &= \dot u(x, t) + u(x, t),\\
+        u(x, 0) &= 0,\\
+        u'(a, t) &= 0,\\
+        u'(b, t) &= 1.
+    \end{aligned}
+\end{equation}
+
+\begin{figure}[tbh]
+    \begin{tikzpicture}[scale=0.8]
+        \coordinate (origin) at (0,0);
+        \coordinate (b0) at (15,0);
+        \coordinate (b1) at ($ (b0) + (0,1.8) $);
+        \coordinate (a0) at (5,0);
+
+        \path [name path=tangent] (origin) -- (b1);
+        \path [name path=avert] (a0) -- +(0,2);
+        \path [name intersections={of=tangent and avert, by=a1}];
+
+        \coordinate (a-1) at ($ (a1)!2!(a0) $);
+        \coordinate (b-1) at ($ (b1)!2!(b0) $);
+
+        \draw [name path=axis, dashed] (origin) -- (b0);
+
+        \draw let \p1 = (a1), \n2 = {0.2*\y1} in
+            (a1) arc (90:270:\n2 and \y1);
+
+        \draw let \p1 = (a1), \n2 = {0.2*\y1} in
+            [densely dotted] (a-1) arc (-90:90:\n2 and \y1);
+
+        \draw let \p1 = (b1), \n2 = {0.2*\y1} in
+            (b0) ellipse (\n2 and \y1);
+
+        \draw (a1) -- (b1) (a-1) -- (b-1);
+        \draw [dotted] (origin) -- (a1) (origin) -- (a-1);
+
+        \pic [draw, "$\theta$", angle radius = 15ex, angle eccentricity=1.2] {angle = b0--origin--b1};
+
+        \coordinate (lright) at ($ (b1)!2!(b0) - (0,40pt) $);
+        \coordinate (lleft) at ($ (origin) + (lright) - (b0) $);
+
+        \path (lleft) node (x0label) { $x=0$ };
+        \path (lright) node (xblabel) { $x=B =b\lambda $ };
+        \path ($ (lleft)!(a0)!(lright) $) node (xalabel) { $x=A =a\lambda$ };
+
+        \draw [dotted] ($ (a-1)-(0,2ex) $) -- ($ (xalabel)+(0,2ex) $);
+        \draw [dotted] ($ (b-1)-(0,2ex) $) -- ($ (xblabel)+(0,2ex) $);
+        \draw [dotted] ($ (origin)-(0,2ex) $) -- ($ (x0label)+(0,2ex) $);
+    \end{tikzpicture}
+    \caption{Tapered cable coordinates.}
+    \label{fig:tapered}
+\end{figure}
+
+Proceeding similarly to the cylinder problem, the solution can be presented as a sum of
+a time-invariant solution $\hat u$ and a series of transients that are solutions to
+a Sturm--Liouville problem.
+
+The time-invariant problem is
+\begin{equation}
+    \begin{gathered}
+        x\hat u''(x) + 2\hat u'(x) - \hat u(x) = 0,\\
+        \hat u'(a) = 0,\quad \hat u'(b) = 1.
+    \end{gathered}
+\end{equation}
+
+Multiplying by $x$ gives a differential equation of the form
+\begin{equation}
+    \label{eq:conicmu}
+    x^2 g''(x) + 2x g'(x) + \mu x g(x) = 0
+\end{equation}
+for some $\mu$. This has solutions in terms of Bessel functions.
+If $\cC_\nu$ is a solution to the Bessel equation
+\begin{equation}
+    z^2 f''(z) + z f'(z) - (z^2 - \nu^2) f(z) = 0,
+\end{equation}
+then $g(x) = x^{-\frac{1}{2}}\cC_1(\alpha x^\frac{1}{2})$ satisfies
+\begin{equation}
+    \label{eq:besg}
+    x^2 g''(x) + 2x g'(x) + \frac{\alpha^2}{4} x g(x) = 0.
+\end{equation}
+Similarly, if $\cZ_\nu$ is a solution to the modified Bessel equation
+\begin{equation}
+    z^2 f''(z) + z f'(z) - (z^2 + \nu^2) f(z) = 0,
+\end{equation}
+then $g(x) = x^{-\frac{1}{2}}\cZ_1(\alpha x^\frac{1}{2})$ satisfies
+\begin{equation}
+    \label{eq:modbesg}
+    x^2 g''(x) + 2x g'(x) - \frac{\alpha^2}{4} x g(x) = 0.
+\end{equation}
+
+The time-invariant solution $\hat u(x)$ satisfies \eqref{eq:modbesg} with $\alpha=2$, and so
+can be written as
+\begin{equation}
+    \hat u(x) = \frac{1}{\sqrt{x}} \left(c_1 I_1(2\sqrt{x}) - c_2 K_1(2\sqrt{x})\right),
+\end{equation}
+with derivative
+\begin{equation}
+    \hat u'(x) = \frac{1}{x} \left(c_1 I_2(2\sqrt{x}) + c_2 K_2(2\sqrt{x})\right).
+\end{equation}
+Solving for the boundary conditions $\hat u'(a)=0$, $\hat u'(b)=1$ then determines the coefficients
+$c_1$ and $c_2$, giving
+\begin{equation}
+    \label{eq:conuhat}
+    \hat u(x) =
+    \frac{b}{\sqrt{x}}\cdot
+    \frac{K_2(2\sqrt{a})I_1(2\sqrt{x}) + I_2(2\sqrt{a})K_1(2\sqrt{x})}
+    {K_2(2\sqrt{a})I_2(2\sqrt{b})-I_2(2\sqrt{a})K_2(2\sqrt{b})}.
+\end{equation}
+
+Separable transient solutions $u(x,t)-\hat u(x)$ will be of the form $e^{-\lambda t}\psi(x)$,
+where $\psi(x)$ satisfies
+\begin{gather}
+    \label{eq:conpsi1}
+    x\psi''(x) + 2\psi'(x) + (\lambda-1)\psi(x) = 0\\
+    \label{eq:conpsi2}
+    \psi'(a) = \psi'(b) = 0.
+\end{gather}
+Multiplying \eqref{eq:conpsi1} by $x$ gives
+\begin{equation}
+    \label{eq:conpsils}
+    x^2\psi''(x) + 2x\psi'(x) + (\lambda-1)x \psi(x) =
+    \frac{d}{dx}(x^2\psi'(x)) + (\lambda-1)x \psi(x) = 0,
+\end{equation}
+which together with \eqref{eq:conpsi2} forms a regular Sturm--Liouville problem with
+non-negative eigenvalues $\lambda_k-1$ and corresponding orthogonal eigenfunctions
+$\psi_k(x)$ with respect to the weight function $\rho(x) = x$.
+
+For $\lambda_0=1$, the eigenfunction is $\psi_0(x)=1$. For $\lambda_k>0$,
+equation \eqref{eq:conpsils} is of the form \eqref{eq:besg}, with $\alpha=2\sqrt{\lambda_k-1}$.
+Writing $\omega_k$ for $\sqrt{\lambda_k-1}$, we have
+\begin{equation}
+    \psi_k(x) = \frac{1}{\sqrt{x}} \left(c_1^{(k)} J_1(2\omega_k \sqrt{x})) + c_2^{(k)} Y_1(2\omega_k \sqrt{x})\right),
+\end{equation}
+with derivative
+\begin{equation}
+    \psi_k'(x) = - \frac{\omega_k}{x} \left(c_1^{(k)} J_2(2\omega_k\sqrt{x})) + c_2^{(k)} Y_2(2\omega_k\sqrt{x})\right).
+\end{equation}
+
+The boundary conditions $\psi_k'(a) = \psi_k'(b) = 0$ imply
+\begin{gather}
+    f_2\left(\sqrt{\tfrac{b}{a}}, 2\omega_k\sqrt{a}\right) = 0,
+    \intertext{where $f_\nu$ is the cross-product Bessel function}
+    f_\nu(q, x) = J_\nu(qx)Y_\nu(x)-J_\nu(x)Y_\nu(qx).
+\end{gather}
+The eigenvalues for $k>0$ are then $\lambda_k=1+\omega_k^2$ with
+\begin{equation}
+    \label{eq:conomegak}
+    \omega_k =  \frac{\chi_{2, k}\left(\sqrt{\scriptstyle \frac{b}{a}}\right)}{2\sqrt{a}}.
+\end{equation}
+where $\chi_{\nu, k}(q)$ is the $k$th positive root of $f_\nu(q, x)=0$. These solutions
+can be computed efficiently via a Newton--Raphson scheme \autocite{sorolla2013}.
+$c_1^{(k)}$ and $c_2^{(k)}$ are then determined up to scale; for example, we can
+set
+\begin{equation}
+    c_1^{(k)} = J_2(2\omega_k\sqrt{x})^{-1},\qquad
+    c_2^{(k)} = -Y_2(2\omega_k\sqrt{x})^{-1},
+\end{equation}
+giving
+\begin{equation}
+    \label{eq:conpsik}
+    \psi_k(x) =\frac{1}{\sqrt{x}}\left(
+    \frac{J_1(2\omega_k\sqrt{x})}{J_2(2\omega_k\sqrt{a})}
+    -
+    \frac{Y_1(2\omega_k\sqrt{x})}{Y_2(2\omega_k\sqrt{a})}
+    \right).
+\end{equation}
+
+The complete solution then will be
+\begin{equation}
+    \label{eq:conuxt}
+    u(x,t) = \hat u(x) - \sum_{k=0}^{\infty} a_k e^{-\lambda t} \psi_k(x),
+\end{equation}
+where the $a_k$ will be determined by the initial conditions.
+
+As the $\psi_k$ are
+eigenfunctions of the Sturm-Liouville problem (\ref{eq:conpsils}, \ref{eq:conpsi2}), they
+are orthogonal with respect to the inner product with weight function $\rho(x)=x$,
+\begin{equation}
+    \langle f, g \rangle = \Int{a}{b} xf(x)g(x) \D x.
+\end{equation}
+Taking the inner product of $\langle u(x,0), \psi_k(x)\rangle$ then gives for each $k$,
+\begin{equation}
+    a_k = \frac{\Int{a}{b} x \hat u(x)\psi_k(x) \D x}{\Int{a}{b} x \psi_k(x)^2 \D x}.
+\end{equation}
+
+For $k=0$,
+\begin{equation}
+    \begin{aligned}
+        a_0 &= \frac{1}{b-a} \cdot \Int{a}{b} x \hat u(x) \D x\\
+        &= \frac{1}{b-a} \cdot \left. x^2\hat u'(x)\right|_a^b\\
+        &= \frac{b^2}{b-a}.
+    \end{aligned}
+\end{equation}
+(See Appendix \ref{ap:conid} for the derivation of this and the following integrals.)
+
+For $k>0$,
+\begin{equation}
+    \begin{aligned}
+        \Int{a}{b} x \hat u(x)\psi_k(x) \D x
+        &= \frac{1}{1+\omega_k^2} \cdot \left. (x^2\hat u'(x)\psi_k(x)-x^2\hat u(x)\psi'_k(x))\right|_a^b\\
+        &= \lambda_k^{-1} b^2\psi_k(b),
+    \end{aligned}
+\end{equation}
+and
+\begin{equation}
+    \begin{aligned}
+        \Int{a}{b} x \psi_k(x)^2 \D x
+        &= \left.
+        x^2\psi_k(x)^2 + \frac{x^2}{\omega_k^2}\psi_k(x)\psi'_k(x)
+        + \frac{x^3}{\omega_k^2}\psi'_k(x)^2
+        \right|_a^b\\
+        &= b^2\psi_k(b)^2 - a^2\psi_k(a)^2,
+    \end{aligned}
+\end{equation}
+giving
+\begin{equation}
+    a_k = \frac{1}{\lambda_k}\cdot \frac{b^2\psi_k(b)}{b^2\psi_k(b)^2-a^2\psi_k(a)^2}.
+\end{equation}
+
+This gives the time-dependent solution
+\begin{equation}
+    u(x,t) = \hat u(x) - \frac{b^2}{b-a}e^{-t} -
+    \sum_{k=1}^{\infty} \lambda_k^{-1} e^{-\lambda_k t}
+    \frac{b^2\psi_k(b)\psi_k(x)}{b^2\psi_k(b)^2-a^2\psi_k(a)^2}
+\end{equation}
+where $\hat u$ is given by \eqref{eq:conuhat}, $\psi_k$ is given by \eqref{eq:conpsik},
+and $\lambda_k = 1 + \omega_k^2$ is given by \eqref{eq:conomegak}.
+
+
+\section{Approximating the gradient}
+
+In a finite volume discretization there
+will be a computation of the approximation to the gradient
+$v'(x)$ as a linear combination of two (or potentially more)
+values of the discrete approximation of the voltage for points
+near $x$.
+
+If the coefficients in this approximation are constant in
+time, they cannot account for source terms or values of the
+voltage outside the points in question. The most faithful
+such approximation should then reproduce the exact gradient $v'(x)$
+in the source-free steady state form of the cable equation,
+\begin{equation}
+    (\sigma v')' = 0.
+\end{equation}
+
+The voltage $v$ is then determined by its values at
+any two distinct points $a$ and $b$,
+\begin{equation}
+    \frac{v(x) - v(a)}{v(b) - v(a)} =
+    \frac{\displaystyle \Int{a}{x} \sigma(z)^{-1} \D z}
+         {\displaystyle \Int{a}{b} \sigma(z)^{-1} \D z},
+\end{equation}
+and correspondingly, the gradient is
+\begin{equation}
+    \sigma(x) v'(x) =
+    \frac{v(b) - v(a)}
+         {\displaystyle\Int{a}{b} \sigma(z)^{-1} \D z}.
+\end{equation}
+If $x$ is a point of discontinuity in $\sigma$,
+the flux $\sigma(x)v'(x)$ is nonetheless well defined.
+
+For a finite volume approximation, however, the available
+estimates may correspond instead to surface area-weighted means
+over control volumes. Consider two adjacent control volumes
+on $[a, m]$ and $[m, b]$, with mean voltages $\bar v_a$ and
+$\bar v_b$ given in terms of a weight function $w(x)$:
+\begin{align}
+    \bar v_a &= w_a^{-1} \Int{a}{m} w(x)v(x) \D x,\\
+    \bar v_b &= w_b^{-1} \Int{m}{b} w(x)v(x) \D x,
+\end{align}
+with normalizing constants $w_a$ and $w_b$. With $\sigma v' = \kappa$
+constant,
+\begin{equation}
+    \begin{aligned}
+        \bar v_a
+        &= w_a^{-1} \Int{a}{m} w(x) \left( v(m) + \kappa \Int{m}{x} \sigma(y)^{-1} \D y \right) \D x \\
+        &= v(m) + \kappa \cdot w_a^{-1} \Int{a}{m} w(x) \Int{m}{x} \sigma(y)^{-1} \D y \D x,
+    \end{aligned}
+\end{equation}
+and similarly, 
+\begin{equation}
+    \begin{aligned}
+        \bar v_b
+        &= v(m) + \kappa \cdot w_b^{-1} \Int{m}{b} w(x) \Int{m}{x} \sigma(y)^{-1} \D y \D x,
+    \end{aligned}
+\end{equation}
+Taking the difference then gives a solution for $\kappa$,
+\begin{equation}
+    \kappa = (\bar v_b - \bar v_a) \cdot \left[
+        w_b^{-1}\!\Int{m}{b} w(x)\! \Int{m}{x} \sigma(y)^{-1} \D y \D x\,-\,
+        w_a^{-1}\!\Int{a}{m} w(x)\! \Int{m}{x} \sigma(y)^{-1} \D y \D x
+        \right]^{\mathrlap{-1}}.
+\end{equation}
+
+\newpage
+\appendix
+\section{The Rallpack 1 model}
+\label{ap:rallpack}
+
+The Rallpack suite \autocite{bhalla1992} is a set of three models for the
+validation and benchmarking of cable-based neuron simulators. The first of
+these comprises a cylindrical passive cable of length $L$ and diameter $d$,
+with a constant current applied at $x=0$. The membrane resistive current
+density is given by $J = R_M^{-1}(v-E)$, where $E$ is a fixed reversal potential.
+The values of the electrial and geometric parameters are given in
+Table~\ref{tbl:rallpack1}.
+
+The membrane potential can then be expressed as
+\begin{equation}
+    v(x, t) = E - \lambda I \cdot u(\frac{b-x}{\lambda}, \frac{t}{\tau}),
+\end{equation}
+where
+\begin{equation}
+    \lambda = \sqrt{\frac{R_M d}{4 R_L}},\quad
+    \tau = R_M C_M,\quad
+    b = L/\lambda,
+\end{equation}
+and $u(x,t)$ is given by \eqref{eq:cylseries}.
+
+\begin{table}[htb]
+    \centering
+    \begin{tabular}{lll}
+        \toprule
+	Parameter & {} & {Value} \\
+        \midrule
+	Cable diameter                    & $d$    & \SI{1.0}{\um} \\
+	Cable length                      & $L$    & \SI{1.0}{\mm} \\
+	Bulk (axial) resistivity          & $R_L$  & \SI{1.0}{\ohm\m} \\
+	Membrane resistivity              & $R_M$  & \SI{4.0}{\ohm\m\squared} \\
+	Membrane specific capacitance     & $C_M$  & \SI{0.01}{\F\per\m\squared} \\
+	Membrane reversal potential       & $E_M$  & \SI{-65.0}{\mV} \\
+	Injected current                  & $I$    & \SI{0.1}{\nA} \\
+        \bottomrule
+    \end{tabular}
+    \caption{Rallpack 1 parameters}
+    \label{tbl:rallpack1}
+\end{table}
+
+
+\section{Series computation for the voltage on the uniform cylinder}
+\label{ap:cylcomp}
+
+How many terms in the series \eqref{eq:cylseries} need to be computed for a given tolerance?
+We obtain bounds on the remainders in the series as follows.
+
+Let $u_n(x,t)$ be the partial sum up to $k=n-1$ in \eqref{eq:cylseries}, and $r_n(x,t) = u(x,t) - u_n(x,t)$
+be the remainder,
+\begin{equation}
+    r_n(x,t) = -\frac{2}{b}\sum_{k=n}^{\infty} (-1)^k\lambda_k^{-1}e^{-\lambda_k t}\cos\frac{k\pi x}{b}.
+\end{equation}
+Then $|r_n(x,t)|\leq R_n(t)$ for all $x$, with
+\begin{equation}
+    R_n(t) = \frac{2}{b}\sum_{k=n}^{\infty} \lambda_k^{-1}e^{-\lambda_k t}\\
+    = \frac{2t}{b}\sum_{k=n}^{\infty} l(k)^{-1}e^{-l(k)}\\
+\end{equation}
+where
+\[
+    l(z) = \left(1+\left(\frac{z\pi}{b}\right)^2\right)\cdot t.
+\]
+The summand is monotonically decreasing in $k$, so
+\begin{equation}
+    \label{eq:rnigamma}
+    R_n(t)
+    \leq \frac{2t}{b}\Int{n}{\infty} l(k)^{-1}e^{-l(k)}\D k
+    = \frac{\sqrt{t}}{b} \Int{\lambda_n t}{\infty} e^{-z}z^{-\frac{3}{2}} \D z.
+\end{equation}
+
+For $a>0$, $s>0$, integrating by parts gives
+$\displaystyle \Int{a}{\infty} e^{-z}z^{-s}\D z < e^{-a}a^{-s}$,
+and so
+\begin{equation}
+    R_n(t) < \frac{1}{b t}\lambda_n^{-3/2}e^{-\lambda_n t}.
+\end{equation}
+
+\section{Identities and integrals for tapered cable eigenfunctions}
+\label{ap:conid}
+
+\newcommand{\Fp}{\smash{F^{\mathrlap{\prime}}}\mkern-2.0mu}
+\newcommand{\Gp}{\smash{G^{\mathrlap{\prime}}}\mkern-2.0mu}
+
+Let
+\begin{equation}
+    \begin{aligned}
+        F_k(z) &= z^{-\frac{k}{2}}\cZ_k(2z^{1/2}),&
+        G_k(z) &= z^{-\frac{k}{2}}\cC_k(2\omega z^{1/2}),
+    \end{aligned}
+\end{equation}
+where $C_k$ and $Z_k$ are solutions to the Bessel and modified Bessel equations
+respectively,
+\begin{equation}
+    \begin{aligned}
+        \cZ_k(z) &= c_1 I_k(z) + c_2 e^{i\pi k} K_k(z),&
+        \cC_k(z) &= c_3 J_k(z) + c_4 Y_k(z)
+    \end{aligned}
+\end{equation}
+for some constants $c_1$, $c_2$, $c_3$, $c_4$.
+
+$F_k$ and $G_k$ are solutions of the differential equations
+\begin{gather}
+    z f''(z) + (k+1) f'(z) - z = 0\\
+    \intertext{and}
+    z g''(z) + (k+1) g'(z) + \omega^2 z = 0
+\end{gather}
+respectively, following the identity (2) of \autocite[p.~512]{lommel1879}.
+
+Applying the derivative and recurrence relations for Bessel functions
+\autocite[\S 10.6, \S10.29]{nistdlmf} then gives the following identities,
+\begin{align}
+    \label{eq:fgrel}
+    \Fp_k(z) &= F_{k+1}(z), & z \Fp_k(z) &= -k F_k(z) + F_{k-1}(z),\\
+    \Gp_k(z) &= -\omega G_{k+1}(z), & z \Gp_k(z) &= -k G_k(z) + \omega G_{k-1}(z).
+\end{align}
+When $k$ is an integer, these imply
+\begin{equation}
+    \begin{aligned}
+        F_{-k}(z) &= z^k F_k(z),\\
+        G_{-k}(z) &= (-z)^k G_k(z),
+    \end{aligned}
+    \qquad \text{for $k\in\mathbb{Z}$}.
+\end{equation}
+
+In the following, derivatives will be respect to $z$, and the argument $z$ to
+$F_k$ and $G_k$ will be omitted where it is unambiguous.
+
+\subsection*{The integral $\Int{}{} z F_1 G_1 \D z$}
+
+The relations \eqref{eq:fgrel} give the identities
+\begin{equation}
+    \begin{aligned}
+        \label{eq:fgrel2}
+        (z^kF_k)' &= z^{k-1}F_{k-1}\\
+        (z^kG_k)' &= \omega z^{k-1}G_{k-1}.
+    \end{aligned}
+\end{equation}
+Consequently
+\begin{equation}
+    \begin{aligned}
+        (z^2 F_2 G_1)' &= z F_1 G_1 + z^2 F_2 \Gp_1 = z F_1 G_1 - \omega z^2 F_2 G_2 \\
+        (z^2 F_1 G_2)' &= z F_1 G_1 + z^2 \Fp_1 G_2 = z F_1 G_1 + z^2 F_2 G_2.
+    \end{aligned}
+\end{equation}
+Cancelling the $F_2G_2$ term gives
+\begin{equation}
+    (z^2 F_2 G_1 + \omega z^2 F_1 G_2)' =
+    (1+\omega^2) z F_1 G_1,
+\end{equation}
+and so
+\begin{equation}
+    (1+\omega^2) \Int{}{} z F_1 G_1 dz = z^2 F_2 G_1 + \omega z^2 F_1 G_2 = z^2 \Fp_1 G_1 - z^2 F_1 \Gp_1.
+\end{equation}
+
+\subsection*{The integral $\Int{}{} z F_1 \D z$}
+
+As observed in \eqref{eq:fgrel2}, $(z^2 F_2)' = z F_1$, and so
+\begin{equation}
+    \Int{}{} z F_1 dz = z^2 F_2 = z^2 \Fp_1.
+\end{equation}
+
+\subsection*{The integral $\Int{}{} z G_1^2 \D z$}
+
+Following the technique of \autocite[p.533]{lommel1879},
+\begin{equation}
+    \begin{aligned}
+        (z^a G_k G_{m+1})'
+        &= a z^{a-1} G_k G_{m+1} + z^a(\Gp_k G_{m+1} + G_k \Gp_{m+1})\\
+        &= a z^{a-1} G_k G_{m+1} - \omega z^a G_{k+1}G_{m+1} + z^{a-1} G_k(\omega G_m - (m+1) G_{m+1})\\
+        &= (a-m-1) z^{a-1} G_k G_{m+1} + \omega z^{a-1} (zG_{k+1}G_{m+1} - G_k G_m).
+    \end{aligned}
+\end{equation}
+The last term is symmetric in $k$ and $m$, so
+\begin{equation}
+    (z^a G_k G_{m+1} - z^a G_m G_{k+1})' =
+    (a-m-1) z^{a-1} G_k G_{m+1} - (a-k-1) z^{a-1} G_m G_{k+1}.
+\end{equation}
+Letting $a=k+1$ then gives,
+\begin{equation}
+    (z^{k+1} G_k G_{m+1} - z^{k+1} G_{k+1} G_m)' = (k-m) z^k G_k G_{m+1},
+\end{equation}
+and in particular,
+\begin{equation}
+    \label{eq:intzg1sq}
+    \Int{}{} z G_1^2 \D z = z^2 G_1^2 - z^2 G_0 G_2.
+\end{equation}
+
+The result \eqref{eq:intzg1sq} can be expressed in terms of just $G_1$ and $\Gp_1$ by
+application of the recurrence relationships:
+\begin{equation}
+    \label{eq:intzg1sqbis}
+    \Int{}{} z G_1^2 \D z = z^2 G_1^2 + \frac{z^2}{\omega^2} G_1 \Gp_1 + \frac{z^3}{\omega^2} G^{\prime 2}_1.
+\end{equation}
+
+
+\printbibliography
+\end{document}

doc/math/passive_cable/cable_computation.tex

-\documentclass[parskip=half]{scrartcl}
-
-\usepackage{amsmath}
-\usepackage{amssymb}
-\usepackage[citestyle=authoryear-icomp,bibstyle=authoryear]{biblatex}
-\usepackage{booktabs}
-\usepackage[style=british]{csquotes}
-\usepackage{fontspec}
-\usepackage{siunitx}
-
-\setmainfont[Numbers=Lowercase]{TeX Gyre Pagella}
-\newfontfamily\dispfamily{TeX Gyre Pagella}
-\addtokomafont{disposition}{\dispfamily}
-
-\MakeAutoQuote{«}{»}
-\DeclareMathOperator{\Res}{Res}
-
-\addbibresource{cable.bib}
-
-\title{Analytic solution to passive cable model}
-\author{Sam Yates}
-\date{October 20, 2016}
-
-\begin{document}
-
-\maketitle
-
-\section{Background}
-
-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}.}
-
-The first two of the tests are based on passive cable models which admit analytic
-solutions, and the Rallpack suite includes code that can generate the reference curves
-for these models. The reference code, however, requires hand-tuning of the iterations;
-the authors state in the internal documentation that
-«the use of 10 terms gives at least six figure accuracy», but this is for the
-fixed time increments used in their dataset and is not validated within the code
-or by the supplied material.
-
-For flexibility of testing and cross-checking of the Rallpack reference data,
-it would be useful to have at hand a robust calculator of the passive cable
-potential solution.
-
-\section{The Rallpack 1 model}
-
-The model comprises a cylindrical cable segment with the physical and electrical
-properties as described in Table~\ref{tbl:rallpack1}. A current of \SI{0.1}{\nA}
-is injected at the left end of the cable from $t=0$, and the initial potential
-is the reversal potential.
-
-\begin{table}[ht]
-    \centering
-    \begin{tabular}{lSl}
-        \toprule
-        Term & {Value} & Property\\
-        \midrule
-        $d$    & \SI{1.0}{\um}                 & cable diameter \\
-        $L$    & \SI{1.0}{\mm}                 & cable length \\
-        $R_A$  & \SI{1.0}{\ohm\m}              & bulk axial resistivity \\
-        $R_M$  & \SI{4.0}{\ohm\m\squared}      & areal membrane resistivity \\
-        $C_M$  & \SI{0.01}{\F\per\m\squared}   & areal membrane capacitance \\
-        $E_M$  & \SI{-65.0}{\mV}               & membrane reversal potential \\
-        \bottomrule
-    \end{tabular}
-    \caption{Cable properties for the Rallpack 1 model.}
-    \label{tbl:rallpack1}
-\end{table}
-
-The potential on a constant-radius passive cable is governed by the PDE
-\begin{equation}
-    \lambda^2 \frac{\partial^2 v}{\partial x^2} =
-    \tau\frac{\partial v}{\partial t} + v - E,
-\end{equation}
-where $E$ is the membrane reversal potential, and $\lambda$ and $\tau$ are the
-electrotonic length and time constants for the cable. It is convenient
-to consider the electrical properties of the cable per unit-length, in terms
-of the linear axial resistivity $r$, the linear membrane capacitance $r$,
-and the linear membrane capacitance $c$. These determine $\lambda$ and $\tau$ by
-\begin{gather*}
-    \lambda = \frac{1}{r\cdot g},\\
-    \tau = r\cdot c.
-\end{gather*}
-
-With the model boundary conditions,
-\begin{subequations}
-    \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,
-    \end{align}
-\end{subequations}
-where $I$ is the injected current and $L$ is the cable length.
-
-The solution $v(x, t)$ can be expressed in terms of the solution $g(x, t; L)$
-to a normalized version of the cable equation,
-\begin{subequations}
-    \begin{align}
-        \label{eq:normcable}
-        \frac{\partial^2 g}{\partial x^2} & =
-        \frac{\partial g}{\partial t} + g,
-        \\
-        \label{eq:normcableinitial}
-        g(x, 0) &= 0,
-        \\
-        \label{eq:normcableleft}
-        \left.\frac{\partial g}{\partial x}\right\vert_{x=0} & = 1,
-        \\
-        \label{eq:normcableright}
-        \left.\frac{\partial g}{\partial x}\right\vert_{x=L} & = 0
-    \end{align}
-\end{subequations}
-by
-\begin{equation}
-    v(x, t)= E - Ir\lambda \cdot g(\frac{x}{\lambda}, \frac{t}{\tau};  \frac{L}{\lambda}).
-\end{equation}
-
-\section{Solution to the normalized cable equation}
-
-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}
-    \frac{\partial^2 G}{\partial x^2} = sG + G.
-\end{equation}
-The boundary conditions \eqref{eq:normcableleft} and \eqref{eq:normcableright} give
-\begin{align}
-    \label{eq:lapleft}
-    \frac{\partial G}{\partial x}(0, s) & = \frac{1}{s}
-    \\
-    \label{eq:lapright}
-    \frac{\partial G}{\partial x}(L, s) & = 0.
-\end{align}
-
-Solutions to \eqref{eq:lap} must be of the form
-\begin{equation}
-    G(x, s) = A(s) e^{mx} +B(s) e^{-mx}
-\end{equation}
-where $m=\sqrt{1+s}$. From \eqref{eq:lapleft} and \eqref{eq:lapright},
-\begin{align*}
-    mA(s) - mB(s) & = \frac{1}{s}
-    \\
-    mA(s)e^{mL} -mB(s)e^{-mL} & = 0
-\end{align*}
-and thus
-\begin{align*}
-    A(s) & = \frac{1}{sm (1-e^{2mL})}
-    \\
-    B(s) & = \frac{e^{2mL}}{sm (1-e^{2mL})}.
-\end{align*}
-
-Consequently,
-\begin{equation}
-    \begin{aligned}
-        G(x, s) &= \frac{1}{ms}\cdot\frac{e^{mx}+e^{2mL-mx}}{1-e^{2mL}}\\
-                &= - \frac{1}{ms}\cdot\frac{\cosh m(L-x)}{\sinh mL}.
-    \end{aligned}
-\end{equation}
-
-\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 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}\\
-                   & = - \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=\sqrt{v/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^{-tv} 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}