** Solve elliptic partial differential equations numerically in Python. ** Since Laplace's equation is widely applied not only in electromagnetism but also in physics, it will be of great learning value. In addition to physics, I think we can learn from the point of numerical solution of boundary value problems of PDEs. However, the method used in this calculation is slow to converge and is not practical. There is a faster numerical solution [2].
** In this paper, the two-dimensional Laplace equation and Poisson equation are differentiated and the potential is determined by the iterative method (Jacobi method. % B3% E3% 83% 93% E6% B3% 95)). The electric field is also obtained from the obtained potential. The derivation of the equation and the outline of the algorithm are summarized in [Addendum]. If you are interested, please refer to it. ** **
● [Addition] August 14, 2017 ** Posted the code to solve Poisson's equation by SOR method (lower part of reference). It is more than 50 times faster than the Jacobi method. </ font> **
Example (1) Solve Laplace's equation. As shown in the figure, the boundary condition is \ phi = V volt at y = 0. It is grounded at the other boundary and is 0 volt. The electrostatic position at this time is determined and visualized.
Example (2) Solve Poisson's equation. Put the external charge Q in the area considered in Example (1). At this time, the electrostatic potential $ \ phi $ and the electric field $ E =-\ nabla \ phi $ are obtained and visualized.
Animation is also generated in the code. This is to investigate how the solution \ phi (x, y) converges.
Postscript December 17, 2017: It is supposed to work on jupyter. If you don't use jupyter, comment out "% matplotlib nbagg" in your code. From masa_stone22's point.
laplace.py3
"""
Laplace equation:Numerical solution,Jacobi(Jacobi)Law:
12 Aug. 2017
"""
%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation #Import methods for creating animations
fig = plt.figure()
anim = []
#
delta_L=0.1 #Grid width
LL = 10 #Square width
L = int(LL/delta_L)
V = 5.0 #Boundary potential on one side
convegence_criterion = 10**-3 #Convergence condition. Decrease this value to improve accuracy.
phi = np.zeros([L+1,L+1])
phi[0,:] = V #boundary condition
phi2 = np.empty ([L+1,L+1])
#for plot
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#Main
delta = 1.0
n_iter=0
conv_check=[]
while delta > convegence_criterion:
if n_iter % 100 ==0: #Monitoring of convergence status
print("iteration No =", n_iter, "delta=",delta)
conv_check.append([n_iter, delta])
for i in range(L+1):
for j in range(L+1):
if i ==0 or i ==L or j==0 or j==L:
phi2[i,j] = phi[i,j]
else:
phi2[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4 #Addendum:formula(11)See
delta = np.max(abs(phi-phi2))
phi, phi2 = phi2, phi
n_iter+=1
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#for plot
plt.colorbar () #Color bar display
plt.xlabel('X')
plt.ylabel('Y')
ani = ArtistAnimation(fig, anim, interval=100, blit=True,repeat_delay=1000)
ani.save("t.gif", writer='imagemagick')
plt.show()
The boundary condition $ \ phi = 5 $ for y = 0 is satisfied. The change in potential is large near the boundary.
The external charge density is $ c = 1000 $, and the boundary potential is $ V = 5 $.
Poisson.py3
"""
Poisson equation:Jacobi(Jacobi)Law
Also displays the electric field:
12 Aug. 2017
"""
%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation #Import methods for creating animations
anim = []
#
delta_L=0.01 #Grid width
LL = 1 #Square width
L = int(LL/delta_L)
V = 5.0 #Boundary y=0 potential(boundary condition)
convegence_criterion = 10**-5
phi = np.zeros([L+1,L+1])
phi[0,:] = V #boundary condition
phi2 = np.empty ([L+1,L+1])
#Charge density setting
eps0=1
charge= np.zeros([L+1,L+1])
c=1000 #Charge density
CC=int(L/4)
CC2=int(L/2)
#Charge density x=[L/4, L/2], y=[L/4,L/2]Embed in the area of
for i in range(CC, CC2+1):
for j in range(CC,CC2+1):
charge[i,j] = c
#
#for plot
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#Main
delta = 1.0
n_iter=0
conv_check=[]
while delta > convegence_criterion:
if n_iter % 100 ==0: #Monitoring of convergence status
print("iteration No =", n_iter, "delta=",delta)
conv_check.append([n_iter, delta])
for i in range(L+1):
for j in range(L+1):
if i ==0 or i ==L or j==0 or j==L:
phi2[i,j] = phi[i,j]
else:
phi2[i,j] = (phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4 + (delta_L**2/(4*eps0))*charge[i,j]
delta = np.max(abs(phi-phi2))
phi, phi2 = phi2, phi
n_iter+=1
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#for plot
pp=plt.colorbar (orientation="vertical") #Color bar display
pp.set_label("phi", fontname="Ariel", fontsize=24)
plt.xlabel('X',fontname="Ariel", fontsize=24)
plt.ylabel('Y',fontname="Ariel", fontsize=24)
plt.title('Solution of the Poisson equation for electrostatic potential ')
#ani = ArtistAnimation(fig, anim, interval=10, blit=True,repeat_delay=1000)
#ani.save("t.gif", writer='imagemagick')
plt.show()
#Electric field calculation
Ex = np.zeros([L,L])
Ey = np.zeros([L,L])
for i in range(L):
for j in range(L):
Ex[i,j] = -(phi[i+1,j]-phi[i,j])/delta_L
Ey[i,j] = -(phi[i,j+1]-phi[i,j])/delta_L
plt.figure()
X, Y= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1)) #Mesh generation
plt.quiver(X,Y,Ex,Ey,color='red',angles='xy',scale_units='xy', scale=40) #Plot vector field
plt.xlim([0,L]) #Range of X to draw
plt.ylim([0,L]) #Range of y to draw
#Graph drawing
plt.grid()
plt.draw()
plt.show()
Result (2): Solution of Poisson's equation of the system with embedded external charge
The boundary condition at y = 0 is $ \ phi = 5 $. Due to the embedded external charge, the potential differs greatly compared to Example (1).
The state of the electric field. (The #y axis has been reversed). This figure was drawn as delta_L = 0.1
.
As expected from the fluctuation of the potential, the electric field fluctuates greatly near the external charge.
The relational expression between the electrostatic field $ E $ and the potential $ \ phi $
Is the starting point.
Substituting (1) into (2)
To get. This is the Poisson equation. Where $ \ epsilon_0 $ is the permittivity in vacuum.
If there is no true electric charge everywhere in space ($ \ rho (x, y, z) = 0 $), the Poisson equation is
Think two-dimensional. Poisson's equation
$ \frac{\partial^2 \phi(x,y)}{\partial x^2} + \frac{\partial^2 \phi(x,y)}{\partial y^2} =-\frac{\rho (x,y,z)} {\epsilon_0} \tag{5}$
Will be. Consider a two-dimensional grid (see figure below) separated by a width of $ \ Delta $ in the $ x $ and $ y $ directions. Taylor expansion of the potential $ \ phi (x, y) $ at the point ($ x + \ Delta, y $)
When Eqs. (6) and Eq. (7) are added together,
$\frac{\partial^2 \phi(x,y)}{\partial x^2} = \frac{\phi(x+\Delta,y)+ \phi(x-\Delta,y)-2\phi(x,y)}{\Delta^2} +\mathcal{O}(\Delta^2) \tag{8} $
Taylor expansion is done for y in the same way, $\frac{\partial^2 \phi(x,y)}{\partial y^2} = \frac{\phi(x,y+\Delta)+ \phi(x,y-\Delta)-2\phi(x,y)}{\Delta^2} +\mathcal{O}(\Delta^2) \tag{9} $
To get. Substitute equations (8) and (9) into (5) and sort them out a little.
$ \phi(x,y) = \frac{1}{4} [\phi(x+\Delta,y)+\phi(x-\Delta,y)+\phi(x,y+\Delta)+\phi(x,y-\Delta)]+\frac{\rho(x,y)}{4 \epsilon_0}\Delta^2+\mathcal{O}(\Delta^4) \tag{10} $
To get. This is a ** differential representation of Poisson's equation **. The error is $ \ mathcal {O} (\ Delta ^ 4) $.
When $ \ rho (x, y) = 0 $
$\phi(x,y) = \frac{1}{4} [\phi(x+\Delta,y)+\phi(x-\Delta,y)+\phi(x,y+\Delta)+\phi(x,y-\Delta)]+\mathcal{O}(\Delta^4) \tag{11} $
Will be. This is a ** differential representation of Laplace's equation **.
In numerical calculation, if the subscript of the array is $ (x_i, y_j) $, equation (10) is
$ \phi(x_i,y_j) = \frac{1}{4} [\phi(x_{i+1},y_{j})+\phi(x_{i-1},y_{j})+\phi(x_{i},y_{j+1})+\phi(x_{i},y_{j-1})]+\frac{\rho(x_{i},y_{j})}{4\epsilon_0}\Delta^2 \tag{12} $ Will be.
The potential $ \ phi (i, j) $ at the point $ (x_i, y_j) $ is approximated by adding the amount proportional to the charge density $ \ rho (x, y) $ to the average value of the values of the surrounding four points. Is done (see the figure below).
** First, give the boundary conditions and the estimated value of the potential in the region, solve Eq. (12) in all regions, and repeat until the solutions in the region converge. This is the Jacobi method. ** The number of repeating steps for convergence is approximately proportional to the square of the number of data points [1].
This method is very slow to converge and impractical, but it provides a basis for understanding modern methods.
The Jacobi method and a slightly modified good calculation method [Gauss-Seidel method](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82% B9% EF% BC% 9D% E3% 82% B6% E3% 82% A4% E3% 83% 87% E3% 83% AB% E6% B3% 95) is slow to converge and impractical. Faster schemes include Sequential Over-relaxation (SOR) [1,2].
[1] ["Numerical Recipe in Sea Japanese Version-Recipe for Numerical Calculation in C Language"](https://www.amazon.co.jp/%E3%83%8B%E3%83%A5 % E3% 83% BC% E3% 83% A1% E3% 83% AA% E3% 82% AB% E3% 83% AB% E3% 83% AC% E3% 82% B7% E3% 83% 94-% E3% 82% A4% E3% 83% B3-% E3% 82% B7% E3% 83% BC-% E6% 97% A5% E6% 9C% AC% E8% AA% 9E% E7% 89% 88- C% E8% A8% 80% E8% AA% 9E% E3% 81% AB% E3% 82% 88% E3% 82% 8B% E6% 95% B0% E5% 80% A4% E8% A8% 88% E7% AE% 97% E3% 81% AE% E3% 83% AC% E3% 82% B7% E3% 83% 94-Press-William-H / dp / 4874085601 / ref = dp_return_1? _ Encoding = UTF8 & n = 465392 & s = books) Technical Review, 1993.
As something that can be viewed on the net
http://www.na.scitec.kobe-u.ac.jp/~yamamoto/lectures/appliedmathematics2006/chapter8.PDF
When
http://web.econ.keio.ac.jp/staff/ito/pdf02/pde2.pdf
I will list.
I posted an article about SOR method and Gauss-Seidel method.
[2] Qiita article, [[Science and technology calculation by Python] Comparison of convergence speeds of Jacobi method, Gauss-Seidel method, and SOR method for Laplace equation, partial differential equation, boundary value problem, numerical calculation](http: // qiita.com/sci_Haru/items/b98791f232b93690a6c3)
The code to solve the Poisson's equation in Example (2) by the SOR method explained in the article [2] is posted. Converges more than 10 times faster than the Jacobi method </ font>, so use this.
"""
Poisson equation:Numerical solution,SOR method
Also displays the electric field:
14 Aug. 2017
"""
%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation #Import methods for creating animations
import matplotlib.animation as animation
anim = []
#
delta_Lx=0.01
delta_Ly=0.01
LLx = 1 #Square width
LLy= 1
Lx = int(LLx/delta_Lx)
Ly = int(LLy/delta_Ly)
V = 5.0 #Voltage
convegence_criterion = 10**-5
phi = np.zeros([Lx+1,Ly+1])
phi_Bound= np.zeros([Lx+1,Ly+1])
phi_Bound[0,:] = V #boundary condition
#Charge density setting
eps0=1
charge= np.zeros([Lx+1,Ly+1])
c=1000
CC=int(Lx/4)
CC2=int(Ly/2)
for i in range(CC, CC2+1):
for j in range(CC,CC2+1):
charge[i,j] = c
#
# for SOR method
aa_recta=0.5*(np.cos(np.pi/Lx)+np.cos(np.pi/Ly))
omega_SOR_recta = 2/(1+np.sqrt(1-aa_recta**2)) #Acceleration parameters in SOR method
print("omega_SOR_rect=",omega_SOR_recta)
#for plot
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#Main
delta = 1.0
n_iter=0
conv_check=[]
while delta > convegence_criterion:
phi_in=phi.copy()
if n_iter % 100 ==0: #Monitoring of convergence status
print("iteration No =", n_iter, "delta=",delta)
conv_check.append([n_iter, delta])
for i in range(Lx+1):
for j in range(Ly+1):
if i ==0 or i ==Lx or j==0 or j==Ly:
phi[i,j] = phi_Bound[i,j]
else:
phi[i,j] = phi[i,j]+omega_SOR_recta *((phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1])/4-phi[i,j]+ (delta_Lx*delta_Ly/(4*eps0))*charge[i,j]) # SOR =Gauss-Seidel+ OR
delta = np.max(abs(phi-phi_in))
n_iter+=1
im=plt.imshow(phi,cmap='hsv')
anim.append([im])
#for plot
pp=plt.colorbar (orientation="vertical") #Color bar display
pp.set_label("phi", fontname="Ariel", fontsize=24)
plt.xlabel('X',fontname="Ariel", fontsize=24)
plt.ylabel('Y',fontname="Ariel", fontsize=24)
plt.title('Solution of the Poisson equation for electrostatic potential ')
#ani = ArtistAnimation(fig, anim, interval=10, blit=True,repeat_delay=1000)
#ani.save("t.gif", writer='imagemagick')
plt.show()
#Electric field calculation
Ex = np.zeros([Lx,Ly])
Ey = np.zeros([Lx,Ly])
for i in range(Lx):
for j in range(Ly):
Ex[i,j] = -(phi[i+1,j]-phi[i,j])/delta_Lx
Ey[i,j] = -(phi[i,j+1]-phi[i,j])/delta_Ly
plt.figure()
X, Y= np.meshgrid(np.arange(0, Lx,1), np.arange(0, Ly,1)) #Mesh generation
plt.quiver(X,Y,Ex,Ey,color='red',angles='xy',scale_units='xy', scale=40) #Plot vector field
plt.xlim([0,Lx]) #Range of X to draw
plt.ylim([0,Ly]) #Range of y to draw
#Graph drawing
plt.grid()
plt.draw()
plt.show()