HFT

Initialization

Print

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

Print["HFT7.m, version 7.01, 1 January 2012; for use with Mathematica versions 7 and higher."]

Print["The HFT7.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 Mathematica notebook called ComputingWithHarmonicFunctions7.nb."]

Print["The most recent version of this HFT7.m package and its documentation notebook ComputingWithHarmonicFunctions7.nb 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 by electronic mail to axler@sfsu.edu."]

Print["* The computer is now unpacking the HFT7.m package."]

BeginPackage["HFT`"]

Usage

Δ::"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[f, x] gives a polynomial "<>"whose Laplacian equals f.  "<>"Here f must be a polynomial function of x."

Ball::"usage"="Ball is one of the values that may be assigned to the option "<>"Orthonormal."

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"="BergmanKernel[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[f , x] solves the "<>"BiDirichlet problem with boundary data f, "<>"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."

DimensionH::"usage"="DimensionH[m, n] gives the vector space dimension of the "<>"space of spherical harmonics of degree m in n variables."

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

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

DoubleBracketingBar::"usage"="DoubleBracketingBar is equal to the Norm."

ExpandNorm::"usage"="ExpandNorm[f] gives f with all terms of the form |x + y| "<>"replaced by Sqrt[ |x|^2 + |y|^2 + 2 x.y ]."

ExteriorSphere::"usage"="ExteriorSphere is an option for Dirichlet, "<>"specifying that the region should be the exterior " <>
"of the unit sphere."

ExteriorNeumann::"usage"="ExteriorNeumann[f , x] solves the "<>"exterior Neumann problem with boundary data f, "<>"as a function of x."

Grad::"usage"="Grad[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[u, x] "<>"gives the harmonic decomposition of u with respect to "<>"the variable x."

HilbertSchmidt::"usage"="HilbertSchmidt[A] gives the "<>"HilbertSchmidt norm of a matrix A."

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

Hyperplane::"usage"="Hyperplane[b, c] denotes the hyperplane b.x = c."

IdMatrix::"usage"="IdMatrix[n] denotes the n-by-n identity matrix. Unlike IdentityMatrix[n], IdMatrix[n] stays unevaluated even when n is a specific positive integer, although it will have the correct properties. For example, IdMatrix[3].x will evaluate to x."

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

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."

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

J::"usage"="J[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."

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

Laplacian::"usage"="Laplicaian 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 "<>"Norm[x]^2, where x is the variable."

Neumann::"usage"="Neumann[f, g, x] solves the "<>"Neumann problem of finding a function of x whose "<>"outward normal derivative (on the unit sphere) is f "<>"and whose Laplacian is g."

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

Orthonormal::"usage"="Orthonormal is an option for BasisH.  "<>"The default value is None.  The value Ball produces a "<>"basis that is orthonormal with respect to volume measure "<>"on the ball.  The value Sphere produces a "<>"basis that is orthonormal with respect to surface area "<>"measure on the sphere."

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."

Region::"usage"="Region is an option for Dirichlet."

S::"usage"="S 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."

Sphere::"usage"="Sphere is one of the values that may be assigned to the option "<>"Orthonormal."

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."

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,Integer,Norm, Power,Times, Tr]

$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[ togetherness , Message[Togetherness::AlreadyOn],$Post=If[zerotozero,Composition[zeroToZero, postTogether, usersPost],Composition[postTogether,usersPost]];togetherness=True; Message[Togetherness::On]]

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

$Post=Composition[postTogether, usersPost];

togetherness=True; zerotozero = False;

TurnOff[ZeroToZero]^:= If[ zerotozero ,$Post=If[togetherness,Composition[postTogether,usersPost], usersPost];zerotozero=False;Message[ZeroToZero::Off],  Message[ZeroToZero::AlreadyOff]]

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

Calculus on (R^n)

IdMatrix

matrix[IdMatrix[_]]^=True;

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

IdMatrix/:IdMatrix[_].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[IdMatrix[n_]]:=IdMatrix[n]

