ベルヌーイ方程式

April 21, 2024

ベルヌーイ方程式

ベルヌーイ方程式は、流体力学における基本的な方程式の一つであり、流れの速さと圧力の関係を表現します。

ベルヌーイ方程式は以下の形で表されます：

$$
P + \frac{1}{2} \rho v^2 + \rho gh = \text{constant}
$$

ここで、各項の意味は次の通りです：

$(P)$ は流体の圧力
$(\rho)$ は流体の密度
$(v)$ は流体の速度
$(g)$ は重力加速度
$(h)$ は流体の高さ

この方程式は、流体のエネルギー保存則を表現しており、特定のポイントでの圧力、速度、高さの変化について関係を示しています。

ベルヌーイ方程式は、流体の流れに関するさまざまな問題で応用されます。

ベルヌーイ方程式をPythonでグラフ化するためには、方程式の関係を適切に表現し、それを可視化する必要があります。

以下に、ベルヌーイ方程式を用いて流体の速度と圧力の関係をグラフ化する例を示します。

import numpy as np
import matplotlib.pyplot as plt

# 定数と流体の特性パラメータ
rho = 1000  # 流体の密度 (kg/m^3)
g = 9.81    # 重力加速度 (m/s^2)
h = 10      # 流体の高さ (m)
constant_term = 100000  # 定数項

# 速度の範囲を定義
v_values = np.linspace(0, 10, 100)  # 0から10までの速度 (m/s)

# 圧力を計算する関数
def calculate_pressure(v):
    return constant_term - 0.5 * rho * v**2 - rho * g * h

# 速度と圧力の関係を計算
pressure_values = calculate_pressure(v_values)

# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')
plt.xlabel('Velocity (m/s)', fontsize=12)
plt.ylabel('Pressure (Pa)', fontsize=12)
plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)
plt.grid(True)
plt.legend()
plt.show()

このコードでは、与えられた流体の密度 $(\rho)$、重力加速度 $(g)$、流体の高さ $(h)$、および定数項を使用して、ベルヌーイ方程式に基づいて速度と圧力の関係を計算しています。

そして、計算結果をグラフ化しています。

得られるグラフは、速度と圧力の関係を示しており、流体力学におけるベルヌーイ方程式の応用を視覚的に理解するのに役立ちます。

[実行結果]

ソースコード解説

ソースコードの詳細な説明を行います。

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

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

numpyは数値計算を効率的に行うためのライブラリであり、npとしてインポートされています。
matplotlib.pyplotはデータの可視化（グラフ描画）を行うためのライブラリであり、pltとしてインポートされています。

2. 定数と流体の特性パラメータの定義

rho = 1000  # 流体の密度 (kg/m^3)
g = 9.81    # 重力加速度 (m/s^2)
h = 10      # 流体の高さ (m)
constant_term = 100000  # 定数項

rhoは流体の密度 $（kg/m^3）$を表します。
gは重力加速度 $（m/s^2）$を表します。
hは流体の高さ $（m）$を表します。
constant_termはベルヌーイ方程式の定数項であり、圧力を計算する際に使用されます。

3. 速度の範囲を定義

1	v_values = np.linspace(0, 10, 100) # 0から10までの速度 (m/s)

np.linspace(0, 10, 100)は$0$から$10$までの速度を等間隔で$100$個のデータポイントに区切って生成します。

4. 圧力を計算する関数の定義

1 2	def calculate_pressure(v): return constant_term - 0.5 * rho * v*2 - rho g * h

calculate_pressure(v)は、速度 vに対する圧力を計算する関数です。
ベルヌーイ方程式に基づいて圧力を計算しています。

5. 速度と圧力の関係を計算

1	pressure_values = calculate_pressure(v_values)

v_valuesに対してcalculate_pressure関数を適用して、速度に対する圧力の値を計算します。

6. グラフのプロット

plt.figure(figsize=(8, 6))
plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')
plt.xlabel('Velocity (m/s)', fontsize=12)
plt.ylabel('Pressure (Pa)', fontsize=12)
plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)
plt.grid(True)
plt.legend()
plt.show()

plt.figure(figsize=(8, 6))で新しい図を作成し、サイズを指定します。
plt.plot(v_values, pressure_values, label='Pressure vs. Velocity')で速度と圧力の関係をプロットします。
plt.xlabel('Velocity (m/s)', fontsize=12)で$x$軸（横軸）のラベルを設定します。
plt.ylabel('Pressure (Pa)', fontsize=12)で$y$軸（縦軸）のラベルを設定します。
plt.title('Bernoulli Equation: Pressure vs. Velocity', fontsize=14)でグラフのタイトルを設定します。
plt.grid(True)でグリッド線を表示します。
plt.legend()で凡例を表示します。
plt.show()でグラフを表示します。

このプログラムは、ベルヌーイ方程式に基づいて速度と圧力の関係を計算し、matplotlibを使用してグラフ化しています。

計算された圧力値は、速度が増加するにつれてどのように変化するかを視覚化することができます。

結果解説

[実行結果]

以下は、グラフに表示される内容の詳細な説明です：

1. x軸（横軸）：速度 (Velocity)

グラフの横軸は速度 $（m/s）$を表しています。
速度の範囲は$0$から$10$までの値で、$100$個の等間隔なデータポイントがプロットされます。
この範囲内で流体の速度がどのように変化するかを示します。

2. y軸（縦軸）：圧力 (Pressure)

グラフの縦軸は圧力（$Pa$：パスカル）を表しています。
圧力はベルヌーイ方程式に基づいて計算され、速度の関数としてプロットされます。
速度が変化するとともに、流体の圧力がどのように変動するかを示します。

3. グラフのラベルとタイトル

