Plotting vectors using derivatives

Question

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_] := 
      -GravitationalConstant*ρ* 
      Integrate[
         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:

Manipulate[
 Graphics3D[
  Table[
   {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!

Answer 1

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.

Answer 2

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.