食事計画(線形計画問題)

食事計画(線形計画問題)

CVXPYを使って線形計画問題を解きます。

ここでは、食事計画の問題を考えます。

この問題では、複数の食品から適切な組み合わせを選んで、特定の栄養要件を満たしつつ、コストを最小限に抑えることを目的とします。

問題設定

食品の種類:

  1. チキン (Chicken)
  2. 牛肉 (Beef)
  3. 豆 (Beans)

各食品のコストと栄養素含有量は以下の通り:

  • チキン: コスト 2ドル、タンパク質 20g、脂肪 5g、炭水化物 0g
  • 牛肉: コスト 3ドル、タンパク質 15g、脂肪 10g、炭水化物 0g
  • 豆: コスト 1ドル、タンパク質 10g、脂肪 2g、炭水化物 20g

栄養要件:

  • タンパク質: 最低 30g
  • 脂肪: 最低 10g、最大 20g
  • 炭水化物: 最低 20g

解法

この問題をCVXPYを使って解きます。

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
import cvxpy as cp
import numpy as np

# 各食品のコストと栄養素含有量
cost = np.array([2, 3, 1])
protein = np.array([20, 15, 10])
fat = np.array([5, 10, 2])
carbs = np.array([0, 0, 20])

# 栄養要件
min_protein = 30
min_fat = 10
max_fat = 20
min_carbs = 20

# 食品の量 (変数)
x = cp.Variable(3, nonneg=True)

# コスト最小化
objective = cp.Minimize(cost @ x)

# 制約条件
constraints = [
protein @ x >= min_protein,
fat @ x >= min_fat,
fat @ x <= max_fat,
carbs @ x >= min_carbs
]

# 問題の定義と解決
problem = cp.Problem(objective, constraints)
problem.solve()

# 結果の表示
print(f"Status: {problem.status}")
print(f"Optimal value: {problem.value:.2f} dollars")
print(f"Optimal solution: {x.value}")
print(f"Chicken: {x.value[0]:.2f} units")
print(f"Beef: {x.value[1]:.2f} units")
print(f"Beans: {x.value[2]:.2f} units")

このコードでは、以下の手順で食事計画問題を解きます:

  1. 各食品のコスト栄養素含有量を定義します。
  2. 栄養要件を設定します。
  3. 変数 x を定義し、各食品の量を表します(非負制約)。
  4. 目的関数コストの最小化として定義します。
  5. 栄養要件に基づく制約を設定します。
  6. 問題を定義し、解決します。
  7. 最適なコスト食品の量を表示します。

このコードを実行すると、栄養要件を満たしつつコストを最小限に抑える最適な食品の組み合わせが得られます。

[実行結果]

Status: optimal
Optimal value: 3.72 dollars
Optimal solution: [0.64 0.48 1.  ]
Chicken: 0.64 units
Beef: 0.48 units
Beans: 1.00 units

ソースコード解説

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

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

1
2
import cvxpy as cp
import numpy as np

ここでは、最適化問題を解くためのCVXPYライブラリと、数値計算を行うためのNumPyライブラリをインポートしています。

2. 各食品のコストと栄養素含有量の定義

1
2
3
4
cost = np.array([2, 3, 1])
protein = np.array([20, 15, 10])
fat = np.array([5, 10, 2])
carbs = np.array([0, 0, 20])

ここでは、各食品(チキン、ビーフ、豆)のコストと栄養素(タンパク質、脂肪、炭水化物)の含有量を配列として定義しています。

例えば、チキンのコストは2ドルで、タンパク質が20グラム、脂肪が5グラム含まれています。

3. 栄養要件の定義

1
2
3
4
min_protein = 30
min_fat = 10
max_fat = 20
min_carbs = 20

ここでは、最低限必要なタンパク質、脂肪の最小および最大量、最低限必要な炭水化物の量を定義しています。

4. 変数の定義

1
x = cp.Variable(3, nonneg=True)

ここでは、最適化問題の変数である各食品の量を表すベクトル x を定義しています。
この変数は非負(負の値を取らない)であることが指定されています。

5. コスト最小化の目的関数

1
objective = cp.Minimize(cost @ x)

ここでは、総コストを最小化する目的関数を定義しています。
cost @ x はコストと食品の量の内積を計算し、これを最小化します。

6. 制約条件の定義

1
2
3
4
5
6
constraints = [
protein @ x >= min_protein,
fat @ x >= min_fat,
fat @ x <= max_fat,
carbs @ x >= min_carbs
]

ここでは、栄養要件に基づく制約条件を定義しています。
具体的には、次の制約を課しています:

  • 総タンパク質が30グラム以上
  • 総脂肪が10グラム以上20グラム以下
  • 総炭水化物が20グラム以上

7. 最適化問題の定義と解決

1
2
problem = cp.Problem(objective, constraints)
problem.solve()

ここでは、目的関数制約条件をもとに最適化問題を定義し、これを解決しています。
problem.solve() によって、最適解が計算されます。

8. 結果の表示

1
2
3
4
5
6
print(f"Status: {problem.status}")
print(f"Optimal value: {problem.value:.2f} dollars")
print(f"Optimal solution: {x.value}")
print(f"Chicken: {x.value[0]:.2f} units")
print(f"Beef: {x.value[1]:.2f} units")
print(f"Beans: {x.value[2]:.2f} units")

ここでは、最適化のステータス最適なコスト、各食品の最適な量を表示しています。

結果解説

[実行結果]

Status: optimal
Optimal value: 3.72 dollars
Optimal solution: [0.64 0.48 1.  ]
Chicken: 0.64 units
Beef: 0.48 units
Beans: 1.00 units
  • Status: optimal:
    問題が最適に解決されたことを示しています。
  • Optimal value: 3.72 dollars:
    最適な食品の組み合わせのコストが3.72ドルであることを示しています。
  • Optimal solution:
    各食品の最適な量を示しています。
    具体的には、チキン0.64単位、ビーフ0.48単位、豆1.00単位です。

これにより、指定された栄養要件を満たしつつ、コストを最小限に抑えるための最適な食事計画が得られます。

巡回セールスマン問題(Traveling Salesman Problem, TSP)

巡回セールスマン問題(Traveling Salesman Problem, TSP)

最適化問題の一例として、「巡回セールスマン問題(Traveling Salesman Problem, TSP)」を挙げます。

TSPは、いくつかの都市を訪問する最短経路を求める問題で、非常に計算困難です。

