# Numerical simulation of magnetic dipole in a inhomogeneous magnetic field

Goal:

I want to use Python to illustrate how a magnetic dipole with magnetic moment m2 moves in a non-homogeneous magnetic field in a 2D-Plane. This field is generated by another magnetic dipole with moment m1. However, I get results that are not plausible.

For the calculation of the trajectory of m2 I use the equations from wiki - force between magnets

Theoretical Background:

Given is a magnetic dipole with constant moment m1 and constant position r1!. m1 generates a magnetic flux density

$$\mathbf {B}_1 (\mathbf {r} )={\frac {\mu _{0}}{4\pi }}{\frac {3\mathbf {\hat {r}} (\mathbf {\hat {r}} \cdot \mathbf {m} )-\mathbf {m} }{|\mathbf {r} |^{3}}}\qquad (1)}$$

where r is a vector from the center of the magnetic dipole to the location where the magnetic field is measured. In this non-homogeneous magnetic field a magnetic dipole with moment m2 and position r2 is placed. Its magnetic field can be neglected B2 = 0.

Only the dynamics of m2 shall be considered. According to the linked Wikipedia page, m2 will be subjected to a torque D and a force F. The torque is neglected D=0. The magnetic moment m2 is therefore constant.

Due to the acting force, the position vector r_2(t) depends on time. Its time evolution is obtained from the solution of the equation of motions, which is my objective:

$$\frac{d^2}{dt^2}\mathbf{r}_2 (t) = \binom{\dot{r}_{2x}}{\dot{r}_{2y}} = \mathbf{F}$$

According to Wikipedia, m1 exerts the following force on m2: (r = r2 - r1, r = ||r||_2):

\begin{align} \mathbf{F} &= \nabla |\mathbf{m}_2\cdot \mathbf{B}_1| \qquad(3) \\\\ &= {\frac {1}{ r^{5}}}\left[(\mathbf {m} _{1}\cdot \mathbf {r} )\mathbf {m} _{2}+(\mathbf {m} _{2}\cdot \mathbf {r} )\mathbf {m} _{1}+(\mathbf {m} _{1}\cdot \mathbf {m} _{2})\mathbf {r} -{\frac {5(\mathbf {m} _{1}\cdot \mathbf {r} )(\mathbf {m} _{2}\cdot \mathbf {r} )}{r^{2}}}\mathbf {r} \right] \qquad (4) \\ &= \frac{\mathbf{m}_1 \mathbf{r}}{r^5} \mathbf{m}_2 + \frac{\mathbf{m}_2\mathbf{r}}{r^5}+\frac{\mathbf{m}_1\mathbf{m}_2}{r^5}\mathbf{r} - \frac{5(\mathbf{m}_1 \mathbf{r})(\mathbf{m}_1\mathbf{r})}{r^7}\mathbf{r}\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + c_3 (\mathbf{r}_2 - \mathbf{r}_1) - c_4 (\mathbf{r}_2 - \mathbf{r}_1) \qquad (5)\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{r}_1 (c_4 - c_3) + \mathbf{r}_2 (c_3 - c_4)\qquad (6)\\\\ &= \mathbf{c}_1 + \mathbf{c}_2 + \mathbf{c} + d\cdot\mathbf{r}_2, \quad \mathbf{c}_1, \mathbf{c}_2, \mathbf{c}\in \mathbb{R}^2, d\in\mathbb{R} \qquad (7) \end{align}

The introduced variables c1,c2, c, d should be self-explanatory. Now we can define 4 ODEs and solve them numerically with the Runge-Kutta method:

\begin{align} \dot{r}_{2x} &= v_{2x} \qquad (8)\\\\ \dot{r}_{2y} &= v_{2y} \qquad (9)\\\\ \dot{v}_{2x} &= c_{1x} + c_{2x} + c_x + d\cdot r_{2x} \qquad (10)\\\\ \dot{v}_{2y} &= c_{1y} + c_{2y} + c_y + d\cdot r_{2y} \qquad (11)\\\\ \end{align}

What to expect?

An attractive force should act on the magnet m2, so that

$$t \gg 1: \mathbf{r}_2 = \mathbf{r}_1$$

