Fields and a simulation
The field of 20,000 random charges by Barnes-Hut tree walk, compared with the exact direct sum; the opening angle theta trades accuracy for speed.
import numpy as np
import tree_d_3dfrom tree_d_3d import FlatQuadtree, Quad, direct_field
rng = np.random.default_rng(0)n = 20_000pos = rng.uniform(-1.0, 1.0, (2, n)) # component-major: pos[0] is xcharge = rng.choice([-1.0, 1.0], n)others = np.arange(n) # each particle excludes its own chargereference = direct_field(pos, pos, charge, softening=1e-3, exclude=others)[1]tree = FlatQuadtree(pos, charge, Quad(0.0, 0.0, 1.0))for theta in (1.0, 0.5, 0.25): field = tree.evaluate(pos, exclude=others, theta=theta, softening=1e-3)[1] error = np.abs(field - reference).max() / np.abs(reference).max() print(f"theta = {theta}: max. relative error {error:.1e}")theta = 1.0: max. relative error 8.6e-03theta = 0.5: max. relative error 1.5e-04theta = 0.25: max. relative error 6.5e-06For line charges (kernel="log", the 2D Poisson kernel) the fast multipole method computes the fields of all particles at once in O(N).
lines = FlatQuadtree(pos, charge, Quad(0.0, 0.0, 1.0), kernel="log", order=16)reference = direct_field(pos, pos, charge, exclude=others, kernel="log")[1]field = lines.fmm(theta=0.5)[1]print( f"FMM: max. relative error {np.abs(field - reference).max() / np.abs(reference).max():.1e}")FMM: max. relative error 3.1e-10An ElectrostaticSimulation advances ChargedBody particles; a cloud of like charges expands, converting potential into kinetic energy.
bodies = [tree_d_3d.ChargedBody(x, y) for x, y in rng.normal(0.0, 0.1, (500, 2))]sim = tree_d_3d.ElectrostaticSimulation( bodies, Quad(0.0, 0.0, 1.0), k_e=1.0, softening=1e-3, dt=1e-4)energy = [sim.get_total_energy()]for _ in range(200): sim.step() energy.append(sim.get_total_energy())print(f"relative energy change {abs(energy[-1] - energy[0]) / abs(energy[0]):.1e}")relative energy change 1.3e-03