(* ::Package:: *) (************************************************************************) (* This file was generated automatically by the Mathematica front end. *) (* It contains Initialization cells from a Notebook file, which *) (* typically will have the same name as this file except ending in *) (* ".nb" instead of ".m". *) (* *) (* This file is intended to be loaded into the Mathematica kernel using *) (* the package loading commands Get or Needs. Doing so is equivalent *) (* to using the Evaluate Initialization Cells menu command in the front *) (* end. *) (* *) (* DO NOT EDIT THIS FILE. This entire file is regenerated *) (* automatically each time the parent Notebook file is saved in the *) (* Mathematica front end. Any changes you make to this file will be *) (* overwritten. *) (************************************************************************) If[$VersionNumber<7,Print["This package may not work correctly with versions of \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) earlier than \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) 7.0. Versions of this package that will work correctly with earlier versions of \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) are available at http://www.axler.net."]] Print["HFT7.m, version 7.00, 1 January 2009; for use with \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) version 7"] Print["The HFT7.m \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) 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 \!\(\* StyleBox[\"Mathematica\",\nFontSlant->\"Italic\"]\) 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 \!\(\* StyleBox[\"Harmonic\",\nFontSlant->\"Italic\"]\)\!\(\* StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\* StyleBox[\"Function\",\nFontSlant->\"Italic\"]\)\!\(\* StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\* StyleBox[\"Theory\",\nFontSlant->\"Italic\"]\) (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`"] \[CapitalDelta]::"usage"="\!\(\*SubscriptBox[\"\[CapitalDelta]\", \"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." \[Delta]::"usage"="\!\(\*SubscriptBox[\"\[Delta]\", \"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 \[CapitalDelta], 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, \!\(\*SubscriptBox[\"x\", \"j\"]\)] gives the partial derivative "<>"of f with respect to \!\(\*SubscriptBox[\"x\", \"j\"]\)." \[CapitalPhi]::"usage"="\[CapitalPhi][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." Begin["`private`"] Unprotect[IntegerQ,Limit,Transpose] protectedWords=Unprotect[Abs,Dot,Integer,Norm, Power,Times, Tr] 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]] matrix[IdMatrix[_]]^=True; IdMatrix/:a_ .IdMatrix[_]:=a IdMatrix/:IdMatrix[_].v_:=v 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[IdMatrix[n_]]:=n Tr[Transpose[A_]]:=Tr[A] Tr[c_ A_]:=c Tr[A]/;matrix[A] 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 = {}; realize[c_] := (real = real\[Union]{c}) ExpandNorm[f_]:=f//.Norm[g_+h_]:>Sqrt[Norm[g]^2+Norm[h]^2+2 g.h] 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[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 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[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] 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] 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[\[Delta], j_].Subscript[\[Delta], k_]:=Subscript[\[Delta], j,k] Subscript[\[Delta], j_] . x_Symbol:=(makeVector[x];Subscript[x, j])/;matrix[x]=!=True a_Symbol.Subscript[\[Delta], j_]:=(makeVector[a];Subscript[a, j])/;matrix[a]=!=True Subscript[\[Delta], j_] . f_Symbol[x_]:=(makeVector[f[_]];f[x][j])/;matrix[f[x]]=!=True/;f=!=List b_Symbol[x_].Subscript[\[Delta], j_]:=(makeVector[b[_]];b[x][j])/;matrix[b[x]]=!=True/;b=!=List Subscript[\[Delta], j_] . Grad[u_][x_]:= \!\(\*SubscriptBox["Partial", RowBox[{"{", "j", "}"}]]\)[u][x] Grad[u_][x_]. Subscript[\[Delta], j_]:= \!\(\*SubscriptBox["Partial", RowBox[{"{", "j", "}"}]]\)[u][x] J[b_][x_].Subscript[\[Delta], j_]:= \!\(\*SubscriptBox["Partial", RowBox[{"{", "j", "}"}]]\)[b][x] Subscript[\[Delta], j_] . J[b_][x_]:=Grad[b.Subscript[\[Delta], j]][x] Subscript[\[Delta], j_] . Y_?matrix:=Subscript[Y, j,"."] A_?matrix.Subscript[\[Delta], j_]:= \!\(\*SubscriptBox["A", RowBox[{"\"\<.\>\"", ",", "j"}]]\) Subscript[\[Delta], k_] . Subscript[Y, j_,"."]:=Subscript[Y, j,k] Subscript[A, j_,"."] . Subscript[\[Delta], k_]:=Subscript[A, j,k] Subscript[\[Delta], j_] . \!\(\*SubscriptBox["Y", RowBox[{"\"\<.\>\"", ",", "k_"}]]\):=Subscript[Y, j,k] \!\(\*SubscriptBox["A", RowBox[{"\"\<.\>\"", ",", "k_"}]]\) . Subscript[\[Delta], j_]:=Subscript[A, j,k] a_?VectorQ . Subscript[\[Delta], j_]:=a[[j]] Subscript[\[Delta], j_] . v_?VectorQ:=v[[j]] A_?MatrixQ . Subscript[\[Delta], j_]:=Transpose[A][[j]] Subscript[\[Delta], 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[\[Delta][_]]=True degree[f_,z_]:=degree1[Expand[f/.Norm[z]->Subscript[z, normReplacement]],z] degree1[0,_]:=-\[Infinity] degree1[HoldPattern[+a__],z_]:=Max[Thread[degree1[{a},z]]] degree1[f_,z_]:=Exponent[f/.Subscript[z, _]->z,z] 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[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] vector[Subscript[\[Delta], _]]^:=True Subscript[\[Delta], j_,k_]:=0/;j!=k Subscript[\[Delta], j_,k_]:=1/;j==k Subscript[\[Delta], j_][k_]:=Subscript[\[Delta], j,k] Norm[Subscript[\[Delta], _]]^=1; 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[\[Delta], j,k]) Partial[y_,y_]:=1 Partial[x_,Subscript[x_, j_]]:=(makeVector[x];Subscript[\[Delta], 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@@( \!\(\*SubscriptBox["Partial", RowBox[{"{", "#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]]], \!\(\*SuperscriptBox["f", "\[Prime]", MultilineFunction->None]\)[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]], \!\(\*SuperscriptBox["f", "\[Prime]", MultilineFunction->None]\)[g] Partial[g,y]]) Subscript[Partial, j_List][Subscript[Partial, k_List][f_]]:=Subscript[Partial, Sort[Join[k,j]]][f] 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[f_,x_Symbol]:=(makeVector[x];0)/;FreeQ[f,x] Grad[Subscript[x_, j_],x_Symbol]:=(makeVector[x];Subscript[\[Delta], 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], \!\(\*SuperscriptBox["f", "\[Prime]", MultilineFunction->None]\)[y] Grad[y,x]]) Grad[Norm]:=#1/Norm[#1]& Grad[Dimension]:=0& Grad[f_,v_List]:=(unmakeVector[v];(Partial[f,#1]&)/@v) 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[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[{ \!\(\*SuperscriptBox["f", "\[Prime]", MultilineFunction->None]\)[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[A_?MatrixQ]:=Plus@@Plus@@(A^2) HilbertSchmidt[IdMatrix[n_]]:=n Subscript[\[CapitalDelta], v_List][f_]:=(unmakeVector[v];Plus@@(Partial[f,{#1,2}]&)/@v) Subscript[\[CapitalDelta], x_Symbol][f_]:=(makeVector[x];0)/;FreeQ[f,x] Subscript[\[CapitalDelta], x_Symbol][Subscript[x_, j_]]:=(makeVector[x];0) Subscript[\[CapitalDelta], x_Symbol][HoldPattern[+a__]]:=Plus@@Thread[Subscript[\[CapitalDelta], x][{a}]] Subscript[\[CapitalDelta], x_Symbol][f_ g_]:=f Subscript[\[CapitalDelta], x][g]+Subscript[\[CapitalDelta], x][f] g+2 Grad[f,x].Grad[g,x] Subscript[\[CapitalDelta], x_Symbol][f_^p_]:=p f^(p-1) Subscript[\[CapitalDelta], 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[\[CapitalDelta], x][p]+2 f^(p-1) Grad[f,x].Grad[p,x]+f^p Log[f]^2 Grad[p,x] . Grad[p,x] Subscript[\[CapitalDelta], x_Symbol][f_[a_. x_+b_.]]:=(makeVector[x];If[b==0,Null,makeVector[b]];a^2 \[CapitalDelta][f][a x+b])/;FreeQ[a,x]&&FreeQ[b,x] Subscript[\[CapitalDelta], x_Symbol][a_ .x_]:=(makeVector[x];0)/;FreeQ[a,x] Subscript[\[CapitalDelta], x_Symbol][x_ .y_,x]:=(makeVector[x];0)/;FreeQ[y,x] Subscript[\[CapitalDelta], x_Symbol][(A_ .x_).y_]:=(makeVector[x];0)/;FreeQ[A,x]&&FreeQ[y,x] Subscript[\[CapitalDelta], x_Symbol][(x_ .Y_).z_]:=(makeVector[x];0)/;FreeQ[Y,x]&&FreeQ[z,x] Subscript[\[CapitalDelta], x_Symbol][a_ .(B_ .x_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[B,x] Subscript[\[CapitalDelta], x_Symbol][a_ .(x_ .Y_)]:=(makeVector[x];0)/;FreeQ[a,x]&&FreeQ[Y,x] Subscript[\[CapitalDelta], x_Symbol][(A_ .x_).x_]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x] Subscript[\[CapitalDelta], x_Symbol][(x_ .A_).x_]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x] Subscript[\[CapitalDelta], x_Symbol][x_ .(A_ .x_)]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x] Subscript[\[CapitalDelta], x_Symbol][x_ .(x_ .A_)]:=(makeVector[x];2 Tr[A])/;FreeQ[A,x] Subscript[\[CapitalDelta], b_Symbol][b_ .Grad[u_][b_]]:=(makeVector[b];2 Subscript[\[CapitalDelta], b][u[b]]+b.Grad[Subscript[\[CapitalDelta], b][u[b]],b]) Subscript[\[CapitalDelta], x_Symbol][Grad[u_][x_].x_]:=(makeVector[x];2 Subscript[\[CapitalDelta], x][u[x]]+x.Grad[Subscript[\[CapitalDelta], x][u[x]],x]) Subscript[\[CapitalDelta], x_Symbol][f_[g_]]:= \!\(\*SuperscriptBox["f", "\[DoublePrime]", MultilineFunction->None]\)[g] Grad[g,x] . Grad[g,x] + \!\(\*SuperscriptBox["f", "\[Prime]", MultilineFunction->None]\)[g] Subscript[\[CapitalDelta], x][g]/;(makeVector[x];vector[g]=!=True) \[CapitalDelta][Norm]:=(Dimension[#1]-1)/Norm[#1]& Subscript[\[CapitalDelta], 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[\[CapitalDelta], 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] \[CapitalDelta][Dimension]:=0& Subscript[\[CapitalDelta], x_Symbol,y_Symbol][f_]:=(unmakeVector[y];Subscript[\[CapitalDelta], x][f]+Partial[f,{y,2}]) (Subscript[\[CapitalDelta], x__]^n_)[f_]:=(Subscript[\[CapitalDelta], x]^(n-1))[Subscript[\[CapitalDelta], x][f]]/;Head[n]==Integer&&Positive[n] Laplacian = \[CapitalDelta] 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]]]{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[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 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[f_,d_,vc___]:=Plus@@(Homogeneous[f,#1,vc]&)/@Range[0,d]/;Head[d]==Integer Volume[n_Integer]:=\[Pi]^(n/2)/Gamma[n/2+1]/;Positive[n] SurfaceArea[n_Integer]:=(2 \[Pi]^(n/2))/Gamma[n/2]/;Positive[n] 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]:=\!\( \*SubsuperscriptBox[\(\[Integral]\), \(-1\), \(1\)]\(\((f /. { \*SuperscriptBox[\(Norm[x]\), \(2\)] -> \*SuperscriptBox[\(x\), \(2\)], Norm[x] -> Abs[x], \*SubscriptBox[\(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]] \!\( \*SubsuperscriptBox[\(\[Integral]\), \(0\), \(1\)]\(Expand[ \*SuperscriptBox[\(r\), \(Dimension[x] - 1\)]\ IntegrateSphere[f /. {Norm[x] -> r, \*SubscriptBox[\(x\), \(j_\)] -> r\ \*SubscriptBox[\(x\), \(j\)], x -> r\ x}, x]] \[DifferentialD]r\)\)]) 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]] Options[Dirichlet]={Region->Sphere, \[CapitalDelta]->0} Dirichlet[p_, x_Symbol,options__] := With[ {t = AntiLaplacian[(\[CapitalDelta] /. {options} /. Options[Dirichlet]),x, Singularity -> 0]},Factor[t+Dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, \[CapitalDelta]->_]]]]/; And[(Region /. {options} /. Options[Dirichlet]) === ExteriorSphere, MemberQ[{options}, \[CapitalDelta]-> _]] Dirichlet[p_, x_Symbol,options__] := With[ {t = AntiLaplacian[(\[CapitalDelta] /. {options} /. Options[Dirichlet]),x]}, t+Dirichlet[p-t,x, Sequence @@ DeleteCases[{options}, \[CapitalDelta]->_]]] /; MemberQ[{options}, \[CapitalDelta]-> _] Dirichlet[p_, x_Symbol,options___] := With[ {region = (Region /. {options} /. Options[Dirichlet])}, Which[region === Sphere, dirichlet[p,x], Head[region]===Quadratic, p + AntiLaplacian[-Subscript[\[CapitalDelta], 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, _],\[Infinity]]] 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] d[0,0]=1 lap[f_,v_]:=Plus@@(\!\( \*SubscriptBox[\(\[PartialD]\), \({#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[{\!\( \*UnderoverscriptBox[\(\[Sum]\), \(j = i\), \(k\)]\(d[i, j]\ \*SuperscriptBox[\(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\[Union]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] 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[\[CapitalDelta], 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[\[CapitalDelta], 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,\[Integral]t^(1-2 m-n) (\[Integral](t^(-1+2 m+n+j) f/.Norm[z]->t)\[DifferentialD]t)\[DifferentialD]t,With[{int=\!\( \*SubsuperscriptBox[\(\[Integral]\), \(0\), \(t\)]\(\(( \*SuperscriptBox[\(s\), \(\(-1\) + 2\ m + n + j\)]\ f /. Norm[z] -> s)\) \[DifferentialD]s\)\)},\[Integral]t^(1-2 m-n) If[FreeQ[int,Indeterminate]&&FreeQ[int,\[Infinity]]&&FreeQ[int,Limit]&&FreeQ[int,If]&&FreeQ[int,Integrate]&&(Singularity/.{options}/.Options[AntiLaplacian])===None,int,\[Integral](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[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[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[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 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[0,2]=1; DimensionH[m_,n_]:=Binomial[n+m-1,n-1]-Binomial[n+m-3,n-1] 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[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]]] Reflection[\[Infinity]]:=0 Reflection[x_]:=(If[Head[x]==Symbol,makeVector[x]];If[Norm[x]===0,\[Infinity],x/Norm[x]^2]) Reflection[x_List, Sphere[r_, c_]]:= Reflection[x, Sphere[r, Table[c.Subscript[\[Delta], j],{j,1,Length[x]}]]] /; Head[c] =!= List Reflection[x_, Sphere[r_, c_List]]:= Reflection[Table[x.Subscript[\[Delta], 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[\[Delta], j],{j,1,Length[x]}],c]] /; Head[b] =!= List Reflection[x_, Hyperplane[b_List,c_]]:= Reflection[Table[x.Subscript[\[Delta], j],{j,1,Length[b]}], Hyperplane[b, c]] /; Head[x] =!= List 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})) \[CapitalPhi][z_List]:=(unmakeVector[z];If[Last[z]===-1&&Norm[Drop[z,-1]]===0,\[Infinity],2 Reflection[z+Append[Table[0,{Length[z]-1}],1]]-Append[Table[0,{Length[z]-1}],1]]) \[CapitalPhi][z_]:=(If[!(z===0||z===\[Infinity]||MatchQ[z,_. S]),makeVector[z];If[Dimension[S]=!= Dimension[z],SetDimension[S,Dimension[z]]]];2 Reflection[z-S]+S) \[CapitalPhi][0,-1]=\[Infinity]; \[CapitalPhi][x_,y_]:= (makeVector[x];unmakeVector[y];{(2 x)/(Norm[x]^2+(y+1)^2),(2 (y+1))/(Norm[x]^2+(y+1)^2)-1}) 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\[CapitalPhi][x]])) KelvinM[u_,z_Symbol]:=(makeVector[z];2^(1/2 (Dimension[z]-2)) Norm[z-S]^(2-Dimension[z]) (u/.{z->\[CapitalPhi][z],Subscript[z, j_]:>\[CapitalPhi][z].Subscript[\[Delta], 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)})) 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[z_,w_]:=(makeVector[z];makeVector[w];(2 (z.Subscript[\[Delta], Dimension[z]]+w.Subscript[\[Delta], Dimension[z]]))/((Dimension[z] Volume[Dimension[z]]) (Norm[z]^2+Norm[w]^2-2 z.w+4 z.Subscript[\[Delta], Dimension[z]] w.Subscript[\[Delta], 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[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[z_,w_]:=(makeVector[z];makeVector[w];With[{n=Dimension[z+w]},(4 (n ((z.Subscript[\[Delta], n])^2+(w.Subscript[\[Delta], n])^2)-Norm[z]^2-Norm[w]^2+2 z.w+(2 n-4) z.Subscript[\[Delta], n] w.Subscript[\[Delta], n]))/(n Volume[n] (Norm[z]^2+Norm[w]^2-2 z.w+4 z.Subscript[\[Delta], n]w.Subscript[\[Delta], 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))]) HarmonicConjugate[u_,x_Symbol,y_Symbol]:=(unmakeVector[x];unmakeVector[y];\!\( \*SubsuperscriptBox[\(\[Integral]\), \(0\), \(y\)]\( \*SubscriptBox[\(\[PartialD]\), \(x\)]u \[DifferentialD]y\)\)-\!\( \*SubsuperscriptBox[\(\[Integral]\), \(0\), \(x\)]\(\(( \*SubscriptBox[\(\[PartialD]\), \(y\)]u /. y -> 0)\) \[DifferentialD]x\)\)) 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[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[m_?IntegerQ+n_?IntegerQ]=True; IntegerQ[m_?IntegerQ n_?IntegerQ]=True; 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]]) Protect[Release[protectedWords]]; End[] *) EndPackage[ ] (* EndPackage[] Print["* You can now use the functions in this package."]