GradPack v1.3 Notebook

Discussion

The built in Mathematica language and standard packages are capable of calculating various differential operators of scalar/vector functions, in various coordinate systems, but only coordinate-wise. That is to say, Mathematica is able to perform calculations such as:

Curl[{x+(y^2),z+Sin[x],xyz+1},Cartesian[x,y,z]]={-1,0,-2 y+Cos[x]}

(with the `VectorCalculus` package), but unable to give results such as:

(Subscript[∇, x])[α((((x^\[Transpose]).x))^n)] = 2 n x α ((((x^\[Transpose]).x))^(-1+n))+(Subscript[∇, x])[α] ((((x^\[Transpose]).x))^n)

Some packages come close to our needs (Symbolic Vector Analysis in Plasma Physics,
VECTAN 1.1),  but none (to my knowledge) is able to get input and produce output at an abstract coordinate-free level, in a readily readable format, for (R^n) vectors (and not just (R^3)) .

The GradPack version 1.3 currently contains implementations of Gradient of scalar fields, and Jacobian of vector fields. To this end, the two main routines that were implemented are:

(Subscript[∇, x])[α]-the gradient of the scalar field α  with respect to the vector variable x.

(Subscript[J, x])[v] -  the Jacobian matrix of the vector field v  with respect to the vector variable x.

The routines can handle dot and cross products, and most standard scalar functions. Vectors involved directly or implicitly in a cross product are internally assumed to be (R^3) vectors (that is, are assigned to replacement rules relevant for (R^3) vectors).
The employed notations are implemented using the standard package Utilities'Notation', and are designed to resemble and sometimes overload the standard ones:

({{Input, Alias Form, Full Form, Meaning}, {(A^ᵀ), A(Ctrl6) EsctrEsc, Transpose[A], The transpose of either a matrix or a vector A}, {A.B, -, Dot[A,B], Standard matrix/matrix or matrix/vector product.
Vector/vector dot product is accordingly inputted as:
		 (v^ᵀ).w}, {v×w, vEsccrossEscw, Cross[v,w], Cross product of the (R^3) vectors v and w}, {(Subscript[∇, x])[f], EscgradEsc, Grad[f,x], Gradient of the scalar function f   w.r.t. the vector x}, {(Subscript[J, x])[v] , EscjacEsc, Jac[v,x], Jacobian of the vector function v   w.r.t. the vector x}, {‖v‖, EscnormEsc, norm[v], Norm of the vector v (euclidean 2-norm)}, {(SuperStar[⌈v⌉]), EscskewEsc, Skew[v], Skew symetric matrix representing cross-by-v:
(SuperStar[⌈v⌉])=(({{0, -(Subscript[v, 3]), (Subscript[v, 2])}, {(Subscript[v, 3]), 0, -(Subscript[v, 1])}, {-(Subscript[v, 2]), (Subscript[v, 1]), 0}})) }})

It is highly unrecommended to make use the Full-Forms listed above: the package interprets internally (Subscript[∇, x])[f] (for example) to the intermediate form Gradd[f,x] to facilitate type checking and error messages, and only when the arguments passed the checks the function Grad[f,x] is invoked.  Thus, direct assignment such as -

Grad[f,x]=0

will bypass the entire argument-check mechanism, and may result in unexpected and even erroneous expressions.

Additional function implemented are:
SetConstant[x], which sets the relevant derivatives (gradients, jacobians) of its argument x to zero, with respect to any variable.
SetIndependent[x], which declares the vectors x[i] to be independent (i.e., (Subscript[J, x[i]])[x[j]]=(Piecewise[{{1, i==j}, {0, i!=j}}])).

The few core formulas are mostly readily-verifyable versions of the Leibnitz formula, such as:

(Subscript[∇, x])[α β]= α (Subscript[∇, x])[β] +β  (Subscript[∇, x])[α]
(Subscript[J, x])[α v]=v.((((Subscript[∇, x])[α]))^ᵀ)+α (Subscript[J, x])[v]
(Subscript[∇, x])[(v^ᵀ).w]=((((Subscript[J, x])[v]))^ᵀ).w + ((((Subscript[J, x])[w]))^ᵀ).v
(Subscript[J, x])[v×w]=(SuperStar[⌈v⌉]).(Subscript[J, x])[w]-(SuperStar[⌈w⌉]) .(Subscript[J, x])[v]
(Subscript[∇, x])[f[α]]=f'[α](Subscript[∇, x])[α]

The current design decision is to require a type declaration, and check only type-operator consistency (as opposed to, say, arithmetic consistency) .  The declaration and check is done by a "DeclareType" syntax, which internally translates to direct "TypeQ" values. Specifically, the statements:

DeclareScalar[{α,β,γ}]
DeclareVector[{v,u,w}]
DeclareMatrix[{A,B,C}]

(multiple declerations must be enclosed in a list) will cause internally results such as:

ScalarQ[α]=True
VectQ[v]=True
MatQ[A]=True

This syntax may be performence-wise inferior to head-assignment (as implemented, e.g. in the VecTan Mathematica package by Brend Fiedler), but is simpler both to implement and interpret. 
Another current design decision is to consider only column vectors as vectors (that is, make VectQ true only for such). This, among other benefits, enforces the possibility to take a Jacobian only of a column vector.

One comment is in order about a seemingly technical limitation: The routines included here aim only to facilitate generation of output in a traditional compact vector/matrix notation. Derivative of a matrix with respect to a vector would create a 3rd  rank tensor, which is unexpressable in such compact notation (there exists a considerable wealth of excellent tensor-calculus Mathematica packages for such needs, both commercial and public). Thus, in view of the objectives of this package, all symbolic matrices that are encountered in a differentiation routine are treated as constants. A warning message is generated whenever this is the case.

Future versions will supply some generic simplification TransformationFunctions (e.g., commuting vectors in an inner product, transforming among various forms of triple products, replacing (v/(‖v‖)) by a symbol (OverHat[v]), etc. ).
It is recommended at first to evaluate the entire Imlementation section, and then skim through the examples and experiment with variations as you go - to get a feeling as to how to use the supplied functions.

Implementation - evaluate this entire section to use.

General, Notations

Unprotect[Power,Dot,Transpose,Cross];

Off[General::spell, General::spell1]

Utilities`Notation`AutoUpdateNotationStyleSheet = False;

Utilities`Notation`AutoLoadNotationPalette = False;
(* Utilities`Notation`AutoUpdateNotationStyleSheet=False; *) 
<<Utilities`Notation`

Notation[((v_)^\[Transpose]) ⟺ Transpose[v_]]

Format[((((((x_?VectQ))^\[Transpose]).(x_?VectQ) ))^(a_))]:=((norm[x])^(2 a))
Format[((((x_?VectQ))^\[Transpose]).(x_?VectQ))]:=(((norm[x]))^2)

Notation[\[LeftDoubleBracketingBar]x_\[RightDoubleBracketingBar] ⟸ norm[x_]]
Notation[\[LeftDoubleBracketingBar]x_\[RightDoubleBracketingBar] ⟹ (((((x_)^\[Transpose]).x_))^(1/2))]

Notation[(v_\[Cross]w_) ⟸ Cross[v_?VectQ,w_?VectQ]]

Notation[(SuperStar[⌈v_⌉]) ⟸ Skew[v_?VectQ]]

Notation[(Subscript[∇, x_])[f_] ⟹ Gradd[f_,x_]]
Notation[(Subscript[∇, x_])[f_] ⟸ Grad[f_,x_]]

Notation[(Subscript[J, x_])[f_] ⟹ Jacc[f_,x_]]
Notation[(Subscript[J, x_])[f_] ⟸ Jac[f_,x_]]

AddInputAlias[\[LeftDoubleBracketingBar]\[Placeholder]\[RightDoubleBracketingBar],"norm"]
AddInputAlias[(SuperStar[⌈\[Placeholder]⌉]),"skew"]
AddInputAlias[(Subscript[∇, \[Placeholder]])[\[Placeholder]],"grad"]
AddInputAlias[(Subscript[J, \[Placeholder]])[\[Placeholder]],"jac"]

