Fourier method for the Neumann problem on a torus

The Fourier method approach to the Neumann problem for the Laplacian operator in the case of a solid torus contrasts in many respects with the much more straight forward situation of a ball in 3-space. Although the Dirichlet-to-Neumann map can be readily expressed in terms of series expansions with toroidal harmonics, we show that the resulting equations contain undetermined parameters which cannot be calculated algebraically. A method for rapidly computing numerical solutions of the Neumann problem is presented with numerical illustrations. The results for interior and exterior domains combine to provide a solution for the Neumann problem for the case of a shell between two tori.


Introduction
The Dirichlet-to-Neumann map for the Laplacian is important in many areas of basic analysis (elliptic boundary value problems [10,23,27,32], inverse problems [5,20]), as well as physics (e.g., fluid mechanics [9], electromagnetic theory [6], electrical impedance tomography [19,22,35], electrical transmission [11]).The map relates the values of a function f : ∂ → R (Dirichlet datum) on the boundary ∂ ( ⊆ R 3 ) to the normal derivative (Neumann datum) of the unique harmonic function u : → R whose boundary values are given by f .The Neumann problem is to recover f (or u) from the normal derivative.
When is a sphere, solving the Dirichlet problem amounts to copying the boundary function's Fourier coefficients in terms of spherical harmonics to obtain the coefficients of the interior harmonic function in terms of solid spherical harmonics.Thus it would make no sense to apply a general-purpose method (such as finite elements) for the Dirichlet problem on the sphere; the same holds for the Neumann problem (cf.Remark 3.7 below).
When is a torus, harmonic functions defined on are naturally expressed as series based on a doubly-indexed collection of toroidal harmonics (involving associated Legendre functions of the first and second kinds of half-integer degree) [17], which are orthogonal with respect to a certain weighted L 2 -inner product.As for the sphere, the coefficients of this expression provide a solution of the Dirichlet problem for the Laplacian quite directly.However, unlike the case for a sphere, the Dirichletto-Neumann mapping for a torus turns out to be far more complicated, and the solution of the Neumann problem involves solving an underdetermined infinite upper-triangular system of linear equations.We express the well-known necessary and sufficient condition for the solvability of the Neumann problem (compatibility condition), as well as the normalization condition, in terms of the Fourier coefficients.The solution of the Neumann problem turns out to involve a special twist in that the free (i.e., undetermined) parameter which appears in the linear system cannot be found algebraically, as far as we know (it is not connected to the normalization condition).However, it can be expressed as a limit of easily calculated algebraic expressions.The solution procedure is illustrated through numerical examples.
A similar phenomenon was found by Hicks [18] in calculating the velocity potential for the fluid motion of a torus moving parallel to a fluid at rest at infinity, in which the unique free parameter leading to convergence of the series was determined in terms of another infinite series.Love [26] approached an inhomogeneous Neumann problem by means of Green's functions for second-order difference equations.In both cases, the construction was based on an expansion of the potential functions in terms of toroidal harmonics.
Of course, one can solve the numerical problem by approximating the normal derivative numerically, but this approach provides little insight into the structure of the solutions (cf.Proposition 5.4 and Corollary 5.5 below.)As far as we know, our explicit formula for the normal derivative in toroidal coordinates has not appeared previously.We have seen applications (e.g., [6]) where the potential in a toroidal conductor is modeled in spherical coordinates, which are not ideally suited for such problems.A method for solving the Neumann problem in a more general class of regions termed "toroidal" appeared in [28], but the notion of "toroidal coordinates" which we use here is not applicable to that method.
The paper is organized as follows.Some properties of toroidal harmonic functions are summarized in Sect. 2. In Sect.3, a computation of the toroidal Neumann derivative leads to a formula for the Dirichlet-to-Neumann mapping.The expansion coefficients of the normal derivative are linear expressions in the Fourier coefficients on the surface of the torus.While the mapping exists in the context of the appropriate Sobolev function spaces associated with the Laplace equation and the boundary data, as well as the trace operator, we only justify the derivation of our formulas for the normal derivative under stronger smoothness assumptions.Sections 4-6 present the algorithm and numerical examples to show the accuracy of the procedure.The computational complexity is shown to be O(N M) in a specific sense, far faster than a general-purpose method such as finite elements.Section 7 explains how to solve the Neumann problem in a toroidal shell using the results for the interior and exterior.

