ヘストネス方程式

4月 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次元プロットによってよく可視化されています。

双曲面の方程式

4月 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次元プロットで表示されます。

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

マクスウェルの方程式

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

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

偏心楕円体

4月 11, 2024

偏心楕円体

偏心楕円体の方程式を考えます：

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

ここで、$ ( a, b, c ) $はそれぞれ楕円体の軸の長さを表します。

この方程式は、中心が原点にあり、軸が $x$, $y$, $z$ 軸に沿った偏心楕円体を表します。

これをPythonで解いて3Dグラフを作成する手順を示します。

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

# 楕円体のパラメータ設定
a = 1.0  # x軸方向の長さ
b = 2.0  # y軸方向の長さ
c = 3.0  # z軸方向の長さ

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

# 楕円体の方程式をパラメータ表示で定義
x = a * np.sqrt(1 - U**2) * np.cos(V)
y = b * np.sqrt(1 - U**2) * np.sin(V)
z = c * 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('Eccentric Ellipsoid')

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

このコードでは、楕円体の方程式をパラメータ表示で定義し、その表面を plot_surface 関数を用いて3Dプロットしています。

ここでは楕円体の軸の長さをそれぞれ$ (a = 1)$、$(b = 2)$、$(c = 3) $と設定しましたが、これらの値を変更することで異なる形状の偏心楕円体を描画することができます。

[実行結果]

このようにして、数学的な方程式をPythonで解いて3Dグラフ化することができます。

ソースコード解説

このソースコードは、Pythonを使用して楕円体を3Dで描画するためのプログラムです。

以下に、コードの詳細な説明を示します。

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

import numpy as np: 数値計算や配列操作のためのライブラリnumpyをnpとしてインポートします。
import matplotlib.pyplot as plt: グラフ描画のためのライブラリmatplotlibからpyplotモジュールをpltとしてインポートします。
from mpl_toolkits.mplot3d import Axes3D: 3Dグラフ描画のためのAxes3Dクラスをインポートします。

2. 楕円体のパラメータ設定:

a = 1.0: $x$軸方向の長さを表すパラメータです。
b = 2.0: $y$軸方向の長さを表すパラメータです。
c = 3.0: $z$軸方向の長さを表すパラメータです。

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

u = np.linspace(-1, 1, 100): パラメータuの値を$-1$から$1$まで$100$個の等間隔の値で生成します。
v = np.linspace(0, 2*np.pi, 100): パラメータ$v$の値を$0$から$2π$まで$100$個の等間隔の値で生成します。
U, V = np.meshgrid(u, v): $u$と$v$の値の組み合わせに基づいて$U$と$V$のメッシュグリッドを生成します。

4. 楕円体の方程式をパラメータ表示で定義:

x = a * np.sqrt(1 - U**2) * np.cos(V): $x$座標を楕円体のパラメータに基づいて計算します。
y = b * np.sqrt(1 - U**2) * np.sin(V): $y$座標を楕円体のパラメータに基づいて計算します。
z = c * U: $z$座標を楕円体のパラメータに基づいて計算します。

5. 3Dプロットの準備:

fig = plt.figure(figsize=(8, 8)): グラフの描画領域（Figureオブジェクト）を作成し、サイズを指定します。
ax = fig.add_subplot(111, projection='3d'): 3D描画用のサブプロット（Axes3Dオブジェクト）を作成します。

6. 楕円体をプロット:

ax.plot_surface(x, y, z, cmap='viridis', edgecolor='k'): 生成した楕円体の座標データを用いて、3D曲面をプロットします。

カラーマップは’viridis’を使用し、エッジカラーは黒色で指定します。

7. 軸ラベルの設定:

ax.set_xlabel('X'): $x$軸のラベルを設定します。
ax.set_ylabel('Y'): $y$軸のラベルを設定します。
ax.set_zlabel('Z'): $z$軸のラベルを設定します。

8. グラフのタイトル設定:

plt.title('Eccentric Ellipsoid'): グラフにタイトルを設定します。

9. グラフの表示:

plt.show(): これまでに設定したグラフを表示します。

これにより、楕円体のパラメータを調整して3Dプロットを行うことができます。

楕円体の長さや形状を変えるには、パラメータa、b、cを変更してみてください。

結果解説

[実行結果]

このプログラムは、楕円体を3Dグラフとして描画するものです。

楕円体の方程式をパラメータ表示で定義し、それを用いて楕円体の表面をプロットしています。

以下に、グラフに表示される内容を詳しく説明します。

1. 楕円体の方程式:

楕円体のパラメータ表示において、楕円体の方程式は次のように定義されています。

x = a * np.sqrt(1 - U**2) * np.cos(V)
y = b * np.sqrt(1 - U**2) * np.sin(V)
z = c * U
ここで、UとVはパラメータであり、Uは楕円体の軸方向に沿った位置を示し、Vは楕円体の周囲を回る位置を示します。
a、b、cはそれぞれ楕円体の長軸、中軸、短軸の長さを表します。

