Rechneranwendungen in der Physik - Elektrostatisches Potential¶
Santiago.R
In [2]:
import numpy as np
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
import time
%matplotlib inline
def sub2ind(size,j,i):
return i*size[1]+j
Differentialgleichung des elektrostatischen Potentials¶
Aus den Maxwell Gleichungen folgt, dass für das elektrostatische Potential $\Phi$ einer Ladungsverteilung im Raum unter Abwesenheit magnetischer Felder $B=0$ und bei normierter Permitivität $\epsilon = 1$ die Beziehung $\vec{E}(\vec{r})= -\vec{\nabla} \Phi$ gilt. Aus dieser Beziehung folgt, dass das elektrostatische Potential gegeben durch $\Phi$ ein skalares Feld ist und durch die Poisson- bzw. für ladungsfreie Spezialfälle die Laplace-Gleichung als $\vec{\nabla}^2 \Phi = - \frac{q}{\epsilon _0}$ und $\vec{\nabla}^2 \Phi = 0$ jeweils beschrieben werden kann. Lösen der Poisson-Gleichung nach $\Phi$ liefert also dann das elektrostatische Potential einer beliebigen Ladungsverteilung.
Das mathematische Problem besteht also darin, mithilfe numerischer Integrationsmethoden eine Lösung der Poisson-Gleichung für bestimmte Randbedingungen zu finden. Diese können unterschiedlich formuliert werden und werden je nach Art entweder Dirichlet-Randbedingungen oder Neumann-Randbedingungen gennant. Für letztere gilt hierbei;
$
\begin{cases}
\Delta \Phi&=&0&\text{in}\ \Omega\\
\frac {\partial\Phi} {\partial n}&=&h&\text{auf}\ \partial\Omega\\
\end{cases}
$
d.h. es wird vom ladungsfreien Spezialfall der Laplace-Gleichung innerhalb des integrierten Bereiches $\Omega$ ausgegangen und die Normalenableitung des Potentials am Rand $\partial \Omega $ ist dann konstant.
Ladungsverteilung Generator [Christoph T. Koch | Institut für Physik | HU Berlin]¶
In [3]:
def Generiere_Rho(Nx,Ny,xRange=4e-9,yRange=3e-9):
e = 1.6022e-19 # Elementarladung
epsilon0 = 8.854e-12 # Vakuum Permittivität (Elektr. Feldkonstante) in As/Vm.
Dimension = 2
# Wir definieren zunächst leere Arrays for rho (Ladung) und Phi
# (Potential), und dann auch die x- und y-Koordinaten für die Gridpunkte:
rho = np.zeros((Ny,Nx))
Phi = np.zeros((Ny,Nx))
x1D = np.linspace(-xRange,xRange,Nx)
y1D = np.linspace(-yRange,yRange,Ny)
rmax = max(2*xRange,2*yRange)
x,y = np.meshgrid(x1D,y1D)
dx = x1D[2]-x1D[1]
dy = y1D[2]-y1D[1]
h2 = dx*dy;
h = np.sqrt(h2);
h3 = h**3;
# Positionieren von 5 Ladungen:
# An 35 Zufälligen Positionen positionieren wir nun ladungen von je einer
# Elementarladung, aber mit zufälligem Vorzeichen. Dazu generieren wir uns zunächst
# einen (5,1)-Vektor von zufälligen x- und y-Positionen:
Nq = 5 # Anzahl der Ladungen
Nrand = 5 # Mindestabstand vom Rand in Pixeln
random.seed(111) # macht die Zufallszahlen reproduzierbar
# Generiere für jeden ZufallsPunkt eine zufälliges Ladungsvorzeichen:
xq = np.zeros((Nq,1),dtype=int)
yq = np.zeros((Nq,1),dtype=int)
for iq in range(Nq):
xq[iq] = random.randint(Nrand,Nx-Nrand)
yq[iq] = random.randint(Nrand,Nx-Nrand)
rho[yq[iq],xq[iq]] = -1+2*random.randint(0,1) # rho = +/-1
# Generieren des Potentials
# Wir wissen, dass jede der Punktladungen das Potential
# $$ V(r) = \frac{e}{4*\pi*\epsilon_0*|r|}$$
# generiert. Das Gesamtpotential ist dann einfach die Superposition dieser
# Punktladungen:
for iq in range(Nq):
# Zunächst berechnen wir den Abstand jedes Pixels von der Ladung iq:
ix = xq[iq]
iy = yq[iq]
r = np.sqrt((x-x[iy,ix])**2+(y-y[iy,ix])**2)
if Dimension == 3:
rinv = (r**-1)/(4*np.pi)
if Dimension == 2:
rinv = np.log(rmax/r)/(2*np.pi*h)
if Dimension == 1:
rinv = 1/h2
# Jetzt müssen wir noch den Wert von 1/(r=0) so definieren, dass die Zweite Ableitung
# dieser Funktion im Pixel links daneben gleich Null ist:
rinv[iy,ix] = 4*rinv[iy,ix-1]-(rinv[iy+1,ix-1]+rinv[iy-1,ix-1]+rinv[iy,ix-2])
Phi = Phi+rho[iy,ix]*e/epsilon0*rinv;
return rho,Phi,x,y
def surf(x,y,z):
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z,rstride=1,cstride = 1,cmap=plt.cm.CMRmap)
In [4]:
#Test
Nx = 20
Ny = 20
rho,Phi0,x,y = Generiere_Rho(Nx,Ny)
surf(x,y,rho)
surf(x,y,Phi0)
Neumann-Problemm¶
In [5]:
#Tatsächlich habe ich es auch nach langen Versuchen leider nicht geschafft, eine Lösungs-Funktion gebacken zu kriegen... insbesondere bin ich mir hier unsicher bei der
#Implementierung der Neumann-Randbedingungen in den Algorithmus und in dem Verfahren, denn dieser konkret mit dem Laplace-Kernel(3) folgen soll.
def Neumann_Solver(rho,max_iter=0,tol=1e-5):
Ny,Nx = rho.shape
if max_iter == 0:
max_iter = 20*Nx*Ny
Phi = np.zeros((Ny,Nx))
x1D = np.linspace(-xRange,xRange,Nx)
y1D = np.linspace(-yRange,yRange,Ny)
dx = x1D[2]-x1D[1]
dy = y1D[2]-y1D[1]
h2 = dx*dy;
h = np.sqrt(h2);
Laplace_Kernel=1/h**2*np.array([[1/6,2/3,1/6],[2/3,-10/3,2/3],[1/6,2/3,1/6]])
N = Nx*Ny
D = lil_matrix((N,N))
count= 0;
for ix in range(1,Nx-1):
for iy in range (1,Ny-1):
D[count,sub2ind(size,iy,ix)] = 4
D[count,sub2ind(size,iy,ix+1)] = -1
D[count,sub2ind(size,iy,ix-1)] = -1
D[count,sub2ind(size,iy+1,ix)] = -1
D[count,sub2ind(size,iy-1,ix)] = -1
count += 1
D = D.tocsr()
spsolve(D, b)
Jacobi Iterator-Funktion [U.Wolff | Institut für Physik | HU Berlin]¶
In [6]:
def JacobiIter(Phi0,max_iter=0,tol=1e-5):
Ny,Nx = Phi0.shape
if max_iter == 0:
max_iter = 20*Nx*Ny
w = 1.0
# Hier definieren wir das Flag Feld, welches im Inneren 1 ist, und am Rand 0:
flag = np.zeros((Ny,Nx)) # einhuellende Matrix mit Nullen
flag[2:Ny-1,2:Nx-1] = 1 # innen Einsen
phi = np.zeros((Ny,Nx)) # = phineu
change = np.zeros((max_iter,))
# Übertrage (skaliert) dessen Rand auf unser Problem.
phi[flag == 0] = Phi0[flag == 0];
phi0 = np.max(np.abs(Phi0));
## Jacobi-Iterationen
count = 0
for iter in range(max_iter):
count += 1
phialt = phi # fuer Jacobi-Verfahren
phi = 0.25*(np.roll(phialt,1,axis=0)+np.roll(phialt,-1,axis=0)+np.roll(phialt,1,axis=1)+np.roll(phialt,-1,axis=1))
phi[flag == 0] = Phi0[flag == 0]
# Berechne Konvergenz:
temp = np.sum(abs(phi-phialt)**2)
# Überrelaxation:
phi = phialt+w*(phi-phialt)
w = 1+0.99*(w-1)
phialt = phi;
change[iter] = np.sqrt(temp/((Nx-2)*(Ny-2)-1)); # Standardabweichung
if (iter > 2) & (w != 1):
if change[iter] >= change[iter-1]:
w = 1+0.5*(w-1);
# print('Nach',iter,'Iterationen, Aenderung =',change[iter]);
if( change[iter] < tol ):
change = change[0:count]
return phi,change
return phi,change
In [7]:
#Test
t_start = time.perf_counter()
phiRec,change = JacobiIter(Phi0,tol=1e-5)
elapsed_time = time.perf_counter()-t_start
print('Anzahl der Iterationen =',len(change));
print('Iterationszeit =',round(elapsed_time,4),'s');
fig = plt.figure()
plt.semilogy(change);
ax = fig.gca()
# ax.grid(color='r', linestyle='-', linewidth=2)
plt.xlabel('Iteration'); plt.ylabel('Aenderung');
surf(x,y,phiRec); plt.title('Rekonstruiertes Potential');
Test für verschiedene Pixelzahlen (Jacobi-Iterator-Funktion)¶
In [8]:
#Speicherlisten
N_liste=[]
t_liste=[]
#Iterations-Loop
for i in range(20,300,10):
Nx=i
Ny=i
rho,Phi0,x,y = Generiere_Rho(Nx,Ny)
t_start = time.perf_counter()
phiRec,change = JacobiIter(Phi0,tol=1e-5)
elapsed_time = time.perf_counter()-t_start
N_liste.append(np.log(i))
t_liste.append(np.log(elapsed_time))
print('Anzahl der Iterationen =',len(change), 'für N=', i);
print('Iterationszeit =',round(elapsed_time,4),'s');
fig = plt.plot(N_liste, t_liste)
plt.xlabel('Logarithmus der Pixelzahl NxN'); plt.ylabel('Logarithmus der Iterationszeit');
In [ ]: