On the Telegrapher’s Equation with Three Space Variables in Non-Rectangular Coordinates

Tristan Guillaume

Thema Research Unit, CY Cergy Paris University, Cergy, France.

DOI: 10.4236/jamp.2020.85070                

Abstract

This article provides a closed form solution to the telegrapher’s equation with three space variables defined on a subset of a sphere within two radii, two azimuthal angles and one polar angle. The Dirichlet problem for general boundary conditions is solved in detail, on the basis of which Neumann and Robin conditions are easily handled. The solution to the simpler problem in cylindrical coordinates is also provided. Ways to efficiently implement the formulae are explained. Minor adjustments result in solutions to the wave equation and to the heat equation on the same domain as well, since the latter are particular cases of the more general telegrapher’s equation.

Keywords

Telegrapher’s Equation, Spherical Coordinates, Fourier-Bessel-Legendre Series Expansion, Associated Legendre Functions, Hyperbolic PDE, Initial Boundary Value Problem

1. Introduction

Hyperbolic partial differential equations (PDEs) are among the three main classes of PDEs in mathematical physics. The most extensively studied hyperbolic PDE is the wave equation. Another major hyperbolic PDE is the telegrapher’s equation. In rectangular coordinates, its canonical form is the following:

$\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}+{\kappa }_1\frac{\partial u\left(x,t\right)}{\partial t}={\kappa }_2{\Delta }_{n}u\left(x,t\right)-{\kappa }_{3}u\left(x,t\right)+\Lambda \left(x,t\right)$ (1)

where $x$ is a vector of n space variables, ${\Delta }_n$ is the n-dimensional Laplacian operator, $\Lambda$ is a “sufficiently” smooth function and ${\kappa }_1,{\kappa }_2,{\kappa }_3$ are three positive constants.

A good introduction to the fundamental properties and main applications of this PDE can be found in Zauderer [1] or Ockendon et al. [2].

Among all the major second order linear PDEs of mathematical physics, the telegrapher’s equation is noticeable for its generality. Indeed, all the most important non-stationary linear equations can be retrieved from (1). In particular, the wave equation is a special case of the telegrapher’s equation obtained by taking ${\kappa }_1={\kappa }_3=\Lambda \left(x,t\right)=0$. The Klein-Gordon equation for a function $w\left(x,t\right)$ is obtained through the following simple relation:

$w\left(x,t\right)=\exp \left(-\frac{{\kappa }_{1}t}{2}\right)u\left(x,t\right)$ (2)

An obvious parameterization of the second-order time derivative in (1) would also enable to recover the parabolic heat equation. For more general background on PDEs related to Equation (1), one can refer to [3].

A number of initial boundary value problems (IBVPs) involving PDE (1) have already been analytically solved, especially in rectangular coordinates. Certainly the most comprehensive compilation of known analytical formulae is the handbook by Polyanin and Nazaikinskii [4], which collects and updates all the results from mathematical physics disseminated in previous handbooks such as [5] and [6]. It appears that not all IBVPs associated with the three-dimensional telegrapher’s equation have had their exact solutions derived in the literature. This can often be explained by the fact that there already exist solution techniques applied to other non-stationary equations that can readily be extended to the multidimensional telegrapher’s equation. However, it ought to be noticed that, in some cases, such an extension entails non-trivial computational issues. This is especially true when it comes to non-rectangular coordinates. On spherical domains, the most general known formulae can cope with boundary conditions with respect to the radial variable, but they cannot handle functional dependency on the boundaries of both the radial and the angular variables. In other words, if we denote the radius by r, the azimuthal (or longitudinal) angle by $\varphi$ and the polar (or latitudinal) angle by $\theta$, then the only exact solutions to the telegrapher’s equation currently available are those that apply to the following domain [4]:

$\left\{0<r_1\le r\le r_2<\infty ,0<\varphi \le 2\pi ,0<\theta \le \pi ,t\ge 0\right\}$ (3)

Such a restriction rules out problems that require functional dependency on the boundaries of the variables $\varphi$ or $\theta$ as well as on the boundary of r. Of course, one can always resort to the usual numerical techniques of approximation, but, in this case, it is not easy to design efficient and accurate finite difference schemes, due to the dimension of the problem and to the complexity of the telegrapher’s equation in a spherical domain. The main contribution of this paper is to provide a closed form solution for a class of IBVPs involving the linear telegrapher’s equation with three space variables and general boundary conditions on the following domain:

$D=\left\{0<r_1\le r\le r_2<\infty ,0<\theta \le {\theta }_1\le \pi ,0<{\varphi }_1\le \varphi \le {\varphi }_2\le 2\pi ,t\ge 0\right\}$ (4)

The domain in (4) allows for partial ranges in both angular variables. As such, it includes subsets of the sphere with many more shapes than the domain in (3), so that more flexible models can be considered, with a greater number of physical specifications.

The linear and separable nature of the PDE under consideration enables the application of classical methods of separation of variables and eigenfunction expansions. One of the main advantages of this approach is that it easily allows to handle non-homogeneity in the equation or in the boundary conditions, once the homogeneous solution has been obtained. However, the solution process involves a non-standard generalized Fourier expansion of the first order initial condition, as well as computational difficulties associated with the integration of associated Legendre functions with non-integral parameters, which may explain why there is currently no available formula in the literature.

It must be emphasized that the telegrapher’s equation is the focus of this article due to its high degree of generality, not because a specific physical application is studied. Thus, it is expected that the analytical results in this article can be fruitfully applied to parabolic problems—after adjustment of the initial conditions, as well as to hyperbolic problems that do not model the same physical phenomena as those initially described by the telegrapher’s equation.

This paper is organized as follows. In Section 2, the main result is stated under the title “Result 1”, i.e. the solution of the three-dimensional telegrapher’s equation on a spherical domain with boundary conditions on all the space variables. The rest of Section 2 is dedicated to the proof of Result 1, which derives from the combination of two lemmas, along with comments on the numerical implementation of the formula. Section 3 presents two corollaries that handle non-homogeneity of the equation and of the boundary conditions, and then states, under the title “Result 2”, the solution of the three-dimensional telegrapher’s equation on a cylindrical domain with boundary conditions on all the space variables.

2. Main Result

The main result of this article is now stated as “Result 1”.

Result 1

Solution to the Dirichlet problem for the telegrapher’s equation with three space variables on a subset of a sphere located within two radii, two azimuthal angles and one polar angle.

Let $u\left(r,\varphi ,\theta ,t\right)$ be the solution to the following partial differential equation:

$\frac{{\partial }^{2}u}{\partial t^2}+{\kappa }_1\frac{\partial u}{\partial t}={\kappa }_2\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{2}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\left(\frac{{\partial }^{2}u}{\partial {\theta }^2}+\cot \left(\theta \right)\frac{\partial u}{\partial \theta }+{\csc }^2\left(\theta \right)\frac{{\partial }^{2}u}{\partial {\varphi }^2}\right)\right)-{\kappa }_{3}u$ (5)

on the domain:

$D=\left\{0<r_1\le r\le r_2<\infty ,0<\theta \le {\theta }_1\le \pi ,0<{\varphi }_1\le \varphi \le {\varphi }_2\le 2\pi ,t\ge 0\right\}$ (6)

with the following boundary conditions:

$u\left(r,\theta ,\varphi ,0\right)=0$ (7)