Toroidal coordinates and toroidal harmonics
We begin with the notation and several well-known facts to be used throughout the paper.
Since toroidal coordinates form an orthogonal coordinate system, we have also n = x η /|x η |, and since η increases as points move towards the interior of η 0 , we see that n is in fact the inward pointing normal vector on ∂ η 0 .The normal derivative of a function f defined in a neighborhood V of a point x ∈ ∂ η 0 (or a half-neighborhood • n, and by orthogonality of the coordinate system the normal derivative is also equal to the following, which is often more convenient for calculations: (3)

Toroidal harmonics
The associated Legendre functions (Ferrer's functions) of the first and second kinds for t > 1 are defined for integer values of n, m ≥ 0, respectively, as where P n (t), Q n (t) denote the classical Legendre polynomials of degree n [1,2,4,8,12,14,17,21,34,36].Extended analytically in the complex plane, these functions are branched at t = ±1.However, they are entire functions of n and m regarded as complex variables [17].In this sense, (4) can be taken as a definition of the associated Legendre functions for half-integer values of n as we will need here.We will abbreviate ν n (θ ) = cos nθ for ν = 1, and ν n (θ ) = sin nθ for ν = −1.The (normalized) interior toroidal harmonic functions are for integers n, m with A derivation of the Laplacian equation in toroidal coordinates and the verification that I ν,μ n,m is harmonic in η 0 can be found in [17, p. 434]; standard results on removable sets for harmonic functions [3] cover the behavior at the singular set.
One may define a weighted L 2 inner product on real-valued functions by with the weight function One deduces the following using the fact that the collection of products [31,Sec. 25]).

Proposition 2.1 The interior toroidal harmonics {I
ν,μ n,m } (for n, m, ν, μ as in (6)) form a complete orthogonal system in the harmonic functions in L 2 ( η 0 , w).Their norms are where ε n = 1 + δ n,0 and δ n,m is the Kronecker delta function.Further, the restrictions to the boundary

Background on the Neumann problem
Given a suitable f : ∂ η 0 → R, the Dirichlet-to-Neumann mapping is described by finding the unique harmonic function u ∈ Har( η 0 ) with boundary values f = u| ∂ η 0 , and then taking the normal derivative h = nor u.The mapping is thus f = h.A common setting [25] is for u to be in the Sobolev space H 1 ( η 0 ) and f in the boundary space H 1/2 (∂ η 0 ).Roughly speaking, H 1 ( η 0 ) consists of L 2 functions having L 2 derivatives, and ) is identified with the space of boundary values of elements of H 1 ( η 0 ), where H 1 0 ( ) denotes the closure of the subspace of functions of compact support.The Neumann problem is to find a function f such that f = h, for a given function h defined on ∂ η 0 .The solution of the Neumann problem, in general, is guaranteed by the following result [10,13,29], valid for domains with sufficiently smooth boundary.
Then there exists an f ∈ H 1/2 (∂ ) such that f = h.This solution is unique up to an additive constant.
The solution f can made unique by requiring the normalization condition for a chosen constant c.

Dirichlet-to-Neumann mapping in toroidal coordinates
The coordinate expression of the Dirichlet-to-Neumann mapping will be based on the following calculation.