2. 楕円体の形状:

楕円体の形状はパラメータa、b、cによって決まります。
aは$x$軸方向の長さ、bは$y$軸方向の長さ、cは$z$軸方向の長さを表します。
これらの値を変更すると、楕円体の形状が変わります。

3. 3Dプロット:

ax.plot_surface(x, y, z, cmap='viridis', edgecolor='k')により、楕円体の表面が3Dプロットされます。
x、y、zは楕円体の表面の座標を示し、cmap='viridis'はカラーマップを設定し、edgecolor='k'はエッジの色を黒に設定しています。

4. 軸ラベル:

ax.set_xlabel('X')、ax.set_ylabel('Y')、ax.set_zlabel('Z')によって、各軸のラベルが設定されます。
これにより、各軸が何を表しているかが分かりやすくなります。

5. グラフのタイトル:

plt.title('Eccentric Ellipsoid')によって、グラフのタイトルが設定されます。
ここでは「Eccentric Ellipsoid」というタイトルが付けられていますが、このタイトルは必要に応じて変更できます。

6. グラフの表示:

plt.show()によって、設定されたグラフが表示されます。
これにより、楕円体の形状や表面が3D空間で視覚化されます。

このプログラムを実行することで、楕円体の形状やパラメータを自由に変更して、異なる楕円体を描画することができます。

楕円体の形状や向きを調整することで、さまざまな3Dグラフを生成できるため、楕円体のパラメータを変更して試してみてください。

奇妙な方程式 3Dグラフ

4月 10, 2024

奇妙な方程式 3Dグラフ

奇妙な方程式を立て、Pythonでそれをプロットしてみます。

3次元のグラフを描くためにMatplotlib の plt.plot_surface() 関数を使用します。

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

# 変な方程式を定義
def weird_equation(x, y):
    return np.sin(x**2 + y**3) * np.cos(x * y) + x * np.exp(-y**2)

# グリッドデータを生成
x = np.linspace(-4, 4, 50)
y = np.linspace(-4, 4, 50)
X, Y = np.meshgrid(x, y)
Z = weird_equation(X, Y)

# 3Dプロットを作成
fig = plt.figure(figsize=(10, 8))
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.view_init(elev=30, azim=-45)
plt.show()

このコードでは、まず奇妙な方程式 weird_equation(x, y) を定義しています。
この方程式は sin、cos、exp などの関数を使って作られた複雑な式です。

次に、np.linspace と np.meshgrid を使って、$x$ と $y$ の値からなるグリッドデータを生成します。
そして、この$ x $と$ y $の値を weird_equation に代入して、$z $の値を計算しています。

最後に、Matplotlib の plt.plot_surface 関数を使って、$X$、$Y$、$Z$ のデータから3次元のサーフェスプロットを描いています。

rstride と cstride はグリッドデータの間引き率を指定するパラメータで、値を大きくすると描画が高速になります。

cmap ではカラーマップを指定しています。
view_init で、プロットの視点を調整しています。

このコードを実行すると、以下のような3Dプロットが表示されるはずです。

[実行結果]

グラフは奇妙な形状をしていますが、これは方程式自体が複雑で変則的なためです。

好みのデザインに合わせて、方程式や表示パラメータを調整することができます。

結果解説

[実行結果]

このグラフを詳しく説明します。

まず、プロットの全体的な形状は、不規則な波打つ曲面になっています。

これは、定義した「奇妙な方程式」の形状を反映したものです。
方程式の中には三角関数の sin、cos、そして指数関数 exp が含まれているため、滑らかではあるものの複雑な形状になっています。

次に、プロットの色分けに注目すると、青系の色から緑、黄色、オレンジ、そして赤へと徐々に変化していることがわかります。

これは、プロットの高さ($z$座標の値)を色で表現しているためです。
青は低い値、赤は高い値を表しています。
具体的には、この色分けは cmap='viridis' というカラーマップを使って指定されています。

プロットの底面($x-y$平面)は、白とグレーの格子状のパターンになっています。
これは、$x軸$と$y軸$の目盛り線を表しています。
目盛りの範囲は、$x$、$y$ともに$-4$から$4$までの値となっています。

さらに、3次元の効果を高めるため、プロットの視点は view_init(elev=30, azim=-45) により、上から$30$度、左から$45$度の角度から見る設定になっています。
この視点を変更すると、形状の立体感が変わります。

最後に、それぞれの軸にはラベル$(X、Y、Z)$が付けられており、データの座標を視覚的に把握しやすくなっています。

このように、この3Dプロットには、方程式から得られる曲面の形状、色分けによる高さの表現、目盛りと軸ラベル、そして視点の設定など、様々な情報が盛り込まれています。

適切にパラメータを調整することで、データの可視化がより効果的に行えます。

ディリクレ問題

4月 9, 2024

ディリクレ問題

