HFT

Initialization

Print

If[$VersionNumber<7,Print["This package may not work correctly with versions of Mathematica earlier than Mathematica 7. Versions of this package that will work correctly with earlier versions of Mathematica are available at http://www.axler.net."]]

Print["HFT11.m, version 11.02, 4 July 2017; for use with Mathematica versions 7 through 11 (and probably later versions of Mathematica)."]

Print["The HFT11.m Mathematica package is designed for computing with harmonic functions."]

Print["Documentation for the use of this package and information about the algorithms used in it is available in the document titled Computing with Harmonic Functions, which is available in both nb and pdf formats."]

Print["The most recent version of this HFT11.m package and its documentation Computing with Harmonic Functions are available at http://www.axler.net."]

Print["For addional information about harmonic functions, see the book Harmonic Function Theory (second edition), by Sheldon Axler, Paul Bourdon, and Wade Ramey, published by Springer."]

Print["This package is copyrighted by Sheldon Axler but is distributed without charge."]

Print["Comments, suggestions, and bug reports should be sent to axler@sfsu.edu."]

Usage

BeginPackage["HFT`"]

Δ::"usage"="(Subscript[Δ, x])[f] gives the Laplacian of f with respect to x."

annulus::"usage"="annulus is an option for region in dirichlet."

antiLaplacian::"usage"="antiLaplacian[u, x] gives a function whose Laplacian 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 \[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]."

basisH::"usage"="basisH[m, x] gives a basis for the space of harmonic polynomials homogeneous of degree m in the variable x."

bergmanKernel::"usage"="bergmanKernel[x, y] gives the Bergman reproducing kernel for the unit ball."

bergmanKernelH::"usage"="bergmanKernelH[z, w] gives the Bergman reproducing kernel for the upper half-space."

bergmanProjection::"usage"="bergmanProjection[ f, x] gives the orthogonal projection of a polynomial f, as a function of x, onto the Bergman space of harmonic functions on the ball."

biDirichlet::"usage"="biDirichlet[p , x] solves the biDirichlet problem with boundary data p, as a function of x."

δ::"usage"="(Subscript[δ, j]) is the vector that equals 1 in the j-th coordinate and 0 in the other coordinates."

dimension::"usage"="dimension[x] is the dimension of x, which might be a vector or a matrix."

dimHarmonic::"usage"="dimHarmonic[m, n] gives the vector space dimension of the space of harmonic polynomials on (R^n) homogeneous of degree m."

dirichlet::"usage"="dirichlet[p, x] gives the harmonic function that equals p on the unit sphere. Here p must be a polynomial function of x."

divergence::"usage"="divergence[f, x] gives the divergence of f with respect to x."

expandNorm::"usage"="expandNorm[f] gives f with all terms of the form \[LeftDoubleBracketingBar]x + y\[RightDoubleBracketingBar] replaced by \!\(\*SqrtBox[\(\*SuperscriptBox[\(\[LeftDoubleBracketingBar]x\[RightDoubleBracketingBar]\), \(2\)]\\  + \\ \*SuperscriptBox[\(\[LeftDoubleBracketingBar]y\[RightDoubleBracketingBar]\), \(2\)]\\  + \\ 2\\ x . y\)]\)."

exteriorNeumann::"usage"="exteriorNeumann[p , x] solves the exterior Neumann problem with boundary data p, as a function of x."

exteriorSphere::"usage"="exteriorSphere is an option for dirichlet, specifying that the region should be the exterior of the unit sphere."

gradient::"usage"="gradient[f, x] gives the gradient of f with respect to x."

harmonicConjugate::"usage"="harmonicConjugate[u, x, y] gives the harmonic conjugate of u with respect to the variables x, y."

harmonicDecomposition::"usage"="harmonicDecomposition[p, x] gives the harmonic decomposition of p with respect to the variable x."

hilbertSchmidt::"usage"="hilbertSchmidt[A] gives the Hilbert-Schmidt norm of a matrix A."

homogeneous::"usage"="homogeneous[f, d, x] gives the homogenous part of f of degree d, with respect to x."

If[ $VersionNumber< 10.2, Hyperplane::"usage"="Hyperplane[n, c] represents the hyperplane with normal n given by the points x that satisfy n.x = c."]

integrateBall::"usage"="integrateBall[f, x] gives the integral of f, as a function of x, over the unit ball."

integrateEllipsoidArea::"usage"="integrateEllipsoidArea[f, b, c, d, x] gives the integral of f/\[LeftDoubleBracketingBar]∇q\[RightDoubleBracketingBar] as a function of x with repsect to surface area measure over the ellipsoid {x: q(x) = 0}, where q(x) = b.(x^2)+ c.x + d. Here b and c must be lists whose length equals dimension[x] and f must be a polynomial function of x. If c is the 0 vector, then c can be omitted. If d = -1, then d can be omitted."

integrateEllipsoidVolume::"usage"="integrateEllipsoidVolume[f, b, c, d, x] gives the integral of f as a function of x with repsect to volume measure over the region {x: q(x) < 0}, where q(x) = b.(x^2)+ c.x + d. Here b and c must be lists whose length equals dimension[x] and f must be a polynomial function of x. If c is the 0 vector, then c can be omitted. If d = -1, then d can be omitted."

integrateSphere::"usage"="integrateSphere[f, x] gives the integral of f, as a function of x, over the unit sphere with respect to normalized surface area measure."

jacobian::"usage"="jacobian[f, x] gives the Jacobian of f with respect to x."

kelvin::"usage"="kelvin[u, x] gives the Kelvin transform of u, thought of as a function of x."

kelvinH::"usage"="kelvinH is the modified Kelvin transform, as defined in Chapter 7 of Harmonic Function Theory."

laplacian::"usage"="laplacian equals Δ, the Laplacian."

multiple::"usage"="multiple is an option for antiLaplacian. The default value is None.  The value norm^2 produces the unique antiLaplacian that is a polynomial multiple of \[LeftDoubleBracketingBar]x(\[RightDoubleBracketingBar]^2), where x is the variable."

neumann::"usage"="neumann[p, x] solves the Neumann problem of finding a harmonic function of x whose outward normal derivative (on the unit sphere) is p."

norm::"usage"="norm[x] gives the norm of a vector x."

normalD::"usage"="normalD[f, z] gives the outward normal derivative of f, as a function of z, with respect to the unit ball."

partial::"usage" = "partial[f, (Subscript[x, j])] gives the partial derivative of f with respect to (Subscript[x, j])."

Φ::"usage"="Φ[z] is the modified reflection defined in Chapter 7 of Harmonic Function Theory."

poissonKernel::"usage"="poissonKernel[x, z] gives the Poisson kernel for the unit ball."

poissonKernelH::"usage"="poissonKernelH[x, y, t] gives the Poisson kernel for the upper half-space."

quadratic::"usage"="quadratic is an option for multiple in antiLaplacian and for region in dirichlet."

reflection::"usage"="reflection[x] gives the reflection of x in the unit sphere."

region::"usage"="region is an option for dirichlet."

southPole::"usage"="southPole is the south pole."

schwarz::"usage"="schwarz[x] gives the maximum of |u[x]|, where u ranges over all harmonic functions on the unit ball with u[0] = 0 and |u| < 1."

setDimension::"usage"="setDimension[x, n] sets the dimension of x to n."

singularity::"usage"="singularity is an option for antiLaplacian. The default value is None."

surfaceArea::"usage"="surfaceArea[n] gives the surface area of the unit ball in n-dimensional real Euclidean space."

taylor::"usage"="taylor[f, d, x] gives the Taylor expansion of f to degree d, with respect to x."

togetherness::"usage"="turnOff[ togetherness ] turns off the feature that simplifies some output. turnOn[ togetherness ] turns this feature back on."

trace::"usage"="trace[A] gives the trace of the matrix A."

turnOff::"usage"="turnOff[F] turns feature F off."

turnOn::"usage"="turnOn[F] turns feature F on."

volume::"usage"="volume[n] gives the volume of the unit ball in n-dimensional real Euclidean space."

zeroToZero::"usage"="turnOn[ zeroToZero ] turns on the feature that changes 0. to 0 via $Post.  turnOff[ zeroToZero ] turns this feature back off."

zonalHarmonic::"usage"="zonalHarmonic[m, x, z] gives the zonal harmonic of degree m with pole x."

Protection

Begin["`private`"]

Unprotect[IntegerQ,Limit,Transpose];

protectedWords=Unprotect[Abs,Dot,IdentityMatrix, Integer, Power,Times]

$Post

SetAttributes[postTogether,HoldAll]

postTogether[Apart[f_]]:=Apart[f]

postTogether[Together[f_]]:=Together[f]

postTogether[Simplify[f_]]:=Simplify[f]

postTogether[Expand[f_]]:=Expand[f]

postTogether[a_=b_]:=a=postTogether[b]

shorter[f_] := With[ {g = Together[f]},
	If[ StringLength[ ToString[ Format[ g, InputForm ] ] ] < 
	StringLength[ ToString[ Format[ f, InputForm ] ] ], g, f ] ]

shorter/:f_Symbol[a___,HoldPattern[shorter][x_],b___]:=f[a,x,b]/;MemberQ[Attributes[f],HoldAll]||MemberQ[Attributes[f],HoldFirst]||MemberQ[Attributes[f],HoldRest]

postTogether[f_]:=Map[shorter,f,{0,-2}]//.HoldPattern[shorter][j_]->j

zeroToZero[f_] := (f /. 0. -> 0);

usersPost=If[ToString[$Post]==="$Post",Identity,$Post];

togetherness::Off="togetherness has been turned off.  "<>"The command  turnOn[togetherness]  will turn it back on."

togetherness::On="togetherness has been turned on.  "<>"The command  turnOff[togetherness]  will turn it back off."

togetherness::AlreadyOn="togetherness is already turned on."

togetherness::AlreadyOff="togetherness is already turned off."

zeroToZero::Off="zeroToZero has been turned off.  "<>"The command  turnOn[zeroToZero]  will turn it back on."

zeroToZero::On="zeroToZero has been turned on.  "<>"The command  turnOff[zeroToZero]  will turn it back off."

zeroToZero::AlreadyOn="zeroToZero is already turned on."

zeroToZero::AlreadyOff="zeroToZero is already turned off."

turnOn[togetherness]^:= If[ togetherness2 , Message[togetherness::AlreadyOn],$Post=If[zerotozero,Composition[zeroToZero, postTogether, usersPost],Composition[postTogether,usersPost]];togetherness2=True; Message[togetherness::On]]