Lemma 3.2
The normal derivatives of the interior toroidal harmonics are given by the formula Proof Apply (3)-( 5) and obtain a formula for nor I Finally, regroup the terms as multiples of the trigonometric expressions of the form ν n .
The boundary function f may be represented conveniently for our purposes via as per Proposition 2.1, with a ν,μ n,m ∈ R. Unless otherwise specified, the indices of summation will always be as in (6) but excluding the cases of (n, m, ν, μ) being (0, m, −1, μ) or (n, 0, ν, −1), taking into account that − 0 = 0 identically.For convenience, we will often write superscripts as "+" in place of 1 and "−" in place of −1.
With this notation, the solution of the Dirichlet problem on η 0 is trivial: Then by Lemma 3.2, we have h = f is in turn given by assuming that f is sufficiently well-behaved to justify the exchange of summation and differentiation.For example, since the trace operator tr : Alternatively, the convergence is valid under the assumption that the sum in (12) together with the sums of the partial derivatives of the terms converge uniformly on compact subsets of η 0 .
We will show how nor u can be expressed in terms of the coefficients of f and certain constants defined in terms of Legendre functions.We will use the abbreviations and τ n,m = τ n,m (η 0 ) associated to η 0 are defined as follows: for all m ≥ 0.
We will use the following well known facts about the growth of Fourier coefficients and the convergence of Fourier series.
for some constant C > 0. Conversely, if (16) The main result of this section is that the coefficients (15) are entries in an infinite upper-triangular matrix which defines the Dirichlet-to-Neumann mapping .Theorem 3.5 For a fixed η 0 , let f : ∂ η 0 → R given by (11) and suppose that the Fourier coefficients satisfy (16) for some constant C where r ≥ 4. Then the coefficients in the absolutely convergent series for the Dirichlet-to-Neumann mapping h = f ∈ C r −3 are given by b for μ = ±1 and all m ≥ 0.
Proof From ( 5) it follows that I ν,μ n,m (η 0 , θ, ϕ) ≤ √ t 0 + 1, so the expansion (11) converges absolutely.In [17, p. 305] it is shown that for fixed η 0 and m, as n → ∞, where ∼ means that the ratio of the two expressions tends to 1. From this it is seen that independently of the value of m.From (20) and Lemma 3.2 one verifies that for some constant C 1 (which depends on η 0 ), and hence by (16), This is enough to guarantee that (13), a double series in θ and ϕ, also converges absolutely.(Our hypothesis on the value of r takes into account that 1/(m+n+1) 2 = ∞.)This in turn permits us to substitute the formula for nor I ν,μ n,m into (13) and then reindex which is (17).By Definition 3.3 and with observation that t 0 − s 0 = e −η 0 it also follows from (20) that Thus by (22), the Fourier coefficients b ν,μ n,m defining h are of order no greater than 1/((m + n + 1) r −1 ), so by Proposition 3.4 we are done.Remark 3. 6 While the differentiability requirement in Theorem 3.5 may appear quite restrictive, we will see that numerically it is often not necessary.It is not difficult to verify that when h in Proposition 3.1 is real analytic, the solution f to the Neumann problem is also real analytic.Since we are mainly interested in the formal relationships among the coefficients we will not go deeper into relaxing the differentiability conditions.Remark 3.7 When one applies the same principle underlying Theorem 3.5 to the Neumann problem on the sphere |x| < 1 in spherical coordinates y 0 = ρ cos θ , y 1 = ρ sin θ cos φ, y 2 = ρ sin θ sin φ, and uses the standard solid spherical harmonics as the basis for the harmonic functions, it is natural to represent the Dirichlet and Neumann functions f (θ, ϕ) and h(θ, ϕ) with the well-known basis {Y ± n,m (1, θ, ϕ)} for L 2 functions on the sphere [32].This yields certain coefficients a ± n,m and b ± n,m for f and h respectively.Since the normal derivative in this situation is equal to the radial derivative, i.e. nor u = (∂u/∂ρ)| ρ=1 , one sees immediately that b ± n,m = na ± n,m , so the spherical analogue of the system (18) is rather trivial.(In the text [24, p. 218], this approach is suggested in a remark after expressing the solution of the Neumann problem for a spheroid, but only for zonal harmonics: functions constant with respect to the angular coordinate, i.e., involving P n but not general P m n .Apart from this, we have not seen this type of solution presented in the literature.)

Algebraic solutions for the Fourier coefficients in the Neumann problem
We now focus on the problem of solving for the coefficients a n,m of h.The upper triangular system (18) can be solved almost trivially by back-substitution.It separates naturally into a separate system for each fixed set of parameters (m, ν, μ), corresponding to the separate modes of the problem.For ν = +1, the solution is Thus one takes the values a 1,m , respectively, to be the only free parameters.Of course, the following fact is necessary for the justification of these statements.Proof In [17, p. 195] it is shown that (in fact, this rather than ( 4) is used to give a definition of Q m n (t)).From this it follows that for all n, m, and η ∈ R + .From (15), we need to show that the value 16s 0 q n+1,m τ n,m = 4 (2n + 3)t 0 q n+1,m + 2(n − m) − 3 q n+2,m (24) does not vanish.Consider the recursion formulas from [4, pp.161-162] and [17, p. 108] Applying (25) to (24), we find Now add and subtract (2n + 3)(2(m + n) + 1)q n+2,m and use (25), yielding 16 s 0 q n+1,m τ n,m = −2(2n + 3)(2n + 2m + 1)q n+1,m−1 + 8 m q n+2,m , which by ( 23) is never zero.An algebraic solution will give rise to a solution of the Neumann problem if it defines a convergent series (11).