My Results

m2 is repelled instead of being attracted.

Probably I made a mistake in the implementation. But I have already checked it several times and found no error... Does anyone of you know what went wrong

Python 3 Code

import numpy as np
import matplotlib.pyplot as plt

"""
Important quantites that will influence the dynamics

m1:  Magnetic Moment of the dipole that creates a constant magnetic field
m2:  Magnetic Moment of a dust particle
r1:  Position of m1: r1(t) = const
r2:  Position of m2: r2(t) will depend on time
"""

# Constants
m1 = np.array([0, 1], float)
m2 = np.array([0, -1], float)
r1 = np.array([0, 0], float)

# Initial conditions of the of m2 R = [r2_x, r2_y, v2_x, v2_y]
R = np.array([0, 3, 0, 0], float)

def distance(r1, r2, n):
"""
computes the distance ||.||_2 between position r1 and r2
"""

return np.linalg.norm(r1 - r2) ** (n)

def pre_factors(r1, r2, m1, m2):
"""
Computes the prefactors in the equation of motion
d^2 r2 / dt^2 = F = c1 + c2 + c + r2 * d   r2,c1,c2,c € R^2, d € R
"""

r = r2 - r1
l = distance(r1, r2, 5)

c1 = ( np.dot(m1, r) * m2 ) / l
c2 = ( np.dot(m2, r) * m1 ) / l

c3 = np.dot(m1, m2) / l
c4 = ( 5 * np.dot(m1, r) * np.dot(m2, r) ) / distance(r1, r2, 7)

c = (c4 - c3) * r2
d = c3 - c4

return c1, c2, c, d

def f(r1, R, m1, m2):
"""
now we can compute the righthand side of the four 1st order ODEs:
d/dt r2_x = v2_x
d/dt r2_y = v2_y
d/dt v2_x = c1_x + c2_x + c_x + d*r2_x
d/dt v2_y = c1_y + c2_y + c_y + d*r2_y
"""

r2x = R
r2y = R
v2x = R
v2y = R

c1, c2, c, d = pre_factors(r1, R[:2], m1, m2)

f_r2x = v2x
f_r2y = v2y
f_v2x = c1 + c2 + c + d * r2x
f_v2y = c1 + c2 + c + d * r2y

return np.array([f_r2x, f_r2y, f_v2x, f_v2y], float)

# Constants for the Simulation
a = 0.0  # t_0
b = 20.0  # t_end
N = 1000 # number of iterations
h = (b-a) / N # time step, e.g. h = (100 - 0) / 1000 = 0.1s

# Create lists to store the computed values
tpoints = np.arange(a,b+h,h)

xpoints, ypoints = [], []
vxpoints, vypoints = [], []

for dt in tpoints:
#print(R)
xpoints.append(R)
ypoints.append(R)
vxpoints.append(R)
vypoints.append(R)

# Runge-Kutta-4th Order Method
k1 = dt * f(r1, R, m1, m2)
k2 = dt * f(r1, R + 0.5 * k1, m1, m2)
k3 = dt * f(r1, R + 0.5 * k2, m1, m2)
k4 = dt * f(r1, R + k3, m1, m2)
R += (k1 + 2*k2 * 2*k3 + k4)

plt.style.use('seaborn')
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))
axes.plot(tpoints, xpoints, label = '$$r_{x_2}(t)$$')
axes.plot(tpoints, ypoints, label = '$$r_{y_2}(t)$$')

axes.legend()

axes.plot(tpoints, vxpoints,'r', label = '$$v_{x_2}(t)$$', alpha = 0.5)
axes.plot(tpoints, vypoints, 'k',label = '$$v_{y_2}(t)$$', alpha = 0.5)
axes.legend()

fig.suptitle("Results for: $$r_1 = (0,0), r_2 = (0,3), m_1 = (0,1), m_2 = (0,-1)$$")
#plt.savefig("2.pdf" , format='pdf',  bbox_inches="tight")
fig.tight_layout()

• Why do you think the two dipoles should be attracted to each other? It depends on their orientation, in general the force can be attraction, repulsion, or zero. Feb 5 at 0:04