今回は、Pythonを使ってTSPを解き、結果をグラフ化します。

巡回セールスマン問題の説明

巡回セールスマン問題は、以下のように定義されます:

  • あるセールスマンが$N$個の都市を訪れ、それぞれの都市を一度だけ訪問し、出発地点に戻る必要がある。
  • すべての都市間の距離が与えられている。
  • セールスマンが移動する総距離を最小化する経路を見つけることが目的。

解法

TSPNP困難問題であり、大規模な問題を厳密に解くのは現実的ではありません。
ここでは、近似解法の一つである遺伝的アルゴリズム(Genetic Algorithm)を用いて解を求めます。

実装

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

1
pip install deap matplotlib numpy

次に、Pythonコードを以下に示します:

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
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import numpy as np
import matplotlib.pyplot as plt
from deap import base, creator, tools, algorithms

# 都市数
num_cities = 20

# ランダムな都市の座標を生成
np.random.seed(42)
cities = np.random.rand(num_cities, 2)

# 距離行列を計算
distance_matrix = np.sqrt(((cities[:, np.newaxis] - cities)**2).sum(axis=2))

# 遺伝的アルゴリズムの設定
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", np.ndarray, fitness=creator.FitnessMin)

toolbox = base.Toolbox()
toolbox.register("indices", np.random.permutation, num_cities)
toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.indices)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

def evalTSP(individual):
distance = 0
for i in range(num_cities):
distance += distance_matrix[individual[i-1], individual[i]]
return distance,

toolbox.register("mate", tools.cxOrdered)
toolbox.register("mutate", tools.mutShuffleIndexes, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)
toolbox.register("evaluate", evalTSP)

def main():
population = toolbox.population(n=300)
ngen = 400
cxpb = 0.7
mutpb = 0.2

result, log = algorithms.eaSimple(population, toolbox, cxpb, mutpb, ngen,
stats=None, halloffame=None, verbose=False)

best_ind = tools.selBest(population, 1)[0]
return best_ind

best_route = main()

# 結果をプロット
plt.figure(figsize=(10, 5))
plt.scatter(cities[:, 0], cities[:, 1], color='red')
for i in range(num_cities):
plt.annotate(i, (cities[i, 0], cities[i, 1]))
route = np.append(best_route, best_route[0])
plt.plot(cities[route, 0], cities[route, 1], 'b-')
plt.title('Traveling Salesman Problem - Best Route')
plt.xlabel('X Coordinate')
plt.ylabel('Y Coordinate')
plt.grid(True)
plt.show()

コードの説明

  1. 都市の生成:

    • num_citiesを設定し、ランダムな都市の座標を生成します。
    • 都市間の距離行列を計算します。
  2. 遺伝的アルゴリズムの設定:

    • DEAPライブラリを使用して遺伝的アルゴリズムを設定します。
    • creator個体(経路)とその適合度を定義します。
    • toolboxで個体の生成交叉突然変異選択、および評価関数を登録します。
  3. 適合度関数:

    • evalTSP関数で個体(経路)の総距離を計算します。
  4. 主関数:

    • 遺伝的アルゴリズムの主要なパラメータを設定し、eaSimple関数を使用して進化を実行します。
    • 最良の経路を選択し、結果を返します。
  5. 結果のプロット:

    • 最良の経路をプロットし、都市の座標経路を視覚化します。

結果

このコードを実行すると、TSP最良の経路がプロットされます。

[実行結果]

都市間の最短経路が青線で表示され、各都市の座標が赤点で示されます。

双曲放物面(hyperbolic paraboloid)

双曲放物面(hyperbolic paraboloid)

双曲放物面(hyperbolic paraboloid)は、特定の二次曲面で、鞍型(saddle shape)として知られる形状を持っています。

これは幾何学的に興味深い構造を持つため、建築土木工学などさまざまな分野で応用されています。

双曲放物面の方程式

双曲放物面は次の一般形の二次方程式で表されます:

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

ここで、$(a) $と$ (b) $は正の実数であり、曲面の形状を定義します。
特に、双曲放物面放物線双曲線の特性を組み合わせたものです。

特徴と性質

  1. 鞍型:

    • 双曲放物面の形状は「鞍型」とも呼ばれ、中心点の周りで異なる曲率を持ちます。
      片方の軸に沿っては上に開いた放物線の形を持ち、もう片方の軸に沿っては下に開いた放物線の形を持ちます。
    • 中心点(原点)では、曲面は凹型凸型が交互に現れ、典型的なのような形状を形成します。
  2. 双曲線と放物線の交差:

    • 双曲放物面は、$x$軸に沿っては放物線$ ( z = \frac{x^2}{a^2} ) $の形状を持ち、$y$軸に沿っては放物線$ ( z = -\frac{y^2}{b^2} ) $の形状を持ちます。
    • 任意の平行断面($z$ = 定数)の形状は双曲線となります。
  3. 無限に広がる曲面:

    • 双曲放物面は無限に広がる曲面であり、$x$軸や$y$軸の方向に制限なく広がります。
  4. 曲率:

    • 双曲放物面の曲率は負であり、これは曲面が中心点の周りで異なる方向に曲がっていることを意味します。

実生活での応用

双曲放物面はその構造的な強度と美しさから、建築土木工学で広く使用されています。

以下はそのいくつかの応用例です:

  1. 建築:

    • モダンな建築物や屋根構造で、双曲放物面はしばしば使用されます。
      例えば、シドニーオペラハウスの屋根の形状は双曲放物面を基にしています。
  2. 橋梁:

    • 双曲放物面の形状は橋のデザインにも応用され、構造的な効率性と美しさを兼ね備えたデザインが可能となります。
  3. インフラストラクチャー:

    • 双曲放物面はパビリオンシェル構造パラボラアンテナなどにも利用されます。
      特に、シェル構造は軽量でありながら高い強度を持ちます。

以下のコードでは、双曲放物面の形状を描き、カラフルなカラーマップを適用します:

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

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

# 放物面の方程式
x = u
y = v
z = u**2 - v**2

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

# 放物面の描画
surf = ax.plot_surface(x, y, z, cmap='coolwarm', edgecolor='none')

# カラーバーの追加
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

# ラベルとタイトルの設定
ax.set_xlabel('X Axis')
ax.set_ylabel('Y Axis')
ax.set_zlabel('Z Axis')
ax.set_title('3D Hyperbolic Paraboloid')

