Computing with Harmonic Functions

Sheldon Axler

24 December 2008		© Sheldon Axler

Getting Started

Harmonic Function Theory

This notebook is the manual for HFT6.m (version 6.22), a Mathematica package for computing with harmonic functions. The HFT6.m package allows the user to make calculations that would take a prohibitive amount of time if done without a computer. For example, the Poisson integral of any polynomial can be computed exactly.

Earlier versions of this software focused on algorithms arising from the material in the book Harmonic Function Theory [ABR] by Sheldon Axler, Paul Bourdon, and Wade Ramey. That book is still the source for many of the algorithms used in the HFT6.m package, but the goal of the package has expanded to include additional symbolic manipulations involving harmonic functions. The HFT6.m package can be used without the Harmonic Function Theory book, but the debt to the book is reflected in the initials chosen for the name of the package.

This documenation and the HFT6.m package are intended for use with Mathematica version 6.

This notebook, the HFT6.m package, and the HFT6.nb notebook that serve as the source file for the HFT6.m package, are available electronically free of charge. To obtain the most recent versions, retrieve them from http://www.axler.net.

Comments, suggestions, and bug reports about this notebook, the HFT6.m package, or the book Harmonic Function Theory are welcome; please send them via e-mail to axler@sfsu.edu.

Loading the Package

To begin a session with the HFT6.m package, first start Mathematica. Then enter the command shown below, making an appropriate modification to indicate the directory or folder in which HFT6.m is located on your computer. After typying <<, you can use the Get File Path command, located in the Input menu, to paste the full pathname of HFT6.m into your input cell:

<<"C:\\MyStuff\\math\\publications\\HFT-Mathematica\\HFT6.m"

After loading the HFT6.m package, you can use the functions described in the rest of this notebook.

This is a live Mathematica notebook, meaning that you can and should evaluate and modify the input cells in this notebook.

On-line Help

This notebook provides the only complete documentation for using the HFT6.m package. However, users can be reminded about the basic use and syntax of each function created by HFT6.m by entering ? followed by the function name:

?Volume

The Mathematica statement Names["HFT`*"] will produce a complete list of the functions created by HFT6.m:

Names["HFT`*"]

Calculus in (R^n)

Norm

Norm[x] is the Euclidean norm of a vector x. When using StandardForm for output (which is the default in Mathematica), this package will use ‖x‖ instead of Norm[x]. When using StandardForm for input (which is the default in Mathematica), you can use either ‖x‖ or Norm[x]:

\[LeftDoubleBracketingBar]5 x\[RightDoubleBracketingBar] -((Norm[2x])^2)

ExpandNorm

ExpandNorm expands norms of sums (and differences) by replacing expressions of the form ‖x+y‖ with (((‖x‖)^2)+((‖y‖)^2)+ ((x.y))^(1/2)). ExpandNorm is often useful for verifying identities. For example, to prove the Symmetry Lemma (1.13 in [ABR]), which states that

‖(y/(‖y‖))-‖y‖ x‖=‖(x/(‖x‖))-‖x‖ y‖

for all nonzero x, y ϵ (R^n), we subtract the left-hand side of the alleged equality from the right-hand side, hoping that the result equals 0:

\[LeftDoubleBracketingBar](x/(\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]))-\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar] y\[RightDoubleBracketingBar]-\[LeftDoubleBracketingBar](y/(\[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]))-\[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar] x\[RightDoubleBracketingBar]

The last output does not look like 0, so we apply ExpandNorm to it. The symbol % used below is the Mathematica abbreviation for the last output:

ExpandNorm[%]

Because the last output equals 0, the Symmetry Lemma has been proved.

Partial

Partial[f, (Subscript[x, j])] is the partial derivative of f with respect to (Subscript[x, j]). Here x denotes a vector in Euclidean space, and (Subscript[x, j]) is the (j^th) coordinate of x. The advantage of using Partial rather than the built-in Mathematica function D or ∂ is that Partial can deal with norms and dot products.  Note that the dot product of b and x is  denoted by b.x:

Partial[((\[LeftDoubleBracketingBar]a+x\[RightDoubleBracketingBar])^(b.x)),(Subscript[x, 5])]

A subscript can be typed in Mathematica by pressing control-dash. Return to normal mode by pressing control-space. A suberscript or exponent can be typed by pressing control-6. Instead of entering (Subscript[x, j]), you can type x[j].

For taking multiple partial derivatives, Partial uses the same syntax as the built-in Mathematica function D. For example, the following command computes the partial derivative of |x| once with respect to the second coordinate, twice with respect to the first coordinate, and three times with respect to the fourth coordinate:

Partial[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar],(Subscript[x, 2]),{(Subscript[x, 1]),2},{(Subscript[x, 4]),3}]

Laplacian

Laplacian[f, x] is the Laplacian of f with respect to x:

Laplacian[(Subscript[x, 3]) \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar],x]

Note that in the last result, the Laplacian depends upon Dimension[x], which is the dimension of the Euclidean space in which x lives.  For work in (R^n), SetDimension[x, n] should be used to set the dimension of x equal to n; here n can be a symbol or a specific integer, say 8. We illustrate this procedure by showing that for each z with |z| = 1, the Poisson kernel P( . , z) is harmonic (in other words, what follows is a proof of Proposition 1.18 of [ABR]):

SetDimension[x,n]

Laplacian[(1-((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^2))/((\[LeftDoubleBracketingBar]x-z\[RightDoubleBracketingBar])^n),x]

ExpandNorm[%] /. Norm[z]->1

Thus the Laplacian of P( . , z) is 0, and hence P( . , z) is harmonic, as claimed.

(Laplacian^m) can be used, with the same syntax as Laplacian, to evaluate the (m^th) power of the Laplacian. For example, here we find the biLaplacian (the Laplacian of the Laplacian) of 1/|x| on (R^n):

(Laplacian^2)[1/\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar],x]

