Partial Differential Equations

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 $ :

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

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

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

  • Poisson equation in 2D

    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.

    Assignment

    We want to solve the 2D Poisson equation for this source term: \begin{equation} \rho(x,y) = sin\left( 6 \pi \frac{x}{L} \right) sin\left( 2 \pi \frac{y}{L} \right) \end{equation} 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: \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)$.

    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: \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)$ :

  • 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)
    print("maxAD = ",maxAD)
    

    Finite difference

    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

    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: \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.

    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?