# 視点の設定
ax.view_init(elev=45, azim=60)

plt.show()

ソースコード解説

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

    • numpyは数値計算を行うために使用します。
    • matplotlib.pyplotはグラフ描画のために使用します。
    • mpl_toolkits.mplot3dは3Dプロットを描画するために使用します。
  2. パラメータの設定:

    • uvは$-2$から$2$までの範囲で等間隔に$100$個の点を生成します。
    • np.meshgridは、2つの1次元配列から2次元の格子を生成します。
  3. 放物面の方程式:

    • x, y, zは放物面の座標を表します。
      放物面の方程式は$z = u^2 - v^2$です。
  4. 3Dプロットの作成:

    • plt.figure(figsize=(10, 8))で新しい図を作成し、figに割り当てます。
      figsizeで図のサイズを指定します。
    • fig.add_subplot(111, projection='3d')で3Dサブプロットを作成し、axに割り当てます。
  5. 放物面の描画:

    • ax.plot_surface(x, y, z, cmap='coolwarm', edgecolor='none')で放物面を描画します。
      cmapはカラーマップを指定し、edgecolor='none'でエッジの色を無効にします。
  6. カラーバーの追加:

    • fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)でカラーバーを追加し、色の範囲を視覚化します。
  7. ラベルとタイトルの設定:

    • ax.set_xlabel, ax.set_ylabel, ax.set_zlabelで各軸のラベルを設定します。
    • ax.set_titleで図のタイトルを設定します。
  8. 視点の設定:

    • ax.view_init(elev=45, azim=60)3Dプロットの視点を設定します。
      elev垂直方向の角度azim水平方向の角度を指定します。
  9. プロットの表示:

    • plt.show()でプロットを表示します。

このコードを実行することで、3Dの放物面がよりカラフルに、美しく描画されます。
カラーマップと視点を工夫することで、視覚的に非常に魅力的なグラフを作成できます。

グラフ解説

[実行結果]

このコードを実行すると、上記のような美しい3Dの放物面が表示されます:

双曲放物面の形状:

$z = u^2 - v^2$に従って計算された双曲放物面が描画されます。
この形状は鞍型で、特定の曲率を持っています。

カラーマップ:

高さに基づいて色が変化するため、放物面の形状が視覚的に強調されます。
coolwarmカラーマップは青から赤へのグラデーションを作成し、見栄えが良くなります。

カラーバー:

放物面の高さに対応する色の範囲が示されるため、特定の高さがどの色に対応するかが一目でわかります。

軸ラベルとタイトル:

$x$軸、$y$軸、$z$軸のラベルとグラフのタイトルが表示され、グラフの内容が明確になります。

  • 視点の設定:
    適切な視点から放物面を観察できるため、形状カラーマップが立体的に見えます。

このグラフは、数学的な関数の視覚化、特に双曲放物面の形状と特性を示すのに非常に有用です。

双曲線2葉面

双曲線2葉面

双曲線2葉面は、特定のパラメータ方程式によって表される3次元曲面です。

この曲面は、次の特徴を持ちます:

  1. 双曲線を含む:
    曲面上の曲線は双曲線であり、無限遠点に向かって広がります。
    これらの双曲線は、2つの独立した葉を形成しています。

  2. 無限遠点に向かって広がる:
    曲面は中心を持たず、無限遠点に向かって広がります。
    これにより、曲面は非常に特殊な形状をしています。

  3. 葉の両側対称性:
    曲面は葉ごとに両側対称的であり、それぞれの葉は双曲線の繰り返しパターンを持ちます。

双曲線2葉面は、数学幾何学の分野で興味深い研究対象とされています。

その特異な形状と幾何学的性質は、数々の興味深い現象や応用を持ちます。

プログラム例

以下は、この双曲線2葉面をPythonで描画するコードです。

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

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

# 双曲線2葉面の方程式
x = np.cosh(u) * np.cos(v)
y = np.cosh(u) * np.sin(v)
z = np.sinh(u)

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

# 双曲線2葉面の描画
ax.plot_surface(x, y, z, cmap='viridis')

# ラベルとタイトルの設定
ax.set_xlabel('X Axis')
ax.set_ylabel('Y Axis')
ax.set_zlabel('Z Axis')
ax.set_title('Hyperbolic Surface')

plt.show()

このコードでは、双曲線2葉面パラメトリック方程式を使用して描画しています。

結果のグラフは、双曲線の特徴を持つ立体的な形状を示します。

[実行結果]

ソースコード解説

このソースコードは、Pythonのnumpymatplotlibを使用して双曲線2葉面を3Dプロットするものです。

それぞれの部分を詳しく見ていきましょう。

  1. パッケージのインポート:
1
2
3
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
  • numpyは数値計算を行うためのPythonパッケージです。
  • matplotlib.pyplotはグラフ描画のためのPythonパッケージで、pltとしてエイリアスされています。
  • mpl_toolkits.mplot3dは3次元の描画を扱うためのサブモジュールです。
  1. パラメータの設定:
1
2
3
u = np.linspace(-2, 2, 100)
v = np.linspace(0, 2*np.pi, 100)
u, v = np.meshgrid(u, v)
  • np.linspaceは指定された範囲内で等間隔の数値を生成します。
    ここでは、uvそれぞれに$100$個の値が生成されます。
  • np.meshgridは、$2$つの1次元配列から2次元の格子を生成します。
    ここでは、uvからそれぞれの点の組み合わせを生成しています。
  1. 双曲線2葉面の方程式:
1
2
3
x = np.cosh(u) * np.cos(v)
y = np.cosh(u) * np.sin(v)
z = np.sinh(u)
  • 双曲線2葉面の方程式を定義しています。
    これは3次元座標系上での曲面です。
  1. 3Dプロットの作成:
1
2
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
  • plt.figure()で新しい図を作成し、figに割り当てます。
  • fig.add_subplot(111, projection='3d')で3Dサブプロットを作成し、axに割り当てます。
  1. 双曲線2葉面の描画:
1
ax.plot_surface(x, y, z, cmap='viridis')
  • ax.plot_surfaceを使用して双曲線2葉面を描画します。
    xyzはそれぞれの座標軸の値であり、cmapはカラーマップを指定します。
  1. ラベルとタイトルの設定:
