

Oh yeah, someone finally asked something about my field.
So, maybe COMSOL is what you’re looking for, but if you want the cool, smooth and amazing stuff of the pros:
Time to learn numerical calculus, like verlet integration, maybe some calculus for turning potential fields into force (maybe some forward differentiation with dual numbers, just like sympy does hahaha, you can use that too). If it’s a hairy equation, that takes some hypergeometric, bessel, etc, scipy has all the special equations.
I think verlet + sympy for turning the potential fields into the force fields at the particle spots will do amazing, if you don’t want paper ready simulations. Like, pos[n, 3], v[n, 3], a[n, 3] (maybe some torque if you want) would be almost all you need.
For visualization, matplotlib has a lot of stuff, you can make animations and videos.
If you want to simulate the fields itself, you can use finite differences or finite element methods.
If you need a lot of performance, I recommend learning Julia or the GODLY Fortran 2008 (fprettify, fortls, gfortran will help to build a modern Fortran program that you can call from python using cdll and ctypes with the ieee c libs for fortran).






Edit: Corrected the integration
This is a simulation I made for a ball falling (the principles can be used everywhere):
import numpy as np import matplotlib.pyplot as plt def calc_acceleration(pos, v): return np.array([0, 0, -9.8]) def integrate(pos, v, a_old, delta_t): pos_new = pos + v * delta_t + 0.5 * a_old * delta_t * delta_t a_new = calc_acceleration(pos_new, v) v_new = v + 0.5 * (a_old + a_new) * delta_t if pos_new[2] < 0.0: v_new[2] = -v_new[2] return pos_new, v_new, a_new # Constants delta_t = 0.01 t_max = 10 # Time evolution times: np.ndarray = np.arange(0, t_max, delta_t) # Initial state v = np.array([15.0, -3.0, 15.0]) pos = np.array([0.0, 0.0, 7.0]) positions = np.zeros([times.size, 3]) a_old = calc_acceleration(pos, v) # Simulating for i in range(times.size): pos, v, a_old = integrate(pos, v, a_old, delta_t) positions[i, :] = pos # Plotting fig = plt.figure() ax = fig.add_subplot(projection="3d") ax.scatter(positions[:, 0], positions[:, 1], positions[:, 2]) plt.show()Very simple, you can see a ball falling and getting up again from that starting point, so you just see at least some starting point on how to do this stuff, the basics are not that difficult. But every simulation is like this, just more complicated hahaha.