${\frac{\partial u\left(r,\theta ,\varphi ,t\right)}{\partial t}|}_{t=0}=f\left(r,\theta ,\varphi \right)$ (8)

$u\left(r_1,\theta ,\varphi ,t\right)=u\left(r_2,\theta ,\varphi ,t\right)=0$ (9)

$u\left(r,{\theta }_1,\varphi ,t\right)=0$ (10)

$u\left(r,\theta ,{\varphi }_1,t\right)=u\left(r,\theta ,{\varphi }_2,t\right)=0$ (11)

where ${\kappa }_1,{\kappa }_2,{\kappa }_3$ are three positive constants and f is a continuous function from D to $ℝ$.

Let us denote:

· $P_{\nu }^{\mu }\left(x\right)$ as the associated Legendre function of the first kind with real variable x, order $\mu$ and degree $\nu$

· $J_{\nu }\left(x\right)$ and $Y_{\nu }\left(x\right)$ as the Bessel functions of the first kind and second kind, respectively, with real order $\nu$

· $j_{\nu }\left(x\right)$ and $y_{\nu }\left(x\right)$ as the spherical Bessel functions of the first kind and second kind, respectively, with real order $\nu$

· $\Gamma \left(x\right)$ as the Gamma function with real variable x

Then, the function $u\left(r,\varphi ,\theta ,t\right)$ is given by:

$u\left(r,\theta ,\varphi ,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{b}_{mnp}T\left(t\right){\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)}}}$ (12)

where the functions $T\left(t\right)$, ${\Phi }_m\left(\varphi \right)$, ${\Theta }_n\left(\theta \right)$ and $R_p\left(r\right)$ are defined as follows:

$T\left(t\right)=\exp \left(-\frac{{\kappa }_1}{2}t\right)\sin \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}\right)$ (13)

${\Phi }_m\left(\varphi \right)=\sin \left({\beta }_m\left(\varphi -{\varphi }_1\right)\right)$ (14)

${\Theta }_n\left(\theta \right)=P_{{\nu }_n}^{{\beta }_m}\left(\cos \theta \right)$ (15)

$\begin{array}{c}{R}_p\left(r\right)=\sqrt{\frac{\pi }{2{\lambda }_p{r}_1}}{J}_{{\nu }_n+1/2}\left({\lambda }_p{r}_1\right)\sqrt{\frac{\pi }{2{\lambda }_{p}r}}{Y}_{{\nu }_n+1/2}\left({\lambda }_{p}r\right)\\ -\sqrt{\frac{\pi }{2{\lambda }_p{r}_1}}{Y}_{{\nu }_n+1/2}\left({\lambda }_p{r}_1\right)\sqrt{\frac{\pi }{2{\lambda }_{p}r}}{J}_{{\nu }_n+1/2}\left({\lambda }_{p}r\right)\end{array}$ (16)

The numbers ${\beta }_m,{\nu }_n$ and ${\lambda }_p$ are defined as follows:

· $\forall m\in \mathbb{N}$, ${\beta }_m$ is the m-th root of the equation

$\sin \left(\beta {\varphi }_2\right)\cos \left(\beta {\varphi }_1\right)-\sin \left(\beta {\varphi }_1\right)\cos \left(\beta {\varphi }_2\right)=0$ (17)

· $\forall n\in \mathbb{N}$, ${\nu }_n$ is the n-th root of the equation

$P_{\nu }^{{\beta }_m}\left(\cos {\theta }_1\right)=0$ (18)

· $\forall p\in \mathbb{N}$, ${\lambda }_p$ is the p-th root of the equation

$y_{{\nu }_n}\left(\lambda r_2\right)j_{{\nu }_n}\left(\lambda r_1\right)-y_{{\nu }_n}\left(\lambda r_1\right)j_{{\nu }_n}\left(\lambda r_2\right)=0$ (19)

The coefficient $b_{mnp}$ is given by:

$\begin{array}{l}{b}_{mnp}=\frac{4{\left({\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)\right)}^{-1/2}{\beta }_m{r}_1{\lambda }_p^4{J}_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)}{\left({\beta }_m\left({\varphi }_2-{\varphi }_1\right)-\sin \left(2{\beta }_m\left({\varphi }_2-{\varphi }_1\right)\right)\right)\left(J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_1\right)-J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)\right){\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2}\\ \times {\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}f\left(\phi ,\vartheta ,\rho \right)\sin \left({\lambda }_p\left(\phi -{\varphi }_1\right)\right)P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)R_p\left(\rho \right){\rho }^2\sin \left(\vartheta \right)\text{d}\phi \text{d}\vartheta \text{d}\rho }}}\end{array}$ (20)

where:

${\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2={\displaystyle \underset{0}{\overset{\pi }{\int }}{\left(\cot \left(\frac{\vartheta }{2}\right)\right)}^{2{\beta }_m}{\left({\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\left(-1\right)}^k\frac{\Gamma \left({\nu }_n+k+1\right)}{\Gamma \left({\nu }_n-k+1\right)}\frac{\sin {\left(\frac{\vartheta }{2}\right)}^{2k}}{\Gamma \left(k+1\right)\Gamma \left(k-{\beta }_m+1\right)}}\right)}^2\sin \left(\vartheta \right)\text{d}\vartheta }$ (21)

which completes the statement of Result 1.

Proof of Result 1:

A solution to the boundary value problem under consideration can be found by the techniques of separation of variables and eigenfunction expansions. The proof consists of two lemmas.

Lemma 1

A solution which verifies both PDE (5) and the boundary conditions associated with the space variables (9)-(11) is given by:

$\begin{array}{c}u\left(r,\theta ,\varphi ,t\right)=\exp \left(-\frac{{\kappa }_1}{2}t\right){\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}(a_{mnp}\cos \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}\right)}}}\\ +b_{mnp}\sin \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}\right)){\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)\end{array}$ (22)

where $a_{mnp}$ and $b_{mnp}$ are coefficients to be determined, and the functions $T\left(t\right)$, ${\Phi }_m\left(\varphi \right)$, ${\Theta }_n\left(\theta \right)$ and $R_p\left(r\right)$ are given in Result 1.

Proof of Lemma 1

The substitution $u\left(r,\varphi ,\theta ,t\right)=R\left(r\right)\Phi \left(\varphi \right)\Theta \left(\theta \right)T\left(t\right)$ carried out into PDE (5) yields:

$\frac{1}{{\kappa }_2}\frac{T^{″}}{T}+\frac{{\kappa }_1}{{\kappa }_2}\frac{T^{\prime }}{T}=\frac{R^{″}}{R}+\frac{2}{r}\frac{R^{\prime }}{R}+\frac{1}{r^2}\left(\frac{{\Theta }^{″}}{\Theta }+\left(\cot \theta \right)\frac{{\Theta }^{\prime }}{\Theta }+\left({\csc }^2\theta \right)\frac{{\Phi }^{″}}{\Phi }\right)-\frac{{\kappa }_3}{{\kappa }_2}$ (23)

Equating both sides of Equation (23) to $-{\lambda }^2$, where ${\lambda }^2$ is a first separation constant, leads to the following two equations:

$\frac{1}{{\kappa }_2}\frac{T^{″}}{T}+\frac{{\kappa }_1}{{\kappa }_2}\frac{T^{\prime }}{T}+{\lambda }^2+\frac{{\kappa }_3}{{\kappa }_2}=0$ (24)

