0
$\begingroup$

I am trying to calculate and visualize the electric field inside and outside of a parallel plate capacitor by assuming a uniform distribution of point charges on each plate and by adding up the Coulomb-fields directly (I want this approach for educational purposes) with the following Python program.

Mathematically speaking, the electric field $\vec{E}(X)$ of a point charge is given by Coulombs law, i.e.

$$ \vec{E}(X) = \frac{1}{4\pi\epsilon_0} \frac{q}{r^3} \cdot \vec{r} $$

where $X$ is the point where you measure the field strength and $\vec{r}$ is the vector pointing from the location $P$ of the point charge to $X$. $r$ denotes the absolute value of $\vec{r}$. $q$ denotes the charge of the point charge and $\epsilon_0$ is a constant.

Now I assume a uniform distribution of two parallel plates with distance $d$ between the plates. On one plate with positive charges on the other with negative ones.

The total electric field at point $X$, is then given by simple superposition of the Coulomb fields, i.e.

$$ E_{\mathrm{total}} = \vec{E}(X) = \frac{q}{4\pi\epsilon_0} \sum_i \frac{1}{r^3} \cdot \vec{r} \propto q \sum_i \frac{1}{r^3} \cdot \vec{r} $$

Since my program should only demonstrate the principle I set $\frac{1}{4\pi\epsilon_0}$ to $1$ and choose $q=1$ for the positive plate and $q=-1$ for the negative one.

To translate the notation to my program I denote the position of a point charge by the vector $\vec{a}$ with coordinates $a_0,a_1,a_2$, and introduce a cartesian coordinate system with coordinates $x$ (to the right), $y$ (to the top) and $z$ (into the plane of view).

Then I can write the field of one charge in componentes as:

$$ \begin{align} E_x &= q\cdot(x-a_0)/r^3 \\ E_y &= q \cdot (y-a_1)/r^3 \end{align} $$

where

$$ r^3 = ((x-a_0)^2 + (x-a_1)^2 + (z-a_2))^{1.5} $$

I look only at the plane $z = 0$ and since I draw in 2D (projection to the plane $z = 0$) I omit $E_z$ (but I am not sure, if this introduces a mistake).

Then I use a loop to equally distribute the charges on my planes and calculate the total field as described above.

# Imports
import matplotlib 
from numpy import *
from pylab import *
import matplotlib.colors as colors


def Cap(N = 10,plus=1,minus=1,app=''):
    # Aspect ratio
    aspekt='auto'

    ## Field between the plates
    xMin,xMax = -0.15,0.15
    yMin,yMax = -10,10
    zMin,zMax = -10,10

    ## Field outside
    # xMin,xMax = -3,3
    # yMin,yMax = -3,3


    ## Plate distance

    d = 0.1

    # Case of plates very far apart
    # d = 50
    # xMin,xMax=-d*1.25,d*1.25
    # yMin,yMax=-d*1.25,d*1.25
    # aspekt = 'equal'

    ## Grids
    xG,yG = meshgrid(linspace(xMin,xMax,18),linspace(yMin,yMax,18))
    xGM,yGM = meshgrid(linspace(xMin,xMax,500),linspace(yMin,yMax,500))

    # Class that represents a point charge
    class charge:
        def __init__(self, q, pos):
            self.q=q
            self.pos=pos

    # Coulomb Law (I only look a the plane z = 0)
    def E_point_charge(q, a, x, y):
        z = 0
        r3 = ((x-a[0])**2+(y-a[1])**2 + (z-a[2])**2)**(1.5)
        return q*(x-a[0])/r3,q*(y-a[1])/r3


    # Total electric field
    def E_total(x, y, charges):
        Ex, Ey=0, 0
        for C in charges:
            E=E_point_charge(C.q, C.pos, x, y)
            Ex=Ex+E[0]
            Ey=Ey+E[1]
            #print(Ex)
        return [Ex, Ey]


    ## Initialize

    close('all')
    fig = figure()

    chargesPlus=[]
    chargesMinus=[]


    ## Generate charges (1/N**2 to have for every N the same total charge on each plate)
    for i in arange(yMin*0.9,yMax*0.9,0.9*(yMax-yMin)/N):
        for j in arange(zMin*0.9,zMax*0.9,0.9*(zMax-zMin)/N):
            if plus != 0:
                chargesPlus.append(charge(1/N**2,[d/2,i,j]))
                print("Charge + at %s,%s" %(i,j))
            if minus != 0:
                chargesMinus.append(charge(-1/N**2,[-d/2,i,j]))
                print("Charge - at %s,%s" %(i,j))
    charges = chargesPlus + chargesMinus



    # 
    E_totalX,E_totalY = E_total(xG,yG,charges)
    EE = sqrt(E_totalX**2 + E_totalY**2)
    E_DirX,E_DirY= E_totalX/EE,E_totalY/EE

    # 
    E_totalXM,E_totalYM = E_total(xGM,yGM,charges)
    EEM = sqrt(E_totalXM**2 + E_totalYM**2)

    # Plots
    ax = subplot(111)
    fig.add_subplot(ax)

    colormap='Spectral_r'

    # Colormap
    I = ax.imshow(EEM,extent=[np.min(xGM),np.max(xGM),np.min(yGM),np.max(yGM)],cmap=colormap, norm=colors.LogNorm(vmin=EEM.min(),vmax=EEM.max()),alpha=0.5,origin='lower',interpolation='bilinear')

    # Vector field
    Q = ax.quiver(xG,yG,E_DirX,E_DirY,color='gray',alpha=0.7,width=0.002,scale_units='inches',scale=4)

    # Colorbar
    fig.colorbar(I,extend='max',orientation='horizontal')


    # Plot point charges
    for C in charges:
        if C.q>0:
            plot(C.pos[0], C.pos[1], 'bo', ms=10*sqrt(C.q))
        if C.q<0:
            plot(C.pos[0], C.pos[1], 'ro', ms=10*sqrt(-C.q))

    # Axis sessings
    xlabel('$x$')
    ylabel('$y$')
    ax.set_aspect(aspekt)

    ax.set_xlim(xMin, xMax)
    ax.set_ylim(yMin, yMax)

    savefig("Capacitor_N_%05d_%app.pdf" %(N,app),dpi=fig.dpi,papersize='a4')
    show()