Scalar / Vector / Matrix Type Properties :

DeclareScalar::usage ="DeclareScalar[x] enables the use of x as a scalar, both in differentiation and in general linearity simplifications. DeclareScalar[{x,y,...}] declares all the contents of the list as scalars in a single statement.";
DeclareVector::usage ="DeclareVector[x] enables the use of x as a vector, both in differentiation and in general linearity simplifications. DeclareVector[{x,y,...}] declares all the contents of the list as vectors in a single statement.";
DeclareMatrix::usage ="DeclareMatrix[x] enables the use of x as a costant matrix, both in differentiation and in general linearity simplifications. DeclareMatrix[{x,y,...}] declares all the contents of the list as matrices in a single statement.";

DeclareScalar[v_List]:=(DeclareScalar/@v;)
DeclareScalar[v_]:=(ScalarQ[v]=True;);
DeclareVector[v_List]:=(DeclareVector/@v;)
DeclareVector[v_]:=(VectQ[v]=True;);
DeclareMatrix[v_List]:=(DeclareMatrix/@v;)
DeclareMatrix[v_]:=(MatQ[v]=True;);

ScalarQ[α_?NumericQ]:=True;
ScalarQ[f_[α__?ScalarQ]]:=True /;MemberQ[Attributes[f], NumericFunction] ;
ScalarQ[Derivative[n_][f_][x_]]:=True /;TrueQ[ScalarQ[f[x]]];
ScalarQ[(((v_?VectQ))^\[Transpose]).w_?VectQ]:=True;
ScalarQ[Dt[α_?ScalarQ,x_]]:=True;
ScalarQ[Sum[α_?ScalarQ,{list__}]]=True;

VectQ[α_?ScalarQ v_?VectQ]=True;
VectQ[v_?VectQ + w_?VectQ]=True;
VectQ[Sum[v_?VectQ,{list__}]]=True;
VectQ[Cross[v_?VectQ,w_?VectQ]]=True;
VectQ[Grad[___]]=True;
VectQ[M_?MatQ . v_?VectQ] =True;
VectQ[Dt[v_?VectQ,x_]]:=True;

MatQ[α_?ScalarQ M_?MatQ]=True;
MatQ[A_?MatQ + B_?MatQ]=True;
MatQ[A_?MatQ.B_?MatQ]=True;
MatQ[Sum[A_?MatQ,{list__}]]=True;
MatQ[Jac[_?VectQ]]=True;
MatQ[v_?VectQ . Transpose[w_?VectQ]]=True;
MatQ[Transpose[A_?MatQ]]=True;
MatQ[Dt[A_?MatQ,x_]]:=True;

Dot Properties:

intended in genral for matrices, applies in particular for vectors, so type checking is omitted.

(x_+ y_) . z_:=x . z + y . z
x_ . (y_+z_):=x . y + x . z

x_ . (α_?ScalarQ y_) := α  (x . y)
(α_?ScalarQ y_) . x_:= α  (y . x)
x_ . (α_?ScalarQ) := α  x
(α_?ScalarQ ) . x_:= α  x
0 . x_ :=0
x_ . 0 := 0

id Properties:

MatQ[id]=True;

id/: id . x_?VectQ:= x
id/:Transpose[x_?VectQ] . id:=Transpose[x]
id/: id . A_?MatQ :=A
id/: A_?MatQ . id :=A
id/:(id^\[Transpose])=id;

Transpose Properties:

Transpose[x_+y_]:=Transpose[x]+Transpose[y]
Transpose[α_?ScalarQ x_]:=α Transpose[x]
Transpose[Transpose[x_]]:=x
Transpose[α_?ScalarQ]:=α
Transpose[A_ .  B_]:=Transpose[B].Transpose[A]/;(Head[A]=!=Transpose || Head[B]=!=Transpose )
Transpose[Sum[f_,{list__}]]:=Sum[Transpose[f],{list}]

