How to calculate the gradient (vector) of a vector field? First question from me. I hope it's clear enough.

I'm trying to make a physics simulator of magnets, magnetic interactions and such.
I've followed the formulas in Wikipedia pages, and some first results are great (same newton force value between math model and reality test).

Now I have to translate everything in a software coding; I'm using c# inside grasshopper inside Rhinoceros.

I can calculate the magnetic field vector "B" of a magnetic dipole vector "m" with this formula: wiki 4 , 
Pic 1

Anyway,
here: wiki 1 , here: wiki 2 and here: wiki 3 it says that the force is this:
Pic 2

That is the gradient of the dot product of "m"(the magnetic dipole moment vector) and "B"(the local magnetic field vector).

I've understood that the gradient is somehow "like" the slope of a function (which you can "manually" calculate by evaluating two points really close and divide their Y difference by their X distance).

For every sample point I'd like to have one gradient vector (of the local resulting magnetic field) to multiply (dot product) with whatever dipole is passing by in that iteration.
(Does this make sense? I'm not native English...)

Can I calculate the gradient by evaluating 2 points near the target location in the same way?

How?

Other ideas?

Thanks in advance.
First of all, since the dipole $m$ on which the force acts is constant, the formula simplifies to
$$
F=\nabla(m\cdot B) = m^TJ_B = J_B^T m,
$$
where $J_B$ is the Jacobian matrix. See also here.

If you want to see the reason why, just work with coordinates and you find
$$
[\nabla(m\cdot B)]_i = \frac{\partial}{\partial x_i} \sum_{j=1}^n m_j B_j
= \sum_{j=1}^n m_j \frac{\partial B_j}{\partial x_i} 
= m^T J_B.
$$

Regarding the question of how to compute $J_B$, there are several approaches:


if $B$ has a specific closed form expression, you can of course use it to compute explicitly its gradient;
you can use finite differences, as you mentioned;
you can use automatic differentiation to compute (a numeric approximation of) the gradient at the same time as you compute the field itself.
No,
$$\operatorname{grad}(\vec{m}\cdot\vec{B}) \neq m \cdot\operatorname{grad}(\vec{B})$$
But rather
$$\begin{align}
\operatorname{grad}(\vec{m}\cdot\vec{B})_i&amp;=\partial_i (m_kB_k)\\
&amp;=(\partial_im_k)B_k+(\partial_iB_k)m_k\\
&amp;=J(\vec{m})_{ik}B_k+J(\vec{B})_{ik}m_k\\
&amp;=(J(\vec{m})\cdot\vec{B})_i+(J(\vec{B})\cdot\vec{m})_i
\end{align}$$
Where $J(\vec{m})$ is the Jacobian of $\vec{m}$ and $J(\vec{B})$ is the jacobian of $\vec{B}$. And if $\vec{m}$ is a constant vector, then you have that
$$\operatorname{grad}(\vec{m}\cdot\vec{B})=J(\vec{B})\cdot\vec{m}$$