Convergent solutions for the toroidal Neumann problem
We observed above that the Dirichlet problem is effectively trivial when data is expressed in terms of the basis of toroidal harmonics.We will now see that for the Neumann problem, some perhaps surprising particularities occur in connection with the convergence properties of algebraic solutions.
The question is when the algebraic solution determined by choice of the initial coefficients a +,μ 0,m and a −,μ 1,m will provide a convergent series in (11).That is to say, one wants to know when the function corresponding to the mode (m, ν, μ) will exist.The principal mode will be considered separately from the remaining modes.

Convergence of the algebraic solution for the principal mode
For the particular indices (m, ν, μ) = (0, +1, +1), the Neumann constants satisfy some special relations which we will need.
Similarly, referring again to Definition 3.3 we find where the last step follows from (25).Finally, ρ n,0 q n−1,0 + σ n,0 q n,0 + τ n,0 q n+1,0 n,0 = 0 for all n (i.e., when h vanishes in the principal mode), the corresponding equations (18) are linear homogeneous, and Lemma 5.1 implies that a +,+ n,0 q n,0 = 2 a +,+ 0,0 q 0,0 (n ≥ 1); i.e., a +,+ n,0 /q n,0 = ε n a +,+ 0,0 /q 0,0 for all n.On comparing the formula of Proposition 5.2 with the exponent determined by α = 1/2, one sees that the resulting solution is indeed equal to a constant function on ∂ η 0 (as must be the case according to Proposition 3.1) with value (π/ √ 2)(a +,+ 0,0 /q 0,0 ).The solution of the Dirichlet problem in η 0 is the same constant.Thus, for every a +,+ 0,0 ∈ R, the algebraic solution gives a convergent series ∞ n=0 a +,+ n,0 I +,+ n,0 (η 0 , θ, ϕ).In the rest of this subsection we give some results which shed some further light on the special nature of the principal mode.First, we show that the compatibility condition and the normalization depend only on the coefficients of the principal mode.We will need the following series expansion.

Proposition 5.2 ([7]
) For all α ∈ C, The compatibility condition (9) applied to h of the form (17) is equivalent to (b) The normalization condition (10) applied to f of the form (11) is equivalent to with the last equality following from Proposition 5.2 with α = 3/2.The proof of (b) follows the same lines as (a).
The following application is of independent interest.

Corollary 5.5 Let f be a particular solution of f = h and set c
Let f be obtained by replacing the coefficients a +,+ n,0 for f with Then f is the unique solution of the Neumann problem which satisfies the normalization condition (10).

Determination of parameter for convergence for non-principal modes
We assume now that (m, ν, μ) = (0, +1, +1).According to Proposition 3.1, there is exactly one algebraic solution of (18) for which (11) converges.The determination of the corresponding unique value a opt of a +,μ 0,m or a −,μ 1,m is not an algebraic question; in fact, we believe that a opt is not the solution of any algebraic or transcendental equation associated with the Neumann data.We will assume that μ = +1 since the case μ = −1 is analogous, the only difference being the initial indexing from n = 1 instead of n = 0.
To simplify the notation, the mode being fixed we will write a n = a ν,μ n,m and let A n (a) denote the value of a n in the solution of equations (18) determined by setting the arbitrary parameter a 0 = a +,μ 0,m (or a −,μ 1,m ) equal to a. Thus A 0 (a) = a, and by a simple induction, one finds recursively defined linear expressions where By construction, every value of a ∈ R gives an algebraic solution {A n (a)} of the system (18).According to (11) and Theorem 3.5, we want to find the unique value a opt provided by Theorem 5.6 for which the series converges absolutely and thus gives a boundary function f ν,μ m (θ, ϕ).In particular, it is necessary that A n (a opt ) → 0 as n → ∞.By (31), this says C n a opt + D n → 0.
Note that the C n do not depend on the data {b n }.It is clear that two consecutive terms C n , C n+1 can never vanish.Under the assumption that the C n do not tend to zero, i.e., C n > > 0 for infinitely many n, we have as n → ∞ on that subsequence.We will look further into the question of the asymptotics of C n in the following section.We summarize our conclusions as follows.
The contrast of the behavior of the principal mode in this theorem is striking, since Eq.(18) show no apparent structural difference between the principal and other modes.The algorithm is essentially as follows:

Calculation of q n,m
From (19) it is clear that q n,m decreases exponentially as n → ∞.In fact, for many values of η 0 in the range we will be considering, |q 100,0 | is less than 10 −300 and thus effectively becomes zero when represented in standard IEEE double-precision format.For many problems values of n as large as 100 would not be necessary, but if one wants to use a large number of terms in a series solution, two ways out of this difficulty present themselves.
The first way is via asymptotic expansions.For large values of n we can estimate q n,m to high accuracy by means of ( 19) and Stirling's formula for the function.However, it turns out there is a gap between small values of n where machine precision is applicable and the large values of n for which Stirling's estimate is sufficiently accurate.For example, for η 0 = 1.5 and m = 2, we find in machine precision that q 131,2 vanishes, while q 130,2 suffers a relative error of 0.0006.The estimate based on (19) gives a relative error of 0.003.The formula (19) is in fact the first term of an asymptotic expansion for Q m n given in [17], and can be made much more accurate by multiplying it by the factor Note that by (15), we only need to know the ratios q n+1,m /q n,m .After cancelling several factors, one derives an approximation which can be substituted into (15) to find, for example, that τ 136,2 vanishes in machine precision, while the approximations of τ n,2 already for n ≥ 20 have relative error less than 0.00001.Similar statements hold for the coefficients ρ n,m and σ n,m .The optimal cutoff value of n at which to switch from a direct calculation to the asymptotic formula would depend on the values of m and η 0 as well as a consideration of the acceptable error in the coefficients.It is likely that other more convenient asymptotic expressions could be developed but this would lead us rather far afield.Multiple-precision arithmetic was used in Mathematica to verify the correctness of the above statements.Multiple precision is also a convenient second alternative for calculating all of the Neumann constants when one is planning to solve many Neumann problems on the same torus: once the coefficients have been obtained and rounded to machine precision, the computation will then easily move forward without further recourse to higher precision arithmetic.

Numerical behavior of C n
We observed that a opt is given by (34) unless it happens that C n → 0. (Recall that C n does not depend on the Neumann data {b ν,μ n,m }.)For small values of n, we have a priori little control over even the sign of the coefficients defined in (15).However, from (22), ρ n,m /τ n,m → 1 and σ n,m /τ n,m → −2t 0 .Therefore if for a single large n we have then by (32) it would follow that i.e., the sequence {C n } grows exponentially.Table 1 lists calculated values of C n corresponding to m = 1 and a range of values of η 0 .Other values of m are shown in Table 2.Even though the initial values can decrease, in all cases that we have examined it appears that C n → ∞ exponentially as n → ∞, and we conjecture that this always holds.Then Algorithm 5.7 is applicable.

Examples of Dirichlet-to-Neumann and Neumann problem calculations
We illustrate the solution of the Neumann problem with numerical examples.For the principal mode (m, ν, μ) = (0, +1, +1) there is not much to be said.Every choice of a 0 = a ++ 0,0 gives an algebraic solution which defines a series which solves the Neumann problem.Corollary 5.5 may be used for normalization.It may be worth noting that when b n = 0 for all n (i.e., the coefficients for a vanishing normal derivative in the mode under consideration), since the sum is a constant, one can dispense with the series altogether for this case.But if one does wish to follow the procedure, the formulas (32) give D n = 0 always, and by (31), we have A n (a) = aC n .Choosing a 0 = 1 without any loss of generality and fixing η 0 , one obtains the values C n by (32) and then the initial coefficients a n by (31).This calculation amounts simply to finding values of the associated Legendre functions of the second kind via the classical recursion formulas, and the only numerical error is that which accumulates due to roundoff.
The numerical examples were calculated in Mathematica on a household laptop computer.All calculations were in machine precision except a few of the coefficients as noted.The calculations took a fraction of a second apart from the graphics.