$x$軸（横軸）のラベルは「Velocity ($m/s$)」として設定されています。
これは流体の速度を表します。
$y$軸（縦軸）のラベルは「Pressure ($Pa$)」として設定されています。
これは流体の圧力を表します。
グラフのタイトルは「Bernoulli Equation: Pressure vs. Velocity」として設定されています。
これはベルヌーイ方程式に基づく圧力と速度の関係を示しています。

4. グリッド線と凡例

グラフにはグリッド線が表示されており、視覚的なガイドとして役立ちます。
凡例（Legend）は「Pressure vs. Velocity」として表示され、プロットされたデータの内容を説明しています。

このグラフは、流体の速度が増加すると圧力がどのように変化するかを示しています。

ベルヌーイ方程式は流体力学で重要な法則であり、このグラフを通じて流体の基本的な特性を理解することができます。

円柱

April 20, 2024

円柱

円柱の方程式は一般的に以下のように表されます：

$$
x^2 + y^2 = r^2
$$

ここで、$( r ) $は円柱の半径を表します。

円柱の方程式をPythonでグラフ化するためには、MatplotlibとNumPyを使用します。
具体的なコード例を示します：

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 円柱の半径
r = 3

# メッシュグリッドの作成
theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(-5, 5, 100)
Theta, Z = np.meshgrid(theta, z)
X = r * np.cos(Theta)
Y = r * np.sin(Theta)

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, color='skyblue', alpha=0.6)

# グラフの設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Cylinder')

# グラフの表示
plt.show()

このコードでは、円柱の方程式 $ ( x^2 + y^2 = r^2 ) $を使用して円柱の表面をプロットしています。

円柱の半径 $( r ) $を変更することで、円柱のサイズを調整することができます。

グラフは3次元で表示され、円柱の構造が視覚化されます。

[実行結果]

ソースコード解説

ソースコードの詳細を説明します。

1. NumPyの利用:

import numpy as np
``
`NumPy`は数値計算やデータ操作に用いられるPythonのライブラリです。
ここでは**円柱の座標**を計算するために使用します。

#### 2. **Matplotlibの利用**:
```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

Matplotlibはデータ可視化のためのPythonライブラリであり、mpl_toolkits.mplot3dモジュールは3次元プロットのための機能を提供します。

3. 円柱の半径を設定:

r = 3

円柱の半径 $ ( r ) $を設定します。
この値を変更すると、描画される円柱のサイズが変わります。

4. メッシュグリッドの作成:

theta = np.linspace(0, 2*np.pi, 100)
z = np.linspace(-5, 5, 100)
Theta, Z = np.meshgrid(theta, z)
X = r * np.cos(Theta)
Y = r * np.sin(Theta)

np.linspace関数を使用して角度 $( \theta ) $と高さ $( z ) $の範囲を定義し、np.meshgrid関数でメッシュグリッドを作成します。

そして、円柱の表面の$ ( x ) $座標と$ ( y ) $座標を計算します。

5. 3Dプロットの作成:

1
2
3

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, color='skyblue', alpha=0.6)

plt.figure()で新しい図を作成し、fig.add_subplot(111, projection='3d')で3次元サブプロットを追加します。

ax.plot_surface()で円柱の表面をプロットします。

6. グラフの設定:

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Cylinder')

軸ラベルやタイトルなど、グラフの表示設定を行います。

7. グラフの表示:

1	plt.show()

最後に plt.show() を呼び出してグラフを表示します。

これにより、円柱の表面が3次元グラフとして描画されます。
円柱の側面は青色で描かれ、透明度は$0.6$に設定されています。

円柱の底面や側面の形状が視覚的に表現され、円柱の構造がわかりやすくなります。

マリーバインッケ方程式

April 19, 2024

マリーバインッケ方程式

マリーバインッケ方程式は、3次元空間における代表的な参照集合の1つです。

この方程式は以下のように定義されます。

$$
z = -0.5 * (x^2 + y^2 - 1)^3 - 0.7 * (x^2 + y^2 - 1)^2 - (x^2 - y^2) / 8
$$

この方程式は、$x$、$y$、$z$の3変数で表された非線形の方程式です。

方程式の右辺は3つの項から構成されています。

1. 第1項:$ -0.5 * (x^2 + y^2 - 1)^3$

この項は、単位球面 $(x^2 + y^2 = 1)$からの距離の3乗に$-0.5$を乗じたものです。
球面から離れるほど値が大きくなります。

2. 第2項: $-0.7 * (x^2 + y^2 - 1)^2 $

この項は、単位球面からの距離の2乗に$-0.7$を乗じたものです。
球面から離れるほど値が大きくなりますが、第1項ほど急激ではありません。

3. 第3項: $-(x^2 - y^2)/8$

この項は、$x^2 - y^2$、つまり$x$軸と$y$軸に沿った楕円の方程式に$-1/8$を乗じたものです。
これにより、$x=y$の直線に沿って値が小さくなる効果があります。

これら3つの項の組み合わせにより、マリーバインッケ曲面の特徴的な形状が生まれます。

具体的には以下のような特徴があります:

中心が最小値$(z=-1.462)$となる深い窪み
8つの放射状に伸びる渦巻き状の突起
$x$軸と$y$軸に関する対称性
複雑な3次元形状

このマリーバインッケ曲面は、フラクタル幾何学や力学系の研究において重要な役割を果たしています。

また、この方程式は最適化問題の目的関数としても使われることがあります。

高度な数学的性質を持つユニークな曲面として知られています。

ソースコード

次のコードは、NumPy と Matplotlib を使用してマリーバインッケ方程式の3D曲面をプロットします。

import numpy as np
import matplotlib.pyplot as plt

# マリーバインッケ方程式
def mariebyinke(x, y):
    return -0.5 * (x**2 + y**2 - 1)**3 - 0.7 * (x**2 + y**2 - 1)**2 - (x**2 - y**2) / 8

# グリッドデータの生成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = mariebyinke(X, Y)

# 3D プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis')
ax.set_title('Mariebyinke Surface')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

mariebyinke 関数は、指定された$ x $と$ y $の値から方程式の値を計算します。
np.linspace を使って$ x $と$ y $の値の範囲を設定し、np.meshgrid でグリッドデータを生成します。
mariebyinke 関数を使って$ Z $値 (方程式の値) を計算します。
plt.figure と add_subplot を使って新しい図とサブプロットを作成します。
plot_surface を使って3D曲面をプロットします。
rstride と cstride は表示の密度を制御します。
軸ラベルとタイトルを設定します。
plt.show() でプロットを表示します。

このコードを実行すると、マリーバインッケ曲面の3Dプロットが表示されます。

曲面のカラーマップは viridis を使用していますが、必要に応じて変更できます。

[実行結果]

グラフ解説

生成されるグラフには、以下の特徴があります:

1. 形状

グラフは、複雑な3次元曲面の形状を持っています。
曲面は、中央部分が窪んでいる形状をしており、外側に向かって持ち上がっています。
曲面には、中心から放射状に伸びる8つの渦巻き状の突起があります。

2. 対称性

曲面は、$x$軸と$y$軸に関して対称な形状をしています。
つまり、$x=0$と$y=0$を回転軸として、曲面を$180$度回転させると同じ形状になります。

3. 極値

曲面の中心$(x=0, y=0)$が最小値となっています。
8つの渦巻き状の突起の先端が局所的な極大値となっています。

4. カラーマップ

デフォルトでは、曲面の高さ($z$値)に応じて、viridisカラーマップが適用されています。
青系の色が低い値、赤系の色が高い値を表しています。

5. 軸

$x$軸、$y$軸、$z$軸がそれぞれ表示されています。
$x$軸と$y$軸の範囲は$-2$から$2$までの値が表示されています。

6. タイトル

グラフのタイトルは「Mariebyinke Surface」と表示されています。

7. 視点

3D空間上の視点は調整可能で、様々な角度から曲面を観察できます。

このグラフは、マリーバインッケ方程式によって生成された特徴的な3次元曲面の形状を可視化しています。

この曲面は数学的に興味深い性質を持っており、様々な分野で研究されています。

楕円錐の方程式

April 18, 2024

楕円錐の方程式

楕円錐の方程式は次のように表されます:

$$
(z - w)^2 / c^2 = (x - u)^2 / a^2 + (y - v)^2 / b^2
$$

$(x, y, z) $は空間内の任意の点の座標
$(u, v, w) $は楕円錐の頂点の座標
$a$, $b$, $c$ はそれぞれ$x$, $y$, $z$軸方向の半径

この方程式は以下のようにも書くことができます:

$$
z = w ± c/a * sqrt(a^2 - (x - u)^2 - (y - v)^2)
$$

符号の$+$ は$z >= w$の上半分、符号の$-$ は$z <= w$の下半分を表します。

この方程式は、$(x, y)$平面上の点$(x-u, y-v)$から頂点$(u, v, w)$を通る直線の長さをcで割ったものが、楕円錐の半径 $a$, $b$の比になることを表しています。

つまり、楕円錐とは、頂点から$x-y$平面に下ろした垂線の長さと、その垂線の足からの距離の比が一定である曲面のことです。

楕円錐は$a=b=c$のとき球錐(単に「錐体」)、$a=b≠c$のとき回転楕円体、$a≠b≠c$のときは一般の楕円錐になります。

このように、楕円錐の方程式は比の関係で表され、その形状は$a$, $b$, $c$の値で決まります。

ソースコード

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータ
a = 2  # x軸の半長径
b = 3  # y軸の半長径
c = 4  # z軸の半長径
u = 1  # x軸の平行移動量
v = 1  # y軸の平行移動量
w = 1  # z軸の平行移動量

# メッシュグリッドデータ
x = np.arange(-6, 6, 0.25)
y = np.arange(-6, 6, 0.25)
x, y = np.meshgrid(x, y)

z1 = (c / a) * np.sqrt(a**2 - (x - u)**2 - (y - v)**2) + w
z2 = -(c / a) * np.sqrt(a**2 - (x - u)**2 - (y - v)**2) + w

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z1, alpha=0.5)
ax.plot_surface(x, y, z2, alpha=0.5)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

このコードは、以下の手順で楕円錐のグラフを描画します。

NumPyとMatplotlibをインポートします。
楕円錐のパラメータ(半長径と平行移動量)を設定します。
メッシュグリッドデータ($x$, $y$座標)を生成します。
楕円錐の方程式から、$z$座標を計算します。
3D空間にプロットし、軸ラベルを設定します。
グラフを表示します。

このコードを実行すると、指定された楕円錐の3Dグラフが表示されます。

パラメータを変更すれば、異なる楕円錐のグラフを描くことができます。

結果解説

[実行結果]

表示されるグラフの内容を詳しく説明します。

1. 表示内容

3次元の楕円錐の表面が、2つの半透明な曲面として描画されています。
1つの曲面は楕円錐の上半分を、もう1つの曲面は下半分を表しています。

2. 座標軸・ラベル

3次元の座標軸($X$軸、$Y$軸、$Z$軸)が表示されています。
各軸にはラベル(X、Y、Z)が付いています。

3. 曲面の形状

曲面の形状は、設定したパラメータ$(a, b, c, u, v, w)$によって決まります。
$a$は$X$軸方向の半長径、$b$は$Y$軸方向の半長径、$c$はZ軸方向の半長径です。
$u$、$v$、$w$は平行移動量で、楕円錐の頂点の位置を決めます。

4. 半透明効果

両方の曲面は$alpha=0.5$と設定されているため、半透明になっています。
これにより、2つの曲面が重なった部分の形状がより分かりやすくなります。

5. 曲面の色

デフォルトでは、曲面の色はオレンジ色とブルーの濃淡で表現されています。
色は曲面の高さ($z$座標の値)に応じて変化します。

このグラフは、与えられた楕円錐の方程式を視覚化したものです。

パラメータを変更すれば、異なる形状や位置の楕円錐を表すことができます。

半透明の重ね描画により、立体的な形状がよりわかりやすくなっています。

カドモフ方程式

April 17, 2024

カドモフ方程式

カドモフ方程式（Kadomtsev-Petviashvili equation）は非線形偏微分方程式の一種であり、主に2次元非線形波動の理論において重要な方程式です。

この方程式は水理学やプラズマ物理学、非線形光学、または海洋学などさまざまな分野で応用されます。

一般的なカドモフ方程式は次の形を取ります：

$$
u_t + u_{xxx} + 6uu_x = 0
$$

ここで各項の意味は以下の通りです：

$( u(x, t) ) $は変数$ ( x ) $と$ ( t ) $の関数であり、波の高さや系の状態を表します。
$( u_t ) $は$ ( u ) $の時間微分であり、時間変化を表します。
$( u_{xxx} ) $は$ ( u ) $の$ ( x ) $に関する3階微分であり、空間的な変化を表します。
$( u_x ) $は$ ( u ) $の$ ( x ) $に関する1階微分であり、空間的な変化率を表します。
$( u_{xxx} ) $は$ ( u_x ) $の$ ( x ) $に関する3階微分であり、空間的な変化率の空間微分を表します。
$( u \cdot u_x ) $は非線形項であり、波の非線形相互作用を表します。

カドモフ方程式はソリトン（孤立波）やその他の非線形波動解を持つことで知られています。

この方程式は安定な非線形波動を記述し、系のエネルギー保存則や非線形な波動相互作用を理解するための基本的なモデルとして広く研究されています。

カドモフ方程式の解析的な解は一般に困難ですが、数値シミュレーションや解析的手法を用いて、方程式の振る舞いや特性を理解することが試みられています。

ソースコード

以下に、Pythonを使用してカドモフ方程式を解き、グラフ化する例を示します。

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

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import solve_ivp

次に、カドモフ方程式の右辺を定義します。

def kadomtsev_petviashvili(t, u):
    u1, u2 = u
    du1_dt = u2
    du2_dt = 6 * u1 * (u1 - u2) * (1 - u1)
    return [du1_dt, du2_dt]

解析に使用する初期値と時間範囲を設定します。

1 2	t_span = [0, 10] u0 = [0.5, 0.5] # 初期値 [u1(0), u2(0)]

solve_ivp() を使用してカドモフ方程式を数値的に解きます。

1 2	sol = solve_ivp(kadomtsev_petviashvili, t_span, u0, t_eval=np.linspace(t_span[0], t_span[1], 1000)) u1, u2 = sol.y

最後に、解をグラフ化します。

plt.figure(figsize=(8, 6))
plt.plot(sol.t, u1, label='u1(t)')
plt.plot(sol.t, u2, label='u2(t)')
plt.xlabel('t')
plt.ylabel('u')
plt.title('Solution of the Kadomtsev-Petviashvili Equation')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、solve_ivp() を使用してカドモフ方程式を解き、その解を時刻 t に対してプロットしています。

適切な初期値や時間範囲、グラフの設定を調整して、カドモフ方程式の数値解をグラフ化することができます。

[実行結果]

結果解説

カドモフ方程式の数値解をグラフ化する際、以下の内容が表示されます：

1. 軸のラベル:

横軸（x軸）: 時間 t
縦軸（y軸）: 解 u1(t) および u2(t) の値

2. 解の変化:

グラフ上には時間 t に対する解 u1(t) と u2(t) の値がプロットされます。
u1(t) と u2(t) は時間の関数として表現され、時間が増加するにつれてそれぞれの解の値がどのように変化するかが示されます。

3. 解の意味:

u1(t) はカドモフ方程式の主要な解であり、系の状態や振る舞いを示します。
u2(t) は u1(t) の時間微分に相当し、系の速度や変化率を表します。

4. タイトルと凡例:

グラフには「Solution of the Kadomtsev-Petviashvili Equation」というタイトルが表示されます。
凡例には u1(t) と u2(t) のラベルが含まれ、どちらがどの曲線を表しているかが示されます。

5. グリッド:

グラフには背景にグリッドが表示され、各目盛りが値の参照を補助します。

このグラフを通じて、時間とともにカドモフ方程式の解がどのように振る舞うかを視覚的に理解することができます。

u1(t) と u2(t) の値の変化を見ながら、方程式が記述する現象や振る舞いを把握することができます。

パラボロイドの方程式

April 16, 2024

パラボロイドの方程式

パラボロイドは数学的な曲面であり、放物線を回転させて得られる回転放物面の一種です。
パラボロイドは二次曲面の一種であり、一般的には次のような形式で表されます：

$$
z = \frac{x^2}{a^2} + \frac{y^2}{b^2}
$$

ここで、$( x )$、$( y )$、$( z ) $は空間内の座標を表し、$( a ) $と$ ( b ) $は定数であり、それぞれ$ x軸$方向と$ y軸$方向の放物線の広がりを制御します。

この方程式における放物面は、$( x ) $と$ ( y ) $の値によって$ ( z ) $座標が決まり、放物線の頂点が原点に位置します。
パラボロイドは、放物線が$ x軸$と$ y軸$の両方向に広がっている形状を持ち、$ ( z ) $の値は放物線の広がりに応じて変化します。

Pythonでパラボロイドをグラフ化する際には、このような方程式を用いて$ ( x )$、$( y ) $の範囲を指定し、それに対応する$ ( z ) $の値を計算します。

得られた$ ( x )$、$( y )$、$( z ) $の値を3次元グラフとしてプロットすることで、パラボロイドの形状を視覚化することができます。

ソースコード

以下のPythonコードは、パラボロイドの方程式を解いてグラフ化するものです。

この例では、NumPyとMatplotlibライブラリを使用しています。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータ設定
a = 1.0
b = 1.0
c = 1.0

# パラメータの範囲設定
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

# パラボロイドの方程式をパラメータ表示で定義
X = a * U
Y = b * V
Z = c * (U**2 + V**2)

# 3Dプロット
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')

# パラボロイドをプロット
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k')

# 軸ラベルの設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# グラフのタイトル
plt.title('Paraboloid')

# グラフを表示
plt.show()

このコードでは、a、b、cはパラメータとしてパラボロイドの形状を制御します。

UとVはパラメータの範囲を定義し、X、Y、Zはパラボロイドの方程式に基づいて表面の座標を計算します。

plot_surface関数を使用してパラボロイドを3Dプロットし、軸ラベルやタイトルを設定してグラフを表示します。

[実行結果]

ソースコード解説

このPythonコードは、パラボロイド曲面を生成し、3次元プロットで視覚化するためのものです。

以下にそれぞれの部分を詳しく説明します。

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

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

numpyは数値計算を行うためのライブラリであり、配列や数学関数を効率的に扱います。
matplotlib.pyplotはグラフ描画ライブラリで、グラフのプロットやカスタマイズに使用します。
Axes3Dは3次元プロット用のクラスを提供します。

2. パラメータ設定:

1
2
3

a = 1.0
b = 1.0
c = 1.0

パラボロイドの方程式における定数$ (a)$、$(b)$、$(c) $を設定します。

3. パラメータの範囲設定:

1
2
3

u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

パラメータ$ (u)$、$(v) $の範囲を設定し、それらの値のメッシュグリッドを作成します。

4. パラボロイドの方程式を定義:

1
2
3

X = a * U
Y = b * V
Z = c * (U**2 + V**2)

パラボロイドの方程式$ (z = c(u^2 + v^2)) $をパラメータ$ (U)$、$(V) $を使って定義します。

5. 3Dプロット:

1
2
3

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k')

plot_surfaceメソッドを使用して、パラボロイド曲面を3次元プロットします。
cmap='viridis'はカラーマップを設定し、’viridis’は色のマッピングを示します。
edgecolor='k'は曲面の枠線の色を黒に設定します。

6. 軸ラベルとタイトルの設定:

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('Paraboloid')

軸ラベルとグラフのタイトルを設定します。

7. グラフの表示:

1	plt.show()

最後に、プロットされたグラフを表示します。

このコードを実行すると、$(u)$、$(v) $の範囲内で定義されたパラボロイド曲面が3次元空間に表示されます。

結果解説

[実行結果]

上記のPythonコードによって生成されるグラフは、パラボロイド（放物体）の3次元プロットです。

このグラフは以下のような特徴を持っています：

1. 形状:

パラボロイドは放物線を回転軸周りに回転させて生成される曲面です。
この例では、パラメータ a、b、c を調整することでパラボロイドの形状を変えることができます。

2. 色:

cmap='viridis' を指定することで、グラフの色をViridisカラーマップに設定しています。
これにより、高さ（Z軸方向）に応じて色が異なります。
Viridisは高さが低いときは青色から、高いときは黄色や緑色に変化するカラーマップです。

3. 座標軸:

グラフには$X軸$、$Y軸$、$Z軸$があります。
各軸はそれぞれ X、Y、Z に対応し、パラボロイドの方程式の座標値を示しています。

4. タイトルとラベル:

グラフには適切なタイトル ('Paraboloid') が付けられています。
$X軸$、$Y軸$、$Z軸$にはそれぞれ 'X'、'Y'、'Z' というラベルが付けられており、各軸の意味を明示しています。

このグラフはパラメータを調整することでパラボロイドの形状を変えることができ、3次元空間で放物体の曲面を視覚化するのに役立ちます。

リャプノフ方程式

April 15, 2024

リャプノフ方程式は、非線形の偏微分方程式の一種で、以下のように表されます。

$$
∂u/∂t = -α u ∂u/∂x + β ∂^2u/∂x^2 + γ ∂(u^2)/∂x
$$

ここで、$u$ は未知関数、$t$ は時間、$x$ は空間座標、$α$、$β$、$γ$ はパラメータです。

この方程式は、様々な分野で現れる非線形現象をモデル化するのに役立ちます。
特に流体力学や プラズマ物理学、光学などの分野で重要な役割を果たします。

第1項　$ -α u ∂u/∂x $は移流項と呼ばれ、流れによる物理量の移動を表します。
第2項　$ β ∂^2u/∂x^2 $は拡散項で、物理量の空間的な広がりを表します。
第3項　$ γ ∂(u^2)/∂x $は非線形項で、物理量同士の相互作用を表します。

$α$、$β$、$γ $のパラメータ値を変えることで、様々な非線形現象をモデル化できます。
例えば、$α$が正の値の場合は前方への伝播、負の値の場合は後方への伝播を表します。

数値的に解くと、リャプノフ方程式の解は初期条件に強く依存しながらも、独特の時空間パターンを形成することがわかります。
このパターンは、対象とする非線形現象の本質的な構造を反映していると考えられています。

従って、リャプノフ方程式を解くことで、複雑な非線形現象の振る舞いを理解し、予測することが可能になります。
特に乱流などの複雑な流れの研究に貢献してきました。

ソースコード

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
alpha = 1.0
beta = 0.5
gamma = -1.0

# 時間と空間の範囲
tmin, tmax, dt = 0, 20, 0.01
xmin, xmax, dx = -10, 10, 0.1

# 初期条件
u0 = np.exp(-(x**2)/2)

# 時間と空間の格子点
t = np.arange(tmin, tmax, dt)
x = np.arange(xmin, xmax, dx)
nt, nx = len(t), len(x)

# 解の格納配列
u = np.zeros((nt, nx))

# 初期条件の設定
u[0] = u0

# Lyapunov 方程式の数値解法
for i in range(nt-1):
    for j in range(1, nx-1):
        u[i+1, j] = u[i, j] - alpha*dt*u[i, j]*(u[i, j+1] - u[i, j-1])/(2*dx) \
                    + beta*dt*(u[i, j+1] - 2*u[i, j] + u[i, j-1])/(dx**2) \
                    + gamma*dt*(u[i, j+1]**2 - u[i, j-1]**2)/(2*dx)

# グラフ描画
plt.figure(figsize=(10, 6))
plt.contourf(x, t, u, levels=np.linspace(np.min(u), np.max(u), 21), cmap='viridis')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('Lyapunov Equation Solution')
plt.show()

上記のコードは、Lyapunov方程式の数値解を計算し、その結果をグラフ化するPythonプログラムです。

まず必要なライブラリ(numpy, matplotlib)をインポートします。
次にLyapunov方程式のパラメータ(alpha, beta, gamma)を設定します。
時間と空間の範囲を設定します。
初期条件$u0$を定義します。
時間と空間の格子点を作成します。
解の格納配列$u$を初期化します。
初期条件を設定します。
Lyapunov方程式の数値解法の本体部分です。
時間ループと$x$方向の空間ループを使って、格子点上の解を計算していきます。
最後に、matplotlibを使って計算結果をコンター図に描画します。

このプログラムを実行すると、Lyapunov方程式の数値解がコンター図として表示されます。
初期条件の影響を受けながら、時間の経過とともに解の様子が変化している様子が確認できます。

[実行結果]

ソースコード解説

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

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

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

Numpyライブラリと、Matplotlibライブラリをインポートしています。
Numpyは数値計算に関する処理を行うのに使用します。
Matplotlibは数値データの可視化を行うのに使用します。

2. パラメータの設定

1
2
3

alpha = 1.0
beta = 0.5 
gamma = -1.0

リャプノフ方程式に現れるパラメータ$α$、$β$、$γ$の値を設定しています。

3. 時間と空間の範囲設定

1 2	tmin, tmax, dt = 0, 20, 0.01 xmin, xmax, dx = -10, 10, 0.1

時間$t$と空間$x$の計算範囲を設定しています。
$t$は$0$から$20$まで、$dt=0.01$間隔で計算します。
$x$は$-10$から$10$まで、$dx=0.1$間隔で計算します。

4. 初期条件の設定

1	u0 = np.exp(-(x**2)/2)

初期条件$u0$をガウス分布$exp(-(x**2)/2)$で与えています。

5. 時間と空間の格子点作成

1
2
3

t = np.arange(tmin, tmax, dt)
x = np.arange(xmin, xmax, dx)
nt, nx = len(t), len(x)

時間$t$と空間$x$の格子点を作成しています。
格子点の個数$nt$、$nx$も求めています。

6. 解の格納配列の初期化

1	u = np.zeros((nt, nx))

計算する解$u$を格納する配列を、$nt×nx$のゼロ行列として初期化しています。

7. 初期条件の代入

u[0] = u0

解の格納配列の最初の時間ステップ$(t=0)$に初期条件$u0$を代入しています。

8. リャプノフ方程式の数値解法

for i in range(nt-1):
    for j in range(1, nx-1):
        u[i+1, j] = u[i, j] - alpha*dt*u[i, j]*(u[i, j+1] - u[i, j-1])/(2*dx) \
                    + beta*dt*(u[i, j+1] - 2*u[i, j] + u[i, j-1])/(dx**2) \
                    + gamma*dt*(u[i, j+1]**2 - u[i, j-1]**2)/(2*dx)

時間と空間の2重ループを使って、リャプノフ方程式の数値解法で次の時間ステップの解を計算しています。
移流項、拡散項、非線形項を陽的に計算しています。

9. グラフ描画

plt.figure(figsize=(10, 6))
plt.contourf(x, t, u, levels=np.linspace(np.min(u), np.max(u), 21), cmap='viridis')
plt.colorbar()
plt.xlabel('x')
plt.ylabel('t')
plt.title('Lyapunov Equation Solution')
plt.show()

最後に、Matplotlibを使って計算結果をコンター図で描画しています。
縦軸が時間$t$、横軸が空間$x$で、コンターの色分けが解$u$の値を表しています。
図のタイトル、軸ラベル、カラーバーも設定しています。

このように、このコードは初期条件から出発し、リャプノフ方程式の数値解法で次の時間ステップの解を計算し、最終的にその結果をコンター図で可視化するプログラムになっています。

結果解説

[実行結果]

このグラフは、時間$t$と空間$x$の2次元のコンター図となっています。
縦軸が$t$で時間の経過を表し、横軸が$x$で空間的な広がりを表しています。

コンター図の色分けは、方程式の解$u$の値を表しています。
コンターの色が赤に近づくほど解の値が大きく、青に近づくほど解の値が小さいことを意味しています。

時間$t=0$における初期条件は、$exp(-(x**2)/2)$によるガウス分布の形状になっていて、グラフの一番下の部分でそれが確認できます。
この初期条件から出発し、時間が経過するにつれて解の様子が変化していきます。

具体的には、初期条件のガウス分布が徐々に変形していき、時間ととともに特徴的なパターンが現れてきます。
このパターンの具体的な形状は、与えられた方程式のタイプや係数値によって異なります。

一般に、このようなコンター図を描くことで、与えられた非線形偏微分方程式の数値解がどのように時間と空間に対して進行するかを視覚的に捉えることができます。
初期条件の影響を受けながらも、独自のダイナミクスに従って解が変化する様子が確認できるのが特徴です。

このように、コンター図から方程式の解の時空間的な振る舞いを詳細に観察でき、対象となる現象の本質的な構造を理解する上で重要な情報を得ることができます。

ヘストネス方程式

April 14, 2024

ヘストネス方程式

ヘストネス方程式は、David Hestenes が1986年に提案した代数的な方程式です。

この方程式は、3次元空間における曲面の形状を記述するものです。

ヘストネス方程式は次のように表されます:

$$
(x^2 + y^2 + z^2 - 1)^3 - x^2y^3 - y^2z^3 - z^2x^3 = 0
$$

この方程式は、非線形の代数方程式で、左辺が$0$になる点の集合が曲面を表しています。

この曲面には以下のような特徴があります:

中心が原点$(0, 0, 0)$にあり、単位球面$(x^2 + y^2 + z^2 = 1)$に内接する曲面です。
曲面は滑らかではなく、いくつかの尖った部分があります。
曲面は$8$回対称性を持っています。
つまり、座標軸周りに$45$度回転させても同じ形状になります。
曲面は自己交差しており、内部に空洞部分があります。

ヘストネス方程式は、幾何学的に興味深い性質を持つ代数方程式として知られており、コンピューター graphicsや数学の視覚化などで利用されています。

この方程式から得られる曲面は、複雑で美しい形状を持つため、芸術的な表現の題材にもなっています。

ソースコード

import numpy as np
import matplotlib.pyplot as plt

def hestenes_equation(x, y, z):
    return (x**2 + y**2 + z**2 - 1)**3 - x**2 * y**3 - y**2 * z**3 - z**2 * x**3

# Meshgridを作成
x = np.linspace(-1.5, 1.5, 50)
y = np.linspace(-1.5, 1.5, 50)
X, Y = np.meshgrid(x, y)

# 方程式を計算
Z = hestenes_equation(X, Y, 0)  # z=0の平面で計算

# 3D プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Hestenes Equation')

plt.show()

上記のコードでは、まずHestenes方程式を定義する関数 hestenes_equation(x, y, z) を作成しています。

この関数は、与えられた x, y, z の値に対するHestenes方程式の値を返します。

次に、numpy を使って x と y の値の範囲を設定し、それらの値の組み合わせからメッシュグリッドを作成しています。

そして、z=0 の平面上でHestenes方程式を計算し、その結果を Z に格納しています。

最後に、matplotlib の plot_surface 関数を使って、計算された Z の値を3次元プロットしています。

実行すると、Hestenes方程式の立体的な曲面がプロットされます。

[実行結果]

ソースコード解説

ソースコードを説明します。

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

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

最初に、NumPyとMatplotlibライブラリをインポートしています。
NumPyは数値計算を、Matplotlibはデータの可視化を行うためのPythonライブラリです。

2. Hestenes方程式の定義

1 2	def hestenes_equation(x, y, z): return (x2 + y2 + z2 - 1)3 - x*2 y3 - y2 * z3 - z2 * x**3

Hestenes方程式を計算する関数hestenes_equationを定義しています。

この関数は、x、y、zの値を引数として受け取り、それらの値に対するHestenes方程式の値を返します。

3. メッシュグリッドの作成

1
2
3

x = np.linspace(-1.5, 1.5, 50)
y = np.linspace(-1.5, 1.5, 50)
X, Y = np.meshgrid(x, y)

NumPyのlinspace関数を使って、xとyの値の範囲を設定しています。
ここでは、$-1.5$から$1.5$までの$50$個の等間隔の値を生成しています。
次に、meshgrid関数を使って、xとyの値の組み合わせからメッシュグリッドを作成しています。

4. Hestenes方程式の計算

1	Z = hestenes_equation(X, Y, 0) # z=0の平面で計算

先程定義したhestenes_equation関数を使って、メッシュグリッド上の各点における方程式の値を計算しています。
ここでは、z=0の平面上での値を計算しています。
計算結果はZに格納されます。

5. 3次元プロットの設定

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Hestenes Equation')

Matplotlibを使って、3次元プロットを設定しています。

まず、新しい図figを作成し、その中に3次元の座標軸axを追加しています。

次に、plot_surface関数を使って、計算されたZの値を3次元プロットしています。
rstrideとcstrideは、プロットする点の間隔を指定するパラメータです。
cmap='viridis'は、曲面の色を指定するためのカラーマップです。

最後に、軸のラベルとタイトルを設定しています。

6. プロットの表示

1	plt.show()

最後に、plt.show()を呼び出すことで、作成したプロットを表示します。

結果解説

[実行結果]

グラフに表示されている内容について説明します。

このグラフは3次元プロットで、Hestenes方程式から得られる曲面を可視化しています。

曲面の形状は以下のような特徴を持っています:

中心が原点$(0, 0, 0)$にあり、単位球面に内接しています。
曲面の外形は球形に近い形をしています。
曲面は滑らかではなく、いくつかの鋭い尖った部分があります。
特に、座標軸の正負方向に$4$つの主要な尖った部分があります。
曲面は$8$回対称性を持っており、座標軸周りに$45$度回転させても同じ形状になります。
曲面は自己交差しており、内部に空洞部分があります。
曲面の内側には複雑な入り組んだ構造があり、一見すると次元が高い幾何学的な形状のようにも見えます。
曲面の色は、viridisカラーマップを使って高さに応じて色分けされています。
原点付近は青紫色、高い部分は黄色がかった緑色になっています。
グラフの背景は白で、曲面の形状が浮き立つようになっています。
$X$軸、$Y$軸、$Z$軸のラベルと、タイトル「Hestenes Equation」が付けられています。

このようにHestenes方程式から得られる曲面は、非常に複雑で美しい形状を持っており、その幾何学的な特徴が3次元プロットによってよく可視化されています。

双曲面の方程式

April 13, 2024

双曲面の方程式

双曲面の方程式をPythonで解いてグラフ化する方法を示します。

ここでは、主軸が$ (z) $軸に平行な双曲面の方程式を例に取ります。

双曲面の方程式:
$$
[ \frac{x^2}{a^2} + \frac{y^2}{b^2} - \frac{z^2}{c^2} = 1 ]
$$

この方程式をグラフ化するために、NumPyとMatplotlibを使用します。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 双曲面のパラメータ
a = 2.0  # x方向の半軸の長さ
b = 3.0  # y方向の半軸の長さ
c = 4.0  # z方向の半軸の長さ

# パラメータの範囲設定
u = np.linspace(-2, 2, 100)
v = np.linspace(-2, 2, 100)
U, V = np.meshgrid(u, v)

# 双曲面の方程式をパラメータ表示で定義
X = a * np.cosh(U) * np.cos(V)
Y = b * np.cosh(U) * np.sin(V)
Z = c * np.sinh(U)

# 3Dプロット
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')

# 双曲面をプロット
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k')

# 軸ラベルの設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# グラフのタイトル
plt.title('Hyperboloid')

# グラフを表示
plt.show()

このコードでは、双曲面のパラメータ$ (a, b, c) $を指定しています。

そして、双曲線の方程式をパラメータ$ (U, V) $の関数として定義し、それを使用して$ (X, Y, Z) $座標を計算しています。

最後に、計算した座標を用いてMatplotlibのplot_surface関数で双曲面をプロットしています。

このコードを実行すると、指定した双曲面が3次元グラフとして表示されます。

[実行結果]

ソースコード解説

ソースコードの詳細な説明は以下の通りです。

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

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

numpy: 数値計算ライブラリ。配列操作や数学関数を提供。
matplotlib.pyplot: グラフ描画ライブラリ。
mpl_toolkits.mplot3d: 3次元プロット用のツールキット。

2. 双曲面のパラメータ設定

1
2
3

a = 2.0  # x方向の半軸の長さ
b = 3.0  # y方向の半軸の長さ
c = 4.0  # z方向の半軸の長さ

双曲面の各軸方向の半軸の長さを設定します。

3. パラメータの範囲設定とメッシュグリッド作成

1
2
3

u = np.linspace(-2, 2, 100)
v = np.linspace(-2, 2, 100)
U, V = np.meshgrid(u, v)

np.linspace: 指定された範囲内で等間隔の数値を生成。
meshgrid: 2つの1次元配列から2次元グリッドを生成。

4. 双曲面の方程式をパラメータ表示で定義

1
2
3

X = a * np.cosh(U) * np.cos(V)
Y = b * np.cosh(U) * np.sin(V)
Z = c * np.sinh(U)

双曲面の方程式をパラメータ$ (U) $と$ (V) $の関数として定義します。
np.cosh: 双曲線余弦関数。
np.sinh: 双曲線正弦関数。

5. 3Dプロットの準備

1 2	fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection='3d')

plt.figure: 新しい図を作成します。
fig.add_subplot: サブプロットを追加します。
ここでは3次元の図を作成しています。

6. 双曲面のプロット

1	ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k')

ax.plot_surface: 双曲面を3Dプロットします。
cmap='viridis': カラーマップをViridisに設定。
edgecolor='k': エッジの色を黒色に設定。

7. 軸ラベルとタイトルの設定

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.title('Hyperboloid')

軸ラベルとグラフのタイトルを設定します。

8. グラフの表示

1	plt.show()

作成した3Dグラフを表示します。

このソースコードを実行すると、定義された双曲面の方程式に基づいて双曲面が描画され、その3次元形状が視覚化されます。

結果解説

[実行結果]

以下はグラフに表示される内容の詳細です。

1. 双曲面の方程式:

双曲面の方程式は次のように表されます。
$$
\frac{x^2}{a^2} + \frac{y^2}{b^2} - \frac{z^2}{c^2} = 1
$$
この方程式は、$ (a, b, c) $を用いて定義される双曲面の形状を示しています。
この例では主軸が$ (z) $軸に平行な双曲面を扱います。

2. パラメータ設定:

(u, v) はそれぞれパラメータ空間の値を表します。
np.linspace(-2, 2, 100) を使用して、それぞれのパラメータの値を生成します。
U, V はこれらのパラメータのメッシュグリッドを作成します。

3. 座標計算:

双曲面の方程式をパラメータ表示で定義します。
$ X = a \cosh(U) \cos(V) $
$ Y = b \cosh(U) \sin(V) $
$ Z = c \sinh(U) $
ここで、$ (\cosh) $は双曲線余弦関数、$ (\sinh) $は双曲線正弦関数です。

4. 3Dプロット:

plot_surface 関数を使用して、双曲面の$ (X, Y, Z) $座標を3次元空間にプロットします。
cmap='viridis' はカラーマップを指定し、edgecolor='k' はエッジの色を黒に設定します。

5. グラフの装飾:

軸ラベルやタイトルを設定して、グラフをわかりやすくします。
$X$軸は ‘X’、$Y$軸は ‘Y’、$Z$軸は ‘Z’ としてラベル付けされます。
グラフのタイトルは ‘Hyperboloid’ となります。

6. グラフの表示:

最後に plt.show() を使用して、作成した3Dグラフを表示します。

このコードを実行すると、$ (a, b, c) $の値に基づいて定義された双曲面が3次元プロットで表示されます。

双曲面の形状や曲線が視覚化され、指定したパラメータによって変化する様子を確認できます。

マクスウェルの方程式

April 12, 2024

マクスウェルの方程式

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

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

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

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

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

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

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

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

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

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

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

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

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

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. シミュレーションパラメータの設定

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. 時間ループ

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. 結果の可視化

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方程式に従って光速度で伝わり、その過程で磁場も生成・伝播する様子がわかります。

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