1
2
3
4
ax.set_xlabel('X Axis')
ax.set_ylabel('Y Axis')
ax.set_zlabel('Z Axis')
ax.set_title('Hyperbolic Surface')
  • ax.set_xlabelax.set_ylabelax.set_zlabelでそれぞれの軸のラベルを設定します。
  • ax.set_titleで図のタイトルを設定します。
  1. プロットの表示:
1
plt.show()
  • plt.show()でプロットを表示します。

これらのステップを組み合わせることで、双曲線2葉面の3Dプロットが生成されます。

グラフ解説

[実行結果]

双曲線2葉面は、3次元空間内に双曲線のような曲線を持つ曲面です。

この曲面は、次のような特徴を持っています:

  1. 曲面は双曲線を繰り返し、2つの異なる双曲線が螺旋状に延びています。
  2. 曲面は中心を持たず、無限遠方に向かって広がっています。
  3. 曲面は両側対称であり、曲線が2つの独立した葉を形成しています。

グラフ上では、双曲線2葉面は双曲線のような複雑なパターンを持ち、立体的な特徴を示しています。

各葉は無限遠方に延びており、立体的な形状を視覚化することができます。

また、色の変化により曲面の表面の凹凸が強調されています。

球面調和関数

球面調和関数

球面調和関数は、球対称の問題に対する解析解を提供するための数学的な関数です。

主に、量子力学天文学などの分野で重要な役割を果たします。

球面調和関数は次のような形式を持ちます:

$$
Y_{lm}(\theta, \phi) = C_{lm} P_l^m(\cos\theta) e^{im\phi}
$$

ここで、$( \theta ) $は極角、$( \phi ) $は方位角です。
$ ( P_l^m ) $は連関勝手関数、$( C_{lm} ) $は規格化定数です。
$ ( l ) $は角運動量量子数、$( m ) $は磁気量子数です。

球面調和関数は、次のような特徴を持ちます:

1. 角度に対する依存性:

$( \theta ) $と$ ( \phi ) $の両方の角度に依存します。
これにより、球面調和関数は球面上の位置に関する情報を提供します。

2. 球対称性:

$( Y_{lm} ) $は球対称関数であり、球面上のすべての方向で同じ値を持ちます。
これは球対称問題に対して解析的な解を提供します。

3. 量子数 ( l ) と ( m ) の影響:

$( l ) $と$ ( m ) $の値によって、球面調和関数の形状と振る舞いが決まります。
$ ( l ) $は角運動量の大きさを示し、$( m ) $は角運動量の $z $成分を示します。

4. 規格化:

$( Y_{lm} ) $は規格化されており、単位球上での値が$1$になります。
これにより、物理的な解釈や計算が容易になります。

球面調和関数は、量子力学で原子や分子の波動関数を記述するために広く使用されます。
また、天文学では、球対称の天体天体の角運動量に関連する問題を解決するためにも使用されます。

プログラム例

球面調和関数を用いてカラフルな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
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# パラメータの設定
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
theta, phi = np.meshgrid(theta, phi)
r = np.abs(np.sin(3*theta) * np.cos(3*phi))

# カラーマップの設定
colors = cm.viridis(r / r.max())

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

# 球面調和関数の描画
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)

ax.plot_surface(x, y, z, facecolors=colors, rstride=1, cstride=1, linewidth=0, antialiased=False)

ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Spherical Harmonics")

plt.show()

このコードは、球面調和関数を使用して3Dグラフを描画します。

球面調和関数には$ ( r = \sin(3\theta) \cdot \cos(3\phi) ) $が使用され、この関数は$ ( \theta ) $と$ ( \phi ) $の角度に応じて球面上の値を定義します。

グラフの色は、球面調和関数の値に応じて変化し、カラーマップが使用されています。

[実行結果]

ソースコード解説

3Dグラフの描画:球面調和関数の可視化

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

  • numpy は数値計算を行うために使用されます。
  • matplotlib.pyplot はグラフ描画のために使用されます。
  • mpl_toolkits.mplot3d は3Dグラフ描画のためのサブモジュールです。
  • matplotlib.cm はカラーマップを操作するために使用されます。

2. パラメータの設定:

  • theta は極角の値を表す$0$から$ ( \pi ) $の範囲を、phi は方位角の値を表す$0$から$ ( 2\pi ) $の範囲をそれぞれ$100$の点に分割します。
  • thetaphi のメッシュグリッドを作成し、球面上の点を表します。
  • r は球面上の各点における球面調和関数の値を計算します。

3. カラーマップの設定:

  • cm.viridis はカラーマップを設定します。
    このカラーマップは、球面調和関数の値に応じて色を変化させます。
  • r の値を正規化して、色を適用します。

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

  • plt.figure() で新しい図を作成します。
  • fig.add_subplot(111, projection='3d') で3Dグラフ用のサブプロットを作成します。
  • ax.plot_surface()球面調和関数を3Dプロットします。
    facecolors パラメータにカラーマップを設定し、球面調和関数の値に基づいて色を変化させます。
  • rstridecstride は、メッシュの行と列のストライド(間隔)を指定します。
    これにより、プロットの解像度が制御されます。
  • linewidthラインの幅を、antialiasedアンチエイリアス処理の有無を制御します。

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

  • ax.set_xlabel()ax.set_ylabel()ax.set_zlabel() を使用して、それぞれ$X$軸、$Y$軸、$Z$軸のラベルを設定します。
  • ax.set_title() で図のタイトルを設定します。

6. グラフの表示:

  • plt.show() を呼び出して、作成したグラフを表示します。

このコードは、球面調和関数を計算し、それを3Dプロットして視覚化します。
球面調和関数の値に応じて色が変化することで、球面の形状や振る舞いが視覚的に理解されます。

グラフ解説

[実行結果]

この3Dグラフは、球面調和関数を可視化しています。
球面調和関数は、球対称のポテンシャル波動関数を表現するための数学的な関数です。
この関数は、量子力学天文学などのさまざまな分野で重要な役割を果たします。

このグラフでは、球面調和関数 $ ( r = \sin(3\theta) \cdot \cos(3\phi) ) $を使用して描画されています。
ここで、$( \theta ) $は極角 $(0から ( \pi ) $まで)、$( \phi ) $は方位角($0$から$ ( 2\pi ) $まで)です。
この関数により、球面上の各点に対する値が定義されます。

