A Simple Simulation of Custom Physical Interaction With Python and Matplotlib
Have fun with physics in this tutorial!
Join the DZone community and get the full member experience.Join For Free
Here, we are going to simulate some vector field (for instance, an electromagnetic field) in an N-d space.
Our plan is:
- To define our theoretical base (like arrays in numpy).
- To define the mechanism of point particles and interaction fields.
- To visualize an electric field.
- To visualize the movement of particles in an electromagnetic field.
You may also like: Vector Based Languages.
The basic element of any physical scene is a Vector. What do we need? Arithmetic operations on a vector, distance, module, and a couple of technical things. The vector we will inherit from List. This is how its initialization looks:
class Vector(list): def __init__(self, *el): for e in el: self.append(e)
Now, we can create a vector with
v = Vector(1, 2, 3)
Let's set the arithmetic operation addition:
class Vector(list): ... def __add__(self, other): if type(other) is Vector: assert len(self) == len(other), "Error 0" r = Vector() for i in range(len(self)): r.append(self[i] + other[i]) return r else: other = Vector.emptyvec(lens=len(self), n=other) return self + other
v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 + v2 >>> [3, 59, 26.2]
We similarly define all the arithmetic operations (the full code of the vector is at the end). Now, we need a distance function. I could simply make dist (v1, v2) — but this is not beautiful, so we will redefine the % operator:
class Vector(list): ... def __mod__(self, other): return sum((self - other) ** 2) ** 0.5
v1 = Vector(1, 2, 3) v2 = Vector(2, 57, 23.2) v1 % v2 >>> 58.60068258988115
We also need a couple of methods for faster vector generation and beautiful output. There is nothing tricky here, so here is the entire code of the Vector class.
Here, in theory, everything is simple — the point has coordinates, speed, and acceleration. In addition, it has mass and a set of custom parameters (for example, for an electromagnetic field, you can set charge).
Initialization will be the following:
class Point: def __init__(self, coords, mass=1.0, q=1.0 speed=None, **properties): self.coords = coords if speed is None: self.speed = Vector(*[0 for i in range(len(coords))]) else: self.speed = speed self.acc = Vector(*[0 for i in range(len(coords))]) self.mass = mass self.__params__ = ["coords", "speed", "acc", "q"] + list(properties.keys()) self.q = q for prop in properties: setattr(self, prop, properties[prop])
To move, immobilize and accelerate our point, we will write the following methods:
class Point: ... def move(self, dt): self.coords = self.coords + self.speed * dt def accelerate(self, dt): self.speed = self.speed + self.acc * dt def accinc(self, force): # Considering exerted force the point obtains acceleration self.acc = self.acc + force / self.mass def clean_acc(self): self.acc = self.acc * 0
Well done, the point itself is done.
We call a field of interaction an object that includes a set of all particles from space and exerts force on them. We will consider a special case of our universe, so we will have one custom interaction (of course, this is easy to expand). Declare the constructor and add a point:
class InteractionField: def __init__(self, F): # F - is a custom force, F(p1, p2, r), p1, p2 - points, r - distance inbetween self.points =  self.F = F def append(self, *args, **kwargs): self.points.append(Point(*args, **kwargs))
Now, the fun part is to declare a function that returns “tension” at that point. Although this concept refers to electromagnetic interaction, in our case, it is some abstract vector, along which we will move the point. In this case, we will have the property of the point q, in the particular case — the charge of the point (in general — whatever we want, even the vector). So, what is tension at point C? Something like this:
Electric intensity in point C is equal to the sum of the forces of all material points acting on some unit point.
class InteractionField: ... def intensity(self, coord): proj = Vector(*[0 for i in range(coord.dim())]) single_point = Point(Vector(), mass=1.0, q=1.0) # That's our "Single point" for p in self.points: if coord % p.coords < 10 ** (-10): # Check whether we compare coord with a point P where P.coords = coord continue d = p.coords % coord fmod = self.F(single_point, p, d) * (-1) proj = proj + (coord - p.coords) / d * fmod return proj
At this point, you can already visualize a vector field, but we will do it at the end. Now, let's take a step in our interaction.
class InteractionField: ... def step(self, dt): self.clean_acc() for p in self.points: p.accinc(self.intensity(p.coords) * p.q) p.accelerate(dt) p.move(dt)
For each point, we determine the intensity in these coordinates and then determine the final force on this particle:
Particle Movement and Vector Field Visualization
We finally reached the most interesting part. Let's start with ...
Modeling the motion of particles in an electromagnetic field
u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 2 + 0.1)) for i in range(3): u.append(Vector.randvec(2) * 10, q=random.random() - 0.5)
Actually, the coefficient k should be equal to some kind of billions (9 * 10 ^ (- 9)), but since it will be quenched by the time t -> 0, it is easier to make both of them normal numbers. Therefore, in our physics k = 300,000.
Next, we add ten points (2-dimensional space) with coordinates from 0 to 10 along each axis. Also, we give each point a charge from -0.25 to 0.25. Then, we run through a loop and draw points according to their coordinates (and traces):
X, Y = ,  for i in range(130): u.step(0.0006) xd, yd = zip(*u.gather_coords()) X.extend(xd) Y.extend(yd) plt.figure(figsize=[8, 8]) plt.scatter(X, Y) plt.scatter(*zip(*u.gather_coords()), color="orange") plt.show()
What should have happened is:
In fact, the plot will be completely random because the trajectory of each point is currently (in 2019) unpredictable.
Vector Field Visualization
We need to go through the coordinates by some step and draw a vector in each of them in the right direction.
fig = plt.figure(figsize=[5, 5]) res =  STEP = 0.3 for x in np.arange(0, 10, STEP): for y in np.arange(0, 10, STEP): inten = u.intensity(Vector(x, y)) F = inten.mod() inten /= inten.mod() * 4 # длина нашей палочки фиксирована res.append(([x - inten / 2, x + inten / 2], [y - inten / 2, y + inten / 2], F)) for r in res: plt.plot(r, r, color=(sigm(r), 0.1, 0.8 * (1 - sigm(r)))) # Цвет по хитрой формуле чтобы добиться градиента plt.show()
You should have gotten something like this:
You can lengthen the vectors themselves, replace * 4 with * 1.5:
Let's create a five-dimensional space with 200 points and an interaction that is not dependent on the square of the distance, but on the 4th degree.
u = InteractionField(lambda p1, p2, r: 300000 * -p1.q * p2.q / (r ** 4 + 0.1)) for i in range(200): u.append(Vector.randvec(5) * 10, q=random.random() - 0.5)
Now, all coordinates, speeds, etc. are defined in five dimensions. Now, let's model something:
velmod = 0 velocities =  for i in range(100): u.step(0.0005) velmod = sum([p.speed.mod() for p in u.points]) # Adding sum of modules of all the velocities velocities.append(velmod) plt.plot(velocities) plt.show()
This is a graph of the sum of all speeds at any given time. As you can see, over time they are slowly accelerating.
Well, that was a short instruction on how to simply simulate elementary physical stuff, thank you for your attention.
Published at DZone with permission of Lev Brovchenkov. See the original article here.
Opinions expressed by DZone contributors are their own.