Implement inner-product commutativity :

Transpose/:(Transpose[x_?VectQ]).( y_?VectQ):= Transpose[y] . x /;({y,x}==Sort[{x,y}] && (x=!=y))

In a 'triple' dot product such as v.(w^ᵀ).u, treat (w^ᵀ).u as a scalar :

Transpose/:v_?VectQ.Transpose[w_?VectQ].u_?VectQ:=(Transpose[w].u) v

The matrix representing cross-product by a vector is skew-symmetric:

Transpose[Skew[v_?VectQ]]:=-Skew[v]

Cross Properties:

Linearity:

Cross[v_?VectQ + w_?VectQ,u_?VectQ]:=Cross[v,u]+Cross[w,u]
 Cross[u_?VectQ,v_?VectQ + w_?VectQ]:=Cross[u,v]+Cross[u,w]
Cross[α_?ScalarQ u_?VectQ,v_?VectQ]:=α Cross[u,v]
Cross[u_?VectQ,α_?ScalarQ v_?VectQ]:=α Cross[u,v]
Cross[0,u_?VectQ]=0;
Cross[u_?VectQ,0]=0;

Geometric identities:

Cross[u_?VectQ,u_?VectQ]=0;
Dot[Transpose[v_?VectQ],Cross[u_?VectQ,v_?VectQ]]=0;
Dot[Transpose[v_?VectQ],Cross[v_?VectQ,u_?VectQ]]=0;
Dot[Transpose[Cross[v_?VectQ,u_?VectQ]],v_?VectQ]=0;
Dot[Transpose[Cross[v_?VectQ,u_?VectQ]],u_?VectQ]=0;
Cross[v_?VectQ ,u_?VectQ]:=-Cross[u,v]/;({u,v}==Sort[{v,u}] && (u=!=v))

MatQ[Skew[v_?VectQ]]=True;
Skew[v_?VectQ].w_?VectQ := Cross[v,w];

Grad / Jac Properties - Main SubSection:

Error messages, results for constants

Jacc::type="Variable `1` not declared as a vector";
Gradd::type="Variable `1` not declared as a scalar";
Jac::ConstMat = "A Jacobian computation encountered the matrix `1` - treating it as constant";
Jac::usage="\!\(J\_x\)[v] gives the Jacobian of the vector field v with respect to the vector variable x";
Grad::usage="\!\(∇\_x\)[α]  gives the Gradient of the scalar field α with respect to the vector variable x";

Jacc[v_,x_]:=If[TrueQ[VectQ[v]],Jac[v,x],(Message[Jacc::type,HoldForm[v]];error)];
Gradd[α_,x_]:=If[TrueQ[ScalarQ[α]],Grad[α,x],(Message[Gradd::type,HoldForm[α]];error)];

Jac[x_,x_]:=id
Grad[α_?NumericQ,x_]=0;

Arithmetic

Grad[f_+g_,x_]:=Gradd[f,x]+Gradd[g,x]
Grad[Sum[f_,{list__}],x]:=Sum[Gradd[f,x],{list}]
Jac[v_+w_,x_]:=Jacc[v,x]+Jacc[w,x]
Jac[Sum[v_,{list__}],x_]:=Sum[Jacc[v,x],{list}]

Leibnitz versions

Grad[α_ β_,x_]:= α Gradd[β,x] +β  Gradd[α,x]
Jac[α_?ScalarQ z_,x_]:=z.(((Gradd[α,x]))^\[Transpose])+α Jacc[z,x]
Grad[Transpose[v_] . w_,x_]:=Transpose[Jacc[v,x]] .  w + Transpose[Jacc[w,x]]. v
Jac[(v_\[Cross]w_),x_]:=Skew[v].Jacc[w,x]-Skew[w] .Jacc[v,x]

Jac[A_?MatQ.v_?VectQ,x_]:=(Message[Jac::ConstMat,HoldForm[A]];A.Jacc[v,x])

General and special scalar derivatives