球面上の各点の$ ( r ) $の値に基づいて、色がグラフに適用されます。
このカラーマップは、色を$ ( r ) $の値に応じて変化させ、球面の表面の色が異なる深さや形状を示すことができます。
こうしたカラーマップを使用することで、球面調和関数の特性を視覚的に理解することができます。

このグラフは、球面調和関数の形状や振る舞いを視覚的に示し、特に量子力学天文学の分野での理解を深めるのに役立ちます。

ラプラス方程式

ラプラス方程式

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

ラプラス方程式は次の形をしています:
$$
\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で示しており、ポテンシャルの変化を視覚的に理解するのに役立ちます。

ディルクレの波動方程式

ディルクレの波動方程式

ディルクレの波動方程式は、数学や物理学で使用される偏微分方程式の一つで、波の伝播をモデル化します。

以下では、ディルクレの波動方程式の概要物理的意味数値解法について詳しく説明します。

ディルクレの波動方程式とは

波動方程式の一般形は次のように表されます:

$$
\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u
$$

  • $( u(x, y, t) ) $は位置 $ ( (x, y) ) $と時間 $ ( t ) $における波の振幅を示します。
  • $( c ) $は波の速度です。
  • $( \nabla^2 ) $はラプラス演算子であり、2次元の場合、次のように書かれます:

$$
\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}
$$

ディルクレ境界条件とは、領域の境界において関数$ ( u ) $がゼロになる条件を指します。
つまり、波動方程式の解$ ( u ) $が領域の端でゼロになるようにします。

物理的意味

波動方程式は、弦の振動電磁波音波など、さまざまな物理現象をモデル化します。
ディルクレ境界条件は、固定された端を持つ弦や膜の振動を表す場合によく使われます。
例えば、ギターの弦の両端が固定されている状態を考えると、弦の端での変位がゼロになるディルクレ条件を適用することができます。

数値解法の概要

波動方程式を数値的に解くためには、有限差分法がよく使用されます。
ここでは、空間と時間を離散化し、差分方程式を用いて波動方程式を解きます。
数値解法のステップは以下の通りです:

  1. 離散化:

    • 空間 $ ( x ) $と$ ( y ) $を格子状に分割し、時間 $ ( t ) $も小さなステップ$ ( dt ) $で分割します。
  2. 初期条件と境界条件の設定:

    • 初期条件として、時間 $( t=0 ) $における波の形状 $( u(x, y, 0) ) $とその時間微分 $( \frac{\partial u}{\partial t}(x, y, 0) ) $を設定します。
    • ディルクレ境界条件として、領域の境界上での波の振幅をゼロに設定します。
  3. 時間ステップごとの更新:

    • 中心差分法を用いて、各時間ステップにおける波の振幅を計算します。
      具体的には、次の差分方程式を使用します:

$$
u_{i,j}^{n+1} = 2u_{i,j}^n - u_{i,j}^{n-1} + \frac{c^2 dt^2}{dx^2} (u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n - 4u_{i,j}^n)
$$

ここで、$( u_{i,j}^n ) $は時間ステップ$ ( n ) $における格子点$ ( (i, j) ) $での波の振幅です。

プログラム例

ディルクレの波動方程式は、典型的には次のような形式をとる二次元の波動方程式です:

$$
\frac{\partial^2 u}{\partial t^2} = c^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)
$$

ここで、$( u(x, y, t) ) $は時間 $( t ) $における位置 $( (x, y) ) $での波の振幅、$ ( c )$ は波の速度です。

ディルクレ境界条件は、領域の境界上で$ ( u ) $がゼロになることを意味します。

この波動方程式を数値的に解いて、3次元グラフとして可視化します。

ここでは、単純な正方形領域内の初期条件と境界条件を使用します。

次の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
44
45
46
47
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータの設定
L = 1.0 # 領域の長さ
T = 0.01 # 時間の長さ
c = 1.0 # 波の速度
dx = 0.01 # 空間ステップ
dt = 0.001 # 時間ステップ

# 空間格子と時間格子の生成
x = np.arange(0, L+dx, dx)
y = np.arange(0, L+dx, dx)
t = np.arange(0, T+dt, dt)
nx = len(x)
ny = len(y)
nt = len(t)

# 初期条件の設定
u = np.zeros((nx, ny, nt))
u[:, :, 0] = np.sin(np.pi * x[:, None]) * np.sin(np.pi * y[None, :])
u[:, :, 1] = u[:, :, 0]

# 波動方程式の数値解法(中心差分法)
for n in range(1, nt-1):
for i in range(1, nx-1):
for j in range(1, ny-1):
u[i, j, n+1] = (2 * u[i, j, n] - u[i, j, n-1] +
c**2 * dt**2 / dx**2 *
(u[i+1, j, n] + u[i-1, j, n] + u[i, j+1, n] + u[i, j-1, n] - 4 * u[i, j, n]))

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