Example 1
To illustrate the calculation of the Dirichlet to Neumann mapping, consider a +,+ n,m = a n where a n = (−1) (n−1)/2 n −2 for n odd, and a n = 0 for n even, 0 ≤ n ≤ 50. Figure 1 shows the function f and the image h = f for a few values of m.While f is smooth, it does not satisfy the regularity condition of Theorem 3.5, and h has jump discontinuities.Except at the jumps, the values of h agree to within 10 −14 to 10 −15 with the normal derivative of the f series obtained by numerical derivation.
As another example, let instead a n = (−1) (n−1)/2 n −1 for n odd.As Fig. 2 shows, now f has a jump discontinuity while the series for the resulting h does not converge, further illustrating the tendency of the Dirichlet to Neumann mapping to reduce the degree of differentiability.
It is readily checked that u is harmonic and (ν, μ) = (+, +).200-digit precision was used to avoid underflow in the calculations We substitute these coefficients into (18) to obtain numerical values for the b ν,μ n,m .Then we compare truncations of the series (21) with the true values of h = nor u according to (36). Figure 3 displays the base-10 logarithm of the absolute error for different combinations of m and η 0 .As is expected, the error is reduced when the number of terms in the series increases.It is also seen that the error increases steadily when larger values of m and η 0 are used.

Example 3
We illustrate our algorithm for solving the Neumann problem using the same function u as in the previous example.The Fourier coefficients b ν,μ n,m are obtained by numerical integration.Then the auxiliary coefficients C n , D n are obtained recursively by (32), and then a opt is approximated by the last value of −C n /D n according to (34).One would expect that the values a ν,μ n,m = A n (a opt ) of (34) provide a convergent series, while for a = a opt , {A n (a)} would not.This is confirmed by Fig. 4, which shows the values of A n (a + ) for small values of .The error in a particular series solution h of the Neumann problem compared to (36) is shown in Fig. 5. Maximum errors for combinations of η 0 , m are shown in Table 3. n,m } respectively.Methods for calculating these coefficients given the values of f or h on ∂ η 0 are well known, depend on a choice of numerical integration procedure, can be found in many standard software packages, and will not be discussed here.In this sense, Algorithm 5.7 can be considered as one part of a larger computation.
The amount of calculation required for Algorithm 5.7 is simple to estimate.Take N terms in the truncated Fourier series and use the values C N , D N to determine a opt .Assume that the Neumann constants ρ n,m , σ n,m , τ n,m have been calculated previously.This is also reasonable when one wishes to solve many Neumann problems on a single surface.Then by (32), approximately 4N multiplications in machine-precision are needed to find a opt , and since we then have all of the C n , D n for n < N , by (31) we need N more multiplications to find the Taylor coefficients for each mode.The 7 Exterior toroidal domain and toroidal shells

Exterior domain
The formula for the normal derivative of an exterior harmonic function u ∈ Har( * η 0 ) = {x : 0 ≤ η < η 0 }, and the solution of the corresponding Neumann problem are quite analogous to that of the interior domain η 0 .The exterior harmonics E ν,μ n,m ∈ Har( * η 0 ), defined by are obtained from the interior harmonics (5) by writing P m n−1/2 (cosh η) in place of Q m n−1/2 (cosh η) and are orthogonal with the same weight function (8) formula but applied in * η 0 .These Legendre functions of the first and second kinds satisfy identical recurrence relationships [17].For this reason, one finds that nor E ν,μ n,m is obtained from the formula of Lemma 3.2 by replacing similarly Q m sn−1/2 (cosh η) with P m n−1/2 (cosh η).Since the solution of the Dirichlet problem in * η 0 with boundary condition f given by ( 11) is one finds that the normal derivative of f will be given by equations (18) when q n,m is replaced in (15) with The method we have described is then applicable with no essential changes for solving the Dirichlet-to-Neumann problem in analogous to the Laurent series for holomorphic functions in an annular domain in the complex plane, converging uniformly in proper closed subdomains.(Note, however, that the inner and outer harmonics together do not form an orthogonal system in η ext ,η int .)A boundary function f : ∂ → R is given collectively by its values for η = η int and η = η ext collectively, let us say For u to be the solution of the Dirichlet problem for f , we combine (42) with (43) to find This might be written symbolically as To solve this system, one needs to verify that it is nonsingular.Instead of a direct verification as in Lemma 4.1, we simply note that if for even one combination of (n, m, ν, μ) there were more than one solution, one could easily construct a Dirichlet problem in the shell with more than one solution.
We see that nor I As in the solution of the Neumann problem for the interior domain, the equations for a fixed value of (m, ν, μ) are independent of those for another value of these indices.They can be solved recursively.The only difference will be that one must solve a pair of equations at each step.

