Possion方程式と静電ポテンシャル② 実装編【Pythonと物理#11】

今回のターゲット：Poisson方程式
Pythonで計算する
実行結果
最後に

今回のターゲット：Poisson方程式

　前回はポワソン方程式を数値計算して静電ポテンシャル $\phi\left(x, y\right)$ を求めるための準備編でした。詳細は以下の記事をご覧ください。

今回のターゲット：Poisson方程式　前回は１次元熱拡散方程式を計算してみました。　今回は電磁気学を扱ってみようかと思います。電磁気学の中でも2次元のPoisson方程式を解いて、静電ポテンシャル$\phi\left(x, y\right...

今回はいよいよ実装してみたいと思います。

勉強の際に参考にさせていただいた書籍もご紹介させていただきます。具体的なPythonのコードも載っており、大変勉強になる１冊です。なお、あくまでも参考にさせていただいている形であり、本記事の内容において数学的・物理学的な誤りがある場合は全て当記事の筆者に責任があります。

リンク

Pythonで計算する

コードは下のようにしました。具体的な説明はコード内にコメントアウトで記載していますので、参考にしていただければと思います。

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csr_matrix as sparse_mx
from scipy import sparse
from scipy.sparse.linalg import spsolve
from scipy.linalg import lu_factor, lu_solve

#x(y)方向の差分行列と単位行列を生成する関数
def generate_diff_mx(q):
    dq = q[1] - q[0]
    nq = q.size
    imx = np.identity(nq, dtype=int)
    imx_1 = np.roll(imx, 1, axis=1)
    imx_1t = imx_1.transpose()
    D2 = sparse_mx(imx_1 - 2*imx + imx_1t)/(dq**2)
    return imx, D2

#電荷分布を生成する関数
def generate_rho(x, y):
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    rho = np.zeros((x.size, y.size), dtype=float)
    X = x.size//2
    Y = y.size//2
    #簡単のために電荷をq=1として、面積dxdyの1つ分のmeshに点電荷が存在すると仮定
    rho[X, Y] = 1.0 / (dx * dy)
    return rho

#境界条件(Boundary Condition:BC)を生成する関数
def generate_BC(x, y):
    #一旦、全てFalseにする
    BC = np.full((x.size, y.size), False)
    #境界をTrueにする
    BC[0, :] = True
    BC[-1, :] = True
    BC[:, 0] = True
    BC[:, -1] = True
    return BC

#境界条件を考慮した連立方程式を解く関数
def solve_simultaneous_eq(D2, rho, BC):
    #rhoとBCを１次元化する
    rho_flat = rho.ravel()
    BC_flat = BC.ravel()
    #境界内部だけで計算する準備
    within_boundary = ~BC_flat #"~"はビット反転演算子でFalseをTrueに、TrueをFalseにする
    #境界内部のみの行列を作る
     #１つ目の[]でwithin_boundaryのTrueの"行"のみ取り出し、２つ目の[]でwithin_boundaryのTrueの"列"のみ取り出す
    D2_wb = D2[within_boundary][:, within_boundary] #2次元であることに注意
    rho_wb = rho_flat[within_boundary] #1次元化していることに注意
    #境界条件を考慮した連立方程式を解く
    sol = spsolve(D2_wb, -rho_wb)
    #静電ポテンシャルを格納するphiを準備、境界内部も一旦0にする
    phi = np.zeros_like(rho_flat)
    #境界内部に計算結果を格納、このとき境界部分は0にしておく
    phi[within_boundary] = sol
    return phi

#Poisson方程式を解く関数
def solve_poisson(x, y, rho, BC):
    nx = x.size
    ny = y.size
    
    #まずは1次元の差分行列と単位行列を生成
    imx_x, d2x = generate_diff_mx(x)
    imx_y, d2y = generate_diff_mx(y)
    #1次元の差分行列と単位行列から2次元の差分行列を生成
    D2x = sparse.kron(d2x, imx_y)
    D2y = sparse.kron(imx_x, d2y)
    #2次元の差分行列からラプラシアンを生成
    D2 = D2x + D2y
    #境界条件を考慮した連立方程式を解く
    phi = solve_simultaneous_eq(D2, rho, BC)
    #2次元配列に戻す
    phi = phi.reshape(rho.shape)
    return phi

#各種パラメータ設定、各種関数呼び出し
Lx, Ly = 2.0, 2.0 #考えている系のx, y方向の長さ
nx, ny = 101, 101 #考えている系のx, y方向の分割数
x = np.linspace(-Lx/2, Lx/2, nx)
y = np.linspace(-Ly/2, Ly/2, ny)
rho = generate_rho(x, y)
BC = generate_BC(x, y)
 #x, yを2次元化
xx, yy = np.meshgrid(x, y, indexing='ij') #indexing='ij'と指定することでxを行、yを列に揃えてくれる
phi = solve_poisson(x, y, rho, BC)

#以下はグラフ作成のためのコード
fig, ax = plt.subplots(dpi=240)
image = ax.pcolormesh(xx, yy, phi, shading='nearest', cmap='rainbow')
fig.colorbar(image, ax=ax)
ax.set_aspect(1)

ax.set_xlabel('x')
ax.set_ylabel('y')

plt.xticks([-1.0, -0.5, 0, 0.5, 1.0])
plt.yticks([-1.0, -0.5, 0, 0.5, 1.0])

ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

plt.show()