中心差分法で電場ベクトルをシミュレーションする

Question

### 

google colab環境下でPythonを用いてポアソン方程式で電位を求めてから, その結果を用いて電場の図を描いています.
**E** = -∇φ
を離散化してベクトル場を表示しようとしたら, 以下のようになりました.
### 発生している問題・エラーメッセージ

![イメージ説明](781075d5452af12c1fdbd3d97bfe49a1.png)

### 該当のソースコード
```python 
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 - 1, j]) / 2*delta_Lx
        Ey[i, j] = -(phi[i, j + 1] - phi[i, j - 1]) / 2*delta_Ly

fig = plt.figure(figsize=(4, 4))

X, Y = np.meshgrid(np.arange(0, Lx, 1), np.arange(0, Ly, 1)) # メッシュ生成
plt.quiver(X, Y, Ex, Ey, color = 'red', angles = 'xy', 
           scale_units = 'xy', scale = 40, 
           label = "Solution of the Poisson equation for electrostatic field") # ベクトル場をプロット

plt.xlim([0, Lx]) # 描くXの範囲
plt.ylim([0, Ly]) # 描くYの範囲
```
### 試したこと

delta_Lxとdelta_Lyを逆にした.

### 補足情報（FW/ツールのバージョンなど）
```python
delta_Lx = 0.01 # グリッドx幅
delta_Ly = 0.01 # グリッドy幅
LLx = 1
LLy = 1

Lx = int(LLx / delta_Lx)
Ly = int(LLy / delta_Ly)

V = 5.0 # 電圧
convegence_criterion = 10**-5

phi = np.zeros([Lx+1, Ly+1])
phi_Bound = np.zeros([Lx+1, Ly+1])
phi_Bound[0, :] = V # 全ての行の0列目を取得
```

# 電荷密度の設定
```python
eps0 = 1
charge = np.zeros([Lx+1, Ly+1])
c = 1500

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
’’’Python
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)) # SOR法における加速パラメータ
print("omega_SOR_rect=", omega_SOR_recta)
```
# メイン
```python
delta = 1.0
n_iter = 0
conv_check = []
while delta > convegence_criterion:
    phi_in = phi.copy()
    if n_iter % 100 == 0: # 収束状況のモニタリング
        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 = ガウス-ザイデル + OR
    delta = np.max(abs(phi - phi_in))

    n_iter += 1

    im = plt.imshow(phi, cmap = 'jet')
    anim.append([im])
```
前進差分法により離散化をしたときの電場のグラフです.
![イメージ説明](156319d1c8da9bb316adae636ae5ced1.png)
このように中心差分法の時も描かれないといけませんが, 上記に貼った通り描かれていないのが現状です,

Accepted Answer

2*delta_Lx → (2*delta_Lx)
2*delta_Ly → (2*delta_Ly)

カッコを付けないと、delta_L*は分母じゃなくて分子に入る

発生している問題・エラーメッセージ

該当のソースコード

試したこと

補足情報（FW/ツールのバージョンなど）

電荷密度の設定

メイン

関連した質問