Exercise_12 : Electric Potentials and Fields: Laplace's Equations

Problem 5.7

冉峰 弘毅班 2014301020064

Abstract

• The partial differential equation is of physicists' great interest for its widely applicaiton in decribing physical systems. This article made a comparsion between Jacobi method and SOR algorithm to solve partial differential equation. Also, it compared the number of iterations and figured out the relation between the number of iterations and the number of grid elements.
  

Background

Introduction

• As we known, the electric potential satisfies Laplace's equation in the space without any electric charges
  $$\begin{equation*} \nabla^2V=0 \end{equation*}$$
  Numerically, we get that:
  $$V(i,j,k) = \frac{1}{6}[V(i+1,j,k)+V(i-1,j,k)+V(i,j+1,k)\\ +V(i,j-1,k)+V(i,j,k+1)+V(i,j,k-1)]$$
  Specially, for the 2-dimensional case:
  $$V(i,j,k) = \frac{1}{4}[V(i+1,j,k)+V(i-1,j,k)+V(i,j+1,k)\\ +V(i,j-1,k)]$$

Jacobi Method

• In numerical linear algebra, the Jacobi method (or Jacobi iterative method) is an algorithm for determining the solutions of a diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges.
  Particularly,in our dicussion in this passage,Jacobi Method can be easily expressed as below:
  $$V_{new}(i,j)=\frac{1}{4}[V_{old}(i+1,j)+V_{old}(i-1,j)+V_{old}(i,j+1)+V_{old}(i,j-1)]$$

Gauss-Sideal Method

• In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and is similar to the Jacobi method.
  Particularly,in our dicussion in this passage,Gauss-Seidel Method can be easily expressed as below:
  $$\begin{eqnarray*} V_{new}(i,j)=\frac{1}{4}[V_{old}(i+1,j)+V_{new}(i-1,j)+V_{old}(i,j+1)+V_{new}(i,j-1)] \end{eqnarray*}$$
  Worth mentioning,whilw the performance if the Gauss-Seidel akgorithm is better than that of the Jacobi method,the imporvement is only modest.It turns out that the number of iterations required for convergence with the Gauss-Siedel method is smaller,but only by a factor of 2.
  But the benefits in Gauss-Seidel method is obvious,cause there is no need for us to save the old values seperately.In other words, the updating can be performed in place with $V(i,j)$ using an algorithm of the form.
  $$\begin{eqnarray*} V(i,j)=\frac{1}{4}[V(i+1,j)+V(i-1,j)+V(i,j+1)+V(i,j-1)] \end{eqnarray*}$$

Simulatneous over-relaxation(SOR)

• With this method ,the problem of relatively low convergence rate can be easily solved.Let $V^*(i,j)$ be the new value of calculated.Then the change of $V$ is represented as:
  $$\begin{eqnarray*} \Delta V(i,j)=V^*(i,j)-V_{old}(i,j) \end{eqnarray*}$$
  However,we have seen in our previous examples that this choice of $\Delta V(i,j)$ is too conservative.So to speed up convergence we will change the potential by a larger amount calculated according to:
  $$\begin{eqnarray*} V_{new}(i,j)=\alpha \Delta V^*(i,j)+V_{old}(i,j) \end{eqnarray*}$$
  where $\alpha$ is a factor that measures how much we "over-relax".Choosing $\alpha=1$ tuekds the Gauss-Seidel method.It turns out that if $\alpha\geq 2$ , the method does bit converge.The best choice of $\alpha$ is:
  $$\begin{eqnarray*} \alpha \approx \frac{2}{1+\pi/L} \end{eqnarray*}$$

The Main Body

• In problem 5.7, we need to write two problems to slove the capacitor problem of figures 5.6 and 5.7, one using the Jacobi method and one using the SOR algorithm. For a fixed accuracy (as set by the convergence test) compare the number of iterations, $N_{iter}$, that each algorithm requires as a function of the number of the grid elelments, $L$. Show that for the Jacobi method $N\sim L^2$, while with SOR $N\sim L$.
  

Code for Jacobi Method

@author: RF
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
class jacobi:
    def __init__(self):
        pass
    def initialization(self,initial_file):
        itxt= open(initial_file)
        self.lattice_in = []
        self.lattice_out = []
        for lines in itxt.readlines():
            lines=lines.replace("\n","").split(",")
            #print (lines[0].split(" ")
            self.lattice_in.append(lines[0].split(" "))
            self.lattice_out.append(lines[0].split(" "))
        itxt.close()
        #print (self.lattice_in)
        #print (self.initial_lattice)
        for i in range(len(self.lattice_in)):
            for j in range(len(self.lattice_in[i])):
                self.lattice_in[i][j] = float(self.lattice_in[i][j])
                self.lattice_out[i][j] = float(self.lattice_out[i][j])
        return 0
    def evolution(self):
        delta = 10
        self.N = 0
        delta_record = []
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
            delta_record.append(delta)
            for i in range(len(self.lattice_in)):
                for j in range(len(self.lattice_in[i])):
                    self.lattice_in[i][j] = self.lattice_out[i][j]
            self.N+= 1
        plt.plot(delta_record,label = 'jacobi method')
        print (self.N)
        print (len(self.lattice_in))
        return 0
    def electric_field(self):
        self.ex = np.zeros([len(self.lattice_in),len(self.lattice_in)])
        self.ey = np.zeros([len(self.lattice_in),len(self.lattice_in)])
        deltax = 0.3333
        deltay = 0.3333
        for i in range(1,len(self.lattice_in) - 1):
            for j in range(1,len(self.lattice_in[0]) - 1):
                self.ey[i][j] = (-(self.lattice_in[i + 1][j] - self.lattice_in[i - 1][j])/(2*deltax))
                self.ex[i][j] = (-(self.lattice_in[i][j + 1] - self.lattice_in[i][j - 1])/(2*deltay))
        return 0
    def plot_field(self):
        self.electric_field()
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        #plt.figure(figsize = (8,8))
        plt.quiver(X, Y, self.ex, self.ey, pivot = 'mid', units='width')
        plt.title('Electric field between two metal plates')
        #plt.show()
        return 0
    def plot_contour(self): 
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        #plt.figure(figsize = (8,8))
        CS = plt.contour(X, Y, self.lattice_in, 13)
        plt.clabel(CS, inline=1, fontsize=10)
        plt.title('Equipotential lines')
        #plt.show()
        return 0
    def plot_surface(self):
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        X, Y = np.meshgrid(np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)), np.arange(-1.00, 1.01, 2./(len(self.lattice_in) - 1)))
        surf = ax.plot_surface(X, Y, self.lattice_in, rstride=1, cstride=1,cmap = cm.coolwarm_r,
                       linewidth=0, antialiased=False)
        ax.set_zlim(-1.01, 1.01)
        ax.zaxis.set_major_locator(LinearLocator(10))
        ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
        fig.colorbar(surf, shrink=0.5, aspect=5)