If x and y are symbols,then Laplacian[f,x,y] gives the Laplacian of f with respect to the vector (x, y), where x is thought of as a vector, y is thought of as a real variable, and (x, y) is thought of as a vector in a Euclidean space whose dimension is one more than the dimension of x. This format is often useful when working with functions defined on the upper half-space, as in Chapter 7 of [ABR]. Note that if we are thinking of (x, y) as an element of (R^n), then x should be defined to have dimension n-1. We illustrate the use of this format of the Laplacian function by showing that for each t ϵ (R^(n-1)), the Poisson kernel (Subscript[P, H])( . ,  t) for the upper half-space is harmonic (see page 145 of [ABR] for another proof of this):

SetDimension[x,n-1]

Laplacian[y/(((((\[LeftDoubleBracketingBar]x-t\[RightDoubleBracketingBar])^2)+(y^2)))^(n/2)),x,y]

ExpandNorm[%]

Thus the Laplacian of (Subscript[P, H])( . ,  t) is 0, and hence (Subscript[P, H])( . ,  t) is harmonic, as claimed.

Grad

Grad[f, x] is the gradient of f with respect to the vector x (Mathematica's convention of spelling out names in full would imply that this function should be called Gradient instead of Grad, but Gradient is already used by Mathematica for an unrelated purpose):

Grad[a.x+(Subscript[x, k])+\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar],x]

In the last output, Delta[k] denotes the vector that equals 1 in the (k^th) coordinate and 0 in the other coordinates.

NormalD

NormalD[f, z] equals the outward normal derivative of f at a point z on the unit sphere.

In the section of Chapter 1 of [ABR] titled The Poisson Kernel for the Ball, the formula for the Poisson kernel was derived by showing that P(x, z) equals the normal derivative (with respect to z) of

(((‖z-x‖)^(2-n))-((‖x‖)^(2-n)) ((‖z-(x/((‖x‖)^2))‖)^(2-n)))/(2-n).

The computation of that normal derivative is a bit complicated and was left to the reader. Now we show how this computation can be done easily by a computer. We begin by taking the normal derivative of the function above:

