# Plotting vectors using derivatives

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},
{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!

• what's ρ in your definition of V? – chris Apr 26 '15 at 7:35
• Even if you define V[x_] := x^2;, D[V,x] is zero, you should write D[V[x],x]=2x. – Mahdi Apr 26 '15 at 7:42
• This demonstration has instructive code: demonstrations.wolfram.com/GravitationalPotentialOfACuboid – ubpdqn Apr 26 '15 at 11:40
• 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! – crockpotveggies Apr 26 '15 at 14:58

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.

• 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...}...? – crockpotveggies Apr 26 '15 at 15:10
• @crockpotveggies: Yes, you call {gx,gy,gz} in Table (or other methods) as a function. – Mahdi Apr 26 '15 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.

• 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? – crockpotveggies Apr 26 '15 at 15:05
• @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. – Christoph Groß Apr 26 '15 at 16:14