A partial differential equation (PDE) of the function $u(x_1,\dots,x_n)$ is an equation of the form: \begin{equation} F\big( \underbrace{x_1,\dots,x_n}_{\text{variables}} , \underbrace{u}_{\text{function}}, \underbrace{ \frac{\partial u}{\partial x_1}, \dots, \frac{\partial u}{\partial x_n}}_{\text{first derivatives}}, \underbrace{ \frac{\partial^2 u}{\partial x_1^2}, \frac{\partial^2 u}{\partial x_1\partial x_2}, \dots }_{\text{second derivatives}}, \underbrace{\dots}_{\text{higher order derivatives}} \big) = 0 \end{equation}
We will focus on linear $2^{nd}$ order PDEs in 2D, which can be written in the form: \begin{equation} a \frac{\partial^2 u}{\partial x^2} + b \frac{\partial^2 u}{\partial x \partial y} + c \frac{\partial^2 u}{\partial y^2} + \underbrace{\dots}_{\text{lower terms}} = 0 \end{equation} Like for conic equations, we can define three classes of linear $2^{nd}$ order PDEs, based on the value of discriminant $b^2 - 4 a c $ :
The Poisson equation is a $2^{nd}$ order PDE of the elliptic type: \begin{equation} \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y) \end{equation}
In a 2D electrostatic problem the Poisson equation can be written as \begin{equation} \nabla^2 V(x,y) = - \rho( x, y ) \end{equation} where $\nabla^2=\sum_i\frac{\partial^2}{\partial x_i^2}$ is the Laplace operator, $V$ is the potential, and $\rho$ is the charge density. We note that the Poisson equation differs from the Lapace equation because of the presence of the source term $\rho$, which makes the equation inhomogeneous.
For the aforementioned source term, it is easy to verify that the solution of the 2D Poisson equation is: \begin{equation} V(x,y)=\frac{L^2}{40\pi^2} sin\left( 6 \pi \frac{x}{L} \right) sin\left( 2 \pi \frac{y}{L} \right) \end{equation} Let's plot $\rho(x,y)$ and $V(x,y)$.
# Define the periodicity
L = 10.0
# Define the number of points used to discretize the source and the potential on a NxN mesh
N = 25
# Create the list that contains the x positions
x = []
for i in range(N) :
x.append(float(i)/float(N)*L)
# Create the list that contains the y positions
y = []
for j in range(N) :
y.append(float(j)/float(N)*L)
# Show the arrays
print("x = ", x)
print("y = ", y)
# Access elements of the array x (Python experts can skip)
print("x[0] = ", x[0])
print("x[1] = ", x[1])
print("x[N-1] = ", x[N-1])
print("x[-1] = ", x[-1])
print("x[-N] = ", x[-N])
# create X and Y grids to easily handle gridpoints
import numpy as np
X, Y = np.meshgrid(x, y)
print("X = ", X)
print("Y = ", Y)
# Define rho and V on the grid
import numpy as np
# rho
rho = np.empty([N,N])
rho = np.sin(X/L*np.pi*6.0)*np.sin(Y/L*np.pi*2.0)
# V analytical solution
Van = np.empty([N,N])
Van = (np.sin(X/L*np.pi*6.0)*np.sin(Y/L*np.pi*2.0))/(40.0*np.pi*np.pi/L/L)
print("rho = ", rho)
print("Van = ", Van)
# method to plot a 3D function
def plot3D( X, Y, Z ) :
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
fig = plt.figure()
ax = fig.gca(projection='3d')
# Plot the surface.
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
# Customize the z axis.
ax.set_zlim(-1.01, 1.01)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Set labels
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
plot3D( X, Y, rho)
plot3D( X, Y, Van)
Both the potential and the source are 2D periodic functions, with periodicity $[L,L]$, and can be therefore expanded in plane waves using the following Fourier expansion: \begin{equation} \rho(\mathbf{r}) = \sum_\mathbf{G} \rho(\mathbf{G}) e^{i\mathbf{G}\cdot\mathbf{r}} \end{equation} where the Fourier components $\rho(\mathbf{G})$ are obtained as: \begin{equation} \rho(\mathbf{G}) = \frac{1}{L^2}\int d\mathbf{r} \rho(\mathbf{r}) e^{-i\mathbf{G}\cdot\mathbf{r}} \end{equation} $\mathbf{r}=(x,y)$, $\mathbf{G}=\frac{2\pi}{L}(n,m)$, with $n=0,\pm 1, \pm 2, \dots $ and $m=0,\pm 1, \pm 2, \dots $.
The Poisson equation in Fourier space becomes a simple multiplication \begin{equation} G^2 V(\mathbf{G}) = \rho(\mathbf{G}) \end{equation} where $G^2 = \mathbf{G}\cdot\mathbf{G}=G_x^2+G_y^2$.
We can therefore follow this protocol to obtain $V(x,y)$ :
# Define G vectors
import numpy as np
gx = np.fft.fftfreq(N,1/N) * 2*np.pi/L
gy = np.fft.fftfreq(N,1/N) * 2*np.pi/L
print("list of n",np.fft.fftfreq(N,1/N))
print("gx = ",gx)
print("gy = ",gy)
Gx, Gy = np.meshgrid(gx, gy)
print("Gx = ", Gx)
print("Gy = ", Gy)
G2 = np.empty([N,N])
G2 = Gx*Gx + Gy*Gy
print("G2 = ", G2)
# Fourier transform rho
import numpy as np
Rrho = rho.astype(complex)
Grho = np.fft.fftn(Rrho)
print(Grho)
print("Grho.real = ",Grho.real)
print("Grho.imag = ",Grho.imag)
# Analize Grho
for j in range(N) :
for i in range(N) :
if (abs(Grho[i,j].real) > 0.00001) :
print('real:',Grho[i,j].real,'n:',Gx[i,j]*L/(2.0*np.pi),'m:',Gy[i,j]*L/(2.0*np.pi))
if (abs(Grho[i,j].imag) > 0.00001) :
print('imag:',Grho[i,j].imag,'n:',Gx[i,j]*L/(2.0*np.pi),'m:',Gy[i,j]*L/(2.0*np.pi))
# Solve Poisson equation in Fourier space
GV = np.empty([N,N],dtype=complex)
GV[0,0] = 0.0 # arbitrary constant
for j in range(N) :
for i in range(N) :
if( i != 0 or j!=0 ) :
GV[i,j] = Grho[i,j] / G2[i,j]
print("GV = ",GV)
# Fourier transform V to real space
RV = np.fft.ifftn(GV)
plot3D(X,Y,RV.real)
# Define the max Abs difference between A[i,j] and B[i,j]
def maxAbsDiff(A,B) :
import numpy as np
C = A-B
Cflat = np.abs(C.flatten())
return max(Cflat)
maxAD = maxAbsDiff(RV.real,Van)
print("maxAD = ",maxAD)
The second derivative of a function $f$ can be computed numerically using the central difference formula \begin{equation} \frac{d^2f}{dx^2} \approx \frac{f(x+h)+f(x-h)-2f(x)}{h^2} \end{equation} where $h$ is a small increment.
In a similar way, we can discretize the Lapace operator on the 2D grid: \begin{equation} \nabla^2 V(x,y) \approx \frac{V(x+\Delta,y)+V(x-\Delta,y)+V(x,y+\Delta)+V(x,y-\Delta)-4V(x,y)}{\Delta^2} \end{equation} This yields the discretized form of the Poisson equation: \begin{equation} \nabla^2V_{i,j} = \frac{V_{i+1,j}+V_{i-1,j}+V_{i,j+1}+V_{i,j-1}-4V_{i,j}}{\Delta^2} =-\rho_{i,j} \end{equation}
We can interpret this as a self-consistent equation: the value of $V$ in every point of the grid depends on the average of its four neighbors plus the source contribution. \begin{equation} V^{n+1}_{i,j} = \frac{V^n_{i+1,j}+V^n_{i-1,j}+V^n_{i,j+1}+V^n_{i,j-1}}{4} + \frac{\Delta^2}{4} \rho_{i,j} \end{equation}
We solve the Possion equation by iterative finite differences
# Define Delta
Delta = x[1]-x[0]
print("Delta = ", Delta )
# Method that plots values in listA and optionally of listB
def plotMaxDiff(listA,listB=[]) :
import matplotlib.pyplot as plt
xB = []
for i in range(len(listB)) :
xB.append(i)
plt.plot(xB,listB,'bo-')
xA = []
for i in range(len(listA)) :
xA.append(i)
plt.plot(xA,listA,'ro-')
plt.xlabel('iteration')
plt.show()
# Define how Vnew is obtained from Vold
def updateV(Vold,source,N,Delta) :
import numpy as np
Vnew = np.empty([N,N])
for j in range(N) :
for i in range(N) :
iplus = i+1
if (iplus==N) :
iplus=0
jplus = j+1
if (jplus==N) :
jplus=0
Vnew[i,j]=0.25*(Vold[i-1,j]+Vold[iplus,j]+Vold[i,j-1]+Vold[i,jplus])+ 0.25*Delta*Delta*source[i,j]
return Vnew
# init
list_max = []
import numpy as np
Vold = np.zeros([N,N])
Vnew = Vold
plot3D(X,Y,Vnew)
list_max.append( maxAbsDiff(Vnew,Van) )
plotMaxDiff(list_max)
#iterate
Vold = Vnew
Vnew = updateV(Vold,rho,N,Delta)
plot3D(X,Y,Vnew)
list_max.append( maxAbsDiff(Vnew,Van) )
plotMaxDiff(list_max)
We introduce a small variation to the previous method. Here self-consistency is obtained by mixing at every iteration the old and new values: \begin{equation} V^{n+1}_{i,j} = (1-\omega)V^{n}_{i,j} +\omega \left[\frac{V^n_{i+1,j}+V^n_{i-1,j}+V^n_{i,j+1}+V^n_{i,j-1}}{4} + \frac{\Delta^2}{4} \rho_{i,j}\right] \end{equation} where $\omega$ is a parameter.
def newUpdateV(Vold,source,N,Delta,omega) :
import numpy as np
Vnew = np.empty([N,N])
for j in range(N) :
for i in range(N) :
iplus = i+1
if (iplus==N) :
iplus=0
jplus = j+1
if (jplus==N) :
jplus=0
Vnew[i,j]=(1.0-omega)*Vold[i,j]+omega*0.25*(Vold[i-1,j]+Vold[iplus,j]+Vold[i,j-1]+Vold[i,jplus])+omega*0.25*Delta*Delta*source[i,j]
return Vnew
#init
list_max_new = []
import numpy as np
Vold = np.zeros([N,N])
Vnew = Vold
plot3D(X,Y,Vnew)
list_max_new.append( maxAbsDiff(Vnew,Van) )
plotMaxDiff(list_max_new,list_max)
# iterate
Vold = Vnew
Vnew = newUpdateV(Vold,rho,N,Delta,2)
plot3D(X,Y,Vnew)
list_max_new.append( maxAbsDiff(Vnew,Van) )
plotMaxDiff(list_max_new,list_max)