(*^

::[paletteColors = 128; 
	fontset = title, "New York", 24, L3, center, bold, nohscroll;
	fontset = subtitle, "New York", 18, L2, center, bold, nohscroll;
	fontset = subsubtitle, "New York", 14, L2, center, bold, nohscroll;
	fontset = section, "New York", 14, L2, bold, nohscroll, grayBox;
	fontset = subsection, "New York", 12, L2, bold, nohscroll, blackBox;
	fontset = subsubsection, "New York", 10, L2, bold, nohscroll, whiteBox;
	fontset = text, "New York", 12, L2;
	fontset = smalltext, "New York", 10, L2;
	fontset = input, "Courier", 12, L2, bold, nowordwrap;
	fontset = output, "Courier", 12, L2, nowordwrap;
	fontset = message, "Courier", 12, L2, R65535, nowordwrap;
	fontset = print, "Courier", 12, L2, nowordwrap;
	fontset = info, "Courier", 12, L2, nowordwrap;
	fontset = postscript, "Courier", 12, L2, nowordwrap;
	fontset = name, "Geneva", 10, L2, italic, B65535, nohscroll;
	fontset = header, "New York", 10, L2;
	fontset = footer, "New York", 12, L2, center;
	fontset = help, "Geneva", 10, L2;
	fontset = clipboard, "New York", 12, L2;
	fontset = completions, "New York", 12, L2;
	fontset = network, "Courier", 10, L2, nowordwrap;
	fontset = graphlabel, "Courier", 12, L2;
	fontset = special1, "New York", 12, L2;
	fontset = special2, "New York", 12, L2, center;
	fontset = special3, "New York", 12, L2, right;
	fontset = special4, "New York", 12, L2;
	fontset = special5, "New York", 12, L2;]
:[font = input; ]
(* Article 573 (27 more) in sci.math.symbolic:
Subject: A mathematica package for calculus in the plane
Date: 31 Jan 89 06:56:02 GMT
Reply-To: dna@emmy.umd.edu (Douglas N. Arnold)
Organization: Math. Dept., Univ. of Maryland, College Park, MD 20742


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


(*  This file tests veccalc2.m by verifying various identities
	
	Remember to uncomment this block

----------- Begin veccalc2.test ---------------------------------------------

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 -----------------------------------------------
*)
^*)