Cap(200)

However the computation time is very long, so is there better approach to do this under the constraint to use the direct summation of Coulomb fields (because of educational purposes) and to use Python as a tool.

The most time consuming step seems to be the calculation of the total field for the colormap plot.

Edit

I tried coarser grid for the colormap and then tried to use the timeit module to find the most time consuming part, but this doesn't seem to report reliable times. Maybe because of asynchronous processes. To show the problem I did a screencast for $N = 70$ here: https://www.pastefile.com/mBh38p

As you can see the plotting itself seems to be the problem. (I used a old i7-3520M CPU in my screencast)

The modified code looks like:

#!/usr/bin/env python3

# Imports
import matplotlib 
from numpy import *
from pylab import *
import matplotlib.colors as colors

from timeit import default_timer as timer
from datetime import timedelta


def Cap(N = 10,plus=1,minus=1,app=''):
    print('### N = %s ###' %N)

    # Aspect ratio
    aspekt='auto'

    ## Field between the plates
    xMin,xMax = -0.15,0.15
    yMin,yMax = -10,10
    zMin,zMax = -10,10

    ## Field outside
    # xMin,xMax = -3,3
    # yMin,yMax = -3,3


    ## Plate distance

    d = 0.1

    # Case of plates very far apart
    # d = 50
    # xMin,xMax=-d*1.25,d*1.25
    # yMin,yMax=-d*1.25,d*1.25
    # aspekt = 'equal'

    ## Grids
    xG,yG = meshgrid(linspace(xMin,xMax,18),linspace(yMin,yMax,18))
    xGM,yGM = meshgrid(linspace(xMin,xMax,40),linspace(yMin,yMax,40))

    # Class that represents a point charge
    class charge:
        def __init__(self, q, pos):
            self.q=q
            self.pos=pos

    # Coulomb Law (I only look a the plane z = 0)
    def E_point_charge(q, a, x, y):
        z = 0
        r3 = ((x-a[0])**2+(y-a[1])**2 + (z-a[2])**2)**(1.5)
        return q*(x-a[0])/r3,q*(y-a[1])/r3


    # Total electric field
    def E_total(x, y, charges):
        Ex, Ey=0, 0
        for C in charges:
            E=E_point_charge(C.q, C.pos, x, y)
            Ex=Ex+E[0]
            Ey=Ey+E[1]
            #print(Ex)
        return [Ex, Ey]


    ## Initialize

    close('all')
    fig = figure()

    charges = []    

    ## Generate charges (1/N**2 to have for every N the same total charge on each plate)
    print("Generate charges")
    start = timer()
    n = 0
    for i in arange(yMin*0.9,yMax*0.9,0.9*(yMax-yMin)/N):
        for j in arange(zMin*0.9,zMax*0.9,0.9*(zMax-zMin)/N):
            if plus != 0:
                charges.append(charge(1/N**2,[d/2,i,j]))
                #print("Charge + at %s,%s" %(i,j))
                #print(n)
                n += 1
            if minus != 0:
                charges.append(charge(-1/N**2,[-d/2,i,j]))
                #print("Charge - at %s,%s" %(i,j))
                #print(n)
                n +=1
    end=timer()
    print(timedelta(seconds=end-start))

    print("Calculating E_total for Quiver Plot")

    start=timer()
    E_totalX,E_totalY = E_total(xG,yG,charges)
    EE = sqrt(E_totalX**2 + E_totalY**2)
    E_DirX,E_DirY= E_totalX/EE,E_totalY/EE
    end=timer()
    print(timedelta(seconds=end-start))

    print("Calculating E_total for colormap")
    start=timer()
    E_totalXM,E_totalYM = E_total(xGM,yGM,charges)
    EEM = sqrt(E_totalXM**2 + E_totalYM**2)
    end=timer()
    print(timedelta(seconds=end-start))

    # Plots
    ax = subplot(111)
    fig.add_subplot(ax)

    #colormap='Spectral_r'
    colormap='coolwarm'

    print("Making colormap")
    start=timer()
    # Colormap
    I = ax.imshow(EEM,extent=[np.min(xGM),np.max(xGM),np.min(yGM),np.max(yGM)],cmap=colormap, norm=colors.LogNorm(vmin=EEM.min(),vmax=EEM.max()),alpha=0.5,origin='lower',interpolation='bilinear')

    # Colorbar
    fig.colorbar(I,extend='max',orientation='horizontal')
    end=timer()
    print(timedelta(seconds=end-start))

    print("Making quiver plot")
    start=timer()
    # # Vector field
    #Q = ax.quiver(xG,yG,E_DirX,E_DirY,color='gray',alpha=0.7,width=0.002,scale_units='inches',scale=4)
    Q = ax.quiver(xG,yG,E_DirX,E_DirY)
    end=timer()
    print(timedelta(seconds=end-start))

    ax.set_title('N = %s' %N)

    # Plot point charges
    for C in charges:
        if C.q>0:
            plot(C.pos[0], C.pos[1], 'bo', ms=10*sqrt(C.q))
        if C.q<0:
            plot(C.pos[0], C.pos[1], 'ro', ms=10*sqrt(-C.q))

    # Axis sessings
    xlabel('$x$')
    ylabel('$y$')
    ax.set_aspect(aspekt)

    ax.set_xlim(xMin, xMax)
    ax.set_ylim(yMin, yMax)


    print("Savefig")
    savefig("Capacitor_N_%05d_%app.pdf" %(N,app),dpi=fig.dpi,papersize='a4')

    print("Rendering")
    show()



Cap(70)

# N = arange(30,300,30)

# for i in N:
#     Cap(i)
$\endgroup$
  • $\begingroup$ Why the downvote? Please explain me, what's wrong with this question and how to improve it to get help. $\endgroup$ – dp21 Mar 2 at 18:26
  • $\begingroup$ Your question is a bit vague I believe, specially I don't see what's the problem here based on your statement: "However the computation time is very long", did you use some sort of systematic inspection of your code to see which function is the most time consuming one in your calculations? "Very long" is itself vague again, how long is "very long" for you? A couple of minutes? A couple of hours? A couple of days? You should not expect that people here run your code and find what's wrong with your script. You should tell us what is exactly the problem and then people be able to help you here. $\endgroup$ – Alone Programmer Mar 2 at 23:24
  • $\begingroup$ It would be useful if you provide what you are doing mathematically. $\endgroup$ – nicoguaro Mar 3 at 16:26
  • $\begingroup$ Also, is this question related to what you want to achieve? $\endgroup$ – nicoguaro Mar 3 at 18:32
  • 1
    $\begingroup$ In that case, I don't understand what your question is. $\endgroup$ – nicoguaro Mar 4 at 22:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.