#A = jacobi()
A = SOR()
A.initialization('initial_capacitor_26.txt')
A.evolution()
#plt.figure(figsize = (8,8))
plt.figure(figsize = (16,8))
plt.subplot(121)
A.plot_contour()
plt.subplot(122)
A.plot_field()
#A.plot_surface()
plt.savefig('chapter5_5.7_26.png', dpi = 256)
plt.show()

Code for SOR Method

@author: RF
class SOR(jacobi):
    def evolution(self):
        delta = 10
        self.N = 0
        delta_record = []
        alpha = 2/(1+np.pi/len(self.lattice_in))
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
                        self.lattice_in[i][j] = alpha * (self.lattice_out[i][j] - self.lattice_in[i][j]) + self.lattice_in[i][j]
            #print (self.lattice_out)
            delta_record.append(delta)
            self.N+= 1
        plt.plot(delta_record, label ='SOR algorithm')
        print (self.N)
        print (len(self.lattice_in))
        #print (self.lattice_in)
        return 0    
class SOR_2(jacobi):
    def evolution(self,alpha):
        delta = 10
        self.N = 0
        delta_record = []
        #alpha = 2/(1+np.pi/len(self.lattice_in))
        while (delta > 0.00001*len(self.lattice_in)*len(self.lattice_in[0])):
            delta = 0
            for i in range(1,len(self.lattice_in) - 1):
                for j in range(1,len(self.lattice_in[i]) - 1):
                    if (self.lattice_in[i][j] != 1 and self.lattice_in[i][j] != -1):
                        self.lattice_out[i][j] = 0.25*(self.lattice_in[i - 1][j] + self.lattice_in[i + 1][j] + self.lattice_in[i][j - 1] + self.lattice_in[i][j + 1])
                        delta+= abs(self.lattice_out[i][j] - self.lattice_in[i][j])
                        self.lattice_in[i][j] = alpha * (self.lattice_out[i][j] - self.lattice_in[i][j]) + self.lattice_in[i][j]
            #print (self.lattice_out)
            delta_record.append(delta)
            self.N+= 1
        #plt.plot(delta_record, label ='SOR algorithm')
        print (self.N)
        #print (len(self.lattice_in))
        #print (self.lattice_in)
        return 0  

Result

Setting initial value and boundary conditions

• The initial value and boundary conditions which have been set in a .txt file
  

  The potential of boundaries is zero. In the inner space, there are two capacitor plates held at $V=1$ (left plate) and $V=-1$ (right plate).

Solving the partial differential equation of two capacitor plates

• Use both of these method, we can easily solve the problem of the potential near two capacitor plates.
  Firatly, solving the problem in the 11*11 grids.
  
  Then, enlarging the grids to 26*26 leads to smoother equipotential curve and a finer structure of electric field.
  

The relationship between the number of iterations and that of grid elements

• In this section, I will make a discussion on the relationship between the number of iterarions and that of grid elements. For convenience, all the grids I used to solve equation is square grids.
  The fitting results are $0.2772L^2$ for the blue curve and $1.701L$ for the red one, which means the number of iteration is proportional to $L^2$ for $Jacobi method$ and is proportional to $L$ for SOR algorithm.
  Therefore, the larger grid we used, the bigger advantage of SOR algorithm would be showed.
  

The Convergence speed

• On the other hand, we can make a comparison of the convergence speed of these methods.
  We can see the convergence speed of SOR algorithm is much larger than that of jacobi method.
  

The $\alpha$ of SOR algorithm

• Additionally, I make a discussion of the relation between the number of iteration and the value of $\alpha$.
  In this case, I set the length of grid is 26, which means the value of $\alpha$ is about 1.783.
  Obviously, this $\alpha$ is not the best choice to decrease the number of iterations.
  

Conclusion

• In the different situations, we can get the same results as the Electromagnetics, so that the method of these two canbe used to solve the problems. Meanwhile, we can easily get the conclusion that the method SOR can be less number of iterations to the Jacobi method.

Thanks For

• The book of Computational Physics, chapter 5
• Chenfeng
• Wikipedia
• Matplotlib Gallery