turnOn[zeroToZero]^:= If[ zerotozero , Message[zeroToZero::AlreadyOn],$Post=If[togetherness2,Composition[zeroToZero, postTogether, usersPost],Composition[zeroToZero,usersPost]];zerotozero=True; Message[zeroToZero::On]]

$Post=Composition[postTogether, usersPost];

togetherness2=True; zerotozero = False;

turnOff[zeroToZero]^:= If[ zerotozero ,$Post=If[togetherness2,Composition[postTogether,usersPost], usersPost];zerotozero=False;Message[zeroToZero::Off],  Message[zeroToZero::AlreadyOff]]

turnOff[togetherness]^:= If[ togetherness2 ,$Post=If[zerotozero,Composition[zeroToZero,usersPost], usersPost];togetherness2=False;Message[togetherness::Off],  Message[togetherness::AlreadyOff]]

Calculus on (R^n)

IdentityMatrix

matrix[IdentityMatrix[_]]^=True;

IdentityMatrix/:a_ .IdentityMatrix[_]:=a

IdentityMatrix/:IdentityMatrix[_].v_:=v

Transpose

Off[Transpose::"nmtx"]

matrix[Transpose[_]]=True

Transpose[Transpose[x_]]:=x

Transpose[{c_ x_?vector}]:=c Transpose[{x}]

Transpose[c_ x_?matrix]:=c Transpose[x]

Transpose[IdentityMatrix[n_]]:=IdentityMatrix[n]

Transpose[{0}]:={0}

trace

trace[IdentityMatrix[n_]]:=n

trace[Transpose[A_]]:=trace[A]

trace[c_ A_]:=c trace[A]/;NumericQ[c]

trace[c_ A_]:=c trace[A]/;matrix[A]

trace[A_]:=Tr[A] /; MatrixQ[A]

norm

Format[norm[x_], StandardForm] := DoubleBracketingBar[x]

MakeExpression[RowBox[{"\[LeftDoubleBracketingBar]", x_, "\[RightDoubleBracketingBar]"}],StandardForm]  := MakeExpression[RowBox[{"norm", "[", x, "]"}], StandardForm]

norm[{x_?vector,y_}]:=(unmakeVector[y];(Sqrt[((norm[x])^2)+(y^2)]))

norm[0] :=0

norm[x_List]:= (Sqrt[Plus @@ (x^2)])

Using norm[x_List]:= Norm[x] leads to (Sqrt[Plus @@ Abs[(x^2)]]), which can lead to differentiation difficulties.

norm[c_ x_?vector]:=(realize[c];Abs[c] norm[x])

norm[x_]:=0/;(makeVector[x];False)

norm[c_?NumericQ x_]:=Abs[c] norm[x]

real and realize

real = {};

realize[c_] := (real = real⋃{c})

expandNorm

expandNorm[f_]:=f//.norm[g_+h_]:>(Sqrt[((norm[g])^2)+((norm[h])^2)+2 g.h])

Abs

Abs[((norm[x_])^(p_.))]:=((norm[x])^p)

Abs/:((Abs[(Subscript[x_, j_])])^(p_?EvenQ)):=((Subscript[x, j])^p)/;vector[x]

Abs/:((Abs[c_])^(p_?EvenQ)):=(c^p)/;MemberQ[real, c] ||MemberQ[real,1/c]||MemberQ[real,(c^(p/2))]||MemberQ[real,(1/(c^(p/2)))]==0

vector

vector[HoldPattern[+a__]]:=If[TrueQ[Or@@Thread[vector[{a}]]],makeVector[+a];True]

vector[HoldPattern[Times[a__]]]:=Or@@Thread[vector[{a}]]

vector[A_?matrix.x_?vector]:=True

vector[x_?vector.Y_?matrix]:=True

vector[a_[".",_]]:=True

vector[a_[_,"."]]:=True

(Subscript[0, _])=0

vector[(Subscript[partial, _])[b_][x_]]:=True/;vector[b[x]]

vector[gradient[_][_]]=True

makeVector

vectors={};

vectorValuedFunctions={};

makeVector[x_Symbol]:=If[MemberQ[vectors,x],Null,x[j_]=(Subscript[x, j]);AppendTo[vectors,x];vector[x]^=True;]

makeVector[c_?NumericQ x_]:=makeVector[x]

makeVector[HoldPattern[+a__]]:=makeVector/@{a}

makeVector[f_Symbol[_]]:=If[MemberQ[vectorValuedFunctions,f],Null,f[x_][j_]=(Subscript[f, j])[x];
		AppendTo[vectorValuedFunctions,f];vector[f[_]]^=True;vector[(f^((_)))[_]]=True;vector[(Subscript[partial, _])[f][_]]=True;]/;f=!=List

makeVector[(Subscript[partial, _])[u_][_]]:=makeVector[u[_]]

Attributes[unmakeVector]={Listable};

unmakeVector[x_Symbol]:= (Unprotect[Abs];Abs/:((Abs[x])^(p_?EvenQ)):=(x^p);If[MemberQ[vectors,x],(x[j$_])=.;vectors=Delete[vectors,Position[vectors,x][[1,1]]];vector[x]^=False;])

matrix

matrix[HoldPattern[Times[a__]]]:=Or@@Thread[matrix[{a}]]

matrix[HoldPattern[+a__]]:=Or@@Thread[matrix[{a}]]

matrix[a_?matrix.x_?matrix]:=True

matrix[jacobian[_][_]]=True;

matrix[{_}]=True;

(Subscript[0, _,_])=0

matrix[Transpose[x_]]:=matrix[x]

makeMatrix

matrices={};

makeMatrix[x_Symbol]:=If[TrueQ[matrix[x]],Null,x[j_,k_]=(Subscript[x, j,k]);AppendTo[matrices,x];matrix[x]^=True;vector[x]^=False]

Dot

Attributes[Dot]={OneIdentity}

b_ . 0 = 0;

0 . b_ = 0;

