This question stems from a previous question I asked here: https://mathematica.stackexchange.com/questions/80877/multiple-integral-of-equation-including-dx-prime?noredirect=1#comment219103_80877

An equation describing a gravitational field of a cuboid is given as:

$V(X,Y,Z) = -G\rho \int_{-D}^{D} \int_{-B}^{B} \int_{-L}^{L} [ \frac{{dx}'{dy}'{dz}'}{\sqrt{(X-{x}')^2+(Y-{y}')^2+(Z-{z}')^2)}} ]$

And in Mathematica:

V[X_, Y_, Z_, Len_, Br_, Dep_] := 
         1 / Sqrt[(X - x)^2 + (Y - y)^2 + (Z - z)^2],
         {x, -Len, Len},
         {y, -Br, Br},
         {z, -Dep, Dep}

My goal is to make a vector plot of the above function, but my problem is that in order to do so each vector is can represented as $\vec{g} = - \bigtriangledown V$, which is a bunch of partial derivatives. I've tried doing the following with no success:

   {Red, Arrowheads[Small], Arrow[{D[V, x], D[V, y], D[V, z]}]
   {X, 0, 10}, {Y, 0, 10}, {Z, 0, 10}]
 "Cuboid Length",
 {Len, 2, 20},
 "Cuboid Breadth",
 {Br, 2, 4},
 "Cuboid Depth",
 {Dep, 2, 4},
 Initialization :> (
   Len = 8; Br = 2; Dep = 2; ρ = 1)

Where I think I'm going wrong is not properly initializing the V function, since when I try to perform a partial derivative on it in a single input line, it returns as 0. I also tried a StreamPlot with no success.