Grad[((α_)^(β_)),x_]:=(α^β)(Log[α]Gradd[β,x]+(β/α)Gradd[α,x])
Grad[f_[y_],x_]:=(D[f[t], t] /.t->y) Gradd[y,x]/;MemberQ[Attributes[f], NumericFunction]

SetConstant , SetIndependent

SetConstant::usage ="SetConstant[x] sets to zero the derivatives of x with respect to any variable. i.e., if x is a scalar \!\(∇\__\)[x]=0, and if x is a vector - \!\(J\__\)[x]=0. SetConstant[{x,y,...}] sets all the contents of the list to constants.";
SetIndependent::usage="SetIndependent[x] for a vector sequence argument x delayed-sets (Subscript[J, x[i]])[x[j]] to 0 if i~=j and to id if i==j";
SetIndependent::type="Variable `1` not declared as a vector";

SetConstant[x_List]:=(SetConstant/@x;)
SetConstant[x_]:= (If[ScalarQ[x],(*then*) (Subscript[∇, y_])[x]^=0;,
								(*else*)If[VectQ[x], (Subscript[J, y_])[x]^=0;],
								(*unevaluated*)If[VectQ[x], (Subscript[J, y_])[x]^=0;]
					] ;       (*      (*either way*)     Dt[x,y_]=0;*)
				    )

SetIndependent[x_]:= If[TrueQ[VectQ[x[1]]],(*then*) (Subscript[J, x[i_]])[x[j_]]^:=DiscreteDelta[i-j] id;,
								(*else*) (Message[SetIndependent::type,HoldForm[x]];error)];

Wrapping it up

Protect[Power,Dot,Transpose,Cross];

Examples

Elementary Examples

Declare some variables to work with:

DeclareVector[{v,u,w}]
DeclareScalar[{α,β,γ}]

(Subscript[∇, v])[(v^\[Transpose]).v]

(Subscript[∇, v])[((((v^\[Transpose]).v))^((-1)/2))]

(Subscript[∇, v])[(1/(\[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]))]

(Subscript[J, v])[%]  (* that is, the Hessian of the original scalar field (1/(\[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar])) *)

%//Simplify

(Subscript[J, v])[((Sin[(α^2)\[LeftDoubleBracketingBar]v\[RightDoubleBracketingBar]])^β) v]//Simplify

(Subscript[J, v])[v\[Cross](α w)]

Clear[v,u,w,α,β,γ]

Spring Forces and Jacobians

Consider a spring connecting the points x[1] and x[2] with rest length L0 and spring constant k. We will form the spring quadratic energy (the usual Hook's law), derive the resulting force on one of the points, and compute the jacobian of this force w.r.t.the same point.

DeclareVector[x[i_]]
SetIndependent[x] (*  practically sets (Subscript[J, x[i_]])[x[j_]] to id\[Cross]delta(i,j)  *)

DeclareScalar[{k,L0}]
SetConstant[{k,L0}]

Δ:=x[2]-x[1]
Energy:=(1/2)k((( \[LeftDoubleBracketingBar]Δ\[RightDoubleBracketingBar]-L0))^2)

Force2:=-(Subscript[∇, x[2]])[Energy]

Force2

which is correct, but too complex. To obtain more readable results, we wish to expand the high-level magnitudes (in this case, say, Δ[i]) to be expanded only when necessary - i.e., only when derivated. To achieve this, we need to clear the explicit definition of the Δ[i], and assign such definitions only to it's derivatives. While this can be done manually, we use the following trick: clear Δ (redeclare it as a vecotr), create auxillary-Δ symbol to contain the explicit definition, and assign auxΔ's derivatives to Δ's derivatives. This way, Δ's explicit definition will be called upon only while trying to evaluate it's derivative.

Clear[Δ]

auxΔ:=x[2]-x[1]

DeclareVector[Δ]

(Subscript[J, v_])[Δ]^:=(Subscript[J, v])[auxΔ]

A much more pleasent expression emerges:

Force2

The same goes for Jacobian computations:

ForceJac2=(Subscript[J, x[2]])[Force2]

Results are subject to ordinary simplification:

%//FullSimplify

Taking the Jacobian of the force w.r.t.the other point, negates the sign of the resut, as expected:

ForceJac1=(Subscript[J, x[1]])[Force2]//FullSimplify

ForceJac1+ForceJac2//Simplify

Clear[x,Δ,k,L0,auxΔ,Energy,Force2,ForceJac2,ForceJac1]

Bending Forces of a Thread

Now for a more 'real life' example: model a thread in 3D as a sequence of keypoints q[i], with interpolating straight segments. We'll setup a 'Bending Energy' to each joint at each keypoint, and derive the resultant justification forces.

DeclareVector[q[i_]]
SetIndependent[q]

We'll define a bending energy at a joint as the minus dot-product of the two adjacent (normalized) segment directions at the joint.

Δ[i_]:=q[i+1]-q[i]
T[i_]:=((Δ[i])/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar]))
Energy[i_]:=-((T[i-1])^\[Transpose]).T[i]