$r^2\frac{R^{″}}{R}+2r\frac{R^{\prime }}{R}+r^2{\lambda }^2=-\left(\frac{{\Theta }^{″}}{\Theta }+\left(\cot \theta \right)\frac{{\Theta }^{\prime }}{\Theta }+\left({\csc }^2\theta \right)\frac{{\Phi }^{″}}{\Phi }\right)$ (25)

The complex roots of the characteristic equation associated with the differential Equation (24) easily yield:

$\begin{array}{c}T\left(t\right)=a_1\exp \left(-\frac{{\kappa }_1}{2}t\right)\cos \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }^2-{\kappa }_3\right)}\right)\\ +a_2\exp \left(-\frac{{\kappa }_1}{2}t\right)\sin \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }^2-{\kappa }_3\right)}\right)\end{array}$ (26)

where $a_1$ and $a_2$ are arbitrary constants.

Equating both sides of Equation (25) with a second separation constant ${\mu }^2$ leads to the following two differential equations:

$r^2{R}^{″}+2r{R}^{\prime }+\left(r^2{\lambda }^2-{\mu }^2\right)R=0$ (27)

$\frac{{\Theta }^{″}}{\Theta }+\left(\cot \theta \right)\frac{{\Theta }^{\prime }}{\Theta }+\left({\csc }^2\theta \right)\frac{{\Phi }^{″}}{\Phi }+{\mu }^2=0$ (28)

Equation (27) is a spherical Bessel equation, while Equation (28) is a spherical harmonics equation.

Rearranging Equation (28) and equating both sides with $-{\beta }^2$, where ${\beta }^2$ is a third separation constant, results in the following two differential equations:

${\Phi }^{″}+{\beta }^2\Phi =0$ (29)

${\Theta }^{″}+\left(\cot \theta \right){\Theta }^{\prime }+\left({\mu }^2-{\beta }^2{\csc }^2\theta \right)\Theta =0$ (30)

The general solution to Equation (29) is given by:

$\Phi \left(\varphi \right)=b_1\cos \left(\beta \varphi \right)+b_2\sin \left(\beta \varphi \right)$ (31)

where $b_1$ and $b_2$ are arbitrary constants.

The boundary condition (11) yields:

$\{\begin{array}{l}{b}_1\cos \left(\beta {\varphi }_1\right)+b_2\sin \left(\beta {\varphi }_1\right)=0\left(\text{32}\right)\\ b_1\cos \left(\beta {\varphi }_2\right)+b_2\sin \left(\beta {\varphi }_2\right)=0(\; 33\; )\end{array}$

Substituting $b_1$ from Equation (32) into Equation (33), we get:

$b_2\left(\sin \left(\beta {\varphi }_2\right)-\frac{\sin \left(\beta {\varphi }_1\right)\cos \left(\beta {\varphi }_2\right)}{\cos \left(\beta {\varphi }_1\right)}\right)=0$ (34)

Since $b_2=0$ entails the trivial solution, $\beta$ must solve:

$\sin \left(\beta {\varphi }_2\right)\cos \left(\beta {\varphi }_1\right)-\sin \left(\beta {\varphi }_1\right)\cos \left(\beta {\varphi }_2\right)=0$ (35)

Classical Sturm-Liouville theory [7] ensures that there is an ordered, countably infinite sequence of real roots ${\beta }_m,m\in \mathbb{N}$ of Equation (35). The evaluation of these roots is an easy task using Newton’s method, as implemented in built-in procedures for equation solving in standard scientific computing packages such as the fsolve command in Maple or the NSolve command in Mathematica.

Eigenfunctions ${\Phi }_m\left(\varphi \right)$ that accommodate the boundary conditions at ${\varphi }_1$ and ${\varphi }_2$ can thus be written as follows:

${\Phi }_m\left(\varphi \right)=\cos \left({\beta }_m{\varphi }_1\right)\sin \left({\beta }_m\varphi \right)-\sin \left({\beta }_m{\varphi }_1\right)\cos \left({\beta }_m\varphi \right)=\sin \left({\beta }_m\left(\varphi -{\varphi }_1\right)\right)$ (36)

Equation (30) is the polar form of the associated Legendre differential equation. It is well-known [3] that a necessary condition for the function $\Theta \left(\theta \right)$ to take bounded values in $\left[0,\pi \right]$ is to have ${\mu }^2=\nu \left(\nu +1\right)$. The general solution of Equation (30) is then given by:

$\Theta \left(\theta \right)=c_1{P}_{\nu }^{{\beta }_m}\left(\cos \theta \right)+c_2{Q}_{\nu }^{{\beta }_m}\left(\cos \theta \right)$ (37)

where $Q_{\nu }^{{\beta }_m}$ is the associated Legendre function of the second kind.

The latter solution has to be eliminated if we want bounded solutions, since $Q_{\nu }^{{\beta }_m}\left(\cos \theta \right)$ goes to infinity as $\theta$ approaches 0 or $\pi$. The boundary condition (10) yields:

$P_{\nu }^{{\beta }_m}\left(\cos {\theta }_1\right)=0$ (38)

Since ${\theta }_1$ is fixed and ${\beta }_m$ has already been determined, Equation (38) has to be solved for $\nu$. A convenient way of doing this is to use the Mehler-Dirichlet representation for the associated Legendre function $P_{\nu }^{\mu }\left(\cos \theta \right)$, which is valid for any $\theta$ between zero and $\pi$ :

$P_{\nu }^{\mu }\left(\cos \theta \right)=\frac{{\left(\sin \theta \right)}^{\mu }}{\sqrt{2\pi }\Gamma \left(\frac{1}{2}-\mu \right)}{\displaystyle \underset{0}{\overset{\theta }{\int }}\frac{\cos \left(\left(\nu +\frac{1}{2}\right)\vartheta \right)}{{\left(\cos \vartheta -\cos \theta \right)}^{\mu +\frac{1}{2}}}\text{d}\vartheta }$ (39)

Some background on the formula (39) can be found in [8]. It is based on the fact that associated Legendre functions can have integral representations by means of contour and definite integrals. Formulae containing contour integrals have the most general character, but, for applications, representations by integrals along some segments of the real axis are more convenient. Of all such representations, the one in (39) is best suited to numerical computations as it is defined along the finite segment $\left[0,\theta \right]$, which only involves a small number of function evaluations. Thus, the zeros of $P_{\nu }^{{\beta }_m}\left(\cos {\theta }_1\right)$ can be determined by the vanishing of the integral in (39) evaluated with $\mu ={\beta }_m$ and $\theta ={\theta }_1$. A simple procedure of numerical evaluation of the integral for a few trial values of $\nu$, accompanied by interpolation of the results, will yield the roots $\left\{{\nu }_n,n\in \mathbb{N}\right\}$ of Equation (38). Since the differences between two successive elements of the sequence ${\left({\nu }_n-{\nu }_{n-1}\right)}_{n\in \mathbb{N}}$ rapidly become small, the distance between the seed of the root-searching algorithm and the exact value of the root at each new iteration is small too, and the algorithm is therefore efficient. Also, although the ${\nu }_n$ ’s can, in principle, take any value, only those values that are greater than −1/2 are actually needed, due to the following identity [8]:

$P_{-\nu -1}^{\mu }\left(z\right)=P_{\nu }^{\mu }\left(z\right)$ (40)

Other interesting properties of the roots of associated Legendre functions with non-integral parameters can be found in [9], where these functions are encountered in a problem of angular confinement in quantum dots.

Each eigenfunction ${\Theta }_n\left(\theta \right)$ can thus be written as:

${\Theta }_n\left(\theta \right)=P_{{\nu }_n}^{{\beta }_m}\left(\cos \theta \right)$ (41)

The general solution to Equation (27) is given by:

$R\left(r\right)=d_1{j}_{{\nu }_n}\left(\lambda r\right)+d_2{y}_{{\nu }_n}\left(\lambda r\right)$ (42)

There are no singularity points associated with the functions in the right hand side of (42), as the variable $r$ is strictly positive. It is then easily seen, applying classical Sturm-Liouville theory [7], that the boundary condition (9) implies that, for each $p\in \mathbb{N}$, there is a real p-th root, denoted as ${\lambda }_p$, of the following equation:

$y_{{\nu }_n}\left(\lambda r_2\right)j_{{\nu }_n}\left(\lambda r_1\right)-y_{{\nu }_n}\left(\lambda r_1\right)j_{{\nu }_n}\left(\lambda r_2\right)=0$ (43)

The zeros of Bessel functions have been extensively studied and this knowledge is embedded in the built-in procedures of standard scientific computing packages such as the BesselJZeros and BesselYZeros commands in Maple.

Eigenfunctions $R_p\left(r\right)$ that accommodate the boundary conditions at $r_1$ and $r_2$ can thus be written as follows:

$R_p\left(r\right)=j_{{\nu }_n}\left({\lambda }_p{r}_1\right)y_{{\nu }_n}\left({\lambda }_{p}r\right)-y_{{\nu }_n}\left({\lambda }_p{r}_1\right)j_{{\nu }_n}\left({\lambda }_{p}r\right)$ (44)

$\begin{array}{c}=\sqrt{\frac{\pi }{2{\lambda }_p{r}_1}}{J}_{{\nu }_n+1/2}\left({\lambda }_p{r}_1\right)\sqrt{\frac{\pi }{2{\lambda }_{p}r}}{Y}_{{\nu }_n+1/2}\left({\lambda }_{p}r\right)\\ -\sqrt{\frac{\pi }{2{\lambda }_p{r}_1}}{Y}_{{\nu }_n+1/2}\left({\lambda }_p{r}_1\right)\sqrt{\frac{\pi }{2{\lambda }_{p}r}}{J}_{{\nu }_n+1/2}\left({\lambda }_{p}r\right)\end{array}$ (45)

which completes the proof of Lemma 1.

The next step consists in finding a way to satisfy the initial conditions (7) and (8).

First, the initial condition (7) immediately yields $a_{mnp}=0$, so that we obtain:

$u\left(r,\theta ,\varphi ,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{b}_{mnp}T\left(t\right){\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)}}}$ (46)

