# Is it possible to do vector calculus in Mathematica?

I am trying to rearrange and manipulate some vector differential equations in Mathematica. As far as I understand you have to tell Mathematica that a variable is a vector by specifying the components of the vector. For example

r = {x, y, z};


If I want to define vector fields I have to do it in the following way

v = {vx[r, t], vy[r, t], vz[r, t]};
n = {nx[r, t], ny[r, t], nz[r, t]};


Now I can use the Div operator to express a differential equation, in this case the continuity equation in hydrodynamics.

D[n, t] + Div[n*v, r] == 0


My problem is the output I get from this. It looks absolutely horrible and I can't do anything with it.

   ({0, 0, 0},1)                                    ({0, 0, 1},0)
{nx             [{x, y, z}, t] + vz[{x, y, z}, t] nz             [{x, y, z}, t] +
({0, 0, 1},0)                                    ({0, 1, 0},0)
nz[{x, y, z}, t] vz             [{x, y, z}, t] + vy[{x, y, z}, t] ny             [{x, y, z}, t] +
({0, 1, 0},0)                                    ({1, 0, 0},0)
ny[{x, y, z}, t] vy             [{x, y, z}, t] + vx[{x, y, z}, t] nx             [{x, y, z}, t] +
({1, 0, 0},0)
nx[{x, y, z}, t] vx             [{x, y, z}, t],
({0, 0, 0},1)                                    ({0, 0, 1},0)
ny             [{x, y, z}, t] + vz[{x, y, z}, t] nz             [{x, y, z}, t] +
({0, 0, 1},0)                                    ({0, 1, 0},0)
nz[{x, y, z}, t] vz             [{x, y, z}, t] + vy[{x, y, z}, t] ny             [{x, y, z}, t] +
({0, 1, 0},0)                                    ({1, 0, 0},0)
ny[{x, y, z}, t] vy             [{x, y, z}, t] + vx[{x, y, z}, t] nx             [{x, y, z}, t] +
({1, 0, 0},0)
nx[{x, y, z}, t] vx             [{x, y, z}, t],
({0, 0, 0},1)                                    ({0, 0, 1},0)
nz             [{x, y, z}, t] + vz[{x, y, z}, t] nz             [{x, y, z}, t] +
({0, 0, 1},0)                                    ({0, 1, 0},0)
nz[{x, y, z}, t] vz             [{x, y, z}, t] + vy[{x, y, z}, t] ny             [{x, y, z}, t] +
({0, 1, 0},0)                                    ({1, 0, 0},0)
ny[{x, y, z}, t] vy             [{x, y, z}, t] + vx[{x, y, z}, t] nx             [{x, y, z}, t] +
({1, 0, 0},0)
nx[{x, y, z}, t] vx             [{x, y, z}, t]} == 0


I would like Mathematica to write the equation in vectorial form so that I can rearrange it and use vector identities to manipulate it.

Is there a way to do this?

-
Maybe this will be of help. –  b.gatessucks Feb 14 '13 at 11:11
Unfortunately that's not a solution to my problem. All the output of the Div, Grad, Curl operators is in component form, i.e. writing out all the components of the vector. This makes it almost impossible to read for more complex equations. I want to manipulate the equations in vector form. For example, I want to use identities such as Div[Cross[v,w]] = w.Curl[v]-v.Curl[w] without having to see it written out into components. –  Holger Schmitz Feb 14 '13 at 11:15
You can format the outputs as different forms. As InputForm, FullForm, OutputForm or TeXForm , ect. –  s.s.o Feb 14 '13 at 11:34
Related/possible duplicate: mathematica.stackexchange.com/questions/3242/… –  nikie Feb 14 '13 at 12:23
Yes, but you are asking to manipulate your vector fields algebraically. The unifying concept is that of the differential operator $d$ which maps a $k$-form to a $k+1$-form, along with the wedge product mapping an ordered pair of a $k$-form and $l$-form to a $k+l$-form. These purely algebraic constructs subsume all your differential operators and most of the identities you might have in mind. (A few more related operations on forms, such as the Hodge $*$, will complete the set of possible algebraic identities.) Abstractly, the smooth vector fields form a differential ring. –  whuber Feb 14 '13 at 18:18

The answer depends a lot on what you mean by "doing" vector calculus. You want results to be displayed without using component notation, and that's in general a difficult thing to achieve.

A prerequisite about doing completely symbolic vector calculus is to define the simplification rules. But even in "non-vector" algebra it's often hard to get Simplify to reach an expression that you would like to obtain. This certainly won't get any easier when going to vector calculus. Trying to go the completely symbolic route may then not be worth your time, especially since I'm guessing you're doing everything in three dimensions, not "n dimensions."

So instead I usually follow a different approach in which Mathematica doesn't derive theorems but instead verifies them. That looks less impressive, but it's a lot easier to do. In your case, it means taking a vector calculus identity lhs == rhs and feeding it to Simplify in hopes of getting the result True. By contrast, starting only with Simplify[lhs], your chances of arriving exactly at rhs will in general by slim to nil.

To verify vector calculus identities, it's typically necessary to define your fields and coordinates in component form, but if you're lucky you won't have to display those components in the end result.

To show some examples, I wasn't able to make up my mind if I should use the VectorAnalysis package or the new version 9 functions. So I'll include my own versions of these functions here, just for the cartesian coordinate case.

Then I'll define something analogous to your vector r, and a general vector field notation. My choice for the coordinate symbols is to define them in a Array using a string "x" as the name for all components, so that there's no potential confusion with existing symbols.

Clear[grad, curl, div, cross, r];

div[f_] := Inner[D, f, Array["x", Length[f]]]

curl[f_?VectorQ] :=
Inner[D, f.Normal[LeviCivitaTensor[Length[f]]], Array["x", Length[f]]]

cross[f_?VectorQ, g_?VectorQ] :=
Dot[g.Normal[LeviCivitaTensor[Length[g]]], f]

r = Array["x", 3];

field[name_, dimension_: 3] :=
Through[Array[name, dimension] @@ Array["x", dimension]]


Now check the vector identity that you mentioned in your comment:

Simplify[
div[cross[field[A], field[B]]] ==
field[B].curl[field[A]] - field[A].curl[field[B]]]


True

Next, check the continuity equation that motivated the question, by defining things in component form but displaying only the truth value of the identity:

Clear[jVec, ρ];

jVec[rVec_?VectorQ, t_] :=
E^(-(rVec.rVec/(1 + t)^2))/(Pi^(3/2) (1 + t)^4) rVec

ρ[rVec_?VectorQ, t_] :=
1/(Sqrt[Pi] Pi (1 + t)^3) Exp[-(rVec.rVec)/(1 + t)^2]

Simplify[
div[jVec[r, t]] + D[ρ[r, t], t] == 0
]


True

When you do look at the coordinate form of the expressions, they will aways look somewhat uglier of course.

-
Thanks Jens, that helped a lot. I have been using this technique in other contexts but thought I could get more readable results by doing the algebra in a fully symbolic way. –  Holger Schmitz Feb 19 '13 at 11:31

Our package VEST (Vector Einstein Summation Tools) is designed to do exactly this. It's described here and can be downloaded from github along with a tutorial.

I realize this question is now very old, but thought this answer might be helpful for others looking for similar capability.

-