ディリクレ問題は、偏微分方程式の境界値問題の一種です。
この問題では、以下の3つの要素が与えられます。

1. 偏微分方程式

通常は、ラプラス方程式やポアソン方程式などの楕円型偏微分方程式が扱われます。

2. 領域($Ω$)

問題が定義される有界な領域です。
例えば、矩形、円、三角形などの2次元または3次元の領域です。

3. 境界条件

領域の境界 ($∂Ω$)上で、未知関数の値が与えられた関数$(f)$に等しくなるという条件です。

ディリクレ問題では、与えられた偏微分方程式と境界条件のもとで、領域内の未知関数 $u(x,y)$の値を求める必要があります。
つまり、領域内部の解を求めることが目的です。

この問題は、数学的にはよく研究されており、様々な解法が知られています。
代表的な解法としては、離散化による数値解法、有限要素法、分離変数法、グリーン関数法などがあります。

ディリクレ問題の応用例としては、電位問題、伝導問題、膜の振動問題、量子力学などがあげられます。
工学や物理学の分野で重要な役割を果たしています。

ソースコード

ディリクレ問題をPythonで解いて、グラフ化します。

ディリクレ問題は、偏微分方程式の境界値問題の一つです。
2次元の場合、方程式は以下のように記述されます。

$$
∂^2u/∂x^2 + ∂^2u/∂y^2 = 0 (on Ω)
$$
$$
u = f(x, y) (on ∂Ω)
$$

$u(x, y) $は未知関数
$Ω $は領域
$∂Ω $は$Ω$の境界
$f(x, y) $は与えられた境界条件

Pythonコードは以下の通りです。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad

# 領域の定義
x_min, x_max = 0, 1
y_min, y_max = 0, 1

# 境界条件
def boundary(x, y):
    return np.sin(np.pi * x) * np.sinh(np.pi * y)

# ディリクレ問題の解
def dirichlet_solution(x, y):
    return np.sin(np.pi * x) * np.sinh(np.pi * y) / np.sinh(np.pi)

# 数値積分
def integrate(f, x_min, x_max, y_min, y_max):
    return dblquad(f, x_min, x_max, lambda x: y_min, lambda x: y_max)[0]

# グリッド点での解の値を計算
x = np.linspace(x_min, x_max, 50)
y = np.linspace(y_min, y_max, 50)
X, Y = np.meshgrid(x, y)
Z = dirichlet_solution(X, Y)

# 3D プロット
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U(x, y)')
ax.set_title('Dirichlet Problem Solution')
plt.show()

このコードでは、最初に領域とその境界条件を定義しています。
境界条件は boundary 関数で与えられています。

次に、ディリクレ問題の解析的な解 dirichlet_solution を定義しています。
この関数は、与えられた境界条件のもとで、ディリクレ問題の解を返します。

integrate 関数は、数値積分を行うための関数です。
ここでは使用していませんが、必要に応じて利用できます。

次に、グリッド点での解の値を計算しています。
np.meshgrid を使って、$X$, $Y$, $Z$ の値を計算しています。

最後に、matplotlib の plot_surface を使って、3D プロットを作成しています。

実行すると、以下のような3Dグラフが表示されます。

[実行結果]

このグラフは、与えられた境界条件のもとでのディリクレ問題の解を視覚化したものです。
境界条件に応じて、解の形状が変化します。

ソースコード解説

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

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

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad

このセクションでは、数値計算や可視化に使用するライブラリ(NumPy, Matplotlib, SciPy)をインポートしています。

2. 問題設定

# 領域の定義
x_min, x_max = 0, 1
y_min, y_max = 0, 1

# 境界条件
def boundary(x, y):
    return np.sin(np.pi * x) * np.sinh(np.pi * y)

ここでは、ディリクレ問題の領域$([0, 1] x [0, 1])$と、境界上での条件を定義しています。
boundary関数は、境界上での未知関数の値を返します。

3. 解析解の定義

1
2
3

# ディリクレ問題の解
def dirichlet_solution(x, y):
    return np.sin(np.pi * x) * np.sinh(np.pi * y) / np.sinh(np.pi)

この関数は、与えられた境界条件のもとでのディリクレ問題の解析解を返します。
つまり、この解が領域内部で満たすべき関数形です。

4. 数値積分関数の定義

1
2
3

# 数値積分
def integrate(f, x_min, x_max, y_min, y_max):
    return dblquad(f, x_min, x_max, lambda x: y_min, lambda x: y_max)[0]

この関数は、2重積分を数値的に計算するためのものです。
ここでは使用していませんが、必要に応じて利用できます。

5. グリッド点での解の値の計算

# グリッド点での解の値を計算
x = np.linspace(x_min, x_max, 50)
y = np.linspace(y_min, y_max, 50)
X, Y = np.meshgrid(x, y)
Z = dirichlet_solution(X, Y)

