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

  • 5
    $\begingroup$ Did you do SetCoordinates[Cartesian[x,y,z]]? It does the trick for me :-) $\endgroup$
    – chris
    Oct 16, 2012 at 20:23
  • $\begingroup$ @chris Boy, do I feel like an idiot. Should I delete this question? $\endgroup$
    – dearN
    Oct 16, 2012 at 20:26
  • $\begingroup$ may be you should check if its already been asked. otherwise its useful idiocy ;-) $\endgroup$
    – chris
    Oct 16, 2012 at 20:27
  • $\begingroup$ @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. $\endgroup$
    – dearN
    Oct 16, 2012 at 20:28
  • 2
    $\begingroup$ Interestingly, 6 of the 8 posts involving SetCoordinates were asked by you... :-) $\endgroup$ Oct 16, 2012 at 22:20

1 Answer 1


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.

  • $\begingroup$ "By default" would suggest that I don't need to specify SetCoordinates[Cartesian[x,y,z]]; , but I had to. $\endgroup$
    – dearN
    Oct 17, 2012 at 18:54
  • $\begingroup$ You don't need to. Div[{Xx,Yy,Zz}] will work without changing the coordinate system. $\endgroup$ Oct 17, 2012 at 18:56
  • $\begingroup$ @MarckThomas but it didn't. I had to specify a cartesian coordinate system before I tried a Div or Grad ... $\endgroup$
    – dearN
    Oct 17, 2012 at 19:00
  • $\begingroup$ 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? $\endgroup$ Oct 17, 2012 at 20:29
  • $\begingroup$ 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. >> $\endgroup$
    – dearN
    Oct 17, 2012 at 20:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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