NormalD[(((\[LeftDoubleBracketingBar]z-x\[RightDoubleBracketingBar])^(2-n))-((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^(2-n)) ((\[LeftDoubleBracketingBar]z-(x/((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^2))\[RightDoubleBracketingBar])^(2-n)))/(2-n),z]

Now we ask the computer to replace ‖z-(x/((‖x‖)^2))‖ in the last output with ‖z-x‖/‖x‖; the Symmetry Lemma (1.13 of [ABR]) implies that  the two expressions are equal:

% /.\[LeftDoubleBracketingBar]z-(x/((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^2))\[RightDoubleBracketingBar]->\[LeftDoubleBracketingBar]z-x\[RightDoubleBracketingBar]/\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]

Applying the Mathematica function PowerExpand to the last output gives the desired result:

PowerExpand[%]

The last output completes the derivation of the formula for the Poisson kernel for the ball.

Matrices

Matrix dimensions

Mathematica uses . to represent matrix multiplication as well as dot product.  When using this package, you must use the SetDimension command to tell Mathematica whenever a symbol is being thought of as a matrix. To think of A as an m-by-n matrix, enter the command SetDimension[A, {m, n}]. For the examples in this section, we want to think of A as a square matrix, so we make the two dimensions of A equal:

SetDimension[A, {n,n}]

Grad[\[LeftDoubleBracketingBar]A.x\[RightDoubleBracketingBar],x]

To interpret input or output involving products of matrices and vectors, think of vectors as either column vectors or row vectors, whichever makes sense in context. Thus in the input above, x is obviously a vector, because we are taking the gradient of some function with respect to x. For A.x to make sense as a vector in (R^n), we must think of x as a column vector (an n-by-1 matrix).  To interpret (A.x).A in the last output, we again think of x as a column vector (an n-by-1 matrix), so that A.x makes sense as an element of (R^n).  Because A.x is then multiplied by the n-by-n matrix A, we must think of A.x as a row vector (an 1-by-n matrix), so that (A.x).A is a vector, as expected here (because it is the gradient of a function).

In the example above we did not enter the command SetDimension[x, n] to tell Mathematica that x is a vector (although it would not have hurt to do so). The HFT6.m package can almost always use the context to distinguish scalars from vectors (if both interpretations make sense, the package assumes that a symbol represents a scalar).  However, you must always explicitly use the SetDimension command to tell Mathematica which objects are matrices.

Rows and columns

If A is a matrix, then (Subscript[A, j,.]) denotes the (j^th) row of A and (Subscript[A, .,j]) denotes the (j^th) column of A:

Partial[\[LeftDoubleBracketingBar]A.x\[RightDoubleBracketingBar]+5\[LeftDoubleBracketingBar]x.A\[RightDoubleBracketingBar],(Subscript[x, 3])]

Note that in the input above, x is thought of as a column vector in the term A.x and as a row vector in the term x.A.  To interpret (A.x).(Subscript[A, .,3]) in the last output, we think of x as a column vector, so that A.x is a vector.  Then (A.x).(Subscript[A, .,3]) is just the dot product of two vectors, so that it is a number, as expected here.

HilbertSchmidt

HilbertSchmidt[A] is the Hilbert-Schmidt norm of a matrix A, which is the square root of the sum of the squares of all entries of A:

Laplacian[\[LeftDoubleBracketingBar]x.A\[RightDoubleBracketingBar],x]

Here Transpose[A] is the matrix obtained from A by interchanging the rows and columns of A.

Divergence

Divergence[f, x] is the divergence of f with respect to the vector x:

SetDimension[A, {n,n}]

Divergence[A.x/\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar],x]

Here Tr[A] is the trace of the square matrix A.

J

J[f, x] is the Jacobian derivative of f with respect to x. Here f should be a function of x taking values in some Euclidean space, and the Jacobian derivative is the usual matrix consisting of partial derivatives of the coordinates of f (Mathematica's convention of spelling out names in full would imply that this function should be called Jacobian instead of J, but Jacobian is already used by Mathematica for an unrelated purpose):

SetDimension[x,n]

J[x/((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^2),x]

Although x can be thought of as a row vector or a column vector, depending on the context, if we were using explicit coordinates we would represent x in Mathematica by something like {(Subscript[x, 1]),(Subscript[x, 2]),... , (Subscript[x, n])}.  Thus {x} would be {{(Subscript[x, 1]),(Subscript[x, 2]),... , (Subscript[x, n])}}, which is a 1-by-n matrix.  Thus Transpose[{x}] is an n-by-1 matrix, and the product Transpose[{x}].{x}, which appears in the last output, is an n-by-n matrix, as expected.

IdentityMatrix[n], which also appears in the last output, denotes the n-by-n matrix with 1's on the diagonal and 0's elsewhere.

Homogeneous

Homogeneous[u, m, x] is the term of degree m in the homogeneous expansion of u at the origin; here u is thought of as a function of the vector x:

SetDimension[x,3]

Homogeneous[(1-((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^2))/((Cos[(Subscript[x, 1])])^(3/2)),4,x]

To find homogeneous expansions about points other than the origin, use Homogeneous[u, m, x, b], which gives the term of degree m in the homogeneous expansion of u, as a function of x, about b:

Homogeneous[(Subscript[x, 1]) (Subscript[x, 2])+((Subscript[x, 3])^4),2,x,b]

The first time in each Mathematica session that you use Homogeneous to find an expansion about a point other than the origin, you may see a message that Togetherness has been turned off.  The Togetherness subsection later in this notebook provides an explanation for this message.

Taylor

Taylor[u, m, x] is the sum of all terms of degree at most m in the Taylor series expansion of u at the origin; here u is thought of as a function of the vector x:

Taylor[((Subscript[x, 1])^2)/(6+(Subscript[x, 2])+Cos[(Subscript[x, 3])]),4,x]

To find Taylor series expansions about points other than the origin, use Taylor[u, m, x, b], which gives the sum of all terms of degree at most m in the Taylor series expansion of u, as a function of x, about b:

Taylor[1+(Subscript[x, 1]) (Subscript[x, 2])+((Subscript[x, 1])^2),2,x,b]

The first time in each Mathematica session that you use Taylor to find an expansion about a point other than the origin, you may see a message that Togetherness has been turned off.  The Togetherness subsection later in this notebook provides an explanation for this message. For now, we turn Togetherness back on so that further output in this notebook will be displayed nicely.

TurnOn[Togetherness]

Nonexplicit functions

Our examples so far have involved only concretely defined functions.  But functions in symbolic form can also be used with each of the differentiation commands in this package. Here, for example, is how to find the gradient of the function whose value at x is f((((g(3x)))^2)):

Grad[f[((g[3x])^2)],x]

Here Grad[g][3x] denotes the gradient of g, evaluated at 3 x. For the example above to make sense, we (and the computer) must think of g as a real-valued function on some Euclidean space (where x lives) and f as a function from R to R.  Let's find the Laplacian of the same function:

Laplacian[f[((g[3x])^2)],x]

Laplacian[g][3x] denotes, of course, the Laplacian of g, evaluated at 3 x.

All our differentiation commands (Partial, Laplacian, Gradient, NormalD, Divergence, Jacobian) and Homogeneous and Taylor can be used with non-explicit functions.  For example, if working in (R^2), to find the term of degree 3 in the homogeneous expansion of f((((g(3x)))^2)) about the origin, you would enter the command SetDimension[x, 2] and then the command Homogeneous[f[((g[3x])^2)], 3, x]. Try it! In the answer you will see terms such as (Subscript[Partial, {1,1,2}])[g][0], which denotes the partial derivative of g, twice with respect to the first variable and once with respect to the second variable, evaluated at 0.

Vector-valued functions

Suppose we compute the partial derivative with respect to (Subscript[x, 1]) of the function that takes x to f(g(((h(x))^2))):

Partial[f[g[((h[x])^2)]],(Subscript[x, 1])]

The example above illustrates a general principle: This package assumes that all functions are real valued, unless told otherwise.  Thus the last output is correct if we are thinking of g as a real-valued function. If we want to think of g as taking values in (R^n), we must first enter the command SetDimension[g[_], n], which instructs the computer that g is a function with range in (R^n):

SetDimension[g[_], n]

Partial[f[g[((h[x])^2)]],(Subscript[x, 1])]

Note that the last output  involves the dot product of two vectors: the gradient of f, evaluated at g(((h(x))^2)), and the derivative of 
g, evaluated at ((h(x))^2). The last two outputs should be compared.

Volume

Volume[n] is the volume (unnormalized) of the unit ball in (R^n):

Volume[4]

Volume is computed using the formula given by Exercise 6 in Appendix A of [ABR].

SurfaceArea

SurfaceArea[n] is the surface area (unnormalized) of the unit sphere in (R^n):

SurfaceArea[57]

SurfaceArea[n] is computed by multiplying the formula for Volume[n] by n; see A.2 in Appendix A of [ABR].

IntegrateSphere

IntegrateSphere[f, x] equals the integral of f, with respect to normalized surface area measure, over the unit sphere in the Euclidean space defined by x. Here f should be a polynomial function of x:

SetDimension[x,n]

IntegrateSphere[((Subscript[x, 1])^2) ((Subscript[x, 2])^4) ((Subscript[x, 3])^6),x]

IntegrateSphere is computed by using the results in Section 3 of [W].

IntegrateBall

IntegrateBall[f, x] equals the integral of f, with respect to (unnormalized) volume measure, over the unit ball in the Euclidean space defined by x. Here f should be a function of x and ||x|| that is a polynomial in x and a function in ||x|| for which Mathematica can find an explicit antiderivative:

SetDimension[x,7]

IntegrateBall[((Subscript[x, 1])^2) ((Subscript[x, 2])^4)/(1+\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]),x]

IntegrateBall is computed by converting to polar coordinates and then using the function IntegrateSphere (see [R], Chapter 8, Exercise 6).

Boundary Value Problems

Dirichlet problems

The Dirichlet problem on the sphere

Dirichlet[p, x] is the solution to the standard Dirichlet problem: find the harmonic function on the unit ball in the Euclidean space defined by x that equals p on the unit sphere.  Thus Dirichlet[p, x] is the Poisson integral of p as a function of x.  Here p must be a polynomial function of x:

SetDimension[x,5]

Dirichlet[((Subscript[x, 1])^4) ((Subscript[x, 2])^2),x]

The last output is a harmonic function on (R^5)that equals ((Subscript[x, 1])^4)((Subscript[x, 2])^2) on the unit sphere. To check this, first we take the Laplacian of the last output, and then we evaluate the last output on the unit sphere:

Laplacian[%, x]

%% /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1

Thus we have indeed solved the specified Dirichlet problem.

The solution to the Dirichlet problem on the sphere is computed by using the algorithm described in [AR].

The Dirichlet problem on the exterior of the sphere

The Region option allows the user to solve Dirichlet problems on regions other than the sphere. Currently supported values for Region are Sphere, ExteriorSphere, Annulus, and Quadratic. The default value of Region is Sphere; thus Dirichlet[p, x, Region->Sphere] is the same as Dirichlet[p, x].

Dirichlet[p,x, Region->ExteriorSphere] gives the solution to the standard exterior Dirichlet problem: find the harmonic function on the exterior (including ∞) of the unit ball in the Euclidean space defined by x that equals p on the unit sphere. This function is the exterior Poisson integral of p as a function of x (see Theorem 4.11 of [ABR]).  Here p must be a polynomial function of x:

Dirichlet[((Subscript[x, 1])^4) ((Subscript[x, 2])^2),x, Region->ExteriorSphere]

The last output is a harmonic function on the exterior of the unit ball in (R^5) (because the dimension of x was set to 5 earlier) that equals ((Subscript[x, 1])^4)((Subscript[x, 2])^2) on the unit sphere. This function is the Kelvin transform of the solution that would have been obtained without the Region->Sphere option; see Chapter 4 of [ABR].

The Dirichlet problem on annular regions

Dirichlet[{p,q},x, Region->Annulus[r,s] ] gives the solution to the annular Dirichlet problem: find the harmonic function on the annular region with inner radius r and outer radius s that equals p on the sphere of radius r and equals q on the sphere of radius s. Here p and q must be polynomial functions of x:

Dirichlet[{((Subscript[x, 1])^3),((Subscript[x, 3])^2)},x, Region->Annulus[1,4]]

Now we check that the last output is harmonic on (R^5) (because the dimension of x was set to 5 earlier), that it equals ((Subscript[x, 1])^3) on the sphere of radius 1, and that it equals ((Subscript[x, 3])^2) on the sphere of radius 4:

Laplacian[%, x]

%% /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1

%%% /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->4

If p=q, then the shortcut Dirichlet[p,x, Region->Annulus[r,s] ] can be used instead of Dirichlet[{p,q},x, Region->Annulus[r,s] ].

The solution to the Dirichlet problem on an annulus is computed by using the techniques described in Chapter 10 of [ABR].

The Dirichlet problem on quadratic surfaces

Dirichlet[p,x, Region->Quadratic[b, c,d] ] gives the harmonic polynomial that equals p on the quadratic surface

b.(x^2)+ c.x + d = 0

Here b and c must be lists whose length equals the dimension of x, and p must be a polynomial function of x. For example, here we find the harmonic polynomial that equals ((Subscript[x, 1])^3)((Subscript[x, 3])^2)on the elliptic paraboloid given by the equation

2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2) = 4(Subscript[x, 3])+ 1

in (R^3):

SetDimension[x,3]

Dirichlet[((Subscript[x, 1])^3)((Subscript[x, 3])^2),x,Region->Quadratic[{2,3,0},{0,0,-4},-1]]

The last output is of the form ((Subscript[x, 1])^3)((Subscript[x, 3])^2)+(2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2)-4(Subscript[x, 3])-1)f, where f is a polynomial. Thus the last output obviously equals ((Subscript[x, 1])^3)((Subscript[x, 3])^2) on the elliptic parabolic given by the equation 2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2) = 4(Subscript[x, 3])+ 1. Thus to verify that the last output is correct, we only need check that it is harmonic:

