1. はじめに

　Python三次元ポアソン方程式による静電場解析に関しては先の投稿において報告した(1)(2)(3)．
　本稿では，二次元ポアソン方程式による静電場解析として十分長い直線状電荷分布を対象に下記を実施する．
(1)二次元ポアソン方程式を差分化
(2)繰り返し法によるプログラムを作成．
①電位を計算
②得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証

2. 例題の問題設定

　各軸のメッシュ数は100とし，中心に一様に帯電した十分長い直線を配置する．

　　　　　図1　問題設定-十分長い直線状電荷の配置

3. ポアソン方程式の離散化

　下記に手順を記す．
　ここで注意が必要な点はρは線密度であることである．

4. プログラム

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import math
from numpy import log as ln

delta_L=1.0

nx = 100
ny = 100
nt = 10000
xmin = 0
xmax = 100*delta_L
ymin = 0
ymax = 100*delta_L

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

p  = np.zeros((ny, nx))
pd = np.zeros((ny, nx))
b  = np.zeros((ny, nx))
x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)

eps0=1
charge= np.zeros((ny, nx))
Q1=1
Q2=Q1/delta_L**2

charge[50,50] = Q2

for it in range(10000):

    pd = p.copy()

    p[1:-1,1:-1] = (pd[1:-1, 2:] + pd[1:-1, :-2] +  
         pd[2:, 1:-1] + pd[:-2, 1:-1] + (delta_L**2/(1*eps0))*charge[1:-1, 1:-1] ) / 4 

    p[0, :] = 0
    p[ny-1, :] = 0
    p[:, 0] = 0
    p[:, nx-1] = 0

def plot2D(x, y, p):
    fig = plt.figure(figsize=(11, 7), dpi=100)
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(x, y)
    surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
    ax.view_init(30, 225)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('p')

plot2D(x, y, p)

fig = plt.figure(figsize=(11,7), dpi=100)
fig

xmin = 0
xmax = 100
ymin = 0
ymax = 100

x  = np.linspace(xmin, xmax, nx)
y  = np.linspace(xmin, xmax, ny)

X, Y = np.meshgrid(x, y)

plt.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plt.ylim([35,65])
plt.xlim([35,65])

fig = plt.figure(figsize=(11,7), dpi=100)
fig
x1  = np.linspace(0, 100, nx)
y1  = np.linspace(xmin, xmax, ny)
for i in range(nx):
    y1[i]=p[i,50]
plt.plot(x1,y1)
plt.title("A simple line graph")
plt.xlim([44,56])

L= 99
Ex = np.zeros([L,L])
Ey = np.zeros([L,L])

for i in range(L):
    for j in range(L):
        Ey[i,j] = -(p[i+1,j]-p[i-1,j])/2/delta_L
        Ex[i,j] = -(p[i,j+1]-p[i,j-1])/2/delta_L

print("Riron Keisan E")
for i in range(1,6):
    print("x=",i*delta_L,"E=",Q1/(2*np.pi)/(i*delta_L))

print("Riron Keisan V")
for i in range(1,6):
    print("x=",i*delta_L,"V=",-Q1/(2*np.pi)*ln(i*delta_L))

print("Keisan E")
for i in range(1,6):
    print("x=",i*delta_L,"E=",Ey[50+i,50]) 

print("Keisan V")
for i in range(1,6):
    print("x=",i*delta_L,"V=",p[50+i,50]-p[50+int(1/delta_L),50])

fig = plt.figure(figsize=(11,7), dpi=100)
fig

X, Y= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1))  # メッシュ生成
plt.quiver(X,Y,Ex,Ey,color='red',angles='xy',scale_units='xy', scale=0.25/delta_L) # ベクトル場をプロット

plt.xlim([45,55])
plt.ylim([45,55])

plt.grid()
plt.draw()

plt.show()