Transpose[{0}]:={0}

Tr

Tr[IdMatrix[n_]]:=n

Tr[Transpose[A_]]:=Tr[A]

Tr[c_ A_]:=c Tr[A]/;matrix[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)]))

If[$VersionNumber < 5, Norm[x_List]:=(Sqrt[x.x]); Norm[0] :=0]

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

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

Norm[c_?NumberQ 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[Grad[_][_]]=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_?NumberQ 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[J[_][_]]=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]||NumberQ[c]

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

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

a_ .{c_ y_}:=c a.{y}/;vector[y]||matrix[y]||NumberQ[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_]) . Grad[u_][x_]:=(Subscript[Partial, {j}])[u][x]

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

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

(Subscript[δ, j_]) . J[b_][x_]:=Grad[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]}]]

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[Grad[_][_]]=True

primitiveVector[δ[_]]=True

degree

degree[f_,z_]:=degree1[Expand[f/.Norm[z]->(Subscript[z, normReplacement])],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[0,x_]:=Null

polynomialDecomposition[p_,x_]:=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, Array[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_] :=
 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]]],J[f][g].Partial[g,(Subscript[x, j])],Grad[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]]],J[f][g].Partial[g,y],Grad[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]:=Grad[f,z].z/.Norm[z]->1

NormalD[f_,v_List]:=(unmakeVector[v];Together[Grad[f,v].v]/.((Norm[v])^2)->1)

Grad

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

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

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

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

Grad[f_ .g_,x_Symbol]:=f.J[g,x]+g.J[f,x]

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

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