Laplacian[%, x]

The solution to the Dirichlet problem on a quadratic surface is computed by using the algorithm developed in [AGV]. That paper also shows that if each coordinate of b is nonnegative and at least one coordinate of b is nonzero, and if the quadratic surface in question is nontrivial in the sense that there exits x ϵ (R^n) such that

b.(x^2)+ c.x + d < 0,

then the polynomial produced by this algorithm is the unique harmonic polynomial that equals p on the quadratic surface in question.

Using the notation above, if c consists of all 0's, then it may be omitted. For example, here we find the harmonic polynomial that equals ((Subscript[x, 1])^4)((Subscript[x, 3])^2)on the ellipsoid given by the equation 2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2)+4((Subscript[x, 3])^2)=1:

Dirichlet[((Subscript[x, 1])^4)((Subscript[x, 3])^2),x,Region->Quadratic[{2,3,4},-1]]

Generalized Dirichlet problems

The generalized Dirichlet problem asks for a function with specified boundary values on some surface and with specified Laplacian inside that surface. Thus if the Laplacian is required to be 0, then this is the same as the usual Dirichlet problem. The generalized Dirichlet problem with specified Laplacian equal to p is solved by adding the option Laplacian-> p to the Dirichlet function; here p must be a polynomial. For example, we can find the function that equals ((Subscript[x, 1])^3)((Subscript[x, 2])^2)on the unit sphere in (R^3) and whose Laplacian equals ((Subscript[x, 2])^2)(Subscript[x, 3]):