このセクションでは、領域内のグリッド点での解の値を計算しています。
np.linspaceを使って$x、y$の座標値を生成し、np.meshgridでそれらの組み合わせを作成しています。
最後に、dirichlet_solution関数を使って各グリッド点での解の値$Z$を計算しています。

6. 3D可視化

# 3D プロット
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U(x, y)')
ax.set_title('Dirichlet Problem Solution')
plt.show()

最後のセクションでは、計算した解の値$Z$を3D曲面プロットで可視化しています。
matplotlib.pyplot.figureでフィギュアを作成し、add_subplotでプロット領域を追加しています。
plot_surfaceを使って$X$、$Y$、$Z$の値から曲面を描画しています。
最後に、軸ラベルとタイトルを設定し、show()でグラフを表示しています。

このコードは、与えられた境界条件のもとでのディリクレ問題の解析解を計算し、3D曲面グラフで可視化しています。
このようにして、解の振る舞いを視覚的に確認することができます。

結果解説

[実行結果]

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

1. グラフの種類

このグラフは3次元の曲面グラフで、matplotlib の plot_surface 関数を使って作成されています。

2. 軸の説明

X軸: これは領域の$x$座標を表しています。
領域は$ [0, 1] $の範囲です。
Y軸: これは領域の$y$座標を表しています。
領域は$ [0, 1] $の範囲です。
Z軸: これは未知関数 $u(x, y)$の値を表しています。

3. 曲面の形状

この曲面は、与えられた境界条件の下でのディリクレ問題の解析解を表しています。
具体的には、境界条件 u(x, y) = sin(πx) * sinh(πy) のもとで、解は u(x, y) = sin(πx) * sinh(πy) / sinh(π) となります。
曲面の形状は、この解析解の式に従って決まります。

4. 境界の値

この曲面は、$x=0$, $x=1$, $y=0$, $y=1$ の境界上で、与えられた境界条件の値をとります。
つまり、曲面の$4$辺は境界条件の値に一致しています。

5. 曲面の振る舞い

$x$方向と$y$方向で、曲面は周期的な振る舞いを示します。
これは境界条件の三角関数とひ過関数の影響です。
中心付近ではほぼ平坦ですが、境界に近づくにつれて値が大きくなります。

このように、このグラフはディリクレ問題の解析解の形状を視覚的に表現しています。

与えられた境界条件の下で、解がどのような振る舞いを示すかがわかります。
境界値問題の解を理解する上で、このようなグラフは非常に有用です。

ヘリコイドの方程式

4月 8, 2024

ヘリコイドの方程式

ヘリコイド（helix）は、螺旋状の曲面を表す数学的な形状です。

ヘリコイドの方程式は以下のように表されます:

$$
x(t) = a \cos(t), \quad y(t) = a \sin(t), \quad z(t) = bt
$$

ここで、$ ( a ) $は螺旋の半径を表し、$ ( b ) $は螺旋の傾斜（ねじれ）を制御するパラメーターです。

$ ( t ) $はパラメーターであり、ヘリコイド上の点の位置を示します。

これをPythonで解いて3Dグラフを描く方法を示します。
必要なライブラリは numpy と matplotlib です。

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

# パラメーターの設定
a = 1.0  # 螺旋の半径
b = 1.0  # 螺旋の傾斜（ねじれ）

# パラメーターの範囲
t = np.linspace(0, 4*np.pi, 1000)  # tを0から4πまで1000点で分割

# ヘリコイドの座標を計算
x = a * np.cos(t)
y = a * np.sin(t)
z = b * t

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

# ヘリコイドをプロット
ax.plot(x, y, z, label='Helix', color='b')

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

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

# 凡例の表示
plt.legend()

# グリッドの表示
ax.grid(True)

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

このコードでは、ヘリコイドのパラメーター$ ( t ) $を$0$から$ ( 4\pi ) $まで$1000$点で分割し、それに基づいて$ ( x(t) $,$ y(t) $,$ z(t) ) $を計算しています。

それらの座標を3Dグラフとして描画し、螺旋状のヘリコイドを表示しています。
a や b の値を変更することで、ヘリコイドの形状を調整することができます。

[実行結果]

ソースコード解説

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

ライブラリのインポート

1
2
3

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

numpy：数値計算を行うためのライブラリです。
ここでは数値計算や配列操作に使用されます。
matplotlib.pyplot：グラフ描画のためのライブラリです。
ここでは2Dグラフおよび3Dグラフの描画に使用されます。
mpl_toolkits.mplot3d：3Dグラフを描画するためのツールキットです。
ここでは3Dの描画機能を使用します。

パラメーターの設定

1 2	a = 1.0 # 螺旋の半径 b = 1.0 # 螺旋の傾斜（ねじれ）

a：螺旋の半径を設定します。
b：螺旋の傾斜（ねじれ）を設定します。

パラメーターの範囲設定

1	t = np.linspace(0, 4*np.pi, 1000) # tを0から4πまで1000点で分割