where $T\left(t\right)$ is given by (13), ${\Phi }_m\left(\varphi \right)$ is given by (36), ${\Theta }_n\left(\theta \right)$ is given by (41), $R_p\left(r\right)$ is given by (45).

Next, the substitution of the initial condition (8) into (46) yields:

$f\left(r,\theta ,\varphi \right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{b}_{mnp}\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}{\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)}}}$ (47)

Lemma 2 shows how to expand $f\left(r,\theta ,\varphi \right)$ into an appropriate Fourier-Legendre-Bessel triple series, thus allowing to determine the coefficients $b_{mnp}$.

Lemma 2

Let $f\left(r,\theta ,\varphi \right)$ be a continuous real function defined on D. Then, $f\left(r,\theta ,\varphi \right)$ admits the following series expansion:

$f\left(r,\theta ,\varphi \right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_{mnp}{\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)}}}$ (48)

where:

$d_{mnp}=\frac{4{\beta }_m{r}_1{\lambda }_p^4{J}_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)}{\left({\beta }_m\left({\varphi }_2-{\varphi }_1\right)-\sin \left(2{\beta }_m\left({\varphi }_2-{\varphi }_1\right)\right)\right)\left(J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_1\right)-J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)\right){\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2}$ (49)

${\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}f\left(\phi ,\vartheta ,\rho \right)\sin \left({\lambda }_p\left(\phi -{\varphi }_1\right)\right)P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)R_p\left(\rho \right){\rho }^2\sin \left(\vartheta \right)\text{d}\phi \text{d}\vartheta \text{d}\rho }}}$ (50)

and ${\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2$ is given in (21).

Proof of Lemma 2:

Based on the separated form of the product $\Phi \left(\varphi \right)\Theta \left(\theta \right)R\left(r\right)$, an appropriate approach is to break down the coefficient $d_{mnp}$ into the following product of three coefficients:

$d_{mnp}=d_m{d}_n{d}_p$ (51)

and to solve successively for $d_m,d_n$ and $d_p.$

It is easy to verify that the set of functions $\left\{{\Phi }_m\left(\varphi \right),m\in \mathbb{N}\right\}$ is orthonormal in $\left[{\varphi }_1,{\varphi }_2\right]$ with respect to the following norm:

${\Vert {\Phi }_m\Vert }^2\triangleq {\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\left[\cos \left({\beta }_m{\varphi }_1\right)\sin \left({\beta }_m\phi \right)-\sin \left({\beta }_m{\varphi }_1\right)\cos \left({\beta }_m\phi \right)\right]}^2\text{d}\phi }$ (52)

$=\frac{1}{2}\left({\varphi }_2-{\varphi }_1\right)-\frac{1}{2{\beta }_m}\sin \left(2{\beta }_m\left({\varphi }_2-{\varphi }_1\right)\right)$ (53)

Therefore, considering r and $\theta$ as constants, (48) can be rewritten as:

$f\left(r,\theta ,\varphi \right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{d}_m{\Phi }_m\left(\varphi \right)}$ (54)

where:

$d_m={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_{np}{\Theta }_n\left(\theta \right)R_p\left(r\right)}}$ (55)

The coefficients $d_m$ in (54) are the Fourier coefficients of the expansion of $f\left(r,\theta ,\varphi \right)$ into a series of sine functions ${\Phi }_m\left(\varphi \right)$ for variable $\varphi$ and fixed $\theta$ and r. Classical Fourier analysis yields:

$d_m=\frac{1}{{\Vert {\Phi }_m\Vert }^2}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}f\left(\phi ,\theta ,r\right)\sin \left({\beta }_m\left(\phi -{\varphi }_1\right)\right)\text{d}\phi }$ (55)

where ${\Vert {\Phi }_m\Vert }^2$ is given by (53), and (54) becomes:

$f\left(r,\theta ,\varphi \right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}\left({\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_{np}{\Theta }_n\left(\theta \right)R_p\left(r\right)}}\right)d_m{\Phi }_m\left(\varphi \right)}$ (56)

where $d_m$ is given by (55) and the coefficient $d_{np}$ remains to be determined.

Given $d_m$ and considering r as a constant, we have:

${\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_{np}{\Theta }_n\left(\theta \right)R_p\left(r\right)}}={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{d}_n{\Theta }_n\left(\theta \right)}$ (57)

where:

$d_n={\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_p{R}_p\left(r\right)}$ (58)