Grad[Norm]:=((#1)/(Norm[#1]))&

Grad[Dimension]:=0&

Grad[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]+Grad[f,x].g/;(makeVector[x];vector[g])

Divergence[A_?matrix.f_,x_Symbol]:=Tr[A.J[f,x]]/;FreeQ[A,x]

Divergence[f_ .T_?matrix,x_Symbol]:=Tr[Transpose[T].J[f,x]]/;FreeQ[T,x]

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

J

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

J[x_,x_Symbol]:=(makeVector[x];IdMatrix[Dimension[x]])

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

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

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

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

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

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

Dimension[J[f_][x_]]:= {Dimension[f[_]],Dimension[x] }

J[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[IdMatrix[n_]]:=n

Δ (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 Grad[f,x].Grad[g,x]

(Subscript[Δ, x_Symbol])[((f_)^(p_))]:=p (f^(p-1)) (Subscript[Δ, x])[f]+p (p-1) (f^(p-2)) Grad[f,x] . Grad[f,x] +2 p (f^(p-1)) Log[f] Grad[f,x].Grad[p,x]+(f^p) Log[f] (Subscript[Δ, x])[p]+2 (f^(p-1)) Grad[f,x].Grad[p,x]+(f^p) ((Log[f])^2) Grad[p,x] . Grad[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 Tr[A])/;FreeQ[A,x]

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

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

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

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

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

(Subscript[Δ, x_Symbol])[f_[g_]]:=f''[g]  Grad[g,x] . Grad[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]

Δ[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]

Laplacian = Δ

Dimension

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

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

Dimension[a_. S+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."

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

SetDimension[x_Symbol,n_]:=(setDimension[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];setDimension[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];setDimension[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];setDimension[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[togetherness&&!(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},setDimension[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,Limit[c_. (r^(Dimension[x]+m_.)),r->0]:=0/;(m>0||m==0)&&FreeQ[c,r]];Dimension[x] Volume[Dimension[x]] (∫_0^1)Expand[(r^(Dimension[x]-1)) IntegrateSphere[f/.{Norm[x]->r,(Subscript[x, j_])->r (Subscript[x, j]),x->r x},x]]\[DifferentialD]r])

IntegrateSphere

IntegrateSphere[f_,v_List]:=(unmakeVector[v];Module[{z},setDimension[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]]

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[ {region = (Region /. {options} /. Options[Dirichlet])}, Which[region === Sphere, dirichlet[p,x],
Head[region]===Quadratic, p + AntiLaplacian[-(Subscript[Δ, x])[p], x, Multiple-> region],
region === ExteriorSphere, exteriorDirichlet[p,x],
Head[region] === Annulus, If[Head[p]===List, annularDirichlet[p[[1]],p[[2]],region[[1]], region[[2]],x], annularDirichlet[p,p,region[[1]], region[[2]],x]]]]

dirichlet[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]->norm,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]->norm,x]]//.{a_.+(norm^(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[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],antiLaplacian[Expand[p],x]]])

antiLaplacianQ[p_, x_, q_ ] := 
With[{ph = blockPolynomialDecomposition[p, x]},If[Length[q] == 1, (q[[1]].((Array[x,Length[q[[1]]]])^2))Sum[antiLaplacianQ2[ph[[j]], x, q[[1]]], {j,1, Length[ph]}],
If[Length[q]==2, (q[[1]].((Array[x,Length[q[[1]]]])^2)+  q[[2]])Sum[antiLaplacianQ4[ph[[j]], x, List @@ q], {j,1, Length[ph]}],(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,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},setDimension[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]]])

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

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

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

antiLaplacian[c_ g_,z_]:=c antiLaplacian[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]] antiLaplacian[((((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]]

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.Grad[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)]])

DimensionH

DimensionH[0,2]=1;

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

BasisH

Options[BasisH]={Orthonormal->None};

BasisH[0,x_,options___?OptionQ]:=With[{orthonormal=Orthonormal/.{options}/.Options[BasisH]},If[orthonormal===None,{1},Which[orthonormal===Sphere,{1},orthonormal===Ball,{(1/(Sqrt[Volume[Dimension[x]]]))}]]]

BasisH[m_Integer,x_Symbol,options___?OptionQ]:=Together[With[{orthonormal=Orthonormal/.{options}/.Options[BasisH]},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},If[orthonormal===None,noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]],orthonormalize[ttt,Which[orthonormal===Sphere,IntegrateSphere,orthonormal===Ball,IntegrateBall][#1 #2,x]&]]]]]]/;Dimension[x]==2

BasisH[m_Integer,x_Symbol,options___?OptionQ]:=Together[With[{orthonormal=Orthonormal/.{options}/.Options[BasisH]},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},If[orthonormal===None,noLeadingMinus/@Transpose[FactorTermsList/@ttt][[2]],orthonormalize[ttt,Which[orthonormal===Sphere,IntegrateSphere,orthonormal===Ball,IntegrateBall][#1 #2,x]&]]]]]]/;Head[Dimension[x]]==Integer

BasisH[m_Integer,v_List,options___?OptionQ]:=(unmakeVector[v];Module[{z},setDimension[z,Length[v]];Expand[BasisH[m,z,options]/.Thread[Table[ (Subscript[z, j]), {j,Length[v]}]->v]]/.Norm[z]->Norm[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[r_, c_]]:= Reflection[x, Sphere[r, Table[c.(Subscript[δ, j]),{j,1,Length[x]}]]] /; Head[c] =!= List

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

Reflection[x_, Sphere[r_, c_]]:= 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];((Norm[x])^(2-Length[x])) (u/.Thread[x->Reflection[x]]))

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,_. S]),makeVector[z];If[Dimension[S]=!= Dimension[z],SetDimension[S,Dimension[z]]]];2 Reflection[z-S]+S)

Φ[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})

S

Norm[S]^=1;

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

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

makeVector[S];

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

(Subscript[S, j_]):=-1/;j==Dimension[S]

(Subscript[S, j_]):=0/;j<Dimension[S]

KelvinM

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

KelvinM[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]]))

KelvinM[u_,z_Symbol]:=(makeVector[z];(2^((1/2) (Dimension[z]-2))) ((Norm[z-S])^(2-Dimension[z])) (u/.{z->Φ[z],(Subscript[z, j_]):>Φ[z].(Subscript[δ, j])}))

KelvinM[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},setDimension[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},setDimension[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."]