ABSTRACT (original article): 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. CITATION (original article): Ashtab, Z., Morais, J. & Porter, R.M. Fourier method for the Neumann problem on a torus. Anal.Math.Phys. 13, 78 (2023). https://doi.org/10.1007/s13324-023-00836-4
Explanation of principal concepts and notation
Explanation of principal concepts and notation
Fourier and Toroidal series and modes
Fourier and Toroidal series and modes
Toroidal coordinates in 3D space are denoted (η,θ,ϕ). The inside of the torus (or solenoid) consists of values for which η<η<∞, where the value η=∞ corresponds to a circle lying in the interior of the torus, going once around the central "hole" . Toroidal coordinates are used because there are convenient formulas for a complete set of harmonic (or potential) functions expressible in these coordinates. The toroidal harmonics are functions of three variables (η,θ,ϕ) and are indexed by two nonnegative integers (n,m) as well as one of four families which we label (Cos,Cos), (Cos,Sin), (Sin, Cos), (Sin,Sin). These families can be written as (ν,μ) where ν, μ are either Cos or Sin. In the programs, for clarity we write trigθ and trigϕ in place of ν and μ.This packages handles Fourier series (real series, in cosines and sines), as well as toroidal series (in toroidal harmonics). Such series are determined by the coefficients of the corresponding terms and thus require an appropriate Mathematica structure to hold the data. Some routines in this package handle lists of coefficients, and others use "labeled" lists, which we now explain.To understand the idea behind labeled lists, we illustrate first for simplicity with a Fourier series of a single angle (although this package always uses series of two angles). A function such asmight be represented by the list of coefficients.Note the zeros in this list used as placeholders for the terms which do not appear in the series. A labeled list might look like ,which might be a more efficient format if there are lots of zero coefficients. However, with this data, one might not know that it is the Cos function which is involved, so perhaps a notation such as {Cos, {{2, 3.1}, {3, 4.2}, {5,5.3}}} or {{{2, Cos}, 3.1}, {{3, Cos} 4.2}, {{5,Cos}, 5.3}} would be useful.In this package, functions of two angles θ, ϕ often require series in θ for trig functions of fixed multiples of ϕ. Therefore the labelling convention is adopted that {{Cos,Cos},m} is the label applicable to a series such as where in this example, all the trig functions of ϕ have the same factor m (but different factors for θ), and also involve the same two functions ν=Cos, μ=Cos.Thus for example, 3.1 Cos[2 θ]Cos[3 ϕ] + 4.2 Cos[3 θ]Cos[3 ϕ] - 5.3 Cos[5 θ]Cos[3 ϕ] will be represented by{ { {Cos,Cos},3}, {0, 0, 3.1, 4.2, 0, -5.3} }.Note that the list of θ-coefficients is not labeled, except for the information that it is a cosine series. Similarly, one has {{{Sin,Cos},3},{0, 0, 3.1, 4.2, 0, -5.3}} as the representation of 3.1 Sin[2 θ]Cos[3 ϕ] + 4.2 Sin[3 θ]Cos[3 ϕ] - 5.3 Sin[5 θ]Cos[3 ϕ]. The first trig function in a pair such as {Sin,Cos} is applied to θ, and the second is applied to ϕ. We will call labels of the form {{trigθ,trigϕ},m} “modes".The reason for labeling the coefficients this way is that the Dirichlet-to-Neumann mapping for a general function is the sum of the Dirichlet-to-Neumann mappings of the parts of the function corresponding to distinct modes.The series we have mentioned are properly called Fourier series. This package mainly works with “toroidal series”, in which the trigonometric functions such as Sin[2θ]Cos[3ϕ] are replaced by toroidal harmonics, sometimes denoted . The main difference is the inclusion of some extra factors: = (coshη)(cosh)For series in toroidal harmonics, the same conventions of labeling coefficients with modes are applied. The formula is called a solid toroidal harmonic because it is a function of three variables (θ,ϕ,η). When is fixed, the formula is a boundary toroidal harmonic (simply called a toroidal harmonic), and differs from a Fourier term by the inclusion of the factor -cosθ.
0
3.1cos2θ+4.2cos3θ-5.3cos5θ
{0,0,3.1,4.2,0,-5.3}
{{2,3.1},{3,4.2},{5,5.3}}
3.1cos2θcosmϕ+4.2cos3θcosmϕ-5.3cos5θcosmϕ
ν,μ
U
n,m
ν,μ
U
n,m
coshη-cosθ
cos(nθ)cos(mϕ)n
Q
m
n
Q
m
η
0
ν,μ
U
n,m
η=
η
0
coshη
0
Dirichlet-to-Neumann mapping and Neumann problem
Dirichlet-to-Neumann mapping and Neumann problem
Let f be a function defined on the boundary of the torus. Let u be the unique harmonic function in the interior of the torus such that the boundary values of u are equal to f. The harmonic function u is called the solution to the Dirichlet problem with boundary values f, and f is called the Dirichlet data of the problem. The function u has a normal derivative at each boundary point of the torus; let h denote the values of these normal derivatives. The operator which sends f to h is called the Dirichlet-to-Neumann mapping. The inverse problem, in which the Neumann data h is given and we want to reconstruct f, is called the Neumann problem. In order for a solution of the Neumann problem to exist, h must satisfy a “compatibility condition” which says that the integral of h over the boundary is zero. When h satisfies the compatibility condition, all solutions of the Neumann problem differ from one another by a constant.
General considerations
General considerations
The symbols have fixed meanings in this package, and are Protected so that the user cannot assign values to them.
x0,x1,x2,θ,ϕ,η
The toroidal harmonics are defined in terms of the Legendre Q functions. The Legendre Q functions tend to zero exponentially, so for large indices they may cause numerical underflow and you may need higher precision arithmetic.
Disclaimer:
While we hope that the formulas are correct, we cannot guarantee their accuracy, so the functions in this packages should not be used for critical work.
References and credits
References and credits
This package contains the results of the paper
[1] “Fourier method for the Neumann problem on a torus” by Z. Ashtab, J. Morais, R. M. Porter, in Analysis and Mathematical Physics 13:78 (2023) DOI: 10.1007/s13324-023-00836-4, research carried out at the Cinvestav and I.T.A.M. (Mexico).
Formulas for toroidal harmonics are from the book
Formulas for toroidal harmonics are from the book
[2] , E. W. Hobson, The Theory of Spherical and Ellipsoidal Harmonics, Cambridge University Press, Cambridge (1931).
This package corrects some errors in [1] (see comments in NP.m).
How to use the package
How to use the package
In[]:=
SetDirectory[NotebookDirectory[]]
In[]:=
<<NP.m
(Type ?NeumannProblemPackage for list of available routines)
After loading NP.m, in order to use most of the functions you must first choose a value of and run InitializeToroidalCoefficients[η0,{nmax,mmax}]. This sets up the necessary constants for the torus defined by η > η. Most functions will fail if this is not done previously.
η
0
0
InitializeToroidalCoefficients[2.5,{20,2}](*useanyvaluesyoulike*)(*TorusandAreaareMathematicafunctions,notfromthispackage:*)tor=Torus[{0,0,0},{Tanh[$η0/2],Coth[$η0/2]}];{Graphics3D[tor,Axes->True],ParametricPlot3D[{x1,x2,x0}/.ToroidalCoordinateChangeRules/.η->$η0,{θ,0,2π},{ϕ,0,2π}]}{Area[tor],TorusArea[20]}
Out[]=
,
Out[]=
{6.61367,6.61367}
If you want to do operations on tori of different shapes, you must run InitializeToroidalCoefficients for each one :
In[]:=
Table[(InitializeToroidalCoefficients[ηval,{20,2}];TorusArea[20]),{ηval,.6,.8,.1}]
Out[]=
{115.463,86.1103,66.9426}
Examples
Examples
Example of an expansion of a boundary function as a Fourier series and as a toroidal series
Example of an expansion of a boundary function as a Fourier series and as a toroidal series
Two illustrate the similarities between Fourier and toroidal series, we consider the function θ + ϕ, which we represent as a Fourier series and then graph it in the (θ , ϕ) plane. You may try using more or fewer coefficients:
In[]:=
InitializeToroidalCoefficients[.9,{20,20}]f=θ+ϕ;coef=Fourier2DCoefficientsLabeled[f,5][10];fs=Fourier2DSeriesFromCoefficients[coef];Plot3D[{fs,f},{θ,0,2Pi},{ϕ,0,2Pi},PlotStyle->{Red,Green},PlotRange->{-1,13}]
Out[]=
In[]:=
Plot3D[{fs,f},{θ,0,2Pi},{ϕ,0,2Pi},PlotStyle->{Red,Green},PlotRange->{-1,10}]
We consider the same function θ + ϕ, which we now represent as a toroidal series and then graph in the (θ , ϕ) plane. You may try using more or fewer coefficients. The function ToroidalNormalizedBoundarySeries combines the two operations of calculating the coefficients and then forming the series with those coefficients:
In[]:=
InitializeToroidalCoefficients[.9,{20,20}]f=θ+ϕ;fs=ToroidalNormalizedBoundarySeries[f,5][10];Plot3D[{fs,f},{θ,0,2Pi},{ϕ,0,2Pi},PlotStyle->{Red,Green},PlotRange->{-1,13}]
Out[]=
Example of D-to-N mapping and its inverse the Neumann problem for a single non-principal mode
Example of D-to-N mapping and its inverse the Neumann problem for a single non-principal mode
Set up for working with an arbitrary value of $η0:
In[]:=
InitializeToroidalCoefficients[.9,{40,40}]
Define a sample harmonic function u:
Choose an arbitrary mode to work with:
We verify that u is harmonic:
Take the boundary values of u (this is called "f" in the article [1])
To get a view of the territory, it is suggested to calculate just a few coefficients:
By inspection we see that this function belongs to a single mode (this also could have been deduced from the formula), so we calculate coefficients for just that mode . You may try more or fewer than 20 coefficients :
Here is the toroidal boundary series corresponding to these coefficients:
We verify that this series is a good approximation of the original boundary function. You can test with any vales of θ,ϕ:
Calculate a formula for the Dirichlet-to-Neumann mapping h of the boundary values. ;(This can only be done because in this example we have an explicit formula for the original harmonic function, and it would not be possible if the given boundary values were, say, from experimental data:
Here we calculate the D-to-N mapping of the boundary values of u, via the coefficients that we have found:
The same could be obtained with a slightly different syntax, by treating it as a single mode:
Here is the toroidal series corresponding to the Dirichlet-to-Neumann image:
(* It takes a simpler form when factored: *)
We check that the series is a good approximation of the theoretical value, assuming enough terms were used. You may try any values of θ,ϕ:
We calculate the inverse operation, that is, we solve the Neumann problem to get the boundary values of the function with given normal derivative h:
The toroidal series corresponding to these coefficients:
We check how good a solution it is. First we compare the coefficients, then the values of the functions:
Finally, the harmonic function in the interior of the torus is obtained by using the same coefficients as for the boundary. We compare this to the function u original. You may use any values you wish (as long as η>$η0):
Example of D-to-N mapping and its inverse the Neumann problem involving many modes:
Example of D-to-N mapping and its inverse the Neumann problem involving many modes:
This example comes from the harmonic function Re[Exp[x1 + I x2]], converted to polar coordinates. The function u is defined, we verify that it is harmonic, initialize the toroidal functions, and define f to be the boundary values of u:
We look at some coefficients, then choose a list of modes which do not give zero coefficients, and calculate a larger number of toroidal coefficients for only those modes:
We form the toroidal series for these boundary coefficients, and compare its values with those of the original boundary function at arbitrarily chosen values of θ and ϕ:
We calculate the Dirichlet - to - Neumann operator applied to the boundary value function f, giving a new function h on the boundary in terms of a series. Also, for this particular example, it is possible to calculate a formula for the normal derivative, which we use to check the values obtained through the series:
Observe that the principal mode is the first one in the list ubdrycoefs and in hcoefsviad2n . Every normal derivative must satisfy the compatibility condition in the principal mode :
Next we solve the NeumannProblem using as input the normal derivative coefficients, and then construct the corresponding series:
There is a snag: the function contains the principal mode {{Cos,Cos},0}, for which the Neumann problem does not have a unique solution. We have found a solution which differs from the original boundary function by a constant, which is the worst that can possibly happen:
The principal mode is the first one in the list ubdrycoefs, and the proper value of "a0" is its first coefficient. With this value we recover the original boundary function:
Finally, the harmonic function in the interior of the torus is obtained by using the same coefficients as for the boundary :
Example of D-to-N and Neumann from the harmonic function Im[Exp[x1 + I x2]]
Example of D-to-N and Neumann from the harmonic function Im[Exp[x1 + I x2]]
Define the harmonic function u, and let f be its boundary restriction:
Inspect a few toroidal coefficients of the boundary function:
We observe that only {Cos,Sin} modes are nonzero, so we calculate coefficients for those modes:
Construct the toroidal series with these coefficients, and compare values with the original function ubdry:
Calculate the normal derivative of the original solid function u, compare with the solution of the D-to-N problem:
Finally, use the above data to show how well the routine NeumannProblem produces the correct result:
Example of D-to-N operator, with less detail
Example of D-to-N operator, with less detail
Example of D-to-N operator, with less detail:
Example of D-to-N operator, with less detail:
Example of several D - to - N calculations on the same torus
Example of several D - to - N calculations on the same torus
It is not necessary to repeat InitializeToroidalCoefficients for solving different boundary problems on the same torus:
Calculation of Example 1 in reference [1]
Calculation of Example 1 in reference [1]
You may try variations, such as InitializeToroidalCoefficients[SetPrecision[1.5,200],{120,5}]
This example shows a list of coefficients which does not converge to zero quite rapidly enough for the general theory to guarantee results. The boundary function f (with toroidal coefficients fcoef) is a smooth function, but its D-to-N image h (with toroidal coefficients hcoef) has discontinuities. This type of phenomenon is typical in the theory of Fourier series.
Calculation of some of the values in the first two columns of reference [1, Table 1]
Calculation of some of the values in the first two columns of reference [1, Table 1]
The function NeumannCnDn produces two sequences (C0,C1,....), (D0,D1,....) by means of a recursion formula. For large n, the value Dn/Cn is a key to solving the Neumann problem in the nonprincipal modes. Some values of these sequences were studied to illustrate the behaviour of these seequences.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Fourier method for the Neumann problem on a torus
by Z. Ashtab, J. Morais & R. Michael Porter
Wolfram Community, STAFF PICKS, December 9, 2025
https://community.wolfram.com/groups/-/m/t/3588303
by Z. Ashtab, J. Morais & R. Michael Porter
Wolfram Community, STAFF PICKS, December 9, 2025
https://community.wolfram.com/groups/-/m/t/3588303