/* Written 12:56 am  Jan 31, 1989 by dna@ballast.cs.umd.edu in uxe.cso.uiuc.edu:sci.math.symbolic */
Mark Phillips writes

>>  ... I'd like to see some Mathematica packages posted to this
>>  newsgroup.  Does anybody have anything they think might be useful to
>>  others?

I agree.  I've been holding off sending in the packages I've been working
on because I'm new to Mathematica, and certainly have a lot to learn.  But
then again, aren't we all?   So here is a package I've been using to work
with differential equations in the plane.  Two files follow. The package,
which defines grad, div, curl, etc., is called  veccalc2.m.
The second file, veccalc2.test uses veccalc2.m to verify 32 identities.
For example, it checks that

               grad[div[s]] - curl[rot[s]] = lap[s]

for an arbtrary function s[x,y].

Comments are welcome.  Even more welcome are other packages.  A particular
request:  does anyone have a package to move between real and complex
forms of differential operators?  For example, how could one use mathematica
to verify Goursat's formula,

           biharmonic[ Real[ Conjugate[z] f[z]+g[z] ] ]  = 0 ,

for all analytic functions f and g?

  -- Douglas N. Arnold

----------- Begin veccalc2.m ------------------------------------------------
BeginPackage["veccalc2`"]

veccalc2::usage =
"The package veccalc2.m defines grad, curl, div, rot, lap, lap2, grad2,
and eps.  These are the gradient, curl, divergence, rotation,
Laplacian, biharmonic, Hessian and infinitesmal strain tensor, for
functions of the 2 variables."

grad::usage =
"If r is a scalar depending of x and y, grad[r,{x,y}] is its gradient,
a 2-vector.  If r is a 2-vector, grad[r,{x,y}] is the gradient by rows,
a 2x2 matrix.  If the second argument is omitted, {x,y} is assumed."

curl::usage =
"If r is a scalar depending on x and y, curl[r,{x,y}] is its curl, a
2-vector.  If r is a 2-vector, curl[r,{x,y}] is the curl by rows, a 2x2
matrix.  If the second argument is omitted, {x,y} is assumed."

div::usage =
"If r is a 2-vector depending on x and y, then div[r] is its
divergence, a scalar.  If r is a 2x2 matrix, div[r] is its divergence
by rows, a 2-vector.  If the second argument is omitted, {x,y} is
assumed."

rot::usage =
"If r is a 2-vector depending on x and y, then rot[r] is its rotation,
a scalar.  If r is a 2x2 matrix, rot[r] is its rotation by rows, a
2-vector.  If the second argument is omitted, {x,y} is assumed."

lap::usage =
"If r depends on x and y, lap[r,{x,y}] is its Laplacian.  If r is
vector- or matrix-valued, lap applies to r componentwise.  If the
second argument is omitted, {x,y} is assumed."

lap2::usage =
"If r depends on x and y, lap2[r,{x,y}] is its biharmonic.  If r is
vector- or matrix-valued, lap applies to r componentwise.  If the
second argument is omitted, {x,y} is assumed."

grad2::usage =
"If r is a scalar depending of x and y, grad2[r,{x,y}] is its Hessian,
a 2x2 matrix.  If r is vector- or matrix-valued, grad2 applies
componentwise, giving a higher order tensor.  If the second argument is
omitted, {x,y} is assumed."

eps::usage =
"If r is a 2-vector depending on x and y, then eps[r,{x,y}] is the
symmetric part of the gradient of r, a 2x2-matrix.  If the second
argument is omitted, {x,y} is assumed."

Begin["`private`"]

(* gradient of a scalar function *)
grad[r_,{x_,y_}] := {D[r,x],D[r,y]}
(* gradient of a vector has gradients of components for rows *)
grad[s_,{x_,y_}] := { grad[s[[1]],{x,y}] , grad[s[[2]],{x,y}] }  /; Dimensions[s] == {2}
grad[s_] := grad[s,{Global`x,Global`y}]

(* (vector) curl of a scalar function *)
curl[r_,{x_,y_}] := {-D[r,y],D[r,x]}
(* curl of a vector has curls of components for rows *)
curl[s_,{x_,y_}] := { curl[s[[1]],{x,y}] , curl[s[[2]],{x,y}] }  /;  Dimensions[s] ==  {2}
curl[s_] := curl[s,{Global`x,Global`y}]

(* divergence of a vector valued function *)
div[s_,{x_,y_}]  := D[s[[1]],x] + D[s[[2]],y]
(* divergence applies to matrices row-wise *)
div[ss_,{x_,y_}]  := { div[ss[[1]],{x,y}], div[ss[[2]],{x,y}] }  /;  Dimensions[ss] == {2,2}
div[s_] := div[s,{Global`x,Global`y}]