For a given m, the set of functions $\left\{{\Theta }_n\left(\theta \right),n\in \mathbb{N}\right\}$ is orthogonal in $\left[0,\pi \right]$ with respect to a weight given by the sine function, i.e.:

${\displaystyle \underset{0}{\overset{\pi }{\int }}{P}_{{\nu }_{n_1}}^{{\beta }_m}\left(\cos \vartheta \right)P_{{\nu }_{n_2}}^{{\beta }_m}\left(\cos \vartheta \right)\sin \left(\vartheta \right)\text{d}\vartheta }=0,n_1\ne n_2$ (59)

The set of functions $\left\{{\Theta }_n\left(\theta \right),n\in \mathbb{N}\right\}$ is orthonormal in $\left[0,\pi \right]$ with respect to the following norm:

${\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2\triangleq {\displaystyle \underset{0}{\overset{\pi }{\int }}{\left(P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)\right)}^2\sin \left(\vartheta \right)\text{d}\vartheta }$ (60)

The coefficients $d_n$ are the Fourier coefficients of the expansion of $f\left(r,\theta ,\varphi \right)$ into a series of functions ${\Theta }_n\left(\theta \right)$ for variable $\theta$ and fixed r and $\varphi$. In practice, one needs to be able to perform the numerical evaluation of the coefficients $d_n$. The function $P_{{\nu }_n}^{{\beta }_m}(.)$ is an associated Legendre function of the first kind with non-integral order ${\beta }_m$ and non-integral degree ${\nu }_n$. Its moment-generating function is undefined and it is therefore less easy to evaluate the norm ${\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2$ than if we were dealing with a standard associated Legendre function with integral parameters. However, any function $P_{{\nu }_n}^{{\beta }_m}\left(z\right)$ is a solution to a second order, linear ordinary differential equation, with three regular singular points and, as such, it can be recast as into a hypergeometric function ${}_{2}F{}_1\left(a,b;c;z\right)$ as follows [10]:

$P_{{\nu }_n}^{{\beta }_m}\left(z\right)=\frac{1}{\Gamma \left(1-{\beta }_m\right)}{\left(\frac{z+1}{z-1}\right)}^{\frac{{\beta }_m}{2}}{}_{2}F{}_1\left(-{\nu }_n,{\nu }_n+1;1-{\beta }_m;\frac{1-z}{2}\right)$ (61)

The function ${}_{2}F{}_1\left(a,b;c;z\right)$, also known as the Gauss hypergeometric function, admits the following series expansion [10]:

${}_{2}F{}_1\left(a,b;c;z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(a\right)}_k{\left(b\right)}_k}{{\left(c\right)}_k}\frac{z^k}{k!}}=\frac{\Gamma \left(c\right)}{\Gamma \left(a\right)\Gamma \left(b\right)}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{\Gamma \left(a+k\right)\Gamma \left(b+k\right)}{\Gamma \left(c+k\right)}\frac{z^k}{k!}}$ (62)

This series representation is well-suited to our problem, since it converges for $\left|z\right|<1$.

Letting $z=\cos \left(\vartheta \right)$ in (61) and performing a little algebra, one can obtain the following formula for the coefficients $d_n$ :

$d_n=\frac{1}{{\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2}{\displaystyle \underset{0}{\overset{\pi }{\int }}f\left(\varphi ,\vartheta ,r\right)P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)\sin \left(\vartheta \right)\text{d}\vartheta }$ (63)

where:

$P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)={\left(\cot \left(\frac{\vartheta }{2}\right)\right)}^{{\beta }_m}{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\left(-1\right)}^k\frac{\Gamma \left({\nu }_n+k+1\right)}{\Gamma \left({\nu }_n-k+1\right)}\frac{\sin {\left(\frac{\vartheta }{2}\right)}^{2k}}{\Gamma \left(k+1\right)\Gamma \left(k-{\beta }_m+1\right)}}$ (64)

and ${\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2$ is given in (21).

A similar trigonometric series is derived by Bremer [11], who also develops a very efficient new algorithm of evaluation. As pointed out in [11], when the order and degree parameters are “small”, i.e. typically close to zero or to the first few natural integers, the coefficients in the series decay rapidly as k increases. The consequence is that only a small number of terms are required to achieve adequate convergence for all practical purposes. Likewise, even when the parameters are of “large” magnitude, i.e. typically greater than ten times some natural integer after the first few integers, the coefficients in the expansion decay rapidly with k provided that $\vartheta$ is sufficiently small. Since $\vartheta \in \left[0,\pi \right]$, the latter condition is guaranteed, so that the series expansion remains efficient in this range of parameters as well.

Therefore, (56) becomes:

$f\left(r,\theta ,\varphi \right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}\left({\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left({\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{d}_p{R}_p\left(r\right)}\right)d_n{\Theta }_n\left(\theta \right)}\right)d_m{\Phi }_m\left(\varphi \right)}$ (65)

where $d_m$ is given by (55), $d_n$ is given by (63), and the coefficient $d_p$ is yet to be determined.

For a given ${\nu }_n$, one can verify that the set of functions $\left\{R_p\left(r\right),p\in \mathbb{N}\right\}$ is orthogonal in $\left[r_1,r_2\right]$ with respect to the weight function $w\left(r\right)=r^2$, i.e. we have:

$\begin{array}{l}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}\left[j_{{\nu }_n}\left({\lambda }_{p_1}{r}_1\right)y_{{\nu }_n}\left({\lambda }_{p_1}\rho \right)-y_{{\nu }_n}\left({\lambda }_{p_1}{r}_1\right)j_{{\nu }_n}\left({\lambda }_{p_1}\rho \right)\right]}{p}_1\ne p_2\\ \left[j_{{\nu }_n}\left({\lambda }_{p_2}{r}_1\right)y_{{\nu }_n}\left({\lambda }_{p_2}\rho \right)-y_{{\nu }_n}\left({\lambda }_{p_2}{r}_1\right)j_{{\nu }_n}\left({\lambda }_{p_2}\rho \right)\right]{\rho }^2\text{d}\rho =0\end{array}$ (66)

The set of functions $\left\{R_p\left(r\right),p\in \mathbb{N}\right\}$ is orthonormal with respect to the following norm:

${\Vert R_p\Vert }^2={\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\left[j_{{\nu }_n}\left({\lambda }_p{r}_1\right)y_{{\nu }_n}\left({\lambda }_p\rho \right)-y_{{\nu }_n}\left({\lambda }_p{r}_1\right)j_{{\nu }_n}\left({\lambda }_p\rho \right)\right]}^2{\rho }^2\text{d}\rho }$ (67)

Therefore, the coefficients $d_p$ are the Fourier coefficients of the expansion of $f\left(r,\theta ,\varphi \right)$ into a series of the combination of spherical Bessel functions defined by the function $R_p\left(r\right)$, for variable r and fixed $\varphi$ and $\theta$, with each $d_p$ given by:

$d_p=\frac{1}{{\Vert R_p\Vert }^2}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}f\left(\varphi ,\theta ,\rho \right)\left[j_{{\nu }_n}\left({\lambda }_p{r}_1\right)y_{{\nu }_n}\left({\lambda }_p\rho \right)-y_{{\nu }_n}\left({\lambda }_p{r}_1\right)j_{{\nu }_n}\left({\lambda }_p\rho \right)\right]{\rho }^2\text{d}\rho }$ (68)

