ラプラス方程式

ラプラス方程式

ラプラス方程式は、電場や重力場、流体力学など多くの物理現象を記述するのに使われる重要な方程式です。

ラプラス方程式は次の形をしています：
$$
\nabla^2 \phi = 0
$$
ここで、$(\phi)$はポテンシャル関数です。

プログラム例

Pythonを用いて、ラプラス方程式の数値解を3Dグラフとして描画するコードを示します。

ここでは、単純な境界条件を用いて解を求めます。

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
43
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータの設定
L = 1.0  # 領域の長さ
N = 50   # 空間分割数
dx = L / N

# 初期条件と境界条件の設定
phi = np.zeros((N, N))
phi[:, 0] = 1.0  # 左側の境界条件 phi = 1

# ラプラス方程式の数値解法（反復法）
def solve_laplace(phi, tol=1e-5, max_iter=10000):
    for _ in range(max_iter):
        phi_new = phi.copy()
        phi_new[1:-1, 1:-1] = 0.25 * (phi[1:-1, :-2] + phi[1:-1, 2:] + phi[:-2, 1:-1] + phi[2:, 1:-1])
        
        # 収束判定
        if np.max(np.abs(phi_new - phi)) < tol:
            break
        phi = phi_new
    return phi

# ラプラス方程式を解く
phi = solve_laplace(phi)

# グラフの作成
x = np.linspace(0, L, N)
y = np.linspace(0, L, N)
X, Y = np.meshgrid(x, y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, phi, cmap='viridis')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Potential $\phi$')
ax.set_title('Solution of Laplace Equation')

plt.show()

上記のソースコードを実行することで、ラプラス方程式の数値解を3Dグラフとして視覚化することができます。

[実行結果]

ソースコード解説

以下に、ソースコードの各部分を詳しく説明します。

1. パラメータの設定

1
2
3
L = 1.0  # 領域の長さ
N = 50   # 空間分割数
dx = L / N
  • L: 計算領域の一辺の長さです。ここでは$1.0$に設定されています。
  • N: 計算領域を分割する数です。ここでは$50$に設定されています。
  • dx: 空間ステップサイズで、$( L ) $を$ ( N ) $で割ったものです。

2. 初期条件と境界条件の設定

1
2
phi = np.zeros((N, N))
phi[:, 0] = 1.0  # 左側の境界条件 phi = 1
  • phi: 計算領域内のポテンシャルを格納する2次元配列です。初期状態ではすべてゼロに設定されています。
  • phi[:, 0] = 1.0: 左側の境界におけるポテンシャルを$1$に設定しています。他の境界条件はゼロのままです。

3. ラプラス方程式の数値解法（反復法）

1
2
3
4
5
6
7
8
9
10
11
def solve_laplace(phi, tol=1e-5, max_iter=10000):
    for _ in range(max_iter):
        phi_new = phi.copy()
        phi_new[1:-1, 1:-1] = 0.25 * (phi[1:-1, :-2] + phi[1:-1, 2:] + phi[:-2, 1:-1] + phi[2:, 1:-1])
        
        if np.max(np.abs(phi_new - phi)) < tol:
            break
        phi = phi_new
    return phi

phi = solve_laplace(phi)
  • solve_laplace: 反復法を用いてラプラス方程式を解く関数です。
    • phi_new = phi.copy(): 現在のポテンシャル分布をコピーして、新しいポテンシャル分布を格納するための配列を作成します。
    • phi_new[1:-1, 1:-1] = 0.25 * (phi[1:-1, :-2] + phi[1:-1, 2:] + phi[:-2, 1:-1] + phi[2:, 1:-1]): 内部の格子点に対して、周囲の格子点のポテンシャルの平均値で更新します。
    • if np.max(np.abs(phi_new - phi)) < tol: 更新前後のポテンシャル分布の差が許容誤差 tol 未満であれば収束とみなし、反復を終了します。

4. グラフの作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
x = np.linspace(0, L, N)
y = np.linspace(0, L, N)
X, Y = np.meshgrid(x, y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, phi, cmap='viridis')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Potential $\phi$')
ax.set_title('Solution of Laplace Equation')

plt.show()
  • x, y: 領域内の格子点の座標を生成するための1次元配列です。
  • X, Y: np.meshgridを使って、2次元グリッドの座標を作成します。
  • fig, ax: Matplotlibを使って3Dプロットを作成します。
  • ax.plot_surface: 3Dグラフの表面をプロットし、viridisカラーマップを使用してポテンシャルの分布を色で表現します。
  • ax.set_xlabel, ax.set_ylabel, ax.set_zlabel: X軸、Y軸、Z軸のラベルを設定します。
  • ax.set_title: グラフのタイトルを設定します。

グラフ解説

[実行結果]

グラフに表示される内容は、以下の通りです。

X軸とY軸:

計算領域内の座標を表します。
ここでは、$(0)$から$(1)$の間を均等に分割しています。

Z軸:

ポテンシャル$ (\phi) $の値を示します。
左側の境界条件として設定された$ ( \phi = 1 ) $から、内部でのポテンシャルの分布が示されています。

ポテンシャルの分布:

色の変化によってポテンシャルの高低が視覚的に表現されています。
左側の境界ではポテンシャルが高く、その他の領域ではポテンシャルが低いことが分かります。

このグラフは、境界条件とラプラス方程式によって決まる2次元領域内の静電ポテンシャルの分布を3Dで示しており、ポテンシャルの変化を視覚的に理解するのに役立ちます。