# Partial Differential Equations

A partial differential equation (PDE) of the function $u(x_1,\dots,x_n)$ is an equation of the form: $$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$$

We will focus on linear $2^{nd}$ order PDEs in 2D, which can be written in the form: $$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$$ 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$ :

• $b^2 - 4 a c < 0$ : elliptic, e.g. the Lapace equation in 2D $$\frac{\partial^2 u}{\partial x^2 } + \frac{\partial^2 u}{\partial y^2} = 0$$

• $b^2 - 4 a c = 0$ : parabolic, e.g. the Diffusion equation in 1D $$\frac{\partial^2 u}{\partial x^2} - \frac{\partial u}{\partial t } = 0$$

• $b^2 - 4 a c > 0$ : hyperbolic, e.g. the Wave equation in 1D $$\frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial t^2} = 0$$

• # Poisson equation in 2D

The Poisson equation is a $2^{nd}$ order PDE of the elliptic type: $$\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f(x,y)$$

In a 2D electrostatic problem the Poisson equation can be written as $$\nabla^2 V(x,y) = - \rho( x, y )$$ 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.

### Assignment

We want to solve the 2D Poisson equation for this source term: $$\rho(x,y) = sin\left( 6 \pi \frac{x}{L} \right) sin\left( 2 \pi \frac{y}{L} \right)$$ whith $L=10$.

### Summary

We will use three methods:

• analytical solution

• spectral method

• finite difference method

• # Analytical solution

For the aforementioned source term, it is easy to verify that the solution of the 2D Poisson equation is: $$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)$$ Let's plot $\rho(x,y)$ and $V(x,y)$.

In [ ]:
# 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)

In [ ]:
# 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])

In [ ]:
# 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)

In [ ]:
# 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)

In [ ]:
# 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()

In [ ]:
plot3D( X, Y, rho)

In [ ]:
plot3D( X, Y, Van)


### Check for understanding

Q: How many periods are along the x and y direction?

# Spectral method

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: $$\rho(\mathbf{r}) = \sum_\mathbf{G} \rho(\mathbf{G}) e^{i\mathbf{G}\cdot\mathbf{r}}$$ where the Fourier components $\rho(\mathbf{G})$ are obtained as: $$\rho(\mathbf{G}) = \frac{1}{L^2}\int d\mathbf{r} \rho(\mathbf{r}) e^{-i\mathbf{G}\cdot\mathbf{r}}$$ $\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 $$G^2 V(\mathbf{G}) = \rho(\mathbf{G})$$ 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)$ :

• Fourier transform $\rho(\mathbf{r}) \rightarrow \rho(\mathbf{G})$

• Solve the Poisson equation in Fourier space $V(\mathbf{G})=\frac{\rho(\mathbf{G})}{G^2}$

• Fourier transform $V(\mathbf{r}) \leftarrow V(\mathbf{G})$

• In [ ]:
# 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)

In [ ]:
# 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)


### Check for understanding

Q: Why do we need to work with complex numbers?

In [ ]:
# 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))


### Check for understanding

Q: Can you explain why we have both positive and negative values for $n$ and $m$?

In [ ]:
# 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)


### Check for understanding

Q: Why is $V(\mathbf{G}=0)$ arbitrary?

In [ ]:
# Fourier transform V to real space
RV = np.fft.ifftn(GV)

In [ ]:
plot3D(X,Y,RV.real)

In [ ]:
# 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)

In [ ]:
maxAD = maxAbsDiff(RV.real,Van)


# Finite difference

The second derivative of a function $f$ can be computed numerically using the central difference formula $$\frac{d^2f}{dx^2} \approx \frac{f(x+h)+f(x-h)-2f(x)}{h^2}$$ where $h$ is a small increment.

In a similar way, we can discretize the Lapace operator on the 2D grid: $$\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}$$ This yields the discretized form of the Poisson 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}$$

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. $$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}$$

We solve the Possion equation by iterative finite differences

In [ ]:
# Define Delta
Delta = x[1]-x[0]
print("Delta = ", Delta )

In [ ]:
# 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()

In [ ]:
# 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

In [ ]:
# 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)

In [ ]:
#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: $$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]$$ where $\omega$ is a parameter.

In [ ]:
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

In [ ]:
#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)

In [ ]:
# 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)


### Check for understanding

Q: Why does the new method improve convergence?