Dirichlet[((Subscript[x, 1])^3)((Subscript[x, 2])^2), x, Laplacian-> ((Subscript[x, 2])^2)(Subscript[x, 3])]

We can check that output above correctly solves the generalized Dirichlet problem by verifying that the Laplacian of the last output indeed equals ((Subscript[x, 2])^2)(Subscript[x, 3]) and that the last output restricted to the unit sphere indeed equals ((Subscript[x, 1])^3)((Subscript[x, 2])^2):

Laplacian[%, x]

%% /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1

The Laplacian option can be used combined with the Region option to solve the generalized Dirichlet problem on any region for which this package can solve the usual Dirichlet problem. The order of the options does not matter. Here we find the function on (R^3) whose Laplacian equals ((Subscript[x, 2])^2)and that equals ((Subscript[x, 1])^4)((Subscript[x, 3])^2)on the ellipsoid 2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2)+4((Subscript[x, 3])^2)- 1 = 0:

Dirichlet[((Subscript[x, 1])^4)((Subscript[x, 3])^2),x,Laplacian->((Subscript[x, 2])^2),Region->Quadratic[{2,3,4},-1]]

The form of the last output shows, without calculation, that it obviously equals ((Subscript[x, 1])^4)((Subscript[x, 3])^2)on the ellipsoid 2((Subscript[x, 1])^2)+3((Subscript[x, 2])^2)+4((Subscript[x, 3])^2)- 1 = 0. Thus to verify that we have correctly solved this generalized Dirichlet problem, we only need to check that the Laplacian equals ((Subscript[x, 2])^2):

Laplacian[%, x]

Dirichlet problems with explicit coordinates

So far we have worked with symbolic vectors, which we have usually called  x or z. Vectors can also be described by giving explicit coordinates in the form of a list. Both formats (symbolic vectors and explicit lists) work with all functions in this package. For example, suppose we are working in (R^3) and want to use x, y, and z instead of (Subscript[x, 1]), (Subscript[x, 2]), and (Subscript[x, 3]). Here is how to find the Poisson integral of (x^3)y(z^2):

Dirichlet[(x^3)y (z^2), {x,y,z}]

To check that the last output is correct, first we verify that it is harmonic by checking that its Laplacian equals 0:

Laplacian[%, {x,y,z}]

Now we check that the difference between the claimed solution to the Dirichlet problem equals (x^3)(yz^2) on the unit sphere by verifying that the difference between these two functions, when restricted to the unit sphere, equals 0:

(x^3)y (z^2)- %% /. z->(Sqrt[1-(x^2)-(y^2)])

Here is another example using explicit coordinates. In this example, we find the function on (R^3) whose Laplacian equals (y^2)z and that equals (x^3)y(z^2) on the ellipsoind 2(x^2)+3(y^2)+4(z^2)- 1 = 0:

Dirichlet[(x^3)y (z^2), {x,y,z}, Region->Quadratic[{2,3,4},-1], Laplacian->(y^2)(z^3)]

Obviously the last output equals (x^3)y(z^2) on the ellipsoind 2(x^2)+3(y^2)+4(z^2)- 1 = 0. We now check that its Laplacian is the desired function:

Laplacian[%, {x,y,z}]

HarmonicDecomposition

HarmonicDecomposition[p, x] gives the decomposition of p into a sum of harmonic polynomials (on the Euclidean space defined by x) times powers of ||x||.  Here p must be a polynomial function of x:

SetDimension[x,n]

HarmonicDecomposition[((Subscript[x, 1])^4),x]

The output given by HarmonicDecomposition[p, x] consists of a list of pairs. The first entry in each pair is a harmonic function, the second entry is the power of ||x|| by which the first entry should be multiplied so that the sum of the resulting terms equals p. For example, in the last output ((3||x((||)^4))/((2+n)(4+n)))-((6||x((||)^2)((Subscript[x, 1])^2))/(4+n))+((Subscript[x, 1])^4) and ((6(-||x((||)^2)+ n ((Subscript[x, 1])^2)))/(n(4+n))) and (3/(n(2+n))) are harmonic functions of x; furthermore

((Subscript[x, 1])^4) = ((3||x((||)^4))/((2+n)(4+n)))-((6||x((||)^2)((Subscript[x, 1])^2))/(4+n))+((Subscript[x, 1])^4) + ((6(-||x((||)^2)+ n ((Subscript[x, 1])^2)))/(n(4+n)))||x((||)^2) + (3/(n(2+n)))||x((||)^4).

For a proof of the existence and uniqueness of the harmonic decomposition, see [ABR], Theorem 5.7. HarmonicDecomposition is computed using the algorithm described in Section 2 of [AR].

AntiLaplacian

AntiLaplacian[u, x] gives an anti-Laplacian of u with respect to x; this is a function whose Laplacian with respect to x equals u.  Here u must  be a polynomial function of x or a sum of terms, each of which is a polynomial in x times a function of ||x||:

SetDimension[x,6]

AntiLaplacian[((Subscript[x, 1])^2) (Subscript[x, 2]) ((\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar])^3) Log[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]],x]