t：パラメーター t の範囲を設定します。
0 から 4π までを 1000 等分した配列を作成します。

ヘリコイドの座標の計算

1
2
3

x = a * np.cos(t)
y = a * np.sin(t)
z = b * t

x：ヘリコイド上の各点の x 座標を計算します。
a * np.cos(t) で計算されます。
y：ヘリコイド上の各点の y 座標を計算します。
a * np.sin(t) で計算されます。
z：ヘリコイド上の各点の z 座標を計算します。
b * t で計算されます。

3Dプロットの準備

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

fig：新しい図を作成します。
ax：3Dサブプロットを追加します。
projection='3d' を指定することで、3次元プロットを有効にします。

ヘリコイドのプロット

1	ax.plot(x, y, z, label='Helix', color='b')

ax.plot：ヘリコイドの座標 (x, y, z) をプロットします。
label='Helix' で凡例のラベルを設定し、color='b' で線の色を青に指定します。

軸ラベルの設定

1
2
3

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

ax.set_xlabel：X軸のラベルを設定します。
ax.set_ylabel：Y軸のラベルを設定します。
ax.set_zlabel：Z軸のラベルを設定します。

グラフのタイトル、凡例、グリッドの表示、グラフの表示

plt.title('Helix')
plt.legend()
ax.grid(True)
plt.show()

plt.title：グラフのタイトルを設定します。
plt.legend：凡例を表示します。
ax.grid(True)：グリッドを表示します。
plt.show()：グラフを表示します。

これにより、螺旋状のヘリコイドが3Dグラフとして描画され、軸ラベルやタイトル、凡例、グリッドが付加された見やすい図が生成されます。

結果解説

[実行結果]

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

1. 螺旋状の形状:

プロットされたグラフは、螺旋状の曲線であるヘリコイドを表します。
この曲線は螺旋状に広がり、中心軸に沿って上方向に立ち上がる形状を持ちます。
螺旋の半径 $ ( a ) $および螺旋の傾斜 $ ( b ) $によって曲線の形状が決まります。

2. 座標軸:

グラフは3次元空間を表現しており、X軸、Y軸、Z軸がそれぞれ描画されています。
X軸とY軸は平面上に、Z軸は垂直方向に伸びる軸を表します。

3. 軸ラベル:

グラフの各軸にはラベルが付けられており、X軸は ‘X’、Y軸は ‘Y’、Z軸は ‘Z’ と表示されています。
これにより、各軸が何を表しているかが分かりやすくなっています。

4. グラフのタイトル:

グラフの上部には ‘Helix’ というタイトルが表示されています。
これは描かれている曲線（ヘリコイド）の名称を示しています。

5. 凡例:

グラフには凡例が表示されており、’Helix’ というラベルが曲線の意味を示しています。

6. グリッド:

グラフ全体にはグリッドが表示されており、3次元空間を視覚的に区切り、座標軸の目盛りを補助します。

このようなグラフは、数学や物理学の学習や可視化に役立ちます。

ヘリコイドの形状を理解するだけでなく、Pythonを使用して数学的な曲線や図形をプロットする方法を学ぶのにも役立ちます。

Tominaga 方程式

4月 7, 2024

Tominaga 方程式

Tominaga 方程式は次のように定義される非線形微分方程式です:

$$
\frac{d^2 y}{dt^2} + \frac{dy}{dt} + y^3 = \cos(t)
$$

この方程式を数値的に解くために、Python を使用してグラフ化する方法を示します。

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

1
2
3

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

次に、Tominaga 方程式の微分方程式を定義します。

def tominaga_equation(t, y):
    dydt = y[1]  # y' = dy/dt
    d2ydt2 = np.cos(t) - y[1] - y[0]**3  # y'' = d^2y/dt^2
    return [dydt, d2ydt2]

ここで tominaga_equation 関数は、Tominaga 方程式を表しています。

この関数は、現在の時刻$ ( t ) $と状態$ ( y ) $を引数として受け取り、微分方程式の右辺の値を計算して返します。

具体的には、$ ( y[0] ) $は$ ( y ) $の値を表し、$ ( y[1] ) $は$ ( \frac{dy}{dt} ) $の値を表します。

次に、微分方程式を数値的に解きます。

ここでは solve_ivp 関数を使用して、初期条件から微分方程式を解きます。

# 初期条件
y0 = [0.0, 0.0]  # 初期位置と初速度

# 解く時間の範囲
t_span = [0, 20]  # 0から20までの時間範囲

# 数値的に微分方程式を解く
sol = solve_ivp(tominaga_equation, t_span, y0, t_eval=np.linspace(t_span[0], t_span[1], 1000))

# 結果をプロット
plt.figure(figsize=(8, 6))
plt.plot(sol.t, sol.y[0], label='y(t)')
plt.plot(sol.t, sol.y[1], label="y'(t)")
plt.xlabel('t')
plt.ylabel('y(t), y\'(t)')
plt.title('Solution of Tominaga Equation')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、solve_ivp 関数を使用して** Tominaga 方程式を数値的に解きます。
初期条件 y0 は [0.0, 0.0] と設定されており、初期位置と初速度を表します。
t_span は解く時間の範囲**を表しています。

