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

Question

google colab環境下でPythonを用いてポアソン方程式で電位を求めてから, その結果を用いて電場の図を描いています.
E = -∇φ
を離散化してベクトル場を表示しようとしたら, 以下のようになりました.

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

該当のソースコード

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/ツールのバージョンなど）

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列目を取得

電荷密度の設定

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)

メイン

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])

前進差分法により離散化をしたときの電場のグラフです.

このように中心差分法の時も描かれないといけませんが, 上記に貼った通り描かれていないのが現状です,

Comments

[jbpb0]

pythonのコードの一番最初の行のすぐ上に
```python
だけの行を追加してください
また、pythonのコードの一番最後の行のすぐ下に
```
だけの行を追加してください

現状、コードがとても読み辛いです

[ikadzuchi]

何が問題なのか伝わりません。
どのような結果になるべきで、現在の結果は求める結果とどのように違うのかを説明してください。

[nagi-kimura]

pythonでしかコードを書いていないんですけど.....

[nagi-kimura]

電場のようにベクトル図が描かれないといけないのですが, ベクトルが描かれていないのでここの点について投稿をしました.

[jbpb0]

現状では、pythonコードのインデントが無くなっているので、コードが解読できません
質問を編集して、私の最初のコメントに書いたように、

```python
(pythonのコード全部)
```

としてくれたら、pythonコードのインデントが他人にも分かり、コードをちゃんと読むことができます

[nagi-kimura]

しました!

[jbpb0]

''''''''''''''''''''''''''''''''''''''''''''''''''Python
と書いてはダメです
```python
と書いてください

同様に
''''''''''''''''''''''''''''''''''''''''''''''''''
と書いてはダメです
```
と書いてください

そう書かないと、他人が質問を見た時にpythonコードのインデントが再現されず、コードの詳細を解読するのが難しいのです

[nagi-kimura]

そうなんですね!
python3ヶ月目なので全然知りませんでした！

[jbpb0]

’’’Python
ではないです
```python
です

’’’
も違います
```
です

自分で手入力するのではなく、この私のコメントからコピペしてください

なお、この書き方は、teratailで質問のコードを他人が読みやすくするためのルールです
実行させる実際のpythonコードには、そんなもの書かないでくださいね
あくまでも、teratailにコードを書く場合だけ、それを追加するのです

[nagi-kimura]

すみません. 申し訳ありません. 心からお詫び申し上げます.

[jbpb0]

うまくいってる前進差分のコードとの違いは、下記だけですか？
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

[nagi-kimura]

前進差分は
Ex[i, j] = -(phi[i + 1, j] - phi[i, j]) / delta_Lx
Ey[i, j] = -(phi[i, j + 1] - phi[i, j]) / delta_Ly
となり, 中心差分との違いはそのコードのみです.
後退差分もうまくいっていて
Ex[i, j] = -(phi[i, j] - phi[i - 1, j]) / delta_Lx
Ey[i, j] = -(phi[i, j] - phi[i, j - 1]) / delta_Ly
となります. 参考までに...
中心差分のみうまく行きません.

Answer 1

2delta_Lx → (2delta_Lx)
2delta_Ly → (2delta_Ly)

カッコを付けないと、delta_L*は分母じゃなくて分子に入る

Comments

[nagi-kimura]

只今実行したら, 求めていた画像が正しく出力されました.
ありごとうございました.