【Python】差分法で電位を求める(1)
- atliqu
- 2025年4月3日
C#でもやってみたけど、Pythonでやると全然簡単にできる。
上辺が10Vで、残り3辺がdΦ/dx=0という初期条件で差分法でやってみた。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
# Python 3.13.2, PyCharm 2024.3.5, numpy 2.2.4, matplotlib 3.10.1 import numpy as np import matplotlib.pyplot as plt fig = plt.figure() dl = 0.1 #微小長さ size = 10 #正方形の大きさ ll = int(size/dl) #分割数 volt = 10.0 #上辺の電位 c_err = 10**-3 #収束判定値 phi1 = np.zeros([ll+1, ll+1]) #電位(保存用） phi2 = np.empty([ll+1, ll+1]) #電位 err = 1.5 #誤差 n = 0 #収束条件モニタ用カウンタ while err > c_err: for i in range(ll+1): for j in range(ll+1): if i==0 : #上辺境界条件 Φ=V phi2[i,j] = volt elif j==0 and i==ll: phi2[i,j] = (phi1[i-1,j] + phi1[i,j+1])/2 #左下頂点境界条件 dΦ/dx=0, dΦ/dy=0 elif j==ll and i==ll: phi2[i, j] = (phi1[i-1,j] + phi1[i,j-1])/2 #左下頂点境界条件 dΦ/dx=0, dΦ/dy=0 elif j==0: phi2[i,j] = (phi1[i-1,j] + phi1[i+1,j] + phi1[i,j+1])/3 #左辺境界条件 dΦ/dx=0 elif i==ll: phi2[i,j] = (phi1[i,j-1] + phi1[i,j+1] + phi1[i-1,j])/3 #下辺境界条件 dΦ/dy=0 elif j==ll: phi2[i,j] = (phi1[i-1,j] + phi1[i+1,j] + phi1[i,j-1])/3 # 右辺境界条件 dΦ/dx=0 else: phi2[i,j] = (phi1[i+1,j] + phi1[i-1,j] + phi1[i,j+1] + phi1[i,j-1])/4 if n % 100 == 0: # 収束状況のモニタリング print(n,err) err = np.max(abs(phi2-phi1)) phi1, phi2 = phi2, phi1 n += 1 im = plt.imshow(phi1, cmap='rainbow') plt.colorbar() plt.xlabel('X') plt.ylabel('Y') plt.show() |
モデルが単純だから絵的に面白くない。C#に比べると、実行速度が相当遅い。高速化の手法があるみたいだけど、まあ、別に急がないので当面このままにする。
実行しているかどうかわからないので、誤差を100回毎に表示するようにした。
C#と比べると、グラデーションの表示とか圧倒的に簡単。 ほぼ何もしなくていい。