Using classical identities verified by Bessel functions, one can obtain the following formula for ${\Vert R_p\Vert }^2$ :

${\Vert R_p\Vert }^2=\frac{J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_1\right)-J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)}{2r_1{\lambda }_p^4{J}_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)}$ (69)

Hence, the coefficients of the expansion of $f\left(r,\theta ,\varphi \right)$ into the Fourier-Legendre-Bessel triple series generated by the product of functions ${\Phi }_m\left(\varphi \right){\Theta }_n\left(\theta \right)R_p\left(r\right)$ are given by $d_{mnp}=d_m{d}_n{d}_p$, where $d_m,d_n$ and $d_p$ are given by (55), (63) and (68), respectively, and Lemma 2 ensues.

The combination of Lemma 1 and Lemma 2 entails Result 1.

Remark: using standard arguments, it can be proven that the continuity of the function $f\left(r,\theta ,\varphi \right)$ ensures uniform convergence of the triple series in (48), just as it ensures the existence of a unique solution in the classical sense of the IBVP (5)-(11).

3. Extensions

Several useful additional results can be inferred from Result 1. First, once can derive the solution to the non-homogeneous telegrapher’s equation, as given in Corollary 1.

Corollary 1

Let $h\left(r,\varphi ,\theta ,t\right)$ be a continuous function from the domain D in (6) to $ℝ$.

The solution of the nonhomogeneous equation:

$\begin{array}{c}\frac{{\partial }^{2}u}{\partial t^2}+{\kappa }_1\frac{\partial u}{\partial t}={\kappa }_2\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{2}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\left(\frac{{\partial }^{2}u}{\partial {\theta }^2}+\cot \theta \frac{\partial u}{\partial \theta }+{\csc }^2\theta \frac{{\partial }^{2}u}{\partial {\varphi }^2}\right)\right)\\ -{\kappa }_{3}u+h\left(r,\varphi ,\theta ,t\right)\end{array}$ (70)

that verifies the homogeneous boundary conditions stated in (7)-(11), is given by:

$\begin{array}{c}u\left(r,\theta ,\varphi ,t\right)={\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}f\left(\rho ,\vartheta ,\phi \right)G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t\right){\rho }^2\sin \left(\vartheta \right)\text{d}\rho \text{d}\vartheta \text{d}\phi }}}\\ +{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}h\left(\rho ,\vartheta ,\phi ,s\right)G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right){\rho }^2\sin \left(\vartheta \right)\text{d}\rho \text{d}\vartheta \text{d}\phi \text{d}s}}}}\end{array}$ (71)

where $G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t\right)$ is the Green’s function, which is immediately derived from Result 1:

$\begin{array}{l}G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t\right)\\ =\exp \left(-\frac{{\kappa }_1}{2}t\right){\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}\sin \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}\right)}}}\\ \times \frac{4{\beta }_m{r}_1{\lambda }_p^4{J}_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)}{{\left({\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)\right)}^{1/2}\left({\beta }_m\left({\varphi }_2-{\varphi }_1\right)-\sin \left(2{\beta }_m\left({\varphi }_2-{\varphi }_1\right)\right)\right)\left(J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_1\right)-J_{{\nu }_n+1/2}^2\left({\lambda }_p{r}_2\right)\right){\Vert P_{{\nu }_n}^{{\beta }_m}\Vert }^2}\\ \times \sin \left({\beta }_m\left(\varphi -{\varphi }_1\right)\right)\sin \left({\beta }_m\left(\phi -{\varphi }_1\right)\right)P_{{\nu }_n}^{{\beta }_m}\left(\cos \theta \right)P_{{\nu }_n}^{{\beta }_m}\left(\cos \vartheta \right)R_p\left(r\right)R_p\left(\rho \right)\end{array}$ (72)

which completes the statement of Corollary 1.

Endowed with the Green’s function, one can also obtain a formula for the Dirichlet problem with non-homogeneous boundary conditions.

Corollary 2

Let $u\left(r,\theta ,\varphi ,t\right)$ be defined on the domain D in (6) and satisfy PDE (70) as well as the following nonhomogeneous boundary conditions:

$u\left(r,\theta ,\varphi ,0\right)=f_1\left(r,\varphi ,\theta \right)$ (73)

${\frac{\partial u\left(r,\theta ,\varphi ,t\right)}{\partial t}|}_{t=0}=f_2\left(r,\theta ,\varphi \right)$ (74)

$u\left(r_1,\theta ,\varphi ,t\right)=f_3\left(\theta ,\varphi ,t\right)$ (75)

$u\left(r_2,\theta ,\varphi ,t\right)=f_4\left(\theta ,\varphi ,t\right)$ (76)

$u\left(r,{\theta }_1,\varphi ,t\right)=f_5\left(r,\varphi ,t\right)$ (77)

$u\left(r,\theta ,{\varphi }_1,t\right)=f_6\left(r,\theta ,t\right)$ (78)

$u\left(r,\theta ,{\varphi }_2,t\right)=f_7\left(r,\theta ,t\right)$ (79)

where $f_1,\cdots ,f_7$ are all continuous real functions such that:

$f_1\left(r,\theta ,\varphi \right)=f_3\left(\theta ,\varphi ,0\right)=f_4\left(\theta ,\varphi ,0\right)=f_5\left(r,\varphi ,0\right)=f_6\left(r,\theta ,0\right)=f_7\left(r,\theta ,0\right)$ (80)

$f_3\left({\theta }_1,\varphi ,t\right)=f_5\left(r_1,\varphi ,t\right),f_4\left({\theta }_1,\varphi ,t\right)=f_5\left(r_2,\varphi ,t\right)$ (81)

$f_3\left(\theta ,{\varphi }_1,t\right)=f_6\left(r_1,\theta ,t\right),f_4\left(\theta ,{\varphi }_1,t\right)=f_6\left(r_2,\theta ,t\right)$ (82)

$f_3\left(\theta ,{\varphi }_2,t\right)=f_7\left(r_1,\theta ,t\right),f_4\left(\theta ,{\varphi }_2,t\right)=f_7\left(r_2,\theta ,t\right)$ (83)

$f_5\left(r_1,{\varphi }_1,t\right)=f_6\left(r_1,{\theta }_1,t\right),f_5\left(r_1,{\varphi }_2,t\right)=f_7\left(r_1,{\theta }_1,t\right)$ (84)

$f_5\left(r_2,{\varphi }_1,t\right)=f_6\left(r_2,{\theta }_1,t\right),f_5\left(r_2,{\varphi }_2,t\right)=f_7\left(r_2,{\theta }_1,t\right)$ (85)

Then, the function $u\left(r,\theta ,\varphi ,t\right)$ is given by:

$u\left(r,\theta ,\varphi ,t\right)=\frac{\partial }{\partial t}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{f}_1\left(\rho ,\vartheta ,\phi \right)G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t\right){\rho }^2\sin \left(\vartheta \right)\text{d}\rho \text{d}\vartheta \text{d}\phi }}}$ (86)

$+{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}\left(f_2\left(\rho ,\vartheta ,\phi \right)+{\kappa }_1{f}_1\left(\rho ,\vartheta ,\phi \right)\right)G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t\right){\rho }^2\sin \left(\vartheta \right)\text{d}\rho \text{d}\vartheta \text{d}\phi }}}$ (87)

