マクスウェルの方程式

マクスウェルの方程式

マクスウェルの方程式は物理学の基礎方程式で、電磁気学の法則を表しています。

これらの方程式を数値的に解くには、偏微分方程式を離散化して計算する必要があります。

ここでは、簡単な例として、1次元の電磁波の伝播をシミュレーションしましょう。

まず、必要なライブラリをインポートします。

1
2
import numpy as np
import matplotlib.pyplot as plt

次に、計算領域と時間ステップを設定します。

1
2
3
4
5
nx = 200 # 空間グリッド数
nt = 100 # 時間ステップ数
dx = 0.01 # 空間グリッドの間隔
dt = 0.001 # 時間ステップ
c = 1.0 # 光速度

電場Ezと磁場Hyを初期化します。

1
2
3
4
5
Ez = np.zeros(nx)
Hy = np.zeros(nx)

# 初期条件を設定
Ez[nx//4] = 1.0

そして、時間ループを回します。

1
2
3
4
5
6
7
8
9
for t in range(nt):
    # 電場の更新
    Ezn = Ez.copy()
    for i in range(1, nx):
        Ez[i] = Ezn[i] + c * dt / dx * (Hy[i] - Hy[i-1])
    
    # 磁場の更新
    for i in range(nx-1):
        Hy[i] = Hy[i] - c * dt / dx * (Ez[i+1] - Ez[i])

最後に、結果をプロットします。

1
2
3
4
5
6
7
8
9
10
plt.subplot(2,1,1)
plt.plot(Ez)
plt.title('Electric Field')

plt.subplot(2,1,2)
plt.plot(Hy)
plt.title('Magnetic Field')

plt.tight_layout()
plt.show()

このコードでは、電場と磁場を離散化された Maxwell 方程式で時間発展させ、その様子をグラフに描画しています。

初期条件として与えた電場のパルスが、時間とともに伝播していく様子が確認できるはずです。

このようにして、Python を使えば Maxwell 方程式を数値的に解くことができ、電磁気現象のシミュレーションが可能になります。

[実行結果]

ソースコード解説

コードを詳しく説明します。

1. ライブラリのインポート

1
2
import numpy as np
import matplotlib.pyplot as plt

このコードではNumPyとMatplotlibの2つのライブラリをインポートしています。

NumPyは数値計算に使われる主要なライブラリで、配列の操作などができます。

Matplotlibはデータの可視化を行うライブラリで、グラフの描画に使われます。

2. シミュレーションパラメータの設定

1
2
3
4
5
nx = 200 # 空間グリッド数
nt = 100 # 時間ステップ数
dx = 0.01 # 空間グリッドの間隔
dt = 0.001 # 時間ステップ
c = 1.0 # 光速度

このコードではシミュレーションに必要なパラメータを設定しています。

  • nxは空間グリッドの総数で、ここでは$20$0と設定されています
  • ntは時間ステップの総数で、$100$と設定されています
  • dxは空間グリッドの間隔で、$0.01$と設定されています
  • dtは時間ステップの間隔で、$0.001$と設定されています
  • cは光速度で、ここでは$1.0$と規格化されています

3. 電場と磁場の初期化

1
2
Ez = np.zeros(nx)
Hy = np.zeros(nx)

このコードでは、電場Ezと磁場Hyを空のNumPy配列として初期化しています。
np.zeros(nx)はサイズがnxの$0$で初期化された1次元配列を作成します。

1
2
# 初期条件を設定
Ez[nx//4] = 1.0

このコードでは、電場Ezの初期条件を設定しています。
具体的には、配列のインデックスnx//4の位置(1/4の位置)に$1.0$を代入しています。
これにより、初期状態で電場の山ができます。

4. 時間ループ

1
2
3
4
5
6
7
8
9
for t in range(nt):
    # 電場の更新
    Ezn = Ez.copy()
    for i in range(1, nx):
        Ez[i] = Ezn[i] + c * dt / dx * (Hy[i] - Hy[i-1])
    
    # 磁場の更新
    for i in range(nx-1):
        Hy[i] = Hy[i] - c * dt / dx * (Ez[i+1] - Ez[i])

このコードが本質的な計算部分です。
時間ステップnt回分のループを回しながら、電場と磁場の値を更新していきます。

  • 電場の更新では、まずEznに現在のEzの値をコピーします。
    そしてi=1からnx-1までの各グリッド点で、Maxwell方程式に従った更新を行います。
  • 磁場の更新では、i=0からnx-2までの各グリッド点で、Maxwell方程式に従った更新を行います。

このようにして、離散化されたMaxwell方程式に従って、電場と磁場の値が次々と更新されていきます。

5. 結果の可視化

1
2
3
4
5
6
7
8
9
10
plt.subplot(2,1,1)
plt.plot(Ez)
plt.title('Electric Field')

plt.subplot(2,1,2)
plt.plot(Hy)
plt.title('Magnetic Field')

plt.tight_layout()
plt.show()

最後のこのコードでは、計算結果の可視化を行っています。

  • plt.subplot(2,1,1)は、グラフ領域を2行1列に分割し、1番目(上側)のサブプロットを選択します
  • plt.plot(Ez)で電場 Ezの値をプロットします
  • plt.title('Electric Field')でそのサブプロットのタイトルを設定します
  • 同様にplt.subplot(2,1,2)で2番目(下側)のサブプロットを選び、plt.plot(Hy)で磁場 Hyをプロットしています

plt.tight_layout()は余白を適切に設定し、plt.show()でグラフを表示します。

このようにして、Maxwell方程式に従った電磁場の時間と伝播をグラフで可視化しています。

以上が全体のコードの説明になります。

Maxwell方程式の離散化、数値計算、可視化がPythonで行われています。

結果解説

[実行結果]

上記のグラフについて詳しく説明します。

グラフは2つの subplot に分かれており、上側が電場(Ez)、下側が磁場(Hy)の時間を示しています。

電場(Ez)のグラフ

  • 横軸は空間グリッドの点($0$から$199$)を表しています。
  • 縦軸は電場の大きさを表しています。
  • 初期条件として、空間グリッドの$1/4$の位置(インデックス$50$付近)に電場の値$1.0$が設定されています。
  • 時間が経つにつれ、この電場の山がグラフの右側(正の方向)に伝播していきます。
  • これは、Maxwell方程式に従い、電場が光速度で伝播する様子を表しています。

磁場(Hy)のグラフ

  • 横軸は空間グリッドの点($0$から$198$)を表しています。
  • 縦軸は磁場の大きさを表しています 。
  • 初期状態では磁場はゼロですが、電場の変化により磁場が生成されます。
  • 時間とともに、電場に伴う磁場の山とその伝播が確認できます。
  • 電場と磁場はお互いに影響を与え合いながら伝播していきます。

つまり、このグラフは1次元の電磁波の伝播を可視化したものです。

初期条件として与えた電場のパルスが、Maxwell方程式に従って光速度で伝わり、その過程で磁場も生成・伝播する様子がわかります。

このようにして、電磁気現象をシミュレーションできます。