HoldPattern[+x__].z_:=Plus@@(#1.z&)/@{x}

a_ .HoldPattern[+x__]:=Plus@@(a.#1&)/@{x}

(c_ x_).y_:=c x.y/;vector[x]||matrix[x]||NumericQ[c]

{c_ x_}.y_:=c {x}.y/;vector[x]||matrix[x]||NumericQ[c]

a_ .(c_ y_):=c a.y/;vector[y]||matrix[y]||NumericQ[c]

a_ .{c_ y_}:=c a.{y}/;vector[y]||matrix[y]||NumericQ[c]

x_ .x_:=((norm[x])^2)/;vector[x]

f_[x_] . f_[x_] := ((norm[f[x]])^2) /; Head[f]==Symbol &&vector[f[x]]

Derivative[s_][f_][x_] . Derivative[s_][f_][x_] := ((norm[Derivative[s][f][x] ])^2)/; Head[f] == Symbol &&vector[f[x]]

(Subscript[partial, j_])[f_][x_] . (Subscript[partial, j_])[f_][x_] := ((norm[ (Subscript[partial, j])[f][x]])^2)/; Head[f] == Symbol &&vector[f[x]]

(Subscript[δ, j_]).(Subscript[δ, k_]):=(Subscript[δ, j,k])

(Subscript[δ, j_]) . x_Symbol:=(makeVector[x];(Subscript[x, j]))/;matrix[x]=!=True

a_Symbol.(Subscript[δ, j_]):=(makeVector[a];(Subscript[a, j]))/;matrix[a]=!=True

(Subscript[δ, j_]) . f_Symbol[x_]:=(makeVector[f[_]];f[x][j])/;matrix[f[x]]=!=True/;f=!=List

b_Symbol[x_].(Subscript[δ, j_]):=(makeVector[b[_]];b[x][j])/;matrix[b[x]]=!=True/;b=!=List

(Subscript[δ, j_]) . gradient[u_][x_]:=(Subscript[partial, {j}])[u][x]

gradient[u_][x_]. (Subscript[δ, j_]):=(Subscript[partial, {j}])[u][x]

jacobian[b_][x_].(Subscript[δ, j_]):=(Subscript[partial, {j}])[b][x]

(Subscript[δ, j_]) . jacobian[b_][x_]:=gradient[b.(Subscript[δ, j])][x]

(Subscript[δ, j_]) . Y_?matrix:=(Subscript[Y, j,"."])

A_?matrix.(Subscript[δ, j_]):=(Subscript[A, ".",j])

(Subscript[δ, k_]) . (Subscript[Y, j_,"."]):=(Subscript[Y, j,k])

(Subscript[A, j_,"."]) . (Subscript[δ, k_]):=(Subscript[A, j,k])

(Subscript[δ, j_]) . (Subscript[Y, ".",k_]):=(Subscript[Y, j,k])

(Subscript[A, ".",k_]) . (Subscript[δ, j_]):=(Subscript[A, j,k])

a_?VectorQ . (Subscript[δ, j_]):=a[[j]]

(Subscript[δ, j_]) . v_?VectorQ:=v[[j]]

A_?MatrixQ . (Subscript[δ, j_]):=Transpose[A][[j]]

(Subscript[δ, j_]) . Y_?MatrixQ:=A[[j]]

x_ .A_?identityMatrixMultiple:=A[[1,1]] x/;primitiveVector[x]

A_?identityMatrixMultiple.x_:=A[[1,1]] x/;primitiveVector[x]

identityMatrixMultiple[A_]:=MatrixQ[A]&&Length[A]>1&&Length[A]===Length[A[[1]]]&&And@@(A[[#1,#1]]===A[[1,1]]&)/@Range[Length[A]]&&And@@Flatten[Table[A[[m,n]]===0||m===n,{m,Length[A]},{n,Length[A]}]]

v_ . Transpose[{x_Symbol}].{x_Symbol} := (v.x)x

a_List.x_Symbol:=(makeVector[x];a.Table[(Subscript[x, j]),{j,Length[a]}])/;matrix[x]=!=True

x_Symbol.y_List:=(makeVector[x];Table[(Subscript[x, j]),{j,Length[y]}].y)/;matrix[x]=!=True

a_List.B_Symbol:=a.Table[ (Subscript[B, j,k]), {j, dimension[B][[1]]}, {k,  dimension[B][[2]]}]/;matrix[B]

B_Symbol.c_List:=Table[ (Subscript[B, j,k]), {j, dimension[B][[1]]}, {k,  dimension[B][[2]]}].c/;matrix[B]

x_?primitiveVector.y_?primitiveVector:=y.x/;Order[ToString[y],ToString[x]]==1&&matrix[y]=!=True&&matrix[x]=!=True

primitiveVector[x_Symbol]:=True/;Context[x]=!="System`"

primitiveVector[gradient[_][_]]=True

primitiveVector[δ[_]]=True

degree

degree[f, z] gives the degree of f as a polynomial function of z. The actual degree of f might be less if there are cancellations involving norms or dot products.

degree[f_,z_]:=degree1[Expand[f/.norm[z]->(Subscript[z, normReplacement])],z]

degree[f_,z_]:=degree1[Expand[f/.{norm[z]->(Subscript[z, normReplacement]), z . _ ->(Subscript[z, normReplacement2]), _ . z -> (Subscript[z, normReplacement3])}] ,z]

degree1[0,_]:=-∞

degree1[HoldPattern[+a__],z_]:=Max[Thread[degree1[{a},z]]]

degree1[f_,z_]:=Exponent[f/.(Subscript[z, _])->z,z]

monomialQ

uninomialQ[f_, x_Symbol] := MatchQ[f,(Subscript[x, _])] || MatchQ[f, ((Subscript[x, _])^(_?((IntegerQ[#1] &&Positive[#1])&)))]

monomialQ[f_, x_Symbol] := uninomialQ[f,x] || (Head[f] === Times && (And@@(uninomialQ[f[[#1]],x]&)/@Range[Length[f]]))

polynomialDecomposition and blockPolynomialDecomposition

polynomialDecomposition[p, x], where p is a polynomial in x, is a list of elements of the form {q, m}, where q is the homogeneous part of p of degree m.

polynomialDecomposition[0,x_Symbol]:=Null

polynomialDecomposition[p_,x_Symbol]:=With[{p1=Expand[p/.norm[x]->(Subscript[x, normReplacement])]},With[{t=({#1,degree[#1,x]}&)/@If[Head[p1]===Plus,List@@p1,{p1}]},With[{s=Table[If[MemberQ[t,{_,j}],{Plus@@Transpose[Cases[t,{_,j}]][[1]],j},Null],{j,0,Max[Transpose[t][[2]]]}]},Delete[s,Position[s,Null]]]]]/.(Subscript[x, normReplacement])->norm[x]

exponentsMod2[t_,x_] :=Drop[Mod[Cases[t  (((Times @@ Array[x, dimension[x]+1]))^2), Power[Subscript[x,_],_]][[All,2]],2], -1]

block[p_, x_, ex_] := Select[p, (exponentsMod2[#,x] == ex) &]

blocks[{p_, m_}, x_] := {{p, m, ConstantArray[0, dimension[x]]}} /; m == 0

blocks[{p_, m_}, x_] := {{p, m, exponentsMod2[p,x]}} /;Not[ Head[p] === Plus]

blocks[{p_, m_}, x_] := Module[{q=p, b = {}, ex, bq},While[Head[q]===Plus,ex = exponentsMod2[q[[1]], x]; bq = block[q, x, ex]; AppendTo[b, {bq,m,ex}]; q = q - bq];If[q===0,b, Append[b, {q, m, exponentsMod2[q,x]}]]]

blockPolynomialDecomposition[p, x], where p is a polynomial in x, is a list of elements of the form {q, m, b}, where b is a list of 0’s and 1’s of length dimension[x] and q is the homogeneous part of p of degree m such that each monomial in q has each ((Subscript[x, j])^k) term with k congruent to b[[j]] mod 2.

blockPolynomialDecomposition[p_, x_] :=
 Flatten[blocks[#, x]& /@ polynomialDecomposition[p,x],1]

δ (delta)

vector[(Subscript[δ, _])]^:=True

(Subscript[δ, j_,k_]):=0/;j≠k

(Subscript[δ, j_,k_]):=1/;j==k

(Subscript[δ, j_])[k_]:=(Subscript[δ, j,k])

norm[(Subscript[δ, _])]^=1;

partial

partial[f_,{(Subscript[x_, _]),0}]:=(makeVector[x];f)

partial[f_,{y_Symbol,0}]:=(unmakeVector[y];f)

partial[f_,{z_,n_Integer}]:=Module[{g = f}, Do[g = Together[partial[g,z]],{n}];g]/;Positive[n]

partial[f_,z1_,z__]:=Module[{g =partial[f,z1]}, Do[g = partial[g, {z}[[j]]], {j, Length[{z}]}]; g]

partial[f_,(Subscript[x_, _])]:=(makeVector[x];0)/;FreeQ[f,x]

partial[f_,y_Symbol]:=(unmakeVector[y];0/;FreeQ[f,y])

partial[(Subscript[x_, j_]),(Subscript[x_, k_])]:=(makeVector[x];(Subscript[δ, j,k]))

partial[y_,y_]:=1

partial[x_,(Subscript[x_, j_])]:=(makeVector[x];(Subscript[δ, j]))

partial[HoldPattern[+a__],y_]:=Plus@@Thread[partial[{a},y]]

partial[f_ g_,y_]:=f partial[g,y]+partial[f,y] g

partial[f_ .g_,y_]:=(If[TrueQ[matrix[f]],Null,makeVector[f]];If[TrueQ[matrix[g]],Null,makeVector[g]];f.partial[g,y]+partial[f,y].g)

partial[f_[v_List],t_]:=Plus@@((Subscript[partial, {#1}])[f][v] * partial[v[[#1]],t]&)/@Range[Length[v]]

partial[((f_)^(p_)),y_]:=p (f^(p-1)) partial[f,y]+(f^p) Log[f] partial[p,y]

partial[f_[g_],(Subscript[x_, j_])]:=(makeVector[x];If[TrueQ[vector[g]],If[TrueQ[vector[f[g]]],jacobian[f][g].partial[g,(Subscript[x, j])],gradient[f][g].partial[g,(Subscript[x, j])]],f'[g] partial[g,(Subscript[x, j])]])

partial[f_[g_],y_Symbol]:=(unmakeVector[y];If[TrueQ[vector[g]],If[TrueQ[vector[f[g]]],jacobian[f][g].partial[g,y],gradient[f][g].partial[g,y]],f'[g] partial[g,y]])

(Subscript[partial, j_List])[(Subscript[partial, k_List])[f_]]:=(Subscript[partial, Sort[Join[k,j]]])[f]

normalD

normalD[f_,z_Symbol]:=gradient[f,z].z/.norm[z]->1

normalD[f_,v_List]:=(unmakeVector[v];Together[gradient[f,v].v]/.((norm[v])^2)->1)

normalD[f_, q_, z_Symbol] := gradient[f, z] . gradient[q, z]/(Sqrt[gradient[q, z]. gradient[q, z]])

gradient

gradient[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x]

gradient[(Subscript[x_, j_]),x_Symbol]:=(makeVector[x];(Subscript[δ, j]))

gradient[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[gradient[{a},x]]

gradient[f_ g_,x_Symbol]:=f gradient[g,x]+gradient[f,x] g

gradient[f_ .g_,x_Symbol]:=f.jacobian[g,x]+g.jacobian[f,x]

gradient[((f_)^(p_)),x_Symbol]:=p (f^(p-1)) gradient[f,x]+(f^p) Log[f] gradient[p,x]

gradient[f_[y_],x_Symbol]:=(makeVector[x];If[TrueQ[vector[y]],gradient[f][y].jacobian[y,x],f'[y] gradient[y,x]])

gradient[ Gamma[t_, f_], x_Symbol]:= -(E^(-f))(f^(t-1))gradient[f, x] /; FreeQ[t, x]

gradient[norm]:=((#1)/(norm[#1]))&

gradient[dimension]:=0&

gradient[f_,v_List]:=(unmakeVector[v];(partial[f,#1]&)/@v)

divergence

divergence[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x]

divergence[x_,x_Symbol]:=(makeVector[x];dimension[x])

divergence[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[divergence[{a},x]]

divergence[f_ g_,x_Symbol]:=f divergence[g,x]+gradient[f,x].g/;(makeVector[x];vector[g])

divergence[A_?matrix.f_,x_Symbol]:=trace[A.jacobian[f,x]]/;FreeQ[A,x]

divergence[f_ .T_?matrix,x_Symbol]:=trace[Transpose[T].jacobian[f,x]]/;FreeQ[T,x]

divergence[u_List,v_List]:=(unmakeVector[v];Plus@@Thread[partial2[u,v]]/.partial2->partial)

jacobian

jacobian[f_,x_Symbol]:=(makeVector[f];makeVector[x];0)/;FreeQ[f,x]

jacobian[x_,x_Symbol]:=(makeVector[x];IdentityMatrix[dimension[x]])

jacobian[HoldPattern[+a__],x_Symbol]:=Plus@@Thread[jacobian[{a},x]]

jacobian[f_ g_?vector,x_Symbol]:=f jacobian[g,x]+Transpose[{gradient[f,x]}].{g}

jacobian[f_Symbol[x_],x_Symbol]:=(makeVector[x];makeVector[f[_]];jacobian[f][x])/;f=!=List

jacobian[f_[y_],x_Symbol]:=If[TrueQ[vector[y]],jacobian[f][y].jacobian[y,x],Transpose[{f'[y]}].{gradient[y,x]}]

jacobian[A_?matrix.f_,x_Symbol]:=(makeVector[f];A.jacobian[f,x])/;FreeQ[A,x]

jacobian[f_ .T_?matrix,x_Symbol]:=(makeVector[f];Transpose[T].jacobian[f,x])/;FreeQ[T,x]

dimension[jacobian[f_][x_]]:= {dimension[f[_]],dimension[x] }

jacobian[f_List,v_List]:=(unmakeVector[v];Table[partial[f[[m]],v[[n]]],{m,Length[f]},{n,Length[v]}])

hilbertSchmidt

hilbertSchmidt[A_?MatrixQ]:=Plus@@Plus@@(A^2)

hilbertSchmidt[IdentityMatrix[n_]]:=n

hilbertSchmidt[ c_?NumericQ A_] := Abs[c] hilbertSchmidt[A]

hilbertSchmidt[ c_ A_] := Abs[c] hilbertSchmidt[A] /; matrix[A]

Δ (laplacian)

(Subscript[Δ, v_List])[f_]:=(unmakeVector[v];Plus@@(partial[f,{#1,2}]&)/@v)

(Subscript[Δ, x_Symbol])[f_]:=(makeVector[x];0)/;FreeQ[f,x]

(Subscript[Δ, x_Symbol])[(Subscript[x_, j_])]:=(makeVector[x];0)

(Subscript[Δ, x_Symbol])[HoldPattern[+a__]]:=Plus@@Thread[(Subscript[Δ, x])[{a}]]

(Subscript[Δ, x_Symbol])[f_ g_]:=f (Subscript[Δ, x])[g]+(Subscript[Δ, x])[f] g+2 gradient[f,x].gradient[g,x]

(Subscript[Δ, x_Symbol])[((f_)^(p_))]:=p (f^(p-1)) (Subscript[Δ, x])[f]+p (p-1) (f^(p-2)) gradient[f,x] . gradient[f,x] +2 p (f^(p-1)) Log[f] gradient[f,x].gradient[p,x]+(f^p) Log[f] (Subscript[Δ, x])[p]+2 (f^(p-1)) gradient[f,x].gradient[p,x]+(f^p) ((Log[f])^2) gradient[p,x] . gradient[p,x]

(Subscript[Δ, x_Symbol])[f_[a_. x_+b_.]]:=(makeVector[x];If[b==0,Null,makeVector[b]];(a^2) Δ[f][a x+b])/;FreeQ[a,x]&&FreeQ[b,x]

(Subscript[Δ, x_Symbol])[a_ .x_]:=(makeVector[x];0)/;FreeQ[a,x]

(Subscript[Δ, x_Symbol])[x_ .y_]:=(makeVector[x];0)/;FreeQ[y,x]

(Subscript[Δ, x_Symbol])[(A_ .x_).y_]:=(makeVector[x];0)/;FreeQ[A,x]&&FreeQ[y,x]

(Subscript[Δ, x_Symbol])[(x_ .Y_).z_]:=(makeVector[x];0)/;FreeQ[Y,x]&&FreeQ[z,x]

(Subscript[Δ, x_Symbol])[a_ .(B_ .x_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[B,x]

(Subscript[Δ, x_Symbol])[a_ .(x_ .Y_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[Y,x]

(Subscript[Δ, x_Symbol])[(A_ .x_).x_]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]

(Subscript[Δ, x_Symbol])[(x_ .A_).x_]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]

(Subscript[Δ, x_Symbol])[x_ .(A_ .x_)]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]

(Subscript[Δ, x_Symbol])[x_ .(x_ .A_)]:=(makeVector[x];2 trace[A])/;FreeQ[A,x]

(Subscript[Δ, b_Symbol])[b_ .gradient[u_][b_]]:=(makeVector[b];2 (Subscript[Δ, b])[u[b]]+b.gradient[(Subscript[Δ, b])[u[b]],b])

(Subscript[Δ, x_Symbol])[gradient[u_][x_].x_]:=(makeVector[x];2 (Subscript[Δ, x])[u[x]]+x.gradient[(Subscript[Δ, x])[u[x]],x])

(Subscript[Δ, x_Symbol])[f_[g_]]:=f''[g]  gradient[g,x] . gradient[g,x] + f'[g] (Subscript[Δ, x])[g]/;(makeVector[x];vector[g]=!=True)

Δ[norm]:=((dimension[#1]-1)/(norm[#1]))&

(Subscript[Δ, x_Symbol])[norm[A_ .x_]]:=(makeVector[x];((((norm[A.x])^2) ((hilbertSchmidt[A])^2)-((norm[A.x.A])^2))/((norm[A.x])^3)))/;FreeQ[A,x]

(Subscript[Δ, x_Symbol])[norm[x_ .Z_]]:=(makeVector[x];((((norm[x.Z])^2) ((hilbertSchmidt[Z])^2)-((norm[x.Z.Transpose[Z]])^2))/((norm[x.Z])^3)))/;FreeQ[Z,x]

(Subscript[Δ, x_Symbol])[Gamma[t_, f_]] :=(E^(-f))(f^(t-1))((norm[ gradient[f, x] ])^2)-(E^(-f))(f^(t-1))(Subscript[Δ, x])[f]  - (t-1)(E^(-f))(f^(t-2))((norm[ gradient[f, x] ])^2) /; FreeQ[t, x]

Δ[dimension]:=0&

(Subscript[Δ, x_Symbol,y_Symbol])[f_]:=(unmakeVector[y];(Subscript[Δ, x])[f]+partial[f,{y,2}])

((Subscript[Δ, x__])^(n_))[f_]:=((Subscript[Δ, x])^(n-1))[(Subscript[Δ, x])[f]]/;Head[n]==Integer&&Positive[n]

(Subscript[Δ, x_Symbol])[u_Symbol[((norm[x_])^(m_))x_]] := (makeVector[x];((norm[x])^(m-2))(((norm[x])^(2+m))Δ[u][((norm[x])^m)x] + m(m + dimension[x]) x.gradient[u][((norm[x])^m)x]+ m(m+2)((norm[x])^m)(jacobian[gradient[u]][((norm[x])^m)x].x).x))

laplacian = Δ

dimension

dimension[a_ x_]:=dimension[x]/;vector[x]||matrix[x]

dimension[a_?NumericQ x_]:=(If[matrix[x]=!=True,makeVector[x]];dimension[x])

dimension[a_. southPole+x_]:=dimension[x]

dimension[x_+w_]:=If[StringLength[ToString[dimension[x]]]<StringLength[ToString[dimension[w]]],dimension[x],dimension[w]]

dimension[v_?VectorQ]:=Length[v]

dimension[A_?MatrixQ]:={Length[A],Length[A[[1]]]}

IntegerQ[dimension[_]]^:=True

setDimension

setDimension::Matrix = "`1` will be considered to be a matrix of size `2`-by-`3`."

setDimension[x_Symbol,{m_,n_}]:=(dimension[x]={m,n};makeMatrix[x];If[Not[IntegerQ[m]],Unprotect[IntegerQ];IntegerQ[m]=True];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
			With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::Matrix, x, m, n]; SetOptions[$FrontEnd, messageOptions];])

setDimension::Vector = "`1` will be considered to be a vector in `2`-dimensional real Euclidean space."

setDimension2[x_Symbol,n_]:=(dimension[x]=n;makeVector[x];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];)

setDimension[x_Symbol,n_]:=(setDimension2[x, n];With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::Vector, x, n]; SetOptions[$FrontEnd, messageOptions];])

setDimension::VectorValuedFunction = "`1` will be considered to be a function taking values in `2`-dimensional Euclidean space."

setDimension[f_Symbol[_],n_]:=(dimension[f[_]]=n;makeVector[f[_]];If[Not[IntegerQ[n]],Unprotect[IntegerQ];IntegerQ[n]=True];
		With[{messageOptions = Options[$FrontEnd, MessageOptions][[1]]}, SetOptions[$FrontEnd, MessageOptions->{Global`WarningAction->{}}];Message[setDimension::VectorValuedFunction, f, n]; SetOptions[$FrontEnd, messageOptions];])/;f=!=List

homogeneous

homogeneous[f_,d_,x_Symbol,c_Symbol]:=(makeVector[x];makeVector[c];setDimension2[c,dimension[x]];homogeneous[f,d,Union[Cases[Level[f,-1],(Subscript[x, _])]],Union[Cases[Level[f,-1],(Subscript[x, _])]]/.x->c]/.Table[ (Subscript[c, j]), {j, Length[Union[Cases[Level[f,-1],(Subscript[x, _])]]]}]->c)/;!(MemberQ[Level[f,{-1}],x])&&Head[d]==Integer

homogeneous[f_,d_,x_Symbol,c_Symbol]:=(makeVector[x];makeVector[c];setDimension2[c,dimension[x]];homogeneous[f,d,x,Table[(Subscript[c, j]), {j, dimension[x]}]]/.Table[(Subscript[c, j]), {j, dimension[x]}]->c)/;MemberQ[Level[f,{-1}],x]&&Head[dimension[x]]==Integer&&Head[d]==Integer

homogeneous[f_,d_,x_Symbol,c_List]:=(makeVector[x];unmakeVector[c];Module[{z},homogeneous[f,d,Table[(Subscript[x, j]), {j,Length[c]}],c]/.(Subscript[x, j_])->(Subscript[z, j])/.x->c/.(Subscript[z, j_])->(Subscript[x, j])])/;Head[d]==Integer

homogeneous[f_,d_,v_List,c_Symbol]:=(makeVector[c];unmakeVector[v];setDimension2[c,Length[v]];homogeneous[f,d,v,Table[(Subscript[c, j]), {j, dimension[c]}]])/;Head[d]==Integer

homogeneous[f_,d_,v_List,c_List]:=(unmakeVector[v];unmakeVector[c];If[togetherness2&&!(norm[c]===0),turnOff[togetherness]];With[{mI=multiIndices[d,Length[v]]},Plus@@(((partial@@Prepend[Transpose[{v,mI[[#1]]}],f]/.Thread[v->c]) Times@@(((v-c))^(mI[[#1]])))/multiFactorial[mI[[#1]]]&)/@Range[Length[mI]]])/;Head[d]==Integer

homogeneous[f_,d_,v_List]:=homogeneous[f,d,v,Table[0,{Length[v]}]]/;Head[d]==Integer

homogeneous[f_,d_,x_Symbol]:=(homogeneous[f,d,Union[Cases[Level[f,-1],(Subscript[x, _])]]]/.Table[0,{Length[Union[Cases[Level[f,-1],(Subscript[x, _])]]]}]->0)/;!(MemberQ[Level[f,{-1}],x])&&Head[d]==Integer

homogeneous[f_,d_,x_Symbol]:=(homogeneous[f,d,x,Table[0,{dimension[x]}]]/.Table[0,{dimension[x]}]->0)/;MemberQ[Level[f,{-1}],x]&&Head[dimension[x]]==Integer&&Head[d]==Integer

multiIndices and multiFactorial

multiFactorial[a_List]:=Times@@(a!)

multiIndices[k_,1]:={{k}}

multiIndices[0,n_]:={Table[0,{n}]}

multiIndices[k_,n_]:=multiIndices[k,n]=Join@@Table[(Prepend[#1,k-j]&)/@multiIndices[j,n-1],{j,0,k}]

multiIndices[k_,1, {t_}]:=If[Mod[k-t,2]==0,{{k}},{}]

multiIndices[0,n_, tt_]:=If[And @@ EvenQ[tt],{Table[0,{n}]},{}]

multiIndices[k_,n_, tt_]:=multiIndices[k,n, tt]=Join@@Table[(Prepend[#1,k-j]&)/@multiIndices[j,n-1, Delete[tt,1]],{j,Mod[k-tt[[1]],2],k, 2}]

taylor

taylor[f_,d_,vc___]:=Plus@@(homogeneous[f,#1,vc]&)/@Range[0,d]/;Head[d]==Integer

volume

volume[n_Integer]:=((π^(n/2))/(Gamma[(n/2)+1]))/;Positive[n]

surfaceArea

surfaceArea[n_Integer]:=((2 (π^(n/2)))/(Gamma[(n/2)]))/;Positive[n]

integrateBall

integrateBall[f_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateBall[f/.Thread[v->Table[(Subscript[z, j]), {j, Length[v]}]],z]])

integrateBall[f_,x_Symbol]:=(makeVector[x];
		f volume[dimension[x]])/;FreeQ[f,x]

integrateBall[f_,x_Symbol]:=(∫_(-1)^1)(f/.{((norm[x])^2)->(x^2),norm[x]->Abs[x],(Subscript[x, 1])->x})\[DifferentialD]x/;dimension[x]==1

integrateBall[f_,x_Symbol]:=(makeVector[x];Module[{r},If[Head[dimension[x]]===Integer,Null,Unprotect[Limit]; Limit[c_. (r^(dimension[x]+m_.)),r->0]:=0/;(m>0||m==0)&&FreeQ[c,r]];dimension[x] volume[dimension[x]] Integrate[Expand[(r^(dimension[x]-1)) integrateSphere[f/.{norm[x]->r,(Subscript[x, j_])->r (Subscript[x, j]),x->r x},x]],{r,0,1}, Assumptions->dimension[x]≥1]])

integrateSphere

integrateSphere[f_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateSphere[f/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],z]])

integrateSphere[f_,x_Symbol]:=(makeVector[x];integrateSphere1[Expand[f/.{norm[x]->1,x.z_:>Table[(Subscript[x, j]),{j, dimension[x]}].Table[(Subscript[z, j]),{j, dimension[x]}],z_ .x:>Table[(Subscript[z, j]),{j, dimension[x]}].Table[(Subscript[x, j]),{j, dimension[x]}]}],x])

integrateSphere1[HoldPattern[+a__],x_]:=Plus@@Thread[integrateSphere1[{a},x]]

integrateSphere1[f_,x_]:=f/;FreeQ[f,x]

integrateSphere1[c_. (Subscript[x_, _]),x_]:=0

integrateSphere1[c_ f_,x_]:=c integrateSphere1[f,x]/;FreeQ[c,x]

integrateSphere1[((Subscript[x_, _])^(j_)),x_]:=integrateSphere3[{j},dimension[x]]

integrateSphere1[f_,x_]:=integrateSphere2[f,x]/;monomialSpecialQ[f,x]

integrateSphere2[f_,x_]:=integrateSphere3[List@@f/.((Subscript[x, _])^(b_))->b,dimension[x]]

integrateSphere3[b_,n_]:=0/;Or@@OddQ[b]

integrateSphere3[b_,n_]:=((Times@@((b-1)!!))/(Times@@Range[n,n+Plus@@b-2,2]))

monomialSpecialQ[f_,x_Symbol]:=Head[f]===Times&&And@@(MatchQ[f[[#1]],((Subscript[x, _])^(_))]&)/@Range[Length[f]]

integrateEllipsoidArea

integrateEllipsoidArea[f_,b__,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateEllipsoidArea[f/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],b,z]])

integrateEllipsoidArea[f_,b_,c_,d_,x_Symbol]:= integrateEllipsoidArea[ f /. Table[ (Subscript[x, j])->(Subscript[x, j]) -c[[j]]/(2b[[j]]),{j,dimension[x]}], b, d - Sum[ ((c[[j]])^2)/b[[j]], {j,dimension[x]}]/4, x]

integrateEllipsoidArea[f_,b_,c_List,x_Symbol]:= integrateEllipsoidArea[f, b, c, -1, x]

integrateEllipsoidArea[f_,b_,d_,x_Symbol]:= -d integrateEllipsoidArea[f, -b/d, x]

integrateEllipsoidArea[f_,b_,x_Symbol]:=(makeVector[x];integrateEllipsoid1[Expand[f/.{norm[Sqrt[b] Table[(Subscript[x, j]),{j, dimension[x]}]]->1,x.z_:>Table[(Subscript[x, j]),{j, dimension[x]}].Table[(Subscript[z, j]),{j, dimension[x]}],z_ .x:>Table[(Subscript[z, j]),{j, dimension[x]}].Table[(Subscript[x, j]),{j, dimension[x]}],((norm[x])^(k_))-> ((Sum[((Subscript[x, j])^2),{j,1, dimension[x]}])^(k/2))}],b,x])

integrateEllipsoid1[HoldPattern[+a__],b_,x_]:=Plus@@Table[ integrateEllipsoid1[{a}[[j]],b, x],{j, 1,Length[{a}]}]

integrateEllipsoid1[f_, b_,x_]:= f dimension[x] volume[dimension[x]]/(2 (Sqrt[Product[b[[j]],{j, dimension[x]}]]))/;FreeQ[f,x]

integrateEllipsoid1[f_,b_,x_]:=0/;Or@@OddQ[Table[Exponent[f,(Subscript[x, j ])],{j,dimension[x]}]]

integrateEllipsoid1[c_ f_,b_,x_]:=c integrateEllipsoid1[f,b, x]/;FreeQ[c,x]

integrateEllipsoid1[f_,b_,x_]:=dimension[x] volume[dimension[x]] Product[ (Exponent[f,(Subscript[x, j])] - 1)!!, {j, dimension[x]}]/(2 (Sqrt[Product[ ((b[[j]])^(Exponent[f,(Subscript[x, j])]+1)), {j, dimension[x]}]])Product[ k, {k, dimension[x], dimension[x] + Sum[Exponent[f,(Subscript[x, j])],{j,dimension[x]}]-2, 2}])

integrateEllipsoidVolume

integrateEllipsoidVolume[f_,b__,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];integrateEllipsoidVolume[f/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],b,z]])

integrateEllipsoidVolume[f_,b_,c_,d_,x_Symbol]:= integrateEllipsoidVolume[ f /. Table[ (Subscript[x, j])->(Subscript[x, j]) -c[[j]]/(2b[[j]]),{j,dimension[x]}], b, d - Sum[ ((c[[j]])^2)/b[[j]], {j,dimension[x]}]/4, x]

integrateEllipsoidVolume[f_,b_,c_List,x_Symbol]:= integrateEllipsoidVolume[f, b, c, -1, x]

integrateEllipsoidVolume[f_,b_,d_,x_Symbol]:= integrateEllipsoidVolume[f, -b/d, x]

integrateEllipsoidVolume[f_,b_,x_Symbol]:=(makeVector[x];integrateEllipsoid1V[Expand[f/.{x.z_:>Table[(Subscript[x, j]),{j, dimension[x]}].Table[(Subscript[z, j]),{j, dimension[x]}],z_ .x:>Table[(Subscript[z, j]),{j, dimension[x]}].Table[(Subscript[x, j]),{j, dimension[x]}],((norm[x])^(k_))-> ((Sum[((Subscript[x, j])^2),{j,1, dimension[x]}])^(k/2))}],b,x])

integrateEllipsoid1V[f_, b_, x_] := integrateBall[f /. Table[(Subscript[x, j])->(Subscript[x, j])/(Sqrt[b[[j]]]), {j, dimension[x]}],x]/(Sqrt[Times @@ b])

Boundary Value Problems

dirichlet

Options[dirichlet]={region->Sphere, Δ->0}

dirichlet[p_, x_Symbol,options__] := With[ {t = antiLaplacian[(Δ /. {options} /. Options[dirichlet]),x, singularity -> 0]},Factor[t+dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, Δ->_]]]]/; And[(region /. {options} /. Options[dirichlet]) === exteriorSphere, MemberQ[{options}, Δ-> _]]

dirichlet[p_, x_Symbol,options__] := With[ {t = antiLaplacian[(Δ /. {options} /. Options[dirichlet]),x]}, t+dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, Δ->_]]] /; MemberQ[{options}, Δ-> _]

dirichlet[p_, x_Symbol,options___] := With[ {region2 = (region /. {options} /. Options[dirichlet])}, Which[region2 === Sphere, dirichlet1[p,x],
Head[region2]===quadratic, p + antiLaplacian[-(Subscript[Δ, x])[p], x, multiple-> region2],
region2 === exteriorSphere, exteriorDirichlet[p,x],
Head[region2] === annulus, If[Head[p]===List, annularDirichlet[p[[1]],p[[2]],region2[[1]], region2[[2]],x], annularDirichlet[p,p,region2[[1]], region2[[2]],x]]]]

dirichlet1[p_,x_]:=Plus@@(#1[[1]]&)/@harmonicDecomposition[p/.norm[x]->1,x] /; PolynomialQ[p /. norm[x]->1, Cases[p, (Subscript[x, _]),∞]]

dirichlet[f_,v_List, opt___]:=explicit[dirichlet,f,v, opt]

exteriorDirichlet[f_,z_]:=Factor[kelvin[dirichlet[f,z],z]]

annular[0,r_,s_,z_]=0

annular[f_,r_,s_,z_]:=Plus@@(If[#1[[2]]===#1[[3]],(((Log[norm[z]]-Log[s]) (r^(#1[[2]])) #1[[1]])/(Log[r]-Log[s])),(((((norm[z])^(2 #1[[2]]-2 #1[[3]]))-(s^(2 #1[[2]]-2 #1[[3]]))) (r^(#1[[2]])) #1[[1]])/((r^(2 #1[[2]]-2 #1[[3]]))-(s^(2 #1[[2]]-2 #1[[3]]))))]&)/@Flatten[(append[harmonicDecomposition[#1[[1]],z],#1[[2]]]&)/@polynomialDecomposition[f,z],1]/;dimension[z]==2

annular[f_,r_,s_,z_]:=Plus@@(((((norm[z])^(2-dimension[z]+2 #1[[2]]-2 #1[[3]]))-(s^(2-dimension[z]+2 #1[[2]]-2 #1[[3]]))) (r^(#1[[2]])) #1[[1]])/((r^(2-dimension[z]+2 #1[[2]]-2 #1[[3]]))-(s^(2-dimension[z]+2 #1[[2]]-2 #1[[3]])))&)/@Flatten[(append[harmonicDecomposition[#1[[1]],z],#1[[2]]]&)/@polynomialDecomposition[f,z],1]

append[b_,d_]:=(Append[#1,d]&)/@b

annularDirichlet[f_,g_,r_,s_,z_Symbol]:=annular[f,r,s,z]+annular[g,s,r,z]

harmonicDecomposition

d[0,0]=1

lap[f_,v_]:=Plus@@((∂_({#1,2}))f&)/@v

harmonicDecomposition[0,x_]:={}

harmonicDecomposition[p_,x_Symbol]:=(makeVector[x];With[{deg=degree[p/.norm[x]->normtemp,x],n=dimension[x]},powers={};hd[_]:=0;With[{partition=Partition[Flatten[Expand[((m=#1[[2]];k=Floor[(m/2)];If[m>1,d[1,1]=(1/(2 (2 m+n-4)))];v=Select[Variables[#1[[1]]],!(FreeQ[#1,x])&];laplac[0]=#1[[1]];laplac[1]=lap[#1[[1]],v];Do[d[j,j]=((d[j-1,j-1] (2 m+n-2 j))/(2 j (2 m+n-4 j+2) (2 m+n-4 j)));laplac[j]=lap[laplac[j-1],v],{j,2,k}];Do[d[i,j]=((d[i,j-1])/(2 (i-j) (2 (m-i-j-1)+n))),{i,0,k-1},{j,i+1,k}];Table[{(∑_(j=i)^k)d[i,j] ((norm[x])^(2 (j-i))) laplac[j],2 i},{i,0,k}])&)/@polynomialDecomposition[p/.norm[x]->normtemp,x]]//.{a_.+(normtemp^(e_.)) b_.,j_}->{{a,j},{b,j+e}}],2]},powers=powers⋃Transpose[partition][[2]];((hd[#1[[2]]]=hd[#1[[2]]]+#1[[1]])&)/@partition];aa=Select[({hd[#1],#1}&)/@powers,#1[[1]]=!=0&];Clear[hd];aa])

harmonicDecomposition[f_,v_List]:=explicit[harmonicDecomposition,f,v]

antiLaplacian

Options[antiLaplacian]={multiple->None,singularity->None}

antiLaplacian[0,x_,options___?OptionQ]:=0

antiLaplacian[0.,x_,options___?OptionQ]:=0.

antiLaplacian[p_,x_Symbol,options___?OptionQ]:=(makeVector[x];
If[Head[multiple/.{options}/.Options[antiLaplacian]]===quadratic,
antiLaplacianQ[p, x, multiple /. {options}],
If[(multiple/.{options}/.Options[antiLaplacian])=== (norm^2)||!(FreeQ[p,norm[x]]),antiLaplacianN[Expand[p],x,options],antiLaplacian1[Expand[p],x]]])

antiLaplacianQ[p_, x_, q_ ] := 
Module[{ph = blockPolynomialDecomposition[p, x]},
Which[ Length[q] ==1, antiLaplacianQ[p,x,quadratic[q[[1]],-1]],
Length[q]==2 , If[Head[q[[2]]]===List,antiLaplacianQ[p,x,quadratic[q[[1]],q[[2]],-1]],(q[[1]].((Array[x,Length[q[[1]]]])^2)+  q[[2]])Sum[antiLaplacianQ4[ph[[j]], x, List @@ q], {j,1, Length[ph]}]],
Length[q]==3,(q[[1]].((Array[x,Length[q[[1]]]])^2)+ q[[2]].Array[x,Length[q[[2]]]]+ q[[3]])Sum[antiLaplacianQ3[ph[[j]], x, List @@ q], {j,1, Length[ph]}]]]

antiLaplacianQ4[p_, x_, q_ ] := antiLaplacianQ2[p, x, q[[1]]] /; q[[2]] ===0

antiLaplacianQ4[p_, x_, q_ ] := Module[{f},With[{m = p[[2]], q2 = q[[1]], q0 = q[[2]], p3= p[[3]]},
f[m+1] = 0;
f[m] = antiLaplacianQ2[p, x, q2];
Do[f[m-j] =  antiLaplacianQ2[{-(Subscript[Δ, x])[q0 f[m-j+2]],m-j,p3},x,q2],{j,1,m} ];
Sum[f[j], {j,0, m}]
]]

antiLaplacianQ3[p_, x_, q_ ] := antiLaplacianQ4[p, x, {q[[1]],q[[3]]}] /;  q[[2]] === Array[0,dimension[x]]

antiLaplacianQ3[p_, x_, q_ ] := Module[{f},With[{m = p[[2]], q2 = q[[1]], q1 = q[[2]].Array[x,Length[q[[2]]]], q0 = q[[3]]},
f[m+1] = 0;
f[m] = antiLaplacianQ2[p, x, q2];
Do[f[m-j] =  antiLaplacianQ2a[-(Subscript[Δ, x])[q1 f[m-j+1] + q0 f[m-j+2]],x,q2],{j,1,m} ];
Sum[f[j], {j,0, m}]
]]

antiLaplacianQ2a[0, x_, q_] := 0

antiLaplacianQ2a[p_, x_, q_ ] := 
With[{ph = blockPolynomialDecomposition[p, x]},Sum[antiLaplacianQ2[ph[[j]], x, q], {j,1, Length[ph]}]
]

antiLaplacianQ2[p_, x_, q_] :=  p[[1]]/(2 Plus @@ q) /; p[[2]] == 0

antiLaplacianQ2[p_, x_, q_] :=  With[{p1=p[[1]], twoA = 2(Plus @@ q), n = dimension[x], multiIndicesmn =multiIndices[p[[2]],dimension[x], p[[3]]]},Module[{d,s},((Plus @@ Map[(Subscript[d, #]) (Times @@ ((Array[x,n])^(#)))/multiFactorial[#]&, multiIndicesmn]) /. Solve[Map[(Expand[(twoA + 4  Plus @@ (# q)) (Subscript[d, #]) + ((Plus @@ (# (# - 1) q Table[(Subscript[d, ReplacePart[#, #[[i]]-2,i]]),{i, n}]))/. (Subscript[d, s_ ]):> Sum[(Subscript[d, ReplacePart[s, s[[j]]+2,j]]),{j,n}])] == multiFactorial[#] Coefficient[p1,Times @@ ((Array[x,n])^(#))])&, multiIndicesmn], Map[(Subscript[d, #])&, multiIndicesmn]])[[1]]]]

antiLaplacianN[HoldPattern[+a__],z_,options___?OptionQ]:=Plus@@Thread[Unevaluated[antiLaplacianN[{a},z,options]]]

antiLaplacianN[c_ g_,z_,options___?OptionQ]:=c antiLaplacianN[g,z,options]/;FreeQ[c,z]

antiLaplacianN[g_,z_,options___?OptionQ]:=antiLap1[If[Head[g]===Times,Select[g,FreeQ[#1,norm[z]]&],If[FreeQ[g,norm[z]],g,1]],If[Head[g]===Times,Select[g,!(FreeQ[#1,norm[z]])&],If[FreeQ[g,norm[z]],1,g]],z,options]

antiLap1[p_,f_,z_,options___?OptionQ]:=Plus@@(antiLap2[#1[[1]],#1[[2]],f,z,options]&)/@harmonicDecomposition[p,z]

antiLap2[p_,j_,f_,z_,options___?OptionQ]:=Plus@@(antiLap3[#1[[1]],#1[[2]],j,f,z,options]&)/@polynomialDecomposition[p,z]

antiLap3[p_,m_,j_,f_,z_,options___?OptionQ]:=With[{n=dimension[z]},p If[(singularity/.{options}/.Options[antiLaplacian])=!=None,∫(t^(1-2 m-n)) (∫((t^(-1+2 m+n+j)) f/.norm[z]->t)\[DifferentialD]t)\[DifferentialD]t,With[{int=(∫_0^t)((s^(-1+2 m+n+j)) f/.norm[z]->s)\[DifferentialD]s},∫(t^(1-2 m-n)) If[FreeQ[int,Indeterminate]&&FreeQ[int,∞]&&FreeQ[int,ConditionalExpression] && FreeQ[int,Limit]&&FreeQ[int,If]&&FreeQ[int,Integrate]&&(singularity/.{options}/.Options[antiLaplacian])===None,int,∫((t^(-1+2 m+n+j)) f/.norm[z]->t)\[DifferentialD]t]\[DifferentialD]t]]/.t->norm[z]]

antiLaplacian[f_,v_List,options___?OptionQ]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];antiLaplacian[f/.((norm[v])^2)->((norm[z])^2)/.Thread[v->Table[(Subscript[z, j]), {j, Length[v]}]],z,options]/.Prepend[Thread[Table[(Subscript[z, j]), {j, Length[v]}]->v],norm[z]->norm[v]]])

antiLaplacian1[HoldPattern[+a__],z_]:=Plus@@Thread[antiLaplacian1[{a},z]]

antiLaplacian1[g_,z_]:=(1/2) g ((Subscript[z, 1])^2)/;FreeQ[g,z]

antiLaplacian1[g_,z_]:=antiLaplacian3[Transpose[Sort[({Exponent[g,(Subscript[z, #1])],#1}&)/@(#1[[2]]&)/@Variables[g]]],g,z]/;monomialQ[g,z]

antiLaplacian1[c_ g_,z_]:=c antiLaplacian1[g,z]/;FreeQ[c,z]

antiLaplacian3[{e_,p_},g_,z_]:=(((Subscript[z, p[[-1]]])^(e[[-1]]+2))/((e[[-1]]+2) (e[[-1]]+1)))/;Length[e]==1

antiLaplacian3[{e_,p_},g_,z_]:=((g ((Subscript[z, Last[p]])^2))/((Last[e]+1) (Last[e]+2)))-Plus@@((((e[[#1]]-1) e[[#1]] antiLaplacian1[((((Subscript[z, Last[p]])^2) g)/((Subscript[z, p[[#1]]])^2)),z])/((Last[e]+1) (Last[e]+2)))&)/@Range[Position[Positive[e],True][[1,1]],Length[e]-1]/;Length[e]>1

neumann

neumann[f_,v_List]:=explicit[neumann,f,v]

neumann[f_,g_,v_List]:=explicit[neumann,f,g,v]

neumann[0,x_Symbol]:=0

neumann[f_,x_Symbol]:=Plus@@(If[#1[[2]]===0,0,((#1[[1]])/(#1[[2]]))]&)/@polynomialDecomposition[dirichlet[f,x],x]

neumann[f_,g_,x_Symbol]:=With[{v=antiLaplacian[g,x]},v+neumann[f-normalD[v,x],x]-Limit[v/.{(Subscript[x, _])->0,norm[x]->t},t->0,Direction->-1]]

antiLaplacianQNa[0,x_Symbol, q_ ] := 0

antiLaplacianQNa[p_, x_Symbol, q_ ] := 
With[{ph = blockPolynomialDecomposition[p, x]},Sum[antiLaplacianQN[ph[[j]], x, q], {j,1, Length[ph]}]
]

makePD[l_,x_] :=
MapIndexed[{(Subscript[x, First[#2]]),#1}&,l]

antiLaplacianQN[p_, x_Symbol, q_] :=  With[{p1=p[[1]], n = dimension[x], multiIndicesmn =multiIndices[p[[2]],dimension[x], p[[3]]]},Module[{c,d,r}, 
c := (q / (2(q.#)+4q)(Subscript[d, #]))& /@ multiIndicesmn;
r:=Laplacian[Table[Times @@ ((Array[x,n])^(ReplacePart[#,#[[i]]+2,i])),{i,n}]& /@ multiIndicesmn,Array[x,dimension[x]]];
((Plus @@ Map[(Subscript[d, #]) (Times @@ ((Array[x,n])^(#)))&, multiIndicesmn]) /.
Solve[Map[  Plus @@ Flatten[c  Apply[D,Prepend[ makePD[#,x],r]]]
==(*multiFactorial[#] Coefficient[p1,Times @@ ((Array[x,n])^(#))]&,multiIndicesmn]*)Apply[D,Prepend[ makePD[#,x],p1]]&,multiIndicesmn],Map[(Subscript[d, #])&, multiIndicesmn]
])[[1]]]]

neumannInvMap[p_,x_Symbol,q_] :=
Plus @@(If[PossibleZeroQ[2q.#[[1]]],#[[2]],(#[[2]]/(2q.#[[1]]))](Times @@ ((Array[x,dimension[x]])^(#[[1]])))& /@CoefficientRules[p,Array[x,dimension[x]]])

neumann2[p_,x_Symbol,q_] :=
With[{ph = Association[#[[2]]->#[[1]]&/@polynomialDecomposition[p,x]],n=dimension[x], qx = q.Array[x,dimension[x]]^2-1},
Module[{f,m,pp},
m = Max[Keys[ph]]-2;
f[m+2] = 0;
f[m+1] = 0;
f[-1] = 0;
Do[
pp = Simplify[(Subscript[Δ, x])[neumannInvMap[Lookup[ph,m+2-j,0]+f[m+2-j],x,q]]];
f[m-j] = antiLaplacianQNa[pp,x,q];
,{j,0,m} 
];
Sum[f[j], {j,0, m}]
]]

neumann1[p_, x_, q_] := neumannInvMap[p,x,q] - neumannInvMap[(q.((Array[x, dimension[x]])^2) -1)neumann2[p,x,q],x,q]

neumann[p_, x_Symbol, options__] := Module[{quad =(region /. options), neu}, 
Which[
Length[quad] == 1,neumann1[p, x, quad[[1]]],
Length[quad] ==2,If[Head[quad[[2]]]===List,neumann[p, x, region-> quadratic[quad[[1]], quad[[2]], -1]],neumann1[p, x, -quad[[1]]/quad[[2]]]/(-quad[[2]])],
Length[quad] == 3,  
neu = neumann[p /. Table[ (Subscript[x, j])->(Subscript[x, j]) -quad[[2]][[j]]/(2quad[[1]][[j]]),{j,dimension[x]}], x, region->quadratic[quad[[1]],Simplify[quad[[3]]-Sum[((quad[[2]][[j]])^2)/quad[[1]][[j]],{j,dimension[x]}]/4]]] /. Table[ (Subscript[x, j])->(Subscript[x, j]) +quad[[2]][[j]]/(2quad[[1]][[j]]),{j,dimension[x]}]; neu -( neu/. Table[(Subscript[x, j])-> 0, {j, dimension[x]}])
]]

neumann[f_, g_, x_Symbol, options__] := Module[{v=antiLaplacian[g,x],quad =(region /. options),q},
Which[Length[quad] == 1, q = quad[[1]].((Array[x, dimension[x]])^2),
Length[quad] == 2, If[Head[quad[[2]]]===List, q = quad[[1]].((Array[x, dimension[x]])^2)+ quad[[2]].Array[x, dimension[x]],  q = quad[[1]].((Array[x, dimension[x]])^2)],
Length[quad]==3, q = quad[[1]].((Array[x, dimension[x]])^2)+ quad[[2]].Array[x, dimension[x]]
];
v+neumann[f-gradient[v,x].gradient[q,x],x, options]]

exteriorNeumann

exteriorNeumann[f_,v_List]:=explicit[exteriorNeumann,f,v]

exteriorNeumann[f_,g_,v_List]:=explicit[exteriorNeumann,f,g,v]

exteriorNeumann[0,x_Symbol]:=0

exteriorNeumann[f_,x_Symbol]:=Plus@@(((((norm[x])^(2-dimension[x]-2 #1[[2]])) #1[[1]])/(dimension[x]+#1[[2]]-2))&)/@polynomialDecomposition[dirichlet[f,x],x]/;dimension[x]=!=2

exteriorNeumann[f_,x_Symbol]:=Plus@@(((norm[x])^(-2 #1[[2]])) #1[[1]] If[#1[[2]]===0,0,(1/(#1[[2]]))]&)/@polynomialDecomposition[dirichlet[f,x],x]/;dimension[x]===2

exteriorNeumann[f_,g_,x_Symbol]:=With[{v=antiLaplacian[g,x,singularity->0]},v+exteriorNeumann[f+normalD[v,x],x]]

biDirichlet

biDirichlet[f_,v_List]:=explicit[biDirichlet,f,v]

biDirichlet[f_,g_,v_List]:=explicit[biDirichlet,f,g,v]

biDirichlet[f_,z_Symbol]:=Together[With[{p=dirichlet[f,z]},p+(1/2) (1-((norm[z])^2)) z.gradient[p,z]]]

biDirichlet[f_,g_,z_Symbol]:=normalIntegral[dirichlet[g,z],z]+biDirichlet[f,z]

normalIntegral[p_,z_]:=(1/2) (((norm[z])^2)-1) p

Spherical Harmonics

zonalHarmonic

zonalHarmonic[0,x_,z_]:=1

zonalHarmonic[m_Integer,x_Symbol,z_Symbol]:=(makeVector[x];makeVector[z];Plus@@((1/(#1! (m-2 #1)! (2^(#1))))((((-1))^(#1)) (dimension[x+z]+2 m-2) Times@@Range[dimension[x+z],dimension[x+z]+2 m-2 #1-4,2] (((x.z))^(m-2 #1)) ((norm[x])^(2 #1)) ((norm[z])^(2 #1)))&)/@Range[0,(m/2)])

zonalHarmonic[m_Integer,x_List,z_List]:=(unmakeVector[x];unmakeVector[z];Expand[Plus@@(((((-1))^(#1)) (Length[x]+2 m-2) Times@@Range[Length[x],Length[x]+2 m-2 #1-4,2] (((x.z))^(m-2 #1)) ((norm[x])^(2 #1)) ((norm[z])^(2 #1)))/(#1! (m-2 #1)! (2^(#1)))&)/@Range[0,(m/2)]])

dimHarmonic

dimHarmonic[0,2]=1;

dimHarmonic[m_,n_]:=Binomial[n+m-1,n-1]-Binomial[n+m-3,n-1]

basisH

ball[f_,g_, x_]:= integrateBall[ f g, x]

sphere[f_,g_, x_]:= integrateSphere[ f g, x]

basisH[0,x_]:={1}

basisH[m_Integer,x_Symbol]:=Together[Module[{xx},With[{ttt=Expand[Expand[(kelvin[#1,xx]&)/@(partial@@Prepend[Transpose[{{(Subscript[xx, 1]),(Subscript[xx, 2])},#1}],Log[norm[xx]]]&)/@Join[(Prepend[#1,0]&)/@multiIndices[m,1],(Prepend[#1,1]&)/@multiIndices[m-1,1]]]]/.xx->x},noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]]]]]/;dimension[x]==2

basisH[m_Integer,x_Symbol]:=Together[Module[{xx},With[{ttt=Expand[Expand[(kelvin[#1,xx]&)/@(partial@@Prepend[Transpose[{Table[ (Subscript[xx, j]), {j, dimension[x]}],#1}],norm[xx]^(2-dimension[x])]&)/@Join[(Prepend[#1,0]&)/@multiIndices[m,dimension[x]-1],(Prepend[#1,1]&)/@multiIndices[m-1,dimension[x]-1]]]]/.xx->x},noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]]]]]/;Head[dimension[x]]==Integer

basisH[m_Integer,x_Symbol, ip_]:= (inprod = ip; If [ip === Ball, inprod = ball]; If[ip=== Sphere, inprod = sphere];orthonormalize[ basisH[m, x], inprod[#1, #2, x] & ])/; Head[dimension[x]] == Integer

basisH[m_Integer,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];Expand[basisH[m,z]/.Thread[Table[ (Subscript[z, j]), {j,Length[v]}]->v]]/.norm[z]->norm[v]])

basisH[m_Integer,v_List, ip_]:= (inprod = ip; If [ip === Ball, inprod = ball]; If[ip=== Sphere, inprod = sphere];orthonormalize[ basisH[m, v], inprod[#1, #2, v] & ])

noLeadingMinus[x_]:=If[Head[x]===Plus,If[Head[x[[1]]]===Times,If[x[[1,1]]<0,Expand[-x],x,x],If[x[[1]]<0,Expand[-x],x,x]],x]

orthonormalize

orthonormalize[b_List,innerProduct_]:= Module[{e,j,t},Do[e[j]=b[[j]]-Plus@@(innerProduct[b[[j]],e[#1]] e[#1]/t[#1]&)/@Range[j-1];t[j] = innerProduct[e[j],e[j]],{j,1,Length[b]}];Together[((e[#1])/(Sqrt[t[#1]]))]&/@Range[Length[b]]]

The Kelvin Transform

reflection

reflection[∞]:=0

reflection[x_]:=(If[Head[x]==Symbol,makeVector[x]];If[norm[x]===0,∞,x/((norm[x])^2)])

reflection[x_List, Sphere[c_, r_]]:= reflection[x, Sphere[Table[c.(Subscript[δ, j]),{j,1,Length[x]}],r]] /; Head[c] =!= List

reflection[x_, Sphere[c_List, r_]]:= reflection[Table[x.(Subscript[δ, j]),{j,1,Length[c]}], Sphere[c, r]] /; Head[x] =!= List

reflection[x_, Sphere[c_, r_]]:= c + (r^2)reflection[x-c]

reflection[x_, Hyperplane[b_, c_]]:= x - ((2((x.b)-c)b)/((norm[b])^2))

reflection[x_List, Hyperplane[b_, c_]]:= reflection[x, Hyperplane[Table[b.(Subscript[δ, j]),{j,1,Length[x]}],c]] /; Head[b] =!= List

reflection[x_, Hyperplane[b_List,c_]]:= reflection[Table[x.(Subscript[δ, j]),{j,1,Length[b]}], Hyperplane[b, c]] /; Head[x] =!= List

kelvin

kelvin[HoldPattern[+u__],x_]:=Plus@@Thread[kelvin[{u},x]]

kelvin[u_,x_List]:=(unmakeVector[x];PowerExpand[Together[((norm[x])^(2-Length[x])) (u/.Thread[x->reflection[x]])],Plus @@ (x^2)])

PowerExpand and Together are included above so that powers of norm[x] get properly collected and simplified.

kelvin[u_,x_Symbol]:=(makeVector[x];((norm[x])^(2-dimension[x])) (u/.{x->reflection[x],(Subscript[x, j_]):>((Subscript[x, j])/((norm[x])^2))}))

Φ

Φ[z_List]:=(unmakeVector[z];If[Last[z]===-1&&norm[Drop[z,-1]]===0,∞,2 reflection[z+Append[Table[0,{Length[z]-1}],1]]-Append[Table[0,{Length[z]-1}],1]])

Φ[z_]:=(If[!(z===0||z===∞||MatchQ[z,_. southPole]),makeVector[z];If[dimension[southPole]=!= dimension[z],setDimension[southPole,dimension[z]]]];2 reflection[z-southPole]+southPole)

Φ[0,-1]=∞;

Φ[x_,y_]:= (makeVector[x];unmakeVector[y];{((2 x)/(((norm[x])^2)+(((y+1))^2))),((2 (y+1))/(((norm[x])^2)+(((y+1))^2)))-1})

southPole

norm[southPole]^=1;

southPole.z_Symbol:=-(Subscript[z, dimension[z]])/;matrix[z]=!=True

a_Symbol.southPole:=-(Subscript[a, dimension[a]])/;matrix[z]=!=True

makeVector[southPole];

N[(Subscript[southPole, j_Real])] := (Subscript[southPole, Round[j]])/; Round[j] == j;

(Subscript[southPole, j_]):=-1/;j==dimension[southPole]

(Subscript[southPole, j_]):=0/;j<dimension[southPole]

kelvinH

kelvinH[HoldPattern[+u__],x_]:=Plus@@Thread[kelvinH[{u},x]]

kelvinH[u_,x_List]:=(unmakeVector[x];(2^((1/2) (Length[x]-2))) (((((norm[x])^2)+2 Last[x]+1))^(2-Length[x])) (u/.Thread[x->Φ[x]]))

kelvinH[u_,z_Symbol]:=(makeVector[z];(2^((1/2) (dimension[z]-2))) ((norm[z-southPole])^(2-dimension[z])) (u/.{z->Φ[z],(Subscript[z, j_]):>Φ[z].(Subscript[δ, j])}))

kelvinH[u_,x_Symbol,y_Symbol]:=(makeVector[x];unmakeVector[y];(2^((1/2) (dimension[x]-1))) (((((norm[x])^2)+(((y+1))^2)))^((1/2) (1-dimension[x]))) (u/.{x->((2 x)/(((norm[x])^2)+(((y+1))^2))),y->((2 (y+1))/(((norm[x])^2)+(((y+1))^2)))-1,(Subscript[x, j_]):>((2 (Subscript[x, j]))/(((norm[x])^2)+(((y+1))^2)))}))

Kernels

poissonKernel

poissonKernel[x_,z_]:=(makeVector[x];makeVector[z];((1-((norm[x])^2) ((norm[z])^2))/(((1-2 x.z+((norm[x])^2) ((norm[z])^2)))^((1/2) dimension[x+z]))))

poissonKernelH

poissonKernelH[z_,w_]:=(makeVector[z];makeVector[w];(2 (z.(Subscript[δ, dimension[z]])+w.(Subscript[δ, dimension[z]])))/((dimension[z] volume[dimension[z]]) (((norm[z])^2)+((norm[w])^2)-2 z.w+4 z.(Subscript[δ, dimension[z]]) w.(Subscript[δ, dimension[z]]))^(((dimension[z])/2))))

poissonKernelH[x_,y_,t_,u_]:=(makeVector[x];makeVector[t];unmakeVector[y];unmakeVector[u];(2 (y+u))/(((dimension[x]+1) volume[dimension[x]+1]) (((((norm[x-t])^2)+(((y+u))^2)))^((1/2) (dimension[x]+1)))))

bergmanKernel

bergmanKernel[x_,y_]:=(makeVector[x];makeVector[y];With[{n=dimension[x+y]},((n+(8 x.y-2 n-4) ((norm[x])^2) ((norm[y])^2)+(n-4) ((norm[x])^4) ((norm[y])^4))/(n volume[n] (((1-2 x.y+((norm[x])^2) ((norm[y])^2)))^(1+(n/2)))))])

bergmanKernelH

bergmanKernelH[z_,w_]:=(makeVector[z];makeVector[w];With[{n=dimension[z+w]},(4 (n ((((z.(Subscript[δ, n])))^2)+(((w.(Subscript[δ, n])))^2))-((norm[z])^2)-((norm[w])^2)+2 z.w+(2 n-4) z.(Subscript[δ, n]) w.(Subscript[δ, n])))/(n volume[n] (((((norm[z])^2)+((norm[w])^2)-2 z.w+4 z.(Subscript[δ, n])w.(Subscript[δ, n])))^(1+(n/2))))])

bergmanKernelH[x_,y_,t_,u_]:=(makeVector[x];makeVector[t];unmakeVector[t];unmakeVector[u];With[{n=dimension[x]+1},((4 ((n-1) (((y+u))^2)-((norm[x-t])^2)))/(n volume[n] (((((norm[x-t])^2)+(((y+u))^2)))^(1+(n/2)))))])

Miscellaneous

harmonicConjugate

harmonicConjugate[u_,x_Symbol,y_Symbol]:=(unmakeVector[x];unmakeVector[y];(∫_0^y)(∂_x)u\[DifferentialD]y-(∫_0^x)((∂_y)u/.y->0)\[DifferentialD]x)

bergmanProjection

bergmanProjection[0,x_]:=0

bergmanProjection[u_,x_Symbol]:=(makeVector[x];Plus@@(Plus@@#1&)/@(Function[j,(((dimension[x]+2 #1[[2]]) #1[[1]])/(#1[[2]]+j+dimension[x]))&][#1[[2]]]/@#1[[1]]&)/@({polynomialDecomposition[#1[[1]],x],#1[[2]]}&)/@({dirichlet[#1[[1]],x],#1[[2]]}&)/@polynomialDecomposition[u,x])

bergmanProjection[u_,v_List]:=explicit[bergmanProjection,u,v]

schwarz

schwarz[0]=0;

schwarz[x_]:=PowerExpand[Simplify[(If[Head[x]==List,unmakeVector[x]];Module[{t,u},With[{s=Integrate[(((Cos[u])^(dimension[x]-2))/(((1+(1/(t^2))-((2 Sin[u])/t)))^((dimension[x])/2))), {u,0,Pi/2},GenerateConditions->False]},Together[((1-(1/(t^2))) (dimension[x]-1) volume[dimension[x]-1] (s-(s/.t->-t)))/(dimension[x] volume[dimension[x]])/.t->(1/(norm[x]))]]]/.HoldPattern[+t_]:>Factor[Select[+t,FreeQ[#1,ArcTan]&]]+Factor[Select[+t,!(FreeQ[#1,ArcTan])&]]/.(-1+t_) (1+t_)->-(1-(t^2)))]/. {(((-1+norm[x]))^(p_?EvenQ)):>(((1-norm[x]))^p), Abs[1 + norm[x]]->1 + norm[x], Abs[-1 + norm[x]]->1 - norm[x]}]

IntegerQ

IntegerQ[m_?IntegerQ+n_?IntegerQ]=True;

IntegerQ[m_?IntegerQ n_?IntegerQ]=True;

Times

HoldPattern[Times[a__]]/;Or@@Thread[vector[{a}]]:=With[{t=Position[Thread[vector[{a}]],True][[1,1]]},Row[{If[Head[Times@@Delete[{a},t]]===Plus,Row[{"(",First[Delete[{a},t]],")"}],Times@@Delete[{a},t]]," ",If[Head[{a}[[t]]]===Plus,Row[{"(",{a}[[t]],")"}],{a}[[t]]]}]]

HoldPattern[Times[a__]]/;Or@@Thread[matrix[{a}]]:=With[{t=Position[Thread[matrix[{a}]],True][[1,1]]},Row[{If[Head[Times@@Delete[{a},t]]===Plus,Row[{"(",First[Delete[{a},t]],")"}],Times@@Delete[{a},t]]," ",If[Head[{a}[[t]]]===Plus,Row[{"(",{a}[[t]],")"}],{a}[[t]]]}]]

explicit

explicit[function_,f_,v_List, opt___]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];function[f/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],z, Sequence @@ ({opt} /.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]])]/.Append[Thread[Table[ (Subscript[z, j]), {j,Length[v]}]->v],z->v]])

explicit[function_,f_,g_,v_List]:=(unmakeVector[v];Module[{z},setDimension2[z,Length[v]];function[f/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],g/.Thread[v->Table[ (Subscript[z, j]), {j,Length[v]}]],z]/.Append[Thread[Table[ (Subscript[z, j]), {j,Length[v]}]->v],z->v]])

Closing

Protect[Release[protectedWords]];

End[]

*)
EndPackage[ ]
(*EndPackage[]

Print["* You can now use the functions in this package."]