(* (scalar) rotation of a vector valued function *)
rot[s_,{x_,y_}]  := D[s[[1]],y] - D[s[[2]],x]
(* rotation applies to matrices row-wise *)
rot[ss_,{x_,y_}]  := { rot[ss[[1]],{x,y}], rot[ss[[2]],{x,y}] }  /;  Dimensions[ss] == {2,2}
rot[s_] := rot[s,{Global`x,Global`y}]

(* Laplacian of a scalar, vector, or matrix-valued function *)
lap[r_,{x_,y_}]  := D[r,{x,2}] + D[r,{y,2}]
lap[s_] := lap[s,{Global`x,Global`y}]

(* biharmonic of a scalar, vector, or matrix-valued function of x and y *)
lap2[r__]  := lap[lap[r]]

(* Hessian *)
grad2[r__] := grad[grad[r]]

(* infinitesimal strain tensor *)
eps[s_,{x_,y_}] := ( grad[s,{x,y}] + Transpose[grad[s,{x,y}]] ) /2  /;  Dimensions[s] == {2}
eps[s_] := eps[s,{Global`x,Global`y}]

End[]
EndPackage[]
----------- End veccalc2.m --------------------------------------------------

----------- Begin veccalc2.test ---------------------------------------------
(*  This file tests veccalc2.m by verifying various identities*)

del = IdentityMatrix[2]
chi = {{0,-1},{1,0}}
r= r0[x,y]
s={s1[x,y],s2[x,y]}
ss = {{s11[x,y],s12[x,y]},{s21[x,y],s22[x,y]}}
Print[" 0:  ", curl[r] == chi.grad[r]   // ExpandAll ]
Print[" 1:  ", rot[s] == div[chi.s]    // ExpandAll ]
Print[" 2:  ", div[grad[r]]  ==  lap[r]    // ExpandAll ]
Print[" 3:  ", div[curl[r]] == 0    // ExpandAll ]
Print[" 4:  ", rot[grad[r]] == 0    // ExpandAll ]
Print[" 5:  ", rot[curl[r]] + lap[r] == 0   // ExpandAll ]
Print[" 6:  ", grad[div[s]] - curl[rot[s]]  ==  lap[s]    // ExpandAll ]
Print[" 7:  ", grad[rot[s]] + curl[div[s]]  ==  chi.lap[s]    // ExpandAll ]
Print[" 8:  ", curl[s] + grad[s].chi == {{0,0},{0,0}}    // ExpandAll ]
Print[" 9:  ", rot[ss] + div[ss.chi] == {0,0}    // ExpandAll ]
Print["10:  ", div[grad[s]]  ==  lap[s]    // ExpandAll ]
Print["11:  ", div[Transpose[grad[s]]]  ==  grad[div[s]]    // ExpandAll ]
Print["12:  ", div[curl[s]] == {0,0}    // ExpandAll ]
Print["13:  ", div[Transpose[curl[s]] ] ==  curl[div[s]]    // ExpandAll ]
Print["14:  ", rot[grad[s]] == {0,0}    // ExpandAll ]
Print["15:  ", rot[Transpose[grad[s]]] == grad[rot[s]]    // ExpandAll ]
Print["16:  ", rot[curl[s]] == -lap[s]    // ExpandAll ]
Print["17:  ", rot[Transpose[curl[s]]] == curl[rot[s]]    // ExpandAll ]
Print["18:  ", div[r del] == grad[r]    // ExpandAll ]
Print["19:  ", div[r chi] == curl[r]    // ExpandAll ]
Print["20:  ", rot[r del] == -curl[r]    // ExpandAll ]
Print["21:  ", rot[r chi] == grad[r]    // ExpandAll ]
Print["22:  ", grad[s] - Transpose[grad[s]] == - rot[s]chi    // ExpandAll ]
Print["23:  ", eps[s] == grad[s] + rot[s]chi/2    // ExpandAll ]
Print["24:  ", eps[s] == Transpose[grad[s]] - rot[s]chi/2    // ExpandAll ]
Print["25:  ", div[eps[s]] == lap[s]/2 + grad[div[s]]/2    // ExpandAll ]
Print["26:  ", div[eps[s]] == lap[s] + curl[rot[s]]/2    // ExpandAll ]
Print["27:  ", div[eps[s]] == grad[div[s]] - curl[rot[s]]/2    // ExpandAll ]
Print["28:  ", grad[curl[r]] == chi.grad[grad[r]]    // ExpandAll ]
Print["29:  ", curl[grad[r]] == - grad[grad[r]].chi    // ExpandAll ]
Print["30:  ", curl[curl[r]] == -chi.grad[grad[r]].chi    // ExpandAll ]
Print["31:  ", div[div[grad[grad[r]]]] == rot[rot[curl[curl[r]]]] // ExpandAll ]
Print["32:  ", div[div[grad[grad[r]]]] == lap[lap[r]]    // ExpandAll ]
----------- End veccalc2.test -----------------------------------------------
/* End of text from uxe.cso.uiuc.edu:sci.math.symbolic */