The function given by the last output has Laplacian equal to ((Subscript[x, 1])^2)(Subscript[x, 2])((‖x‖)^3)log(‖x‖)on (R^6), as we can easily verify:

Laplacian[%, x]

To understand how AntiLaplacian is computed, suppose we want to find an anti-Laplacian of p(x)f(||x||), where p is a polynomial on (R^n) and f is a continuous function on (0, ∞).  Applying HarmonicDecomposition to p, using linearity, and replacing f(||x||) with a power of ||x|| times f(||x||), we see that we only need to find anti-Laplacians of functions of the form q(x)f(||x||), where q is a harmonic polynomial on (R^n) homogeneous of degree m.  To find an anti-Laplacian of q(x)f(||x||), suppose h is a twice-differentiable function on (0, ∞).  The Laplacian of q(x)h(||x||) equals

q(x) (((2m + n - 1) h'(||x||) + ||x|| h''(||x||))/(||x||)),

where we have used the product rule for Laplacians (see 1.19 in [ABR]) and Exercise 29, Chapter 1 of [ABR].  We want to make the expression above equal  q(x)f(||x||), so we need only solve the differential equation

(((2m + n - 1) h'(t) + t h''(t))/t) = f(t).

A solution to this equation is given by

h(t) = (∫_a^t)(s^(1-2m-n)) (∫_b^s)(r^(2m+n-1))f(r)drds

for any choice of the constants a and b. AntiLaplacian uses the formula above with b = 0, although if a singularity is encountered at 0 it automatically makes a different choice. If you know there is a singularity at 0, you can save some time by using the option Singularity -> 0, as illustrated by the following example. The two timing tests below show a large speed-up when using the option Singularity -> 0:

Timing[AntiLaplacian[((Subscript[x, 1])^2) (Subscript[x, 2])   ((Log[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]])^10),x]]

Timing[AntiLaplacian[((Subscript[x, 1])^2) (Subscript[x, 2])   ((Log[\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]])^10),x, Singularity-> 0]]

The anti-Laplacian of a given function is never unique.  However, each polynomial has a unique anti-Laplacian that is a polynomial multiple of ((‖x‖)^2), where x is  the variable. The option Multiple -> (Norm^2) produces the unique anti-Laplacian that is a polynomial multiple of ((‖x‖)^2):

AntiLaplacian[((Subscript[x, 1])^2) ((Subscript[x, 2])^5),x,Multiple->(Norm^2)]

Neumann

Neumann[f, x] is the solution to the standard Neumann problem: find the harmonic function on the unit ball in the Euclidean space defined by x whose outward normal derivative on the unit sphere equals f and whose value at the origin equals 0. Here f must be a polynomial function of x.  Furthermore, the integral of f over the unit sphere with respect to surface area measure must equal 0 (Green's identity shows that this condition is necessary for the existence of a solution to the standard Neumann problem):

SetDimension[x,4]

Neumann[((Subscript[x, 1])^6) (Subscript[x, 2]),x]

ExteriorNeumann

ExteriorNeumann[f, x] is the solution to the standard exterior Neumann problem: find the harmonic function on the exterior of the unit ball in theEuclidean space defined by x whose outward normal derivative on the unit sphere equals f and whose limit at ∞ equals 0.   Here f must be a polynomial function of x. If the dimension of the Euclidean space equals 2, then the integral of f over the unit circle with respect to arc length measure must equal 0:

SetDimension[x,5]

ExteriorNeumann[((Subscript[x, 1])^6) (Subscript[x, 2]),x]

BiDirichlet

BiDirichlet[f, x] is the solution to the standard biDirichlet problem: find the biharmonic function on the unit ball in the Euclidean space defined by x that equals f on the unit sphere and whose normal derivative on the unit sphere equals 0.   (A function is called biharmonic if the Laplacian of its Laplacian equals 0.) Here f must be a polynomial function of x:

SetDimension[x,3]

BiDirichlet[((Subscript[x, 1])^4) ((Subscript[x, 2])^3),x]

To check that the last output is correct, first we take its Laplacian squared to check that we have a biharmonic function:

(Laplacian^2)[%, x]

Next we check that the outward normal derivative on the unit sphere indeed equals 0:

NormalD[%%, x]

Finally, we check that we have the correct boundary values:

%%% /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1

Spherical Harmonics

BasisH

BasisH[m, x] gives a vector space basis for (Subscript[H, m])((R^n)), the space of homogeneous harmonic polynomials on (R^n)of degree m, where n is the dimension of x. The user must first set the dimension of x to an positive integer value. Wrapping the Mathematica command TableForm around BasisH[m, x] will produce an output display that is more readable than the usual list format:

SetDimension[x,3]

TableForm[BasisH[4,x]]

The option the option Orthonormal -> Ball, when used with BasisH, will produce a basis for (Subscript[H, m])((R^n))that is orthonormal with respect to the inner product on (L^2)(B, dV) (here B denotes the open unit ball in Euclidean space and V denotes volume measure on B):

TableForm[BasisH[4,x,Orthonormal->Ball]]

The option Orthonormal -> Sphere, when used with BasisH, will produce a basis for (Subscript[H, m])((R^n))that is orthonormal with respect to the inner product on (L^2)(S, dσ) (here S denotes the unit sphere in Euclidean space and σ denotes normalized surface area measure on S):

TableForm[BasisH[4,x,Orthonormal->Sphere]]

To produce an orthonormal basis for (Subscript[H, m])(S), the space of spherical harmonics on the unit sphere S, with respect to the inner product on (L^2)(S, dσ), use the option Orthonormal -> Sphere and in the result replace ‖x‖ with 1:

TableForm[BasisH[5,x,Orthonormal->Sphere]] /. \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1

BasisH is computed by using Theorem 5.25 of [ABR].

ZonalHarmonic

ZonalHarmonic[m, x, y] is the extended zonal harmonic (Subscript[Z, m])(x,y)as defined by 8.7 of [ABR]:

SetDimension[x,n]

ZonalHarmonic[5,x,y]

To obtain the zonal harmonic (Subscript[Z, m])(x,y)as a homogeneous polynomial in x of degree m with pole z on the unit sphere (as in Theorem 5.38 of [ABR]), use ZonalHarmonic[m, x, y] and in the result replace ‖z‖ with 1:

SetDimension[x,7]

ZonalHarmonic[6,x,z] /. \[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]->1

To obtain the zonal harmonic (Subscript[Z, m])(x,z) as an element of (Subscript[H, m])(S) as originally defined in Chapter 5 of [ABR], use ZonalHarmonic[m,x,z] and in the result replace ‖x‖ and ‖z‖ with 1:

ZonalHarmonic[8,x,z] /. {\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]->1,\[LeftDoubleBracketingBar]z\[RightDoubleBracketingBar]->1}

DimensionH

DimensionH[m, n] is the vector space dimension of (Subscript[H, m])((R^n)), the space of harmonic polynomials on (R^n)homogeneous of degree m:

DimensionH[12,100]

DimensionH is computed by using Proposition 5.8 of [ABR].

Inversion and the Kelvin Transform

Reflection

Reflection[x] is the reflection of a vector x with respect to the unit sphere. Thus the reflection of x equals (x/((‖x‖)^2)). The argument of Reflection can be either an expression representing a vector or a list of coordinates for a vector:

Reflection[2x+y]

Reflection[{1,-2,5, 11}]

Reflection[x, Sphere[r, c]] is the reflection of a vector x with respect to the sphere with radius r centered at c. Here each of  x and c can be either an expression representing a vector or a list of coordinates for a vector:

Reflection[{2,4, 5}, Sphere[7, {3, 1, 6}]]

Reflection[2x+y, Sphere[3, z]]

Reflection[x, Sphere[r, c]] is computed using equation 4.13 in [ABR].

Reflection[x, Hyperplane[b, t]] is the reflection of a vector x with respect to the hyperplane consisting of those vectors whose inner product with b equals t. In other words, if n denotes the dimension of x (which must equal the dimension of b), then the reflection is taken with respect to the hyperplane {y ϵ (R^n): b.y=t}. Here each of  x and b can be either an expression representing a vector or a list of coordinates for a vector:

Reflection[x , Hyperplane[{1, 4,5}, 7]]

Reflection[{(Subscript[x, 1]),(Subscript[x, 2])}, Hyperplane[{4,5}, 7]]

Reflection[x, Hyperplane[b, t]] is computed using Exercise 9 on page 71 of [ABR].

Kelvin

Kelvin[u, x] is the Kelvin transform of u, thought of as a function of x:

SetDimension[x,n]

Kelvin[(((x.y))^2)+((Subscript[x, 3])^5) Norm[x],x]

Phi

Φ[z] is is the modified inversion introduced on page 154 of [ABR]. The modified inversion Φ is defined by

Φ(z) = 2((z-S)/((‖z-S‖)^2)) + S,

where S is the south pole (0, 0, ..., 0, -1). The modified inversion Φ is useful for translating questions about harmonic functions on balls to questions about harmonic functions on half-spaces, and vice versa.

Here is a proof of that Φ(Φ(z)) = z, as claimed in Proposition 7.18 (a) of [ABR]:

Φ[z]

Φ[%]

Another, and often more convenient, form of Φ is given by Φ[x, y]. Here x is a vector, y is a real number, and thus (x, y) is a vector in a Euclidean space whose dimension is one more than the dimension of x.

Here is a proof of the assertion of Exercise 6 in Chapter 7 of [ABR], which asserts that

1 - ((‖Φ(x, y)‖)^2) = ((4y)/(((‖x‖)^2)+ (((y + 1))^2))):

Φ[x,y]

1-((\[LeftDoubleBracketingBar]%\[RightDoubleBracketingBar])^2)

KelvinM

KelvinM[u, z] is the modified Kelvin transform defined to be equal to

(2^((n-2)/2)) ((‖z - S‖)^(2-n)) u(Φ(z)).

This modified Kelvin transform is defined on page 155 of [ABR]. Here is a proof of the assertion made there that the modified Kelvin transform is its own inverse:

SetDimension[z,n]

KelvinM[u[z],z]

KelvinM[%, z]

PowerExpand[%]

Another form of KelvinM is given by KelvinM[u, x, y]. Here x is a vector, y is a real variable, and u is a function of (x, y), which is a vector in a Euclidean space whose dimension is one more than the dimension of x:

SetDimension[x,n-1]

KelvinM[x[4],x,y]

KelvinM[z[4],z]

If we had used the input KelvinM[ z[4], z ], the result would have included terms of the form (Subscript[S, 4]), which denotes the fourth coordinate of the south pole S. Of course (Subscript[S, 4]) equals either 0 or -1, but since the dimension of z has been set to equal the symbol n, the computer has no way of determining whether or not n equals 4. With the input KelvinM[ z[4], z ] there is no ambiguity; (Subscript[x, 4]) cannot be the last coordinate, because the last coordinate is denoted by y.

Kernels

PoissonKernel

PoissonKernel[x, y] is the extended Poisson kernel P(x,y) for the unit ball as defined by Equation 6.21 of [ABR]:

SetDimension[x,n]

PoissonKernel[x,y]

The usual Poisson kernel P(x, ζ), where ζ lies on the unit sphere, can be obtained by replacing ‖ζ‖ with 1:

PoissonKernel[x,ζ] /. \[LeftDoubleBracketingBar]ζ\[RightDoubleBracketingBar]-> 1

PoissonKernelH

PoissonKernelH[z, w] is the extended Poisson kernel (Subscript[P, H])(z, w) for the upper half-space H, as defined by Equation 8.22 of [ABR]:

SetDimension[z,5]

PoissonKernelH[z,w]

The usual Poisson kernel (Subscript[P, H])(z, w) for the upper half-space, where we think of w ϵ ∂H, can be obtained by replacing the last coordinate of w with 0:.

% /. (Subscript[w, 5])->0

PoissonKernelH[x, y, t, u] equals (Subscript[P, H])((x, y), (t, u)) and thus is another format for the extended Poisson kernel on the upper half-space. Here x and t denote vectors, and y and u denote nonnegative numbers. Note that if we want to work in the upper half-space of (R^n), then the dimension of x should be set to n-1:

SetDimension[x,n-1]

PoissonKernelH[x,y,t,u]

The usual Poisson kernel  (Subscript[P, H])((x, y), t)  for the upper half-space, where t ϵ ∂H, can be setting the last argument in this format to 0:

PoissonKernelH[x,y,t,0]

BergmanKernel

BergmanKernel[x, y] is the reproducing kernel for the harmonic Bergman space of the unit ball:

SetDimension[x,10]

BergmanKernel[x,y]

BergmanKernel is computed using Theorem 8.13 of [ABR].

BergmanKernelH

BergmanKernelH[z, w] is the reproducing kernel for the harmonic Bergman space of the upper half-space:

SetDimension[z,10]

BergmanKernelH[z,w]

BergmanKernelH[x, y, t, u] is another format for the reproducing kernel for the harmonic Bergman space of the upper half-space. Here x and t denote vectors, y and u denote nonnegative numbers, and the kernel is evaluated at ((x, y), (t, u)). Note that if we want to work in the upper half-space of (R^n), then the dimension of x should be set to n-1:

SetDimension[x,n-1]

BergmanKernelH[x,y,t,u]

BergmanKernelH is computed by using Theorem 8.24 of [ABR].

Miscellaneous

BergmanProjection

BergmanProjection[u, x] is the orthogonal projection of uonto (b^2)(B), the harmonic Bergman space of the unit ball (see Chapter 8 of [ABR]). Here u must be a polynomial function of x:

SetDimension[x,5]

BergmanProjection[((Subscript[x, 1])^5) ((Subscript[x, 2])^3),x]

BergmanProjection is computed by using Corollary 8.15 of [ABR].

HarmonicConjugate

HarmonicConjugate[u, x, y] is the harmonic conjugate of u on (R^2), where the coordinates in (R^2)are denoted by x, y and u is a harmonic function of (x, y):

HarmonicConjugate[15(x^2) y+12(x^3) y-5(y^3)-12x (y^3),x,y]

HarmonicConjugate is computed by using Exercise 11 in Chapter 1 of [ABR].

Schwarz

Schwarz[x] is the maximum of |u(x)|, where u ranges over all harmonic functions on the unit ball (in the Euclidean space whose dimension is determined by x) with u(0) = 0 and |u| < 1 on the ball. Thus, in the notation used in Theorem 6.24 of [ABR] to present the Harmonic Schwarz Lemma, Schwarz[x] equals U(‖x‖ N).  The user must first set the dimension of the Euclidean space to a positive integer value:

SetDimension[x,5]

Schwarz[x]

Schwarz[x] is computed by evaluating the Poisson integral that defines U(‖x‖ N), iterating an appropriate formula from Appendix A of [ABR].

Togetherness

When not using the HFT6.m software package, Mathematica might tell the user that the result of a calculation is (t^2)+(1-t) (1 + t).  Of course the user will easily recognize that the last quantity equals 1, but with more complicated results the user will need to apply appropriate Mathematica commands (such as Simplify, Expand, or Together) to get the result in simplest form.  Applying Simplify to all Mathematica output is poor practice, because Mathematica can take a long time to realize that some expressions cannot be  simplified. For almost all output generated by the HFT6.m package,  applying Together to the result and all its subexpressions puts  it in simplest form.  Thus the HFT6.m package redefines the Mathematica variable $Post so that for each output X, Mathematica will display either X or the result of applying Together to X and all its subexpressions, whichever is shorter (if the user had defined $Post before starting the HFT6.m package, the HFT6.m version of $Post is composed with the user's). This procedure is invisible to the user and usually requires only a short amount of computer time. If you want to turn off this feature, use the command TurnOff[Togetherness] (if the user had defined $Post before starting the HFT6.m package, this restores the user's version).  The command TurnOn[Togetherness] turns this feature back on again.

If Homogeneous or Taylor is used to find an expansion about a point other than the origin, then Togetherness will be automatically turned off. If this were not done then, for example, the last output in the Taylor subsection earlier in this notebook would be simplified to 1 + (Subscript[x, 1])(Subscript[x, 2]) + ((Subscript[x, 1])^2), which is not what the user wanted to see when asking for an expansion about b.

References

[ABR] Sheldon Axler, Paul Bourdon, and Wade Ramey, Harmonic Function Theory, second edition, Springer Graduate Texts in Mathematics, 2001.

[AGV] Sheldon Axler, Pamela Gorkin, and Karl Voss, The Dirichlet problem on quadratic surfaces, Mathematics of Computation 73 (2004), 637-651.

[AR] Sheldon Axler and Wade Ramey, Harmonic polynomials and Dirichlet-type problems, Proceedings of the American Mathematical Society 123 (1995), 3765-3773.

[R] Walter Rudin, Real and Complex Analysis, third edition, McGraw-Hill, 1987.

[W] Hermann Weyl, On the volume of tubes, American Journal of Mathematics 61 (1939), 461-472.