X, Y = np.meshgrid(x, y)
Z = u[:, :, nt//2]

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

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U')
ax.set_title('Wave Equation at t = {:.3f}'.format(t[nt//2]))

plt.show()

このコードでは、波動方程式を数値的に解くために中心差分法を使用しています。

初期条件としては、正弦波の形状を使用しています。

また、ディルクレ境界条件を満たすために、境界上では常にゼロになるようにしています。

このスクリプトを実行すると、時間 $( t = T/2 ) $における波動の振幅 $( u(x, y, t) ) $を3Dプロットとして視覚化することができます。

[実行結果]

ソースコード解説

このPythonスクリプトは、ディルクレの波動方程式を2次元空間内で数値的に解き、時間のある時点での波動の振幅を3Dグラフとして可視化します。

以下では、スクリプトを章立てて詳しく説明します。

1. 必要なライブラリのインポート

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

この部分では、数値計算に必要なNumPyと、グラフ作成のためのMatplotlib、および3Dプロット用のツールキットをインポートしています。

2. パラメータの設定

1
2
3
4
5
6
# パラメータの設定
L = 1.0 # 領域の長さ
T = 0.01 # 時間の長さ
c = 1.0 # 波の速度
dx = 0.01 # 空間ステップ
dt = 0.001 # 時間ステップ

ここでは、計算領域のサイズ$ ( L )$、シミュレーションの時間長$ ( T )$、波の速度$ ( c )$、空間ステップ $( dx )$、および時間ステップ$ ( dt ) $を設定しています。

3. 空間格子と時間格子の生成

1
2
3
4
5
6
7
# 空間格子と時間格子の生成
x = np.arange(0, L+dx, dx)
y = np.arange(0, L+dx, dx)
t = np.arange(0, T+dt, dt)
nx = len(x)
ny = len(y)
nt = len(t)

ここでは、空間と時間の格子点を生成しています。

np.arangeを使って、$0$から$( L )$までの間をステップサイズ$( dx )$で分割した配列$ ( x ) $と$ ( y )$ を作成します。

時間についても同様に配列$ ( t ) $を作成します。

4. 初期条件の設定

1
2
3
4
# 初期条件の設定
u = np.zeros((nx, ny, nt))
u[:, :, 0] = np.sin(np.pi * x[:, None]) * np.sin(np.pi * y[None, :])
u[:, :, 1] = u[:, :, 0]

この部分では、波動の振幅を格納する3次元配列$ ( u ) $を初期化しています。

初期条件として、正弦波の形状を与えています。

x[:, None]y[None, :]を用いて二次元配列を作成し、初期時刻$ ( t=0 ) $での振幅を設定します。

次の時間ステップ$ ( t=dt ) $も同じ初期条件に設定します。

5. 波動方程式の数値解法(中心差分法)

1
2
3
4
5
6
7
# 波動方程式の数値解法(中心差分法)
for n in range(1, nt-1):
for i in range(1, nx-1):
for j in range(1, ny-1):
u[i, j, n+1] = (2 * u[i, j, n] - u[i, j, n-1] +
c**2 * dt**2 / dx**2 *
(u[i+1, j, n] + u[i-1, j, n] + u[i, j+1, n] + u[i, j-1, n] - 4 * u[i, j, n]))

この部分では、中心差分法を用いて波動方程式を数値的に解いています。

各時間ステップ$ ( n ) $に対して、空間格子点$ ( (i, j) ) $での波動の振幅を更新します。

中心差分法を用いることで、時間と空間の二階微分を近似的に計算しています。

6. 3Dプロットの作成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

X, Y = np.meshgrid(x, y)
Z = u[:, :, nt//2]

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

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('U')
ax.set_title('Wave Equation at t = {:.3f}'.format(t[nt//2]))

plt.show()

最後に、時間$ ( t = T/2 ) $における波動の振幅を3Dプロットとして表示します。

np.meshgridを用いて2次元の格子点$ ( X ) $と$ ( Y ) $を作成し、それに対応する波動の振幅 $( Z )$ を取得します。

plot_surfaceを用いて3Dプロットを描画し、軸ラベルとタイトルを設定しています。

総括

このスクリプトは、ディルクレの波動方程式を2次元空間内で数値的に解き、特定の時点での波動の振幅を3Dグラフとして視覚化するものです。

波動方程式の数値解法として中心差分法を使用し、初期条件として正弦波を設定しています。

結果として得られた波動の振幅は、時間の中間点で3Dプロットとして表示されます。

グラフ解説

[実行結果]

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

1. X軸とY軸:

  • グラフの$X$軸と$Y$軸は、空間内の位置を表しています。
    具体的には、正方形領域内の点$ ( (x, y) ) $を示しています。
    この領域は、$0$から$1$の範囲に設定されています。

2. Z軸:

  • $Z$軸は、波動の振幅 $( u(x, y, t) ) $を示しています。
    これは、時間のある特定の時点$ ( t ) $における空間の各点での波の高さを意味します。

3. 波動の振幅の変化:

  • グラフは、時間 $( t = T/2 ) $における波動の振幅を3次元的に表示しています。

波動の振幅は、初期条件として与えた正弦波形に基づいて、時間と共に変化します。
この特定の時間における波の形状が、3Dプロットとして描画されます。

4. 色の変化:

  • プロットの色は振幅の高さに対応しています。
    通常、色のスキーム(ここではviridisカラーマップ)が使用され、低い振幅は寒色(青や緑)で、高い振幅は暖色(黄やオレンジ)で示されます。

三角波

三角波

三角波は、音楽信号処理でよく使われる基本的な波形の一つで、名前の通り三角形の形をしています。

三角波は直線的に上昇し、最高点に達すると直線的に下降するという特徴を持ちます。

これは、次の特徴を持ちます。

  1. 周期的: 三角波は一定の周期で繰り返されます。
  2. 対称性: 上昇部分下降部分が対称で、鋭い頂点を持ちます。
  3. 線形変化: 各周期内で直線的に増加し、減少します。

数学的表現

三角波の数学的表現は以下のように定義されます。

$$
\text{tri}(t) = 2 \left| \frac{t}{T} - \left\lfloor \frac{t}{T} + 0.5 \right\rfloor \right|
$$

ここで、$ ( T ) $は周期を表します。

三角波はこの数式により、時間に対して周期的に増減する波形として描かれます。

ソースコード例

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

# 三角波の定義
def triangle_wave(t, period=1):
return 2 * np.abs(2 * (t / period - np.floor(t / period + 0.5))) - 1

# 時間軸と空間軸の設定
t = np.linspace(0, 10, 500)
x = np.linspace(-5, 5, 500)
T, X = np.meshgrid(t, x)
Z = triangle_wave(T)

# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(T, X, Z, cmap='viridis')

# グラフのラベル
ax.set_xlabel('Time')
ax.set_ylabel('X')
ax.set_zlabel('Amplitude')
ax.set_title('3D Triangle Wave')

plt.show()

[実行結果]

ソースコード解説

このPythonソースコードの説明を行います。

1. 必要なライブラリのインポート

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

ここでは、NumPyMatplotlib、および3Dプロットのためのmpl_toolkits.mplot3dをインポートしています。

NumPy数値計算用Matplotlibデータ可視化用のライブラリです。

2. 三角波の定義

1
2
def triangle_wave(t, period=1):
return 2 * np.abs(2 * (t / period - np.floor(t / period + 0.5))) - 1

この関数は、指定された周期で三角波を生成します。

数式に基づいて、時間 t に対する三角波の値を計算します。

period波の周期を設定する引数です。

3. データの設定

1
2
3
4
t = np.linspace(0, 10, 500)
x = np.linspace(-5, 5, 500)
T, X = np.meshgrid(t, x)
Z = triangle_wave(T)

ここでは、時間軸 t空間軸 x を設定します。

np.linspace 関数は、指定した範囲で等間隔に分割された値を生成します。

np.meshgrid 関数を使用して、2Dグリッドデータ TX を作成し、三角波関数 triangle_wave を用いて Z を計算します。

4. 3Dプロットの作成

1
2
3
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(T, X, Z, cmap='viridis')

ここでは、3Dプロットを作成します。

fig は図を生成し、ax は3Dプロットのための軸を設定します。

ax.plot_surface は、TX、および Z を使って表面プロットを作成し、カラーマップ viridis を適用します。

5. グラフのラベル付け

1
2
3
4
ax.set_xlabel('Time')
ax.set_ylabel('X')
ax.set_zlabel('Amplitude')
ax.set_title('3D Triangle Wave')

ここでは、3Dグラフにラベルとタイトルを設定します。

set_xlabelset_ylabel、および set_zlabel メソッドで各軸にラベルを付け、set_title メソッドでグラフのタイトルを設定します。

6. グラフの表示

1
plt.show()

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

これにより、作成した3D三角波グラフが表示されます。

グラフ解説

[実行結果]

1. 時間軸 (Time):

  • $( t ) $は時間を表し、三角波の周期的な変化を示します。
  • このプロットでは、時間$ ( t ) $が$0$から$10$までの範囲を持ちます。

2. 空間軸 (X):

  • $( x ) $は空間を表し、プロットの幅を示します。
  • このプロットでは、空間$ ( x ) $が$-5$から$5$までの範囲を持ちます。

3. 振幅 (Amplitude):

  • $( z ) $は三角波の振幅を示します。
  • 三角波の特徴として、振幅は周期ごとに上下します。

4. 3Dプロット:

  • plot_surface 関数を使って、時間 $ ( t ) $と空間 $ ( x ) $に対する三角波の振幅を3Dで視覚化します。
  • カラーマップ viridis は、振幅の値に対する色を割り当てます。

この3Dプロットは、時間軸に沿って周期的に変化する三角波を視覚的に表現し、三角波の特性をより直感的に理解するのに役立ちます。

レムニスケート(Lemniscate)

レムニスケート(Lemniscate)

レムニスケート(Lemniscate)は、双曲線のような形状を持つ数学的な曲線です。

上記のグラフは、レムニスケート・オブ・ベルヌーイ(Lemniscate of Bernoulli)をプロットしたものです。

この曲線は、極座標で表現すると次のようになります:

$$
r^2 = 2 \cos(2\theta)
$$

このグラフの特性をいくつか説明します:

  1. 形状: 曲線は数字の8の字のような形をしています。
    これが「レムニスケート」という名前の由来です。
    レムニスケートはラテン語で「リボン」を意味します。

  2. シンメトリー: 曲線は原点に対して対称であり、$x$軸および$y$軸に対しても対称です。

  3. 定義域: このグラフでは、パラメータ$ ( t ) $を$ (-\pi) $から$ (\pi) $まで変化させて曲線を描いています。

  4. 座標の変換: $x$と$y$の座標は次の式に基づいて計算されます。
    $$
    x = \frac{\sqrt{2} \cos(t)}{1 + \sin^2(t)}
    $$
    $$
    y = \frac{\sqrt{2} \cos(t) \sin(t)}{1 + \sin^2(t)}
    $$

このように、レムニスケートは美しい数学的性質を持ち、解析学代数幾何学において重要な役割を果たしています。

ソースコード例

以下にレムニスケートをPythonでグラフ化するためのソースコードを示します。

このコードを実行すると、レムニスケート・オブ・ベルヌーイの形状が描かれます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
import matplotlib.pyplot as plt

# Define the parameter t
t = np.linspace(-np.pi, np.pi, 400)

# Lemniscate of Bernoulli equations
x = np.sqrt(2) * np.cos(t) / (1 + np.sin(t)**2)
y = np.sqrt(2) * np.cos(t) * np.sin(t) / (1 + np.sin(t)**2)

# Plot the Lemniscate
plt.figure(figsize=(8, 8))
plt.plot(x, y)
plt.title("Lemniscate of Bernoulli")
plt.xlabel("x")
plt.ylabel("y")
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

このコードについての簡単な説明:

  1. numpymatplotlib のライブラリをインポートします。
  2. t という変数を$ (-\pi) $から$ (\pi) $までの$400$個の等間隔の点に設定します。
  3. レムニスケート・オブ・ベルヌーイの$ x $座標と$ y $座標を計算します。
  4. plt.plot() 関数を使用して、計算された$ x $と$ y $の値をプロットします。
  5. タイトルや軸ラベル、グリッド線を追加して、プロットを整えます。
  6. plt.show() 関数を使用して、プロットを表示します。

このコードを実行すると、レムニスケートの美しい8の字の形が描かれます。

グラフ解説

[実行結果]

レムニスケート・オブ・ベルヌーイ(Lemniscate of Bernoulli)のグラフについて詳しく説明します。

この曲線は、極座標で次のように表される特定の形状の曲線です:

$$
r^2 = 2 \cos(2\theta)
$$

グラフの内容の説明:

  1. 形状:

    • グラフは数字の「8」の字の形をしています。
      これはレムニスケートの特徴的な形状であり、名前の由来ともなっています(ラテン語で「リボン」を意味します)。
  2. シンメトリー:

    • この曲線は、原点(中心点)に対して対称です。
      また、$x$軸および$y$軸に対しても対称です。
      つまり、グラフの各象限に同じ形状の部分が現れます。
  3. 範囲:

    • $t$パラメータは$(-\pi) $から$ (\pi) $までの範囲を取ります。
      これにより、曲線の全体が描かれます。
  4. 数学的な背景:

    • レムニスケートは、楕円関数複素解析において重要な役割を果たす数学的な曲線です。
      特に、代数幾何学における研究対象となります。
  5. 座標の計算:

    • 曲線の$ x $座標と$ y $座標は以下の式によって計算されます:
      $$
      x = \frac{\sqrt{2} \cos(t)}{1 + \sin^2(t)}
      $$
      $$
      y = \frac{\sqrt{2} \cos(t) \sin(t)}{1 + \sin^2(t)}
      $$
    • これらの式により、曲線の各点の座標が計算され、プロットされます。
  6. プロットの特徴:

    • グラフは8の字を描き、中心に向かって絞られた形状となります。
    • 縦横の長さが等しいため、グラフは正方形の図形に内接するように表示されます。

リサージュ曲線(Lissajous curve)

リサージュ曲線(Lissajous curve)

リサージュ曲線(Lissajous curve)は、2つの調和振動(サイン波やコサイン波)の合成によって得られる曲線で、フランスの物理学者ジュール・アントワーヌ・リサージュによって1860年代に研究されました。

これらの曲線は、特に電子工学信号処理物理学などの分野で、周期的な現象や波動の比較に役立ちます。

基本方程式

リサージュ曲線は次のパラメトリック方程式で表されます:

$$
x(t) = A \sin(a t + \delta)
$$
$$
y(t) = B \sin(b t)
$$

  • $( t ) $は時間またはパラメータです。
  • $( A ) $と $( B ) $は、それぞれ$ ( x ) $軸と$ ( y ) $軸方向の振幅です。
  • $( a ) $と $( b ) $は、それぞれ$ ( x ) $軸と$ ( y ) $軸方向の周波数です。
  • $( \delta ) $は位相差で、$ ( x ) $軸方向の波と$ ( y ) $軸方向の波の間の初期位相の違いを示します。

特性と形状

リサージュ曲線の形状は、パラメータ$ ( a )$、$ ( b )$、$ ( \delta ) $によって決まります。
これらのパラメータの値によって、曲線は非常に多様な形状をとります。

  1. 周波数比$ ( \frac{a}{b} ) $の効果:

    • 整数比の場合(例えば$、( \frac{a}{b} = 1 )$、$( \frac{a}{b} = 2 ) $など)、リサージュ曲線は閉じた形になります。
    • 非整数比の場合、曲線は閉じられない複雑な形状を持ち、長い時間範囲で描かれるとカオス的なパターンになります。
  2. 位相差$ ( \delta ) $の効果:

    • 位相差が$ ( 0 )$ または$ ( \pi ) $の場合、曲線は直線になります(特に、振幅$ ( A ) $と$ ( B ) $が同じ場合)。
    • $( \delta ) $が$ ( \frac{\pi}{2} ) $の場合、典型的なリサージュ曲線(楕円形や八の字形)を生成します。

実用例

リサージュ曲線は、以下のような分野で重要な役割を果たします:

  • オシロスコープ: オシロスコープで2つの信号の周波数位相差を視覚的に比較するために使用されます。
    1つの信号が$ ( x ) $軸、もう1つの信号が$ ( y ) $軸に接続されると、リサージュ図形が表示されます。
  • 信号処理: 2つの波形信号の同期や比較に使用されます。
  • 物理学とエンジニアリング: 振動分析周期現象の研究に役立ちます。

プログラム例

以下は、Pythonを用いてリサージュ曲線をプロットするコードです:

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
import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
A = 1 # x軸の振幅
B = 1 # y軸の振幅
a = 3 # x軸の周波数
b = 2 # y軸の周波数
delta = np.pi / 2 # 位相差

# 時間の配列
t = np.linspace(0, 2 * np.pi, 1000)

# x(t) と y(t) の定義
x = A * np.sin(a * t + delta)
y = B * np.sin(b * t)

# グラフのプロット
plt.figure(figsize=(8, 8))
plt.plot(x, y)
plt.title('Lissajous Curve')
plt.xlabel('x(t)')
plt.ylabel('y(t)')
plt.grid(True)
plt.axis('equal')
plt.show()

結果の解釈

  • 閉じた曲線: $( \frac{a}{b} ) $が有理数であれば、曲線は閉じます。
    例えば、上記のコードでは$ ( \frac{a}{b} = \frac{3}{2} ) $であり、曲線は繰り返しパターンを示します。
  • 対称性: $( \delta ) $の値により対称性が変わります。
    例えば、$ ( \delta = 0 ) $の場合、$ ( x ) $と$ ( y ) $の振動は同期し、直線楕円になります。

リサージュ曲線は、異なる調和振動の相互作用を視覚的に表現するための強力なツールであり、その形状はパラメータの調整により多様に変化します。

これにより、信号の特性相互関係を解析する際に非常に有用です。

グラフ解説

[実行結果]

リサージュ曲線のグラフの内容について説明します。

グラフの内容説明

  1. $X$軸と$Y$軸の振動:

    • $(x(t) = A \sin(a t + \delta))$
    • $(y(t) = B \sin(b t))$
      この2つの振動が組み合わさることで、リサージュ曲線が形成されます。
  2. 振幅$ (A) $と$ (B)$:

    • 振幅$ (A) $と$ (B) $は、それぞれ$ (x) $軸と$ (y) $軸方向の振動の最大値を示します。
    • 振幅が大きいほど、波の高さが高くなります。
  3. 周波数$ (a) $と$ (b)$:

    • 周波数$ (a) $と$ (b) $は、$ (x)$ 軸と$ (y) $軸方向の振動の周期を決定します。
    • 周波数が高いほど、波が速く振動します。
  4. 位相差$ (\delta)$:

    • 位相差$ (\delta) $は、$ (x) $軸方向の振動が$ (y) $軸方向の振動に対してどれだけずれているかを示します。
    • 位相差が異なると、リサージュ曲線の形状が変わります。

プロットされたグラフ

  • 曲線の形状

    • 曲線の形状は、パラメータ$ (a)$、$ (b)$、$ (\delta) $によって決まります。
      例えば、$ (a) $と$ (b) $が整数比であれば、閉じた滑らかな曲線が得られます。

    非整数比の場合、複雑なカオス的な形状になることがあります。

  • 例示されたパラメータ

    • $(A = 1)$、$ (B = 1)$: $(x) $軸と$ (y) $軸方向の振幅は同じです。
    • $(a = 3)$、$ (b = 2)$: $(x) $軸方向の振動が$ (y) $軸方向の振動に対して$3:2の$比率を持ちます。
    • $(\delta = \pi / 2)$: $(x) $軸方向の振動が$ (y) $軸方向の振動に対して$ (\pi / 2) $ラジアン($90$度)ずれています。

結果のプロット

  • 軸のラベル: $(x) $軸と$ (y) $軸がそれぞれラベル付けされており、曲線が2次元平面上に描かれています。
  • グリッド: グリッドが表示され、曲線の形状を確認しやすくしています。
  • 等スケール:$ (x) $軸と$ (y) $軸のスケールを同じにすることで、曲線の形状が正確に表示されます。

これらの特性により、リサージュ曲線は2つの調和振動の位相差振幅比を視覚的に理解するための強力なツールとなります。

リサージュ曲線は、電子工学音響学などで信号の特性を解析する際によく用いられます。