$+{\kappa }_2^2{r}_1^2{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}{f}_3\left(\vartheta ,\phi ,s\right){\frac{\partial G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right)}{\partial \rho }|}_{\rho =r_1}\sin \left(\vartheta \right)\text{d}\vartheta \text{d}\phi \text{d}s}}}$ (88)

$-{\kappa }_2^2{r}_2^2{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}{f}_4\left(\vartheta ,\phi ,s\right){\frac{\partial G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right)}{\partial \rho }|}_{\rho =r_2}\sin \left(\vartheta \right)\text{d}\vartheta \text{d}\phi \text{d}s}}}$ (89)

$-{\kappa }_2^2\sin \left({\theta }_1\right){\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{f}_5\left(s,\phi ,\rho \right){\frac{\partial G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right)}{\partial \vartheta }|}_{\vartheta ={\theta }_1}{\rho }^2\text{d}\vartheta \text{d}\phi \text{d}s}}}$ (90)

$+{\kappa }_2^2{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{f}_6\left(s,\vartheta ,\rho \right){\frac{\partial G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right)}{\partial \phi }|}_{\phi ={\varphi }_1}{\rho }^2\sin \left(\vartheta \right)\text{d}\vartheta \text{d}\phi \text{d}s}}}$ (91)

$-{\kappa }_2^2{\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{f}_7\left(s,\vartheta ,\rho \right){\frac{\partial G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right)}{\partial \phi }|}_{\phi ={\varphi }_2}{\rho }^2\sin \left(\vartheta \right)\text{d}\vartheta \text{d}\phi \text{d}s}}}$ (92)

$+{\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{0}{\overset{{\theta }_1}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{0}{\overset{t}{\int }}h\left(\rho ,\vartheta ,\phi ,s\right)G\left(r,\rho ,\varphi ,\phi ,\theta ,\vartheta ,t-s\right){\rho }^2\sin \left(\vartheta \right)\text{d}\rho \text{d}\vartheta \text{d}\phi \text{d}s}}}}$ (93)

which ends the statement of Corollary 2.

The theory underpinning the derivation of Corollary 1 and Corollary 2 is explained, e.g., in [12]. The role of the Green’s function is elaborated on in [13]. The conditions listed in (80)-(85) are necessary requirements for the solution $u\left(r,\theta ,\varphi ,t\right)$ given in (86)-(93) to be unique.

It should also be pointed out that Neumann or Robin boundary conditions can be easily handled in the same framework as the one used to solve the Dirichlet problem in Section 2. For example, if the radial boundary condition in (9) is replaced with the following:

${\frac{\partial }{\partial r}u\left(r,\theta ,\varphi ,t\right)|}_{r=r_1}={\frac{\partial }{\partial r}u\left(r,\theta ,\varphi ,t\right)|}_{r=r_2}=0$ (94)

Then, each ${\lambda }_p$ becomes the p-th root of the following equation:

${y^{\prime }}_n\left(\lambda r_2\right){j^{\prime }}_n\left(\lambda r_1\right)-{y^{\prime }}_n\left(\lambda r_1\right){j^{\prime }}_n\left(\lambda r_2\right)=0$ (95)

and the eigenfunctions of the Sturm-Liouville problem associated with $R\left(r\right)$ become:

$R_p\left(r\right)={j^{\prime }}_n\left({\lambda }_p{r}_1\right)y_n\left({\lambda }_{p}r\right)-{y^{\prime }}_n\left({\lambda }_p{r}_1\right)j_n\left({\lambda }_{p}r\right)$ (96)

Likewise, if the azimuthal angular boundary condition in (11) is replaced with the following:

${\frac{\partial u\left(r,\theta ,\varphi ,t\right)}{\partial \varphi }|}_{\varphi ={\varphi }_1}={\frac{\partial u\left(r,\theta ,\varphi ,t\right)}{\partial \varphi }|}_{\varphi ={\varphi }_2}=0$ (97)

Then, each ${\beta }_m$ becomes the m-th root of the following equation:

$\cos \left(\beta {\varphi }_2\right)\sin \left(\beta {\varphi }_1\right)-\cos \left(\beta {\varphi }_1\right)\sin \left(\beta {\varphi }_2\right)=0$ (98)

and the eigenfunctions of the Sturm-Liouville problem associated with $\Phi \left(\varphi \right)$ become:

${\Phi }_m\left(\varphi \right)={\beta }_m\left(b_1\sin \left(-{\beta }_m{\varphi }_1\right)\sin \left({\beta }_m\varphi \right)-b_2\cos \left({\beta }_m{\varphi }_1\right)\cos \left({\beta }_m\varphi \right)\right)$ (99)

Finally, the solution to the three-dimensional telegrapher’s problem in cylindrical coordinates can also be computed and is given in the following Result 2.

Result 2

Solution to the Dirichlet problem for the telegrapher’s equation with three space variables on a subset of a cylinder located within two radii, two angles and two heights.

Let $u\left(r,\varphi ,z,t\right)$ be the solution of the following partial differential equation:

$\frac{{\partial }^{2}u}{\partial t^2}+{\kappa }_1\frac{\partial u}{\partial t}={\kappa }_2\left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\varphi }^2}+\frac{{\partial }^{2}u}{\partial z^2}\right)-{\kappa }_{3}u$ (100)

on the domain:

$D=\left\{0\le r_1\le r\le r_2<\infty ,0\le {\varphi }_1\le \varphi \le {\varphi }_2\le 2\pi ,0\le z_1\le z\le z_2<\infty ,t\ge 0\right\}$ (101)

with the following boundary conditions:

$u\left(r,\varphi ,z,0\right)=0$ (102)

${\frac{\partial u\left(r,\varphi ,z,t\right)}{\partial t}|}_{t=0}=f\left(r,\varphi ,z\right)$ (103)

$u\left(r_1,\varphi ,z,t\right)=u\left(r_2,\varphi ,z,t\right)=0$ (104)

$u\left(r,{\varphi }_1,z,t\right)=u\left(r,{\varphi }_2,z,t\right)=0$ (105)

$u\left(r,\varphi ,z_1,t\right)=u\left(r,\varphi ,z_2,t\right)=0$ (106)

where ${\kappa }_1,{\kappa }_2,{\kappa }_3$ are three positive constants and f is a continuous function from D to $ℝ$.

Then, the function $u\left(r,\varphi ,z,t\right)$ is given by:

$u\left(r,\varphi ,z,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\displaystyle \underset{p=1}{\overset{\infty }{\sum }}{b}_{mnp}T\left(t\right){\Phi }_m\left(\varphi \right)R_n\left(r\right)Z_p\left(z\right)}}}$ (107)

The functions $T\left(t\right)$, ${\Phi }_m\left(\varphi \right)$, $R_n\left(r\right)$ and $Z_p\left(z\right)$ are defined as follows:

$T\left(t\right)=\exp \left(-\frac{{\kappa }_1}{2}t\right)\sin \left(t\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}\right)$ (108)

${\Phi }_m\left(\varphi \right)=\sin \left({\beta }_m\left(\varphi -{\varphi }_1\right)\right)$ (109)