What am I doing wrong, and how can I fix it? Is it even necessary to use partial derivatives to define each vector? Thanks to those who can help!

  • $\begingroup$ what's ρ in your definition of V? $\endgroup$
    – chris
    Apr 26, 2015 at 7:35
  • $\begingroup$ Even if you define V[x_] := x^2;, D[V,x] is zero, you should write D[V[x],x]=2x. $\endgroup$
    – Mahdi
    Apr 26, 2015 at 7:42
  • $\begingroup$ This demonstration has instructive code: demonstrations.wolfram.com/GravitationalPotentialOfACuboid $\endgroup$
    – ubpdqn
    Apr 26, 2015 at 11:40
  • $\begingroup$ Sorry for the late reply, just woke up. In the above, ρ = 1 represents density. It can remain 1 since it is considered uniform. Did not realize that I could write D[V[x]... that way! Wow not sure how I missed that demonstration, I actually searched through the demos and must have missed it! $\endgroup$ Apr 26, 2015 at 14:58

2 Answers 2


Assuming $G=\rho=1$, move the gradient operator into the integral:

$$ \begin{align} - \mathbf{g} = \nabla V(X,Y,Z) &= \nabla \int_{-D}^{D} \int_{-B}^{B} \int_{-L}^{L} \frac{{dx}'{dy}'{dz}'}{\sqrt{(X-{x}')^2+(Y-{y}')^2+(Z-{z}')^2)}} \\ &= \int_{-D}^{D} \int_{-B}^{B} \int_{-L}^{L} \nabla \left(\frac{1}{\sqrt{(X-{x}')^2+(Y-{y}')^2+(Z-{z}')^2)}}\right) {dx}'{dy}'{dz}' \\ \end{align} $$

Say you want $\mathbf{g}_x$:

$ \mathbf{g}_{x}(X,Y,Z,D,B,L)=\int_{-D}^{D} \int_{-B}^{B} \int_{-L}^{L} \left(\frac{X-x}{\left((X-x)^2+(Y-y)^2+(Z-z)^2\right)^{3/2}} \right){dx}'{dy}'{dz}' $

then you can define:

gx[X_, Y_, Z_, Len_, Br_, Dep_] := NIntegrate[
       (-x + X)/((-x + X)^2 + (-y + Y)^2 + (-z + Z)^2)^(3/2), 
       {x, -Len, Len}, {y, -Br, Br}, {z, -Dep, Dep}];
gx[1, 7, 3, 8, 2, 2]
(* 0.190816 *)

then find other components gy and gz.

  • $\begingroup$ Interesting, I didn't even realize that I could move the gradient operator inside the integral. Just to make sure I understand this, by using gx I can plug that into the arrow function I'm using in my question like... Arrow[{gx..., gy..., gz...}...? $\endgroup$ Apr 26, 2015 at 15:10
  • $\begingroup$ @crockpotveggies: Yes, you call {gx,gy,gz} in Table (or other methods) as a function. $\endgroup$
    – Mahdi
    Apr 26, 2015 at 17:44

You guessed right that one problem in your code is the missing definition of V. In the snippet you posted above, you defined a function V[X_, Y_, Z_, Len_, Br_, Dep_]. Later, you try to call it by writing D[V, x]. This does not work, because Mathematica distinguishes between V and V[X_, Y_, Z_, Len_, Br_, Dep_]. You can see this by looking into the down values. Here, I have used your definition for V and also made another one in order to illustrate my point:

V[X_, Y_, Z_] := -GravitationalConstant * \[Rho] * Integrate[
  1/Sqrt[(X - x)^2 + (Y - y)^2 + (Z - z)^2],
  {x, -\[Infinity], \[Infinity]},
  {y, -\[Infinity], \[Infinity]},
  {z, -\[Infinity], \[Infinity]}];

V[X_, Y_, Z_, Len_, Br_, Dep_] :=
-GravitationalConstant * \[Rho] * Integrate[
  1/Sqrt[(X - x)^2 + (Y - y)^2 + (Z - z)^2],
  {x, -Len, Len}, {y, -Br, Br}, {z, -Dep, Dep}]

The Mathematica function DownValues lets us inspect all definitions of an object.

DownValues[V] // InputForm

  HoldPattern[V[X_, Y_, Z_]] :>
    (-GravitationalConstant) * \[Rho] * Integrate[
      1/Sqrt[(X - x)^2 + (Y - y)^2 + (Z - z)^2],
      {x, -Infinity, Infinity}, 
      {y, -Infinity, Infinity},
      {z, -Infinity, Infinity}],
  HoldPattern[V[X_, Y_, Z_, Len_, Br_, Dep_]] :>
    (-GravitationalConstant) * \[Rho] * Integrate[
      1/Sqrt[(X - x)^2 + (Y - y)^2 + (Z - z)^2], 
      {x, -Len, Len}, {y, -Br, Br}, {z, -Dep, Dep}]

This means that the same function name with different arguments will be resolved as completely different objects in Mathematica. So when calling D[V, x], Mathematica looks for a definition of V, finds none that matches your call and leaves it unevaluated. Therefore, the differentiation of V with respect to x gives zero.

Next Mahdi has made a valid point: It appears to me that you do not really need the potential function at all. If you are only interested to calculate the vector field from a given density distribution \[Rho], then I would suggest you use his method and calculate the vector field directly (of course, if [\Rho] depends on spacial coordinates, it has to be kept inside the integral). Also, using NIntegrate instead of Integrate should be a lot faster for generating a plot. You could even think about parallelization of the computation as each vector is independent from the other ones.

  • $\begingroup$ DownValues is a function I was not familiar with. Let's say I want to become better at plotting in the future and run into a similar problem, is DownValues a function I can typically use to diagnose my problems? $\endgroup$ Apr 26, 2015 at 15:05
  • $\begingroup$ @crockpotveggies Well, that would depend on the situation, but of course there are issues you can diagnose with DownValues. One possible scenario: You can identify problems originating from old definitions. Even after you have deleted an input cell, definitions from that cell remain in the kernel until you clear them e.g. by calling ClearAll or you restart the kernel itself. $\endgroup$ Apr 26, 2015 at 16:14

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.