The overall justification force on q[i] results from contributions of all the bending energies which are affected by q[i]'s position:

F[i_]:=-(Subscript[∇, q[i]])[Energy[i-1]+Energy[i]+Energy[i+1]]

F[i]

Which is entirely true, and yet utterly useless. Let us re-apply the approach taken in the previous example, which tends to yield much more readable results:

(1) Identify 'high-level' quantities that if used to express the desired expression, will probably give more compact results (in our case, the candidates are the Δ[i]'s and the T[i]'s).
(2) Replace the explicit definitions for the high-level objects with definitions to their derivatives only.
(3) Compute the desired derivative-this time using the high-level building blocks instead of elementary ones.

Clear[Δ,T]

DeclareVector[{Δ[i_],T[i_]}]

auxΔ[i_]:=q[i+1]-q[i]
auxT[i_]:=((Δ[i])/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar])) (* and NOT ((auxΔ[i])/(\[LeftDoubleBracketingBar]auxΔ[i]\[RightDoubleBracketingBar])) ! *)

T/:(Subscript[J, v_])[T[i_]]:=(Subscript[J, v])[auxT[i]]
Δ/:(Subscript[J, v_])[Δ[i_]]:=(Subscript[J, v])[auxΔ[i]]

F[i]

%//Simplify

We can obtain even more compact results - at present, only manually. Recognize the occurence of T[i] on (Subscript[J, q[j]])[T[i]]:

(Subscript[J, q[i]])[T[i]]
(Subscript[J, q[i+1]])[T[i]]

(Subscript[J, q[i]])[T[i]]-(((auxT[i].((auxT[i])^\[Transpose])-id)/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar])))//Simplify
(Subscript[J, q[i+1]])[T[i]]-(-((auxT[i].((auxT[i])^\[Transpose])-id)/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar])))//Simplify

So, we manually assign these jacobian values:

T/:(Subscript[J, v_])[T[i_]]=.

T/:(Subscript[J, q[j_]])[T[i_]]:=(Piecewise[{{((T[i].((T[i])^\[Transpose])-id)/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar])), j-i==0}, {-((T[i].((T[i])^\[Transpose])-id)/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar])), j-i==1}, {0, True}}])

F[i]

which can be further simplified manually with-

a[i_]:=(1/(\[LeftDoubleBracketingBar]Δ[i]\[RightDoubleBracketingBar]))(id-T[i].((T[i])^\[Transpose])).(T[i-1]+T[i+1])

F[i]-(a[i-1]-a[i])//Simplify

Clear[F,a,T,Δ,auxΔ,auxT,Energy]

Package Loading

This notebook is accompanied by a parallel package, with identical functionality .To use it, first copy " GradPack.m" to a directory within theMathematica path (for instance, on a windows platform the user path 'C:\Documents and Settings\All Users\Application Data\Mathematica\Applications' is recommended). To invoke the package, use-

<<GradPack.m

Credits

Copyright © 2005-2006, Ofek Shilon, Simbionix Ltd.

For bug reports, comments and suggestions: contact ofek@simbionix.com.