$R_n\left(r\right)=Y_{{\beta }_m}\left({\mu }_{n}r\right)J_{{\beta }_m}\left({\mu }_n{r}_1\right)-Y_{{\beta }_m}\left({\mu }_n{r}_1\right)J_{{\beta }_m}\left({\mu }_{n}r\right)$ (110)

$Z_p\left(z\right)=\sin \left(\sqrt{{\lambda }_p^2-{\mu }_n^2}\left(z-z_1\right)\right)$ (111)

The numbers ${\beta }_m,{\mu }_n$ and ${\lambda }_p$ are defined as follows:

· $\forall m\in \mathbb{N}$, ${\beta }_m$ is the m-th root of the following equation

$\sin \left(\beta {\varphi }_2\right)\cos \left(\beta {\varphi }_1\right)-\sin \left(\beta {\varphi }_1\right)\cos \left(\beta {\varphi }_2\right)=0$ (112)

· $\forall n\in \mathbb{N}$, ${\mu }_n$ is the n-th root of the following equation

$Y_{{\beta }_m}\left(\mu r_2\right)J_{{\beta }_m}\left(\mu r_1\right)-Y_{{\beta }_m}\left(\mu r_1\right)J_{{\beta }_m}\left(\mu r_2\right)=0$ (113)

· $\forall p\in \mathbb{N}$, ${\lambda }_p$ is the p-th root of the following equation

$\sin \left(\sqrt{{\lambda }^2-{\mu }_n^2}{z}_2\right)\cos \left(\sqrt{{\lambda }^2-{\mu }_n^2}{z}_1\right)-\sin \left(\sqrt{{\lambda }^2-{\mu }_n^2}{z}_1\right)\cos \left(\sqrt{{\lambda }^2-{\mu }_n^2}{z}_2\right)=0$ (114)

The coefficient $b_{mnp}$ is given by:

$\begin{array}{c}{b}_{mnp}=\frac{8{\beta }_m\sqrt{{\lambda }_p^2-{\mu }_n^2}}{\left({\beta }_m\left({\varphi }_2-{\varphi }_1\right)-\sin \left(2{\beta }_m\left({\varphi }_2-{\varphi }_1\right)\right)\right)\left(\sqrt{{\lambda }_p^2-{\mu }_n^2}\left(z_2-z_1\right)-\sin \left(2\sqrt{{\lambda }_p^2-{\mu }_n^2}\left(z_2-z_1\right)\right)\right)}\\ \times \frac{J_{{\beta }_m}^2\left({\mu }_n{r}_1\right)-J_{{\beta }_m}^2\left({\mu }_n{r}_2\right)}{{\pi }^2{\mu }_n^2{J}_{{\beta }_m}^2\left({\mu }_n{r}_2\right)\sqrt{{\kappa }_1^2-4\left({\kappa }_2{\lambda }_p^2-{\kappa }_3\right)}}\\ \times {\displaystyle \underset{{\varphi }_1}{\overset{{\varphi }_2}{\int }}{\displaystyle \underset{r_1}{\overset{r_2}{\int }}{\displaystyle \underset{z_1}{\overset{z_2}{\int }}f\left(\phi ,\rho ,\xi \right){\Phi }_m\left(\phi \right)R_n\left(\rho \right)Z_p\left(\xi \right)\rho \text{d}\xi \text{d}\rho \text{d}\phi }}}\end{array}$ (115)

Proof of Result 2 is omitted for the sake of brevity, as the solution process is quite similar to the one leading to Result 1, while being simpler. In cylindrical coordinates, non-homogeneous equations, non-homogeneous boundary conditions, as well as Neumann and Robin conditions, can be handled in the same way as in spherical coordinates.

4. Conclusions

This article has shown how to obtain a closed form solution to the non-stationary version of one of the most general second order linear PDEs of mathematical physics, namely the telegrapher’s equation, in a three-dimensional spherical domain, under general boundary conditions and flexible specifications of the ranges of the various parameters. It is quite remarkable that such a sophisticated computational problem admits an exact analytical solution as a uniformly convergent series of triple integrals involving functions whose numerical evaluation is relatively easy. Although most of the problems encountered in the physical sciences are in three dimensions in space, it is worth noticing that it would be a simple extension of this work to derive an exact solution to the same problem with four space variables instead of three, i.e. in a spherical domain defined as follows:

$D=\left\{0<r_1\le r\le r_2<\infty ,0<\theta \le {\theta }_1\le \pi ,0<{\theta }^{\prime }\le {\theta }_2\le \pi ,0<{\varphi }_1\le \varphi \le {\varphi }_2\le 2\pi ,t\ge 0\right\}$ (116)

where ${\theta }^{\prime }$ is an additional angle, i.e. the relations between the spherical coordinates and the rectangular ones, denoted by $x_1,x_2,x_3,x_4$ are now given by:

$x_1=r\sin {\theta }^{\prime }\sin \theta \cos \varphi ,x_2=r\sin {\theta }^{\prime }\sin \theta \sin \varphi ,x_3=r\sin {\theta }^{\prime }\cos \theta ,x_4=r\cos {\theta }^{\prime }$ (117)

However, from a practical perspective, one would have to weigh the benefit of having an analytical formula, relative to the cost of increased difficulty of numerical evaluation due to the greater dimension of the integration process.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Zauderer, E. (2006) Partial Differential Equations of Applied Mathematics. 3rd Edition, Wiley-Interscience, Hoboken, New Jersey. https://doi.org/10.1002/9781118033302
[2]	Ockendon, J., Howison, S., Lacey, A. and Movchan, A. (2003) Applied Partial Differential Equations. Oxford University Press, Oxford.
[3]	Myint-U, T. and Debnath, L. (2007) Linear Partial Differential Equations. 4th Edition, Birkhäusen, Boston.
[4]	Polyanin, A.D. and Nazaikinskii, V.E. (2016) Handbook of Linear Partial Differential Equations for Engineers and Scientists. 2nd Edition, CRC Press, Boca Raton. https://doi.org/10.1201/b19056
[5]	Stakgold, I. (2000) Boundary Value Problems of Mathematical Physics. SIAM, Philadelphia. https://doi.org/10.1137/1.9780898719888
[6]	Tikhonov, A.N. and Samarskii, A.A. (1990) Equations of Mathematical Physics. Dover Publications, New York.
[7]	Zettl, A. (2005) Sturm-Liouville Theory. American Mathematical Society, Providence.
[8]	Lebedev, N.N. (1972) Special Functions and Their Applications. Dover Publications, New York.
[9]	Lew Yan Voon, L.C. and Willatzen, M. (2003) Angular Confinement, Non-Integral Quantum Numbers and Controllable Degeneracies. Europhysics Letters, 62, 299-305. https://doi.org/10.1209/epl/i2003-00395-x
[10]	Andrews, G.E., Askey, R. and Roy, R. (2006) Special Functions. Cambridge University Press, Cambridge.
[11]	Bremer, J. (2018) An Algorithm for the Numerical Evaluation of the Associated Legendre Functions That Runs in Time Independent of Degree and Order. Journal of Computational Physics, 360, 15-38. https://doi.org/10.1016/j.jcp.2018.01.014
[12]	Dezin, A.A. (1987) Partial Differential Equations. An Introduction to the General Theory of Linear Boundary Value Problems. Springer-Verlag, Berlin-New York.
[13]	Duffy, D.G. (2001) Green’s Functions with Applications. Chapman & Hall/CRC, Boca Raton.