解が得られた後は、解の時間変化に対する y(t) と y'(t)（速度）の値をプロットしています。
これにより、Tominaga 方程式の解の振る舞いや特性をグラフ上で視覚化することができます。

このコードを実行すると、Tominaga 方程式の解の時間変化が示されたグラフが表示されます。
解の振動や周期性などの特性を観察することができます。

[実行結果]

ソースコード解説

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

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

1
2
3

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

numpy ライブラリは数値計算用の基本的な機能を提供します。
matplotlib.pyplot ライブラリはグラフ描画用の機能を提供します。
scipy.integrate.solve_ivp 関数は微分方程式を数値的に解くための機能を提供します。

2. Tominaga 方程式の定義

def tominaga_equation(t, y):
    dydt = y[1]  # y' = dy/dt
    d2ydt2 = np.cos(t) - y[1] - y[0]**3  # y'' = d^2y/dt^2
    return [dydt, d2ydt2]

tominaga_equation 関数は、Tominaga 方程式を定義します。
引数 t は現在の時間を表し、y は状態ベクトル [y, dy/dt] を表します。
関数内では、y[1] を dy/dt として定義し、np.cos(t) - y[1] - y[0]**3 を d^2y/dt^2 として定義しています。
関数は、[dy/dt, d^2y/dt^2] をリストとして返します。

3. 初期条件の設定

1	y0 = [0.0, 0.0] # 初期位置と初速度

y0 は Tominaga 方程式の初期条件を設定します。
ここでは初期位置を 0.0、初速度を 0.0 としています。

4. 解く時間の範囲の設定

1	t_span = [0, 20] # 0から20までの時間範囲

t_span は解く時間の範囲を指定します。
ここでは 0 から 20 までの時間範囲を指定しています。

5. 数値的に微分方程式を解く

1	sol = solve_ivp(tominaga_equation, t_span, y0, t_eval=np.linspace(t_span[0], t_span[1], 1000))

solve_ivp 関数を使用して、tominaga_equation 関数で定義された微分方程式を数値的に解きます。
tominaga_equation は微分方程式を表す関数です。
t_span は解く時間の範囲を指定します。
y0 は初期条件を指定します。
t_eval=np.linspace(t_span[0], t_span[1], 1000) は、解を評価する時間点を指定します。
ここでは時間範囲を$1000$個の等間隔で分割して解を評価します。

6. 結果のプロット

plt.figure(figsize=(8, 6))
plt.plot(sol.t, sol.y[0], label='y(t)')
plt.plot(sol.t, sol.y[1], label="y'(t)")
plt.xlabel('t')
plt.ylabel('y(t), y\'(t)')
plt.title('Solution of Tominaga Equation')
plt.legend()
plt.grid(True)
plt.show()

plt.figure(figsize=(8, 6)) は図のサイズを指定します。
plt.plot(sol.t, sol.y[0], label='y(t)') と plt.plot(sol.t, sol.y[1], label="y'(t)") で、解 sol の時間変化をグラフ化します。
sol.t は時間、sol.y[0] は位置 y(t)、sol.y[1] は速度 y'(t) を表します。
plt.xlabel('t') と plt.ylabel('y(t), y\'(t)') で軸ラベルを設定します。
plt.title('Solution of Tominaga Equation') でグラフのタイトルを設定します。
plt.legend() で凡例を表示します。
plt.grid(True) でグリッドを表示します。
plt.show() でグラフを表示します。

これにより、Tominaga 方程式の解 y(t) とその導関数 y'(t) の時間変化がグラフで可視化されます。

グラフから、解の振動や特性を観察することができます。

結果解説

[実行結果]

グラフを説明します。

1. 横軸 (t):

時間$ ( t ) $を表します。
この例では、$0$から$20$までの時間範囲を等間隔で分割した点が横軸上に表示されます。

2. 縦軸 (y(t), y’(t)):

$ ( y(t) ) $と$ ( y’(t) ) $の値を表します。
$ ( y(t) ) $は Tominaga 方程式の解であり、時間 $ ( t ) $に対する位置や変位を表します。
$ ( y’(t) ) $は$ ( \frac{dy}{dt} ) $の値であり、時間 $ ( t ) $に対する速度を表します。

3. グラフの表示:

グラフには、時間$ ( t ) $に対する$ ( y(t) ) $と$ ( y’(t) ) $の変化がプロットされます。
$ ( y(t) ) $の曲線は Tominaga 方程式の解の時間変化を示し、振動や非線形な特性が視覚化されます。
$ ( y’(t) ) $の曲線は$ ( y(t) ) $の時間変化に対する速度を示し、解の変化率や速度の振る舞いを表します。

4. タイトル:

