# Simulate / plot electrostatic field

I have something like this: For given voltage (9 V red wire, 0 V blue wire) I want to simulate electrostatic field.

I already know, for example that in the point (10.5;0) there are 10 V and etc. I want to simulate and plot electrostatic field for future comparison and for my report.

Is there any software for doing that? I like to work with Maple, so if there is any toolboxes, I would be glad to hear about them.

For slightly different position plot looks awful now:

• For this sort of thing, I use Comsol Multiphysics. comsol.com
– vnb
May 19, 2014 at 2:33

As I said earlier in my comment, Comsol is pretty useful for these sort of things. I took your problem, and this is the plot that I got. I hope it's useful.

• Can you share some quick guide how did you do this? In which mode and etc? May 19, 2014 at 7:47
• In AC/DC mode I use electrostatics. For materials I used copper, and as boundaries the bottom one as ground ($V=0$) and the top one with electric potential $V=9$. Also, to get the field lines add streamlines to your solution.
– vnb
May 19, 2014 at 8:18
• @vnb: Hi, I have some problems to perform a simulation with COMSOL and actually I don't know many people over here that know how to use this software. Honestly I couldn't find any useful document or example/tutorial model to resolve my problems. I really need this. I noticed you've used this software to answer a questions on SE. Could you please help me resolve my problem? I really appreciate that.Thank you! (andx@chmail.ir) Jun 2, 2014 at 17:17
• "Buy this expensive program" isn't really a good stackexchange answer. How about explaining some of the physics? Also, the spacing of the field lines is a bit funky, they should be fairly evenly spaced away from corners, or have a well-behaved gradient.
– uhoh
Oct 19, 2017 at 16:16

The brine tank is a historical and excellent way to find solutions to the Laplace equation! This was used to design things like cyclotron magnets and electrostatic lenses for charged particles before computers were fast enough.

Image from Tanabe and Yamada 1958; Science reports of the Research Institute, Tohoku University, Series A, 10, 133-174

The simplest way to solve it numerically is to do "Jacobi relaxation". Read more about it here for example. In cartesian coordinates, you make a numerical grid and set the points that correspond to the surface of the conductors to their potential, and set the potential to the boundary to some intermediate value. Below I've chosen -1, +1 and 0. The boundary needs to be far enough away that it doesn't affect your region of interest much. This is an approximate technique.

Then you iterate. For each pass, you make a new potential that is the average of the +/-x and +/-y nearest neighbors, four in total. You update only the points that are not on the conductor or boundary. After thousands of iterations the potential approaches a solution to the Laplace equation.

If you want a better solution, make the spacing smaller and the boundary farther away. If it's too slow (frequently it is in 3D) then you can use

• Multigrids; solve on a coarse (fast) grid, then interpolate to a fine grid and iterate a little longer

• Successive Overrelaxation; "overshoot" the averaging by changing the values by an amount larger than 100% of what you would have with simple averaging. Try values like 110% to 150%, it can go unstable if they are too large.

If you have any questions or need clarification, leave a comment. I used this technique in this question, where I saved intermediate snapshots and turned them into a GIF.

Here is an example in Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve

d = 0.25

x = np.arange(0, 35+d, d)
y = np.arange(0, 48+d, d)

X, Y = np.meshgrid(x, y)

phi = np.zeros_like(X)

isneg = np.zeros_like(X, dtype=bool)
ispos = np.zeros_like(X, dtype=bool)
bound = np.zeros_like(X, dtype=bool)

isneg[(X>=9)*(X<=19)*(Y>=10)*(Y<=12)] = True

ispos[(X>=10)*(X<=20)*(Y>=28)*(Y<=30)] = True
ispos[(X>=18)*(X<=20)*(Y>=20)*(Y<=28)] = True

bound[:,0] = bound[:,-1] = bound[0,:] = bound[-1,:] = True

phi[isneg] = -1.0
phi[ispos] = +1.0
phi[bound] =  0.0

updateme = np.ones_like(X, dtype=bool)
updateme[isneg] = updateme[ispos] = updateme[bound] = False

kernel = 0.25*np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], dtype=float)

phi0 = phi.copy()  # keep the original handy
keepers = 0, 100, 500, 1000, 2000, 9999
keep = []
for i in range(10000):
if i in keepers:
keep.append(phi.copy())
print i
phi2 = convolve(phi, kernel)
phi[updateme] = phi2[updateme]

plt.figure()
for i, (thing, n) in enumerate(zip(keep, keepers)):
plt.subplot(2, 3, i+1)
plt.imshow(thing, origin='lower')
plt.title("n="+str(n))
plt.show()