Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

For the non linear partial differential equation below, why are the Gradient, Laplacian and Divergence being evaluated to zero despite using the VectorAnalysis package.

Bo = 1/300;
\[Delta] = 10^-3;
\[Epsilon] = 10^-6;
r = 0;
m = 0.05;
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(h[x, y, t]\)\) + 
 Div[-h[x, y, t]^3 Bo Grad[h[x, y, t]] + 
   h[x, y, t]^3 Grad[
     Laplacian[h[x, y, t]]] + ((\[Delta] h[x, y, t]^3) Grad[
     h[x, y, t]])/(Bi h[x, y, t] + K1)^3 + 
   m (h[x, y, t]/(K1 + Bi h[x, y, t]))^2 Grad[
     h[x, y, t]]] + \[Epsilon]/(Bi h[x, y, t] + K1) + r \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\((
\*SubscriptBox[\(\[PartialD]\), \(x\)]
SuperscriptBox[\(h[x, y, t]\), \(2\)], \(K1 + Bi\ h[x, y, t]\)]\ 
\*SuperscriptBox[\(h[x, y, t]\), \(3\)])\)\)

I am doing this because I'd like to evaluate the magnitude of each of these terms at a certain time step.

Here is a simple sample code without the Vector Analysis package that works fine:

I solve a PDE using NDSolve and then plot each term (I call each of the two terms Eq1 and Eq2 here) at a certain time step.

sol = u /. 
  NDSolve[{D[u[t, x], t] == 
      0.5 D[u[t, x], x, x] + u[t, x] D[u[t, x], x],  
     u[t, -Pi] == u[t, Pi] == 0 , u[0, x] == Sin[x]}, 
    u, {t, 0, 2}, {x, -Pi, Pi}][[1]];

Eq1 = 0.5 D[u[t, x], x, x];
Magnitude1[t_, x_] = Eq1 /. u -> sol
Plot[Magnitude1[1.5, x], {x, -\[Pi], \[Pi]}]

enter image description here

Eq2 = u[t, x] D[u[t, x], x];
Magnitude2[t_, x_] = Eq2 /. u -> sol
Plot[Magnitude2[1.5, x], {x, -\[Pi], \[Pi]}]

enter image description here

share|improve this question
Did you do SetCoordinates[Cartesian[x,y,z]]? It does the trick for me :-) – chris Oct 16 '12 at 20:23
@chris Boy, do I feel like an idiot. Should I delete this question? – drN Oct 16 '12 at 20:26
may be you should check if its already been asked. otherwise its useful idiocy ;-) – chris Oct 16 '12 at 20:27
@chris I checked. This hasn't been asked before. If you feel it's useful idiocy, I'll let it be. I am thinking however, I should receive negative points for dumb-itude such as this. – drN Oct 16 '12 at 20:28
Interestingly, 6 of the 8 posts involving SetCoordinates were asked by you... :-) – Sjoerd C. de Vries Oct 16 '12 at 22:20
up vote 4 down vote accepted

Operators like Div require the specifications of a coordinate system. By default, this coordinate system is cartesian and in terms of variables Xx, Yy and Zz. If you specify the vector field in terms of x, y and z, its derivatives with respect to the coordinates will be zero and the Div operator will correctly evaluate to zero.

The coordinate system can be changed manually. The following code will yield 3 as intended.

share|improve this answer
"By default" would suggest that I don't need to specify SetCoordinates[Cartesian[x,y,z]]; , but I had to. – drN Oct 17 '12 at 18:54
You don't need to. Div[{Xx,Yy,Zz}] will work without changing the coordinate system. – Marcks Thomas Oct 17 '12 at 18:56
@MarckThomas but it didn't. I had to specify a cartesian coordinate system before I tried a Div or Grad ... – drN Oct 17 '12 at 19:00
Interesting. The VectorAnalysis reference reads: "The initial default coordinate system is Cartesian, with coordinate variables Xx, Yy, and Zz." What does your Coordinates[] initially give? – Marcks Thomas Oct 17 '12 at 20:29
After loading VectorAnalysis, Coordinates[] gives me an error message: Coordinates::shdw: Symbol Coordinates appears in multiple contexts {VectorAnalysis,Global}; definitions in context VectorAnalysis` may shadow or be shadowed by other definitions. >> – drN Oct 17 '12 at 20:31

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.