Conclusions
We have presented an approach to studying the Dirichlet-to-Neumann mapping and solving the Neumann problem for the Laplacian on a torus.It is shown how the Dirichlet-to-Neumann mapping may be expressed by means of certain infinite series based on toroidal harmonics.We describe the well-known necessary and sufficient condition for the solvability of the Neumann problem (compatibility condition), as well as the normalization condition in terms of the Fourier coefficients.These results show that in this context the Neumann problem involves an infinite upper-triangular system of linear equations.The solution of the problem involves a special twist in that the unique value of the free parameter in this underdetermined linear system which truly gives a solution, cannot be found algebraically.Therefore we express it as a limit of easily calculated algebraic expressions.Numerical results are displayed for the accuracy of the algorithm.While the theoretical justification requires a fairly strict smoothness assumption, the method appears applicable in great generality.The Neumann problem is also solved in a toroidal shell.
, Now eliminate all occurrences of the derivative (Q m n+1 ) by means of the recurrence formula[4, pp.161-162]

Lemma 4 . 1
The values τ n,m are never zero.

Definition 4 . 2
We will say that a collection of real numbers {a ν,μ n,m } is an algebraic solution of the Neumann problem posed by {b ν,μ n,m } when all of the equations (18) are satisfied.

Fig. 3
Fig.3Base-10 logarithm indicating number of significant figures of approximation of the Dirichletto-Neumann mapping given by Eq. (18) truncating the series (21) to 0 ≤ n ≤ N for varying values of N .Accuracy is lost as m or η 0 increases.For comparison purposes 100-digit precision was used to obtain the Neumann constants.It was found that machine precision was sufficient for η 0 > 0.5, which is applicable for most "reasonably-shaped" toroidal domains

.
| ∂ int is obtained from the formula of Lemma 3.2 with η 0 replaced with η int , while nor I ν,μ n,m | ∂ ext is obtained by using η ext instead.The boundary values nor E ν,μ n,m | ∂ int and nor E ν,μ n,m | ∂ ext are then obtained by replacing Q m n−1/2 with P m n−1/2 .Once we have the harmonic function u as in (42), we have then nor u When the convergence of the series is absolute, one may apply the same rearranging and reindexing as described in the proof of Theorem 3.5 to obtain the coefficients in the Dirichlet-to-Neumann mapping h = f , h(η int , θ, ϕ) = cosh η 0 − cos θ b int ν,μ n,m ν n (θ ) μ m (ϕ), h(η ext , θ, ϕ) = cosh η 0 − cos θ b ext ν,μ

Table 3
Significant figures in the numerical solution of the Neumann problem on the torus showing the increase in accuracy with the number of terms We now combine our results in the interior and exterior domains for solving the Neumann problem in a toroidal shell.Let η int < η ext .Common to an interior and an exterior domain are the points of the toroidal shell= η int ,η ext = η ext ∩ * η int .A general harmonic function u in η ext ,η int and continuous in the closure can be expressed via an integral of its boundary values over ∂ η ext ,η int using the Poisson kernel for the torus [16, Ch. 1].This integral is the difference of the integrals over ∂ η ext and ∂ η int , which give a decomposition u = u 0 + u 1 with u 0 ∈ Har η int and u 1 ∈ Har * η ext .Consequently, we may express u as the sum of two series u (20)0 .It is worth noting that parallel to(20)we