グラフの上部には「Solution of Tominaga Equation」というタイトルが表示されます。
これは、表示されている曲線が Tominaga 方程式の数値解を表していることを示します。

5. 凡例 (Legend):

グラフには「y(t)」と「y’(t)」という凡例が表示されています。
「y(t)」は Tominaga 方程式の解$ ( y(t) ) $を表し、位置や変位の時間変化を示します。
「y’(t)」は Tominaga 方程式の解の導関数$ ( y’(t) ) $を表し、速度や変化率の時間変化を示します。

グラフを通じて、Tominaga 方程式の解の時間変化や振る舞いが視覚的に理解できます。
特に非線形な微分方程式の場合、解の振動や周期性、または解がどのように非線形項$ ( \cos(t) ) $によって影響を受けるかがグラフから読み取れます。

数値的な解法を用いた場合は、計算の精度や初期条件の影響を考慮する必要があります。
また、解の振動や安定性を観察するために、時間範囲や時間刻みの選択によって解の挙動が異なることもあります。

ブロック波動方程式

4月 6, 2024

ブロック波動方程式

ブロック波動方程式は、1次元の時間依存しない波動方程式です。

一般的には、以下のように表されます：

$$
\frac{\partial^2 u}{\partial x^2} = -k^2 u
$$

ここで、$ ( u ) $は波動の変位、$ ( x ) $は空間座標、$ ( k ) $は波数です。

以下は、ブロック波動方程式を解いてグラフ化する Python コードの例です。

import numpy as np
import matplotlib.pyplot as plt

# 空間座標の範囲
x_values = np.linspace(0, 10, 100)

# 波数
k = 2

# ブロック波動方程式の解の計算
u_values = np.sin(k * x_values)

# グラフの描画
plt.plot(x_values, u_values)
plt.xlabel('Position (x)')
plt.ylabel('Displacement (u)')
plt.title('Solution of the Block Wave Equation')
plt.grid(True)
plt.show()

このコードでは、$ ( x ) $の範囲を$ (0) $から$ (10) $までの間で取り、その範囲における波動の変位 $ ( u ) $を計算しています。

そして、その解をグラフにプロットしています。

[実行結果]

ソースコード解説

このPythonのソースコードは、ブロック波動方程式の解を計算し、グラフ化するものです。

以下はソースコードの詳細です：

1. NumPyとMatplotlibのインポート:

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

NumPyは数値計算を行うためのPythonライブラリであり、ここでは配列や数学関数の使用に役立ちます。
Matplotlibはグラフ描画ライブラリであり、ここではグラフを描画するために使用されます。

2. 空間座標の範囲設定:

1	x_values = np.linspace(0, 10, 100)

np.linspace() 関数は、指定した範囲内で等間隔の値を生成します。
ここでは、$ (0) $から$ (10) $の範囲を$ (100) $個の点に分割しています。

3. 波数の設定:

k = 2

ブロック波動方程式の解における波数$ (k) $を定義しています。

4. ブロック波動方程式の解の計算:

1	u_values = np.sin(k * x_values)

$ (x) $座標に対するブロック波動方程式の解を計算しています。
ここでは単純化のために、sin関数を用いています。

5. グラフの描画:

plt.plot(x_values, u_values)
plt.xlabel('Position (x)')
plt.ylabel('Displacement (u)')
plt.title('Solution of the Block Wave Equation')
plt.grid(True)
plt.show()

plt.plot() 関数で$ (x) $座標と$ (u) $座標の値を使ってグラフを描画します。
plt.xlabel() と plt.ylabel() で軸ラベルを設定し、 plt.title() でグラフのタイトルを設定します。
plt.grid(True) でグリッドを表示し、 plt.show() でグラフを表示します。

これにより、ブロック波動方程式の解のグラフが生成され、 $ (x) $座標と$ (u) $座標の関係が視覚化されます。

結果解説

[実行結果]

このグラフは、ブロック波動方程式の解を示しています。

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

横軸（X軸）：

空間座標$ (x)$。この軸は、波の位置を表します。
範囲は$ (0) $から$ (10) $までです。

縦軸（Y軸）：

変位$ (u)$。この軸は、波の変位を表します。
範囲は$ (-1) $から$ (1) $までです。

波形：

$ (u = \sin(kx)) $の関数を表しており、これはブロック波動方程式の解です。
$ (k) $は波数を示し、この例では$ (k = 2) $としています。
この波形は、$ (x) $の値に応じて振幅が変化し、正弦関数の波形を持ちます。

タイトル：

“Solution of the Block Wave Equation”。グラフの内容を簡潔に示しています。

軸ラベル：
“Position (x)” と “Displacement (u)”。それぞれ横軸と縦軸の内容を説明しています。

このグラフを見ると、波が$ (x) $の値に応じて振動している様子がわかります。

波の変位は$ (x) $に対して正弦関数の形をしており、周期的な振動を示しています。

シュワルツの不等式

4月 5, 2024

シュワルツの不等式

シュワルツの不等式とは、2つのベクトルの内積の絶対値の2乗は、それぞれのベクトルのノルム(長さ)の積以下であるということを示す不等式です。

つまり、ベクトル$a = (a1, a2, …, an)$、ベクトル$b = (b1, b2, …, bn)$に対して、

$ (a1 * b1 + a2 * b2 + … + an * bn)^2 ≦ (a1^2 + a2^2 + … + an^2) * (b1^2 + b2^2 + … + bn^2) $

が成り立ちます。

左辺はベクトル$a$・$b$の内積の絶対値の2乗、右辺は $a$のノルムと$b$のノルムの積を表しています。

この不等式は、ベクトル内積が最大になるのは、2つのベクトルが同じ方向を向いている時であり、ベクトル間の角度が大きくなるほど内積は小さくなることを意味しています。

幾何学的には、2つのベクトルを組み合わせて作った平行四辺形の面積は、底辺と高さの積以下であることを表しています。

このようにシュワルツの不等式は、ベクトルの内積とノルムの関係を不等式の形で表した重要な結果となっています。

サンプルコード

シュワルツの不等式は、以下の式で表されます。

$$
(x^2 + y^2) * (1 + x^2 + y^2) >= 1
$$

この不等式は、2変数 x と y の関数として表すことができます。

次のコードでは、シュワルツの不等式を3D空間上に可視化しています。

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

# シュワルツの不等式を表す関数
def f(x, y):
    return (x**2 + y**2) * (1 + x**2 + y**2)

# グラフ化する範囲
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

# グラフを描画
fig = plt.figure(figsize=(10, 8))
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('Schwarz Inequality Visualization')

plt.show()

まず、f(x, y) 関数でシュワルツの不等式の右辺を定義しています。

次に、x と y の値の範囲を設定しています。
ここでは、xとyがそれぞれ -2 から 2 までの範囲を 100 分割して計算しています。
np.meshgrid関数を使うことで、XとYの格子点上の値を求めています。

Z = f(X, Y) では、XとYの格子点上での関数値を計算しています。

最後に、matplotlibとmpl_toolkits.mplot3dを使って、X、Y、Zの値から3D曲面を描画しています。
plot_surface関数を使うことで、格子点上の値から滑らかな曲面を描くことができます。
cmap='viridis'は、曲面の色付けに使われるカラーマップを指定しています。

実行すると、x、y、z軸に対して、シュワルツの不等式の右辺を表す3D曲面が描画されます。
この曲面の形状から、不等式の性質を視覚的に確認することができます。

このように、任意の$2$変数関数を指定して、その関数値を3D空間上に可視化することができます。

[実行結果]

ソースコード解説

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

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	def f(x, y): return (x2 + y2) * (1 + x2 + y2)

f(x, y) 関数は、2つの引数 x と y を受け取り、シュワルツの不等式を表す関数です。
式 (x**2 + y**2) * (1 + x**2 + y**2) を計算して結果を返します。

3. グラフ化する範囲の設定:

x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

np.linspace(-2, 2, 100) によって、範囲$ ([-2, 2]) $を$100$個の等間隔に区切った数列を生成し、x と y に代入します。
np.meshgrid(x, y) によって、x と y のすべての組み合わせに対する格子状の座標行列 X と Y を生成します。
f(X, Y) によって、X と Y の各座標における関数 f(x, y) の値を計算し、Z に格納します。

4. グラフの描画:

fig = plt.figure(figsize=(10, 8))
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('Schwarz Inequality Visualization')

plt.show()

plt.figure(figsize=(10, 8)) によって、図のサイズを設定します。
fig.add_subplot(111, projection='3d') によって、3Dグラフを描画するためのサブプロットを追加します。
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis') によって、3D 曲面を描画します。
rstride と cstride はメッシュの密度を制御するパラメータです。
ax.set_xlabel('x')、ax.set_ylabel('y')、ax.set_zlabel('z') によって、各軸のラベルを設定します。
ax.set_title('Schwarz Inequality Visualization') によって、図のタイトルを設定します。
plt.show() によって、描画したグラフを表示します。

結果解説

[実行結果]

出力される3D曲面は、シュワルツの不等式 $(x^2 + y^2) * (1 + x^2 + y^2) >= 1 $の右辺の値を、$x-y$平面上の点$(x, y)$について表しています。
カラーマップの濃い色ほど、不等式の右辺の値が大きくなっています。

この曲面の形状から、不等式の性質を視覚的に確認できます。
例えば、原点$(0, 0)$周辺では不等式の両辺が等しくなり、原点から離れるほど右辺の値が大きくなっていきます。
また、曲面は原点を中心として点対称の形状をしており、不等式の性質である「$x^2 + y^2 >= 0$」も確認できます。

このグラフを見ることで、シュワルツの不等式の形状的な特徴や、成り立つ範囲などを直感的に把握することができます。