バーニー龍（Barnsley Dragon）

バーニー龍（Barnsley Dragon）は、数学的なフラクタルであり、フラクタルジオメトリの一例です。

これは、1988年にマイケル・バーンズリー（Michael Barnsley）によって提案されました。
バーニー龍は、確率的な方法で生成されるフラクタル図形であり、自己相似性を示します。

バーニー龍の生成は、$2$つの確率的な変換に基づいています。
これらの変換は、$4\times4$の行列で表され、点を新しい位置に移動させます。

その後、ランダムに変換が選択され、次の位置に点が移動します。
このプロセスを繰り返すことで、複雑なフラクタル図形が生成されます。

バーニー龍は、ジュリア集合の一種であり、幾何学的に美しいパターンを持ちます。
多くの場合、計算機グラフィックスや美術、数学の教育などの分野で使用され、その視覚的な魅力と数学的な興味深さから広く注目されています。

ソースコード

以下は、Pythonを使用してバーニー龍を生成し、matplotlibを使用してグラフ化する例です。

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

import numpy as np
import matplotlib.pyplot as plt

def barnsley_dragon(n_points):
    # 初期点
    x, y = 0, 0
    # バーニー龍の変換確率
    p = [0.7875, 0.2125]
    # 変換行列
    m = np.array([[0.824074, 0.281482, -0.212346, 0.864198, 1.6, -0.4],
                  [-0.212346, 0.864198, 0.281482, 0.824074, 0.0, 1.6]])

    points = [(x, y)]
    for _ in range(n_points):
        # ランダムに変換を選択
        i = np.random.choice(range(len(p)), p=p)
        # 変換を適用
        x, y = m[i, 0]*x + m[i, 1]*y + m[i, 4], m[i, 2]*x + m[i, 3]*y + m[i, 5]
        points.append((x, y))
    
    return points

def plot_barnsley_dragon(points):
    # バーニー龍を描画
    plt.figure(figsize=(8, 8))
    plt.title("Barnsley Dragon")
    plt.scatter(*zip(*points), s=0.2, color='green')
    plt.axis('equal')
    plt.axis('off')
    plt.show()

# ポイントの数
n_points = 100000
# バーニー龍の生成
points = barnsley_dragon(n_points)
# バーニー龍の描画
plot_barnsley_dragon(points)

このコードは、指定された数の点を生成し、それらの点をmatplotlibを使用してプロットします。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

この部分は、NumPyとMatplotlibをインポートしています。
NumPyは数値計算のためのライブラリであり、Matplotlibはグラフの描画や可視化のためのライブラリです。

2. バーニー龍の生成関数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def barnsley_dragon(n_points):
   # 初期点
   x, y = 0, 0
   # バーニー龍の変換確率
   p = [0.7875, 0.2125]
   # 変換行列
   m = np.array([[0.824074, 0.281482, -0.212346, 0.864198, 1.6, -0.4],
                 [-0.212346, 0.864198, 0.281482, 0.824074, 0.0, 1.6]])

   points = [(x, y)]
   for _ in range(n_points):
       # ランダムに変換を選択
       i = np.random.choice(range(len(p)), p=p)
       # 変換を適用
       x, y = m[i, 0]*x + m[i, 1]*y + m[i, 4], m[i, 2]*x + m[i, 3]*y + m[i, 5]
       points.append((x, y))
   
   return points

この部分は、バーニー龍を生成する関数です。
バーニー龍は確率的な方法で生成されるため、初期点から始めてランダムな変換を繰り返し適用します。

3. バーニー龍の描画関数:

1
2
3
4
5
6
7
8

def plot_barnsley_dragon(points):
   # バーニー龍を描画
   plt.figure(figsize=(8, 8))
   plt.title("Barnsley Dragon")
   plt.scatter(*zip(*points), s=0.2, color='green')
   plt.axis('equal')
   plt.axis('off')
   plt.show()

この部分は、バーニー龍を描画する関数です。
matplotlib.pyplot.scatterを使用して、生成された点をプロットし、matplotlib.pyplot.showでグラフを表示します。

4. バーニー龍の生成と描画:

1
2
3
4
5
6

# ポイントの数
n_points = 100000
# バーニー龍の生成
points = barnsley_dragon(n_points)
# バーニー龍の描画
plot_barnsley_dragon(points)

この部分では、指定された数の点を生成し、barnsley_dragon関数でバーニー龍を生成し、plot_barnsley_dragon関数で描画しています。

結果解説

[実行結果]

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

1. バーニー龍の生成:

barnsley_dragon関数は、指定された数の点を生成します。
この関数は、バーニー龍の確率的な変換を実行し、各変換に基づいて点を移動します。
バーニー龍は、$2$つの確率変換からなり、各変換には$6$つの変換行列のパラメータがあります。
これらのパラメータは、バーニー龍の形状と配置に影響します。

2. グラフの描画:

plot_barnsley_dragon関数は、生成された点をmatplotlibを使用してプロットします。
各点は緑色の小さな円で表示され、その位置はバーニー龍の形状を反映しています。
グラフのタイトルは「Barnsley Dragon」となっています。

3. グラフの設定:

plt.figure関数によって、図のサイズが設定されています。
plt.title関数によって、グラフにタイトルが追加されます。
plt.scatter関数は、点を散布図としてプロットします。
plt.axis('equal')によって、$x$軸と$y$軸の比率が等しくなります。
plt.axis('off')は、$x$軸と$y$軸の目盛りや枠線を非表示にします。

4. 表示される内容:

グラフに表示される内容は、生成されたバーニー龍の形状です。
バーニー龍は、幾何学的に複雑なフラクタル構造を持つ図形であり、ランダムな変換を繰り返すことで生成されます。
このグラフは、そのようなフラクタル図形の一例を示しています。

マンデルブロ集合

マンデルブロ集合は、複素数平面上で定義される非常に興味深い集合です。
その定義と性質は以下の通りです。

定義

マンデルブロ集合 Mは、複素平面 C上の点の集合で、以下の漸化式が発散しない点の集まりです。

$$
Z(n+1) = Z(n)^2 + C
$$

ここで、$Z(0) = 0 $で、$C $は複素平面上の任意の点です。

性質

• マンデルブロ集合はフラクタル構造を持ち、自己相似性があります。
  つまり、任意の拡大率で見ても同じような模様が現れ続けます。
• 集合は完全に連結しており、ある一点から出発すれば、集合の任意の点に到達できます。
• 境界上には無数の小さな突起があり、さらにズームインすると、より小さな構造が無限に現れます。
• 主要な構造体の周りには、無数の小さな副次的構造体があります。
• 黒い部分は集合に属し、発散しない点を表します。
  有色の部分は集合に属さず、その色は発散までの反復回数に応じて変化します。

数学的側面

• マンデルブロ集合は、複素数方程式 $ Z(n+1) = Z(n)^2 + C $の解の挙動を視覚化したものです。
• この方程式は、カオスや分岐理論などの分野で重要な役割を果たします。
• マンデルブロ集合の形状は、方程式の解の安定性や分岐を反映しています。

マンデルブロ集合は美しさと複雑さを併せ持ち、コンピューターによる視覚化によって発見された驚くべき数学的オブジェクトです。

単純な式から生み出されながら、無限の構造を含む不思議な性質を持っています。

ソースコード

Python でマンデルブロ集合をプロットすることができます。

以下のコードでは、numpy と matplotlib を使用しています。

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

def mandelbrot(c, max_iter=100):
    z = c
    for n in range(max_iter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return max_iter

width, height = 4, 4
re_start, re_end = -2, 1
im_start, im_end = -1.5, 1.5

re = np.linspace(re_start, re_end, width * 100)
im = np.linspace(im_start, im_end, height * 100)

mandel = np.zeros((height * 100, width * 100))
for i, x in enumerate(re):
    for j, y in enumerate(im):
        mandel[j, i] = mandelbrot(x + 1j * y)

fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(mandel, cmap='hot', extent=(re_start, re_end, im_start, im_end))
ax.set_title("Mandelbrot Set", fontsize=20)

plt.show()

このコードの解説:

1. mandelbrot 関数は、複素数 c に対して、マンデルブロ集合の値を計算します。
   最大反復回数は max_iter に設定されています。

2. width、height でプロット領域のアスペクト比を設定し、re_start、re_end、im_start、im_end で複素平面の表示範囲を設定しています。

3. np.linspace を使って、実数部と虚数部の値の範囲を生成しています。

4. mandel という空の NumPy 配列を作成し、その要素に mandelbrot 関数の出力値を格納していきます。

5. plt.imshow で mandel 配列をプロットし、cmap でカラーマップを設定しています。extent で表示範囲を指定しています。

実行すると、以下のようなマンデルブロ集合のプロットが表示されます。

[実行結果]

このプログラムでは、繰り返し計算を効率的に行うために NumPy の配列計算を活用しています。

計算量が多いため、大きな解像度で描画するとかなり時間がかかります。

必要に応じて max_iter を調整したり、計算範囲を狭めるなどの最適化が必要です。

また、カラーマップを変更したり、ズームして細部を描画したりと、さまざまな表現が可能です。

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

このコードは、NumPyとmatplotlibライブラリをインポートしています。
NumPyは数値計算に使用され、matplotlibはデータの可視化に使用されます。

2. マンデルブロ集合の計算関数

1
2
3
4
5
6
7

def mandelbrot(c, max_iter=100):
    z = c
    for n in range(max_iter):
        if abs(z) > 2:
            return n
        z = z*z + c
    return max_iter

この関数は、複素数cに対するマンデルブロ集合の値を計算します。
最大反復回数max_iterを設定し、その回数までに発散しない場合はmax_iterを返します。
発散した場合は、発散するまでの反復回数を返します。

3. グラフ描画領域の設定

1
2
3

width, height = 4, 4
re_start, re_end = -2, 1
im_start, im_end = -1.5, 1.5

ここでは、グラフの描画領域のアスペクト比(width, height)と、複素平面上の表示範囲(re_start, re_end, im_start, im_end)を設定しています。

4. 複素平面上の点の生成

1
2

re = np.linspace(re_start, re_end, width * 100)
im = np.linspace(im_start, im_end, height * 100)

np.linspaceを使って、実数部と虚数部の値の範囲を生成しています。
これにより、複素平面上の点を表すグリッドが作成されます。

5. マンデルブロ集合の値の計算と格納

1
2
3
4

mandel = np.zeros((height * 100, width * 100))
for i, x in enumerate(re):
    for j, y in enumerate(im):
        mandel[j, i] = mandelbrot(x + 1j * y)

まず、mandel配列を$0$で初期化します。
次に、実数部と虚数部の値を組み合わせて複素数を作り、mandelbrot関数を呼び出します。
その結果を、mandel配列の対応する位置に格納します。

6. グラフの描画

1
2
3
4

fig, ax = plt.subplots(figsize=(10, 10))
ax.imshow(mandel, cmap='hot', extent=(re_start, re_end, im_start, im_end))
ax.set_title("Mandelbrot Set", fontsize=20)
plt.show()

ここでは、matplotlibを使ってマンデルブロ集合をプロットしています。
ax.imshowでmandel配列を画像として表示し、cmap='hot'でカラーマップを設定しています。
extent引数で表示範囲を指定し、ax.set_titleでタイトルを設定しています。
最後にplt.show()でグラフを表示します。

このコードでは、NumPyの効率的な配列計算を活用しながら、マンデルブロ集合の計算と可視化を行っています。
計算量が多いため、大きな解像度で描画するとかなり時間がかかります。

必要に応じてmax_iterを調整したり、計算範囲を狭めるなどの最適化が必要です。

結果解説

[実行結果]

マンデルブロ集合のグラフには、以下のような特徴が表れています。

全体の構造

• 中心付近に主要な構造体があり、その周りに無数の小さな構造体が存在しています。
• 主要な構造体は、見た目が肺のような形をしており、複雑で自己相似的なフラクタル構造を持っています。
• グラフの外側は黒くなっており、これは発散する点を表しています。

主要な構造体の詳細

• 主要な構造体の内側は、黒い部分と輝く黄色の部分が入り組んだ複雑な模様になっています。
• この構造体の境界線上には、無数の小さな突起や飛び出した部分があり、さらに細かい構造が無限に続いています。
• 主要な構造体の外側にも、小さな副次的な構造体が多数存在しています。

カラーマップの意味

• 使用されているカラーマップ (‘hot’)は、青から赤へと徐々に変化する色合いを表しています。
• 青色の部分は発散するまでの反復回数が少ないことを示し、赤色の部分は反復回数が多いことを示します。
• 黒い部分は発散しないで無限に計算が続く点を表しています。

自己相似性

• マンデルブロ集合の最も印象的な特徴は、自己相似性です。
• つまり、構造体の一部を拡大していくと、似たような構造が無限に現れ続けます。
• この自己相似性は、フラクタルの典型的な性質です。

マンデルブロ集合は、単純な複素数方程式から生み出された驚くべき構造を持つフラクタルであり、数学的にも美しさがあります。

このグラフは、その複雑で魅力的な構造を視覚化したものです。

三次元のサーフェスプロット

plotly ライブラリを使って、三次元のサーフェスプロットを描画してみましょう。

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 plotly.graph_objects as go

# データの範囲を設定
x = np.linspace(-3, 3, 50)
y = np.linspace(-3, 3, 50)

# 三次元のサーフェスプロットのデータを生成
X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X**2 + Y**2))

# サーフェスプロットを描画
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])

# グラフのタイトルと軸ラベルを設定
fig.update_layout(
    title='Beautiful 3D Surface Plot',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z'
    ),
    width=800,
    height=600
)

# カラースケールを設定
fig.update_traces(colorscale='Viridis')

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

このコードでは、以下の手順で3Dサーフェスプロットを描画しています。

1. np.linspace 関数を使って、$x$座標と$y$座標の範囲を設定します。
2. np.meshgrid 関数を使って、$x$座標と$y$座標の組み合わせを生成します。
3. 三次元のサーフェスプロットのデータ Z を生成します。
   ここでは、$sin(sqrt(x^2 + y^2))$ という式を使っています。
4. plotly ライブラリの go.Surface 関数を使って、サーフェスプロットを描画します。
5. グラフのタイトルと軸ラベルを設定します。
6. fig.update_traces 関数を使って、サーフェスプロットのカラースケールを設定します。
   ここでは 'Viridis' というカラースケールを使っています。
7. fig.show 関数を使って、グラフを表示します。

実行すると、以下のような美しい3Dサーフェスプロットが表示されます。

[実行結果]

このサーフェスプロットは、sin関数を使った美しい波打つ形状をしています。
カラースケールの 'Viridis' は、緑から紫へとなめらかに変化する色合いで、サーフェスの立体感を強調しています。

このようなサーフェスプロットは、様々な分野のデータを可視化するのに役立ちます。

例えば、地形データ、気象データ、物理現象のシミュレーションデータなどを3次元で表現できます。

Pythonの plotly を使えば、このような美しい3Dグラフを簡単に描画できます。

ソースコード解説

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

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

1
2

import numpy as np
import plotly.graph_objects as go

最初に、NumPyとPlotlyライブラリをインポートしています。
NumPyは数値計算に使われる便利なライブラリで、Plotlyは対話的な可視化ツールです。

2. データの範囲の設定

1
2

x = np.linspace(-3, 3, 50)
y = np.linspace(-3, 3, 50)

ここでは、$x$座標と$y$座標のデータ範囲を設定しています。
np.linspace関数を使って、$-3$から$3$までの範囲を$50$個の等間隔のデータポイントに分割しています。

3. サーフェスプロットのデータ生成

1
2

X, Y = np.meshgrid(x, y)
Z = np.sin(np.sqrt(X**2 + Y**2))

次に、サーフェスプロットのデータを生成しています。
np.meshgrid関数を使って、$x$と$y$の組み合わせを2次元配列として生成しています。
そして、 Z = np.sin(np.sqrt(X**2 + Y**2))という式を使って、$z$座標のデータを計算しています。
この式は、円形の波打つ形状を表しています。

4. サーフェスプロットの描画

1

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])

Plotlyのgo.Surface関数を使って、サーフェスプロットを描画するためのfigオブジェクトを作成しています。
x、y、zには、先ほど生成したデータを渡しています。

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

1
2
3
4
5
6
7
8
9
10

fig.update_layout(
    title='Beautiful 3D Surface Plot',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z'
    ),
    width=800,
    height=600
)

fig.update_layoutメソッドを使って、グラフのタイトルと軸ラベルを設定しています。
また、グラフのサイズもwidthとheightで指定しています。

6. カラースケールの設定

1

fig.update_traces(colorscale='Viridis')

fig.update_tracesメソッドを使って、サーフェスプロットのカラースケールを設定しています。
ここでは'Viridis'というカラースケールを使っています。
これは緑から紫へとなめらかに変化する色合いです。

7. グラフの表示

1

fig.show()

最後に、fig.show()メソッドを実行して、作成したグラフを表示します。

このコードでは、NumPyを使ってデータを生成し、Plotlyを使って対話的な3Dサーフェスプロットを描画しています。

サーフェスプロットは、波打つ形状をしており、カラースケールの設定により、美しい色合いで表現されています。

グラフのタイトルと軸ラベルも付けられているので、データの可視化がより分かりやすくなっています。

ハイゼンベルクの不確定性原理

ハイゼンベルクの不確定性原理は、位置 $ ( x ) $と運動量 $ ( p ) $の間の関係を表現します。

具体的には、位置の不確定性 $ ( \Delta x ) $と運動量の不確定性 $ ( \Delta p ) $の積が、プランク定数 $ ( \hbar ) $の半分以上になることを示します。

数学的には、不確定性原理は次のように表されます：

$$
\Delta x \cdot \Delta p \geq \frac{\hbar}{2}
$$

以下は、ハイゼンベルクの不確定性原理をPythonで解いてグラフ化する例です。

ここでは、位置 $ ( x ) $の不確定性 $ ( \Delta x ) $を横軸に、運動量 $ ( p ) $の不確定性 $ ( \Delta p ) $を縦軸に取り、等式$ ( \Delta x \cdot \Delta p = \frac{\hbar}{2} ) $の直線をプロットします。

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

import numpy as np
import matplotlib.pyplot as plt

# プランク定数
h_bar = 1

# 不確定性原理の関数
def uncertainty_principle(delta_x):
    return h_bar / (2 * delta_x)

# 不確定性原理を満たす直線の式
def uncertainty_line(delta_x):
    return h_bar / (2 * delta_x)

# 位置の不確定性の範囲
delta_x_values = np.linspace(0.1, 2, 100)

# 運動量の不確定性を計算
delta_p_values = uncertainty_principle(delta_x_values)

# 不確定性原理を満たす直線を計算
line_values = uncertainty_line(delta_x_values)

# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(delta_x_values, delta_p_values, label='Uncertainty Principle Curve', color='blue')
plt.plot(delta_x_values, line_values, label='Uncertainty Principle Line', linestyle='--', color='red')
plt.xlabel('Position Uncertainty (Delta x)')
plt.ylabel('Momentum Uncertainty (Delta p)')
plt.title('Heisenberg Uncertainty Principle')
plt.grid(True)
plt.legend()
plt.show()

このコードでは、位置の不確定性 $ ( \Delta x ) $を横軸に取り、運動量の不確定性 $ ( \Delta p ) $を縦軸に取ります。

そして、不確定性原理を満たす直線$ ( \Delta x \cdot \Delta p = \frac{\hbar}{2} ) $を赤色の点線でプロットします。

[実行結果]

ソースコード解説

このソースコードは、ハイゼンベルクの不確定性原理を解説するためのPythonプログラムです。

以下に、各部分の詳細を説明します。

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

• import numpy as np: 数値計算を行うためのNumPyライブラリをインポートし、npという別名で使用します。
• import matplotlib.pyplot as plt: グラフのプロットに使用するMatplotlibライブラリのpyplotモジュールをインポートし、pltという別名で使用します。

2. プランク定数の定義:

• プランク定数$ ( \hbar ) $を$ 1 $として定義します。

3. 不確定性原理の関数:

• uncertainty_principle(delta_x): 位置の不確かさ$ ( \Delta x ) $を引数として受け取り、不確定性原理$ ( \Delta x \cdot \Delta p \geq \frac{\hbar}{2} ) $を用いて運動量の不確かさ$ ( \Delta p ) $を計算します。

4. 不確定性原理を満たす直線の式:

• uncertainty_line(delta_x): 位置の不確かさ$ ( \Delta x ) $を引数として受け取り、不確定性原理を満たす直線$ ( \Delta x \cdot \Delta p = \frac{\hbar}{2} ) $の式を表します。

5. 位置の不確定性の範囲の設定:

• delta_x_values = np.linspace(0.1, 2, 100): 位置の不確かさ$ ( \Delta x ) $の範囲を$0.1$から$2$まで$100$等分した値を生成し、delta_x_valuesに代入します。

6. 運動量の不確定性の計算:

• delta_p_values = uncertainty_principle(delta_x_values): 上記で定義した不確定性原理の関数を使用して、位置の不確かさに対応する運動量の不確かさ$ ( \Delta p ) $を計算します。

7. 不確定性原理を満たす直線の計算:

• line_values = uncertainty_line(delta_x_values): 上記で定義した不確定性原理を満たす直線の式を使用して、与えられた位置の不確かさに対応する運動量の不確かさの値を計算します。

8. グラフのプロット:

• Matplotlibを使用して、位置の不確かさ$ ( \Delta x ) $を横軸、運動量の不確かさ$ ( \Delta p ) $を縦軸としてプロットします。
• 不確定性原理を満たす曲線と直線をプロットし、それぞれ青色と赤色で表示します。
• グリッド線を表示し、軸のラベルとタイトルを設定します。
• 凡例を追加して、曲線と直線の説明を表示します。
• 最後に、グラフを表示します。

結果解説

[実行結果]

このグラフは、ハイゼンベルクの不確定性原理を可視化しています。

不確定性原理は、位置 $ ( x ) $と運動量 $ ( p ) $の間の不確かさの関係を表現します。

具体的には、位置の不確かさ $ ( \Delta x ) $と運動量の不確かさ $ ( \Delta p ) $の積が、プランク定数 $ ( \hbar ) $の半分以上になることを示します。

このグラフには以下の要素が含まれています：

1. 青色の曲線 (Uncertainty Principle Curve):

• 横軸に位置の不確かさ $ ( \Delta x )$、縦軸に運動量の不確かさ $ ( \Delta p ) $を取ります。
• ハイゼンベルクの不確定性原理を満たすように、$( \Delta x \cdot \Delta p = \frac{\hbar}{2} ) $の曲線を表します。
• 曲線上の点は、位置と運動量の不確かさの組み合わせを示します。

2. 赤色の点線 (Uncertainty Principle Line):

• 不確定性原理を満たす直線を表します。
• $ ( \Delta x \cdot \Delta p = \frac{\hbar}{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
35
36
37
38
39
40
41

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

# 連立方程式を定義する関数
def equations(vars):
    x, y, z = vars
    eq1 = x**2 + y**2 + z**2 - 9
    eq2 = x**2 - y**2 + 2*z**2 - 4
    eq3 = 3*x**2 + y**2 - 2*z**2 - 12
    return [eq1, eq2, eq3]

# 初期値を設定
x0 = 1
y0 = 1
z0 = 1

# 連立方程式を解く
solution = fsolve(equations, [x0, y0, z0])
print(f"Solution: x={solution[0]:.3f}, y={solution[1]:.3f}, z={solution[2]:.3f}")

# グラフ化
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)

Z1 = np.sqrt(9 - X**2 - Y**2)
Z2 = np.sqrt((4 + Y**2 - X**2) / 2)
Z3 = np.sqrt((12 - 3*X**2 - Y**2) / 2)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z1, alpha=0.5, color='r', label='eq1')
ax.plot_surface(X, Y, Z2, alpha=0.5, color='g', label='eq2')
ax.plot_surface(X, Y, Z3, alpha=0.5, color='b', label='eq3')
ax.scatter(solution[0], solution[1], solution[2], c='k', marker='o', s=50, label='Solution')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

このコードでは、3つの非線形方程式から成る連立方程式を定義しています。

$$
x^2 + y^2 + z^2 - 9 = 0
$$
$$
x^2 - y^2 + 2z^2 - 4 = 0
$$
$$
3x^2 + y^2 - 2z^2 - 12 = 0
$$

scipy.optimize.fsolve関数を使って、この連立方程式の解を数値的に求めています。

次に、meshgridを使って$x-y$平面上の点を生成し、各点における$z$値を3つの方程式から計算しています。
これらの点をplot_surfaceで3次元プロットしています。

最後に、連立方程式の解をscatterでプロットしています。

実行すると、以下のような3次元グラフが表示されます。
赤、緑、青の曲面が3つの方程式を表し、黒い点が連立方程式の解を表しています。

[実行結果]

このようにPythonを使えば、複雑な連立方程式の解を視覚的に確認することができます。

ソースコード解説

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

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

1
2
3
4

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

• NumPyは数値計算ライブラリ。
• MatplotlibはPythonでグラフを描くライブラリ。
• mpl_toolkits.mplot3dは3次元プロットに必要。
• scipy.optimizeからfsolve関数を使って連立方程式を解きます。

2. 連立方程式を定義する関数

1
2
3
4
5
6

def equations(vars):
    x, y, z = vars
    eq1 = x**2 + y**2 + z**2 - 9
    eq2 = x**2 - y**2 + 2*z**2 - 4
    eq3 = 3*x**2 + y**2 - 2*z**2 - 12
    return [eq1, eq2, eq3]

• この関数は3変数$x$, $y$, $z$の非線形連立方程式を定義しています。
• $eq1$, $eq2$, $eq3$がそれぞれの方程式。
• リストとして[eq1, eq2, eq3]を返します。

3. 初期値の設定

1
2
3

x0 = 1
y0 = 1 
z0 = 1

• 連立方程式を解くための初期値$x0$, $y0$, $z0$を設定しています。

4. 連立方程式の解を求める

1
2

solution = fsolve(equations, [x0, y0, z0])
print(f"Solution: x={solution[0]:.3f}, y={solution[1]:.3f}, z={solution[2]:.3f}")

• fsolve関数に先程定義した方程式とベクトル$[x0, y0, z0]$を渡して解を求めます。
• 解の$x$, $y$, $z$値をprintで出力します。

5. グラフ化の準備

1
2
3
4
5
6
7

x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)

Z1 = np.sqrt(9 - X**2 - Y**2)
Z2 = np.sqrt((4 + Y**2 - X**2) / 2)  
Z3 = np.sqrt((12 - 3*X**2 - Y**2) / 2)

• $x$, $y$の範囲を$-3$から$3$まで$100$点に分けてリストを作成。
• meshgridで$X$, $Y$の格子点を作成。
• 3つの方程式からそれぞれ$Z1$, $Z2$, $Z3$を計算。

6. プロットとグラフの設定

1
2
3
4
5
6
7
8
9
10

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z1, alpha=0.5, color='r', label='eq1')
ax.plot_surface(X, Y, Z2, alpha=0.5, color='g', label='eq2')
ax.plot_surface(X, Y, Z3, alpha=0.5, color='b', label='eq3')
ax.scatter(solution[0], solution[1], solution[2], c='k', marker='o', s=50, label='Solution')
ax.set_xlabel('X')
ax.set_ylabel('Y') 
ax.set_zlabel('Z')
plt.show()

• 新しい図とプロット領域を作成。
• plot_surfaceで3つの曲面をプロット。
• scatterでソリューションの点をプロット。
• 軸ラベルを設定。
• グラフを表示。

このコードでは、3変数の非線形連立方程式を定義し、その解を求めています。
その上で、解と3つの曲面を3次元プロットで可視化しています。

曲面の交点がちょうど解になっていることが一目でわかるグラフになっています。

結果解説

[実行結果]

この3次元グラフには以下の内容が表示されています。

1. 3つの曲面 (赤、緑、青)

• これらは3つの連立方程式をそれぞれ表した曲面です。
• 赤い曲面は $x^2 + y^2 + z^2 - 9 = 0$ を満たす点の集合。
• 緑の曲面は $x^2 - y^2 + 2z^2 - 4 = 0$ を満たす点の集合。
• 青の曲面は $3x^2 + y^2 - 2z^2 - 12 = 0$ を満たす点の集合。

2. 黒い点

• この点が連立方程式の解$ (x, y, z) $を表しています。
• 具体的な値は $x=1.826$, $y=-1.030$, $z=1.714$ です。
• この点は3つの曲面が交わる点に位置しています。

3. 座標軸

• $X$軸 (赤): $x$座標を表します。
• $Y$軸 (緑): $y$座標を表します。
• $Z$軸 (青): $z$座標を表します。

4. 曲面とグラフの範囲

• $X$, $Y$軸方向に$-3$から$3$まで。
• 曲面はこの範囲内で描かれています。

5. 透明度(alpha値)

• 曲面には$0.5$の透明度が設定されています。
• これにより、3つの曲面が重なる部分が見やすくなっています。

6. 凡例

• 曲面の色と方程式の対応関係が示されています。
• 黒い点が「Solution」と表示されています。

このように、この3次元グラフは与えられた3つの連立方程式の様子と、その解をわかりやすく可視化しています。

曲面の交点が連立方程式の解になっていることが一目でわかります。

カルマンフィルタ

カルマンフィルタは、状態推定やシステム制御に広く使用される信号処理手法です。
以下では、Pythonでカルマンフィルタを実装して、状態推定の例をグラフ化します。

まずは、NumPyとMatplotlibを使用してカルマンフィルタを実装します。
以下のコードでは、1次元のシステムに対してカルマンフィルタを適用し、真の状態と推定された状態の推移をグラフ化します。

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

import numpy as np
import matplotlib.pyplot as plt

def kalman_filter(initial_state, initial_estimate_error_covariance, process_variance, measurement_variance, measurements):
    # 初期状態の設定
    x_estimate = initial_state
    P_estimate = initial_estimate_error_covariance
    
    # 推定された状態と真の状態の履歴を保存するリスト
    x_history = []
    true_state_history = []
    
    # カルマンフィルタの更新
    for z in measurements:
        # 予測ステップ
        x_predict = x_estimate
        P_predict = P_estimate + process_variance
        
        # 更新ステップ
        K = P_predict / (P_predict + measurement_variance)
        x_estimate = x_predict + K * (z - x_predict)
        P_estimate = (1 - K) * P_predict
        
        # 結果を保存
        x_history.append(x_estimate)
        true_state_history.append(z)
    
    return x_history, true_state_history

# システムの真の状態を生成
np.random.seed(0)
true_state = np.cumsum(np.random.randn(100))  # ランダムウォークモデル

# ノイズのある観測値を生成
measurements = true_state + 0.5 * np.random.randn(100)

# カルマンフィルタのパラメータ設定
initial_state = 0
initial_estimate_error_covariance = 1
process_variance = 1e-5
measurement_variance = 0.5**2

# カルマンフィルタを適用
estimated_states, true_states = kalman_filter(initial_state, initial_estimate_error_covariance, process_variance, measurement_variance, measurements)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(true_states, label='True State', color='blue', linestyle='--')
plt.plot(estimated_states, label='Estimated State', color='red')
plt.title('Kalman Filter Example')
plt.xlabel('Time Step')
plt.ylabel('State')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、ランダムウォークモデルに従うシステムの真の状態を生成し、観測にノイズが加わった状態をシミュレートします。
その後、カルマンフィルタを用いて真の状態を推定し、推定された状態と真の状態の推移をグラフ化します。

カルマンフィルタは、システムの動的な状態を推定するために、システムモデルと観測モデルを組み合わせて効果的に利用されます。

この例では、1次元のシステムに対するカルマンフィルタの動作を示していますが、より複雑なシステムや高次元の状態推定にも応用することができます。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

NumPyは数値計算や配列操作のためのライブラリであり、Matplotlibはグラフ描画のためのライブラリです。

2. カルマンフィルタ関数の定義

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

def kalman_filter(initial_state, initial_estimate_error_covariance, process_variance, measurement_variance, measurements):
    # 初期状態の設定
    x_estimate = initial_state
    P_estimate = initial_estimate_error_covariance
    
    # 推定された状態と真の状態の履歴を保存するリスト
    x_history = []
    true_state_history = []
    
    # カルマンフィルタの更新
    for z in measurements:
        # 予測ステップ
        x_predict = x_estimate
        P_predict = P_estimate + process_variance
        
        # 更新ステップ
        K = P_predict / (P_predict + measurement_variance)
        x_estimate = x_predict + K * (z - x_predict)
        P_estimate = (1 - K) * P_predict
        
        # 結果を保存
        x_history.append(x_estimate)
        true_state_history.append(z)
    
    return x_history, true_state_history

この関数はカルマンフィルタを実装します。
引数として初期状態、初期推定誤差共分散行列、プロセスの分散、観測の分散、観測値を受け取ります。
関数内ではカルマンフィルタのアルゴリズムに基づいて真の状態を推定します。

3. シミュレーション用データの生成

1
2
3

np.random.seed(0)
true_state = np.cumsum(np.random.randn(100))  # ランダムウォークモデル
measurements = true_state + 0.5 * np.random.randn(100)  # ノイズのある観測値

ランダムウォークモデルに従う真の状態 true_state と、観測ノイズを含んだ観測値 measurements を生成します。

4. カルマンフィルタのパラメータ設定

1
2
3
4

initial_state = 0
initial_estimate_error_covariance = 1
process_variance = 1e-5
measurement_variance = 0.5**2

カルマンフィルタの初期状態やパラメータを設定します。

5. カルマンフィルタの適用と結果のプロット

1
2
3
4
5
6
7
8
9
10
11
12

estimated_states, true_states = kalman_filter(initial_state, initial_estimate_error_covariance, process_variance, measurement_variance, measurements)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(true_states, label='True State', color='blue', linestyle='--')
plt.plot(estimated_states, label='Estimated State', color='red')
plt.title('Kalman Filter Example')
plt.xlabel('Time Step')
plt.ylabel('State')
plt.legend()
plt.grid(True)
plt.show()

カルマンフィルタを適用し、真の状態 true_states と推定された状態 estimated_states の時間変化をグラフで表示します。
青の破線が真の状態を示し、赤の実線がカルマンフィルタによって推定された状態を示します。
推定された状態が真の状態に近づくことが期待されます。

このグラフはカルマンフィルタがシステムの状態を効果的に推定できることを示しています。

結果解説

[実行結果]

上記のカルマンフィルタのグラフは、真の状態とカルマンフィルタによって推定された状態の時間変化を示しています。
以下はグラフに表示される内容の詳細です：

True State (青の破線):

• 真の状態を表します。
• ランダムウォークモデルに従って生成されたデータで、時間とともにランダムに変化します。
• ノイズが加わった観測値を元にカルマンフィルタが推定を行います。

Estimated State (赤の実線):

• カルマンフィルタによって推定された状態を表します。
• 真の状態と観測値を基にして、カルマンフィルタが状態を推定します。
• 推定された状態は観測値と真の状態の間を滑らかに補完したものとなります。

Time Step (時間ステップ):

• $x$軸は時間ステップを表します。
• 各ステップで真の状態とその推定値が示されています。

State (状態):

• $y$軸はシステムの状態を表します。
• 真の状態と推定された状態の値が示されています。

このグラフを通じて、カルマンフィルタがシステムの真の状態を観測ノイズを考慮しながら効果的に推定している様子が確認できます。

推定された状態は真の状態に近い軌道を描いており、ノイズの影響を抑えつつシステムの状態を正確に追跡しています。

カルダノの公式（Cardano's formula）

カルダノの公式（Cardano’s formula）は、2次方程式$ ( ax^2 + bx + c = 0 ) $の解の公式です。

この公式は以下のように表されます：

$$
x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}
$$

この公式をPythonで使用して、2次方程式の解を計算し、グラフ化する方法を示します。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

次に、2次方程式の解を計算する関数を定義します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

def solve_quadratic(a, b, c):
    discriminant = b**2 - 4*a*c  # 判別式を計算

    if discriminant > 0:
        # 2つの異なる実数解を持つ場合
        x1 = (-b + np.sqrt(discriminant)) / (2*a)
        x2 = (-b - np.sqrt(discriminant)) / (2*a)
        return x1, x2
    elif discriminant == 0:
        # 重解を持つ場合
        x = -b / (2*a)
        return x, x
    else:
        # 実数解を持たない場合
        return None

この関数 solve_quadratic は、2次方程式$ ( ax^2 + bx + c = 0 ) $の係数$ ( a, b, c ) $を受け取り、解を計算して返します。

解が実数の場合は$ ( (x_1, x_2) ) $の形で返し、重解の場合は$ ( (x, x) ) $の形で返します。

実数解が存在しない場合は None を返します。

次に、計算した解をグラフ化します。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

def plot_quadratic_solution(a, b, c):
    x_values = np.linspace(-10, 10, 400)  # xの範囲を設定
    y_values = a*x_values**2 + b*x_values + c  # 2次方程式のグラフを作成

    plt.figure(figsize=(8, 6))
    plt.plot(x_values, y_values, label=f'{a}x^2 + {b}x + {c} = 0')  # 2次方程式のグラフをプロット

    # 解の計算
    solutions = solve_quadratic(a, b, c)
    if solutions:
        for solution in solutions:
            plt.plot(solution, 0, 'ro')  # 解を赤い点でプロット

    plt.title('Quadratic Equation')
    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid(True)
    plt.axhline(0, color='black', linewidth=0.5)
    plt.axvline(0, color='black', linewidth=0.5)
    plt.show()

この関数 plot_quadratic_solution は、2次方程式のグラフを描画します。
2次方程式の解が存在する場合は、解を赤い点でプロットします。

これらの関数を使って、具体的な2次方程式を解いてグラフ化する例を示します。

1
2
3
4
5
6
7

# 例: 2x^2 - 5x + 2 = 0 の解とグラフ化
a = 2
b = -5
c = 2

# 2次方程式の解を計算し、グラフ化
plot_quadratic_solution(a, b, c)

この例では、2次方程式$ ( 2x^2 - 5x + 2 = 0 ) $の解を計算し、その解を含むグラフを描画しています。

他の2次方程式に対しても同様に関数を使用できます。

グラフ解説

[実行結果]

2次方程式$ ( ax^2 + bx + c = 0 ) $の解とそのグラフを詳しく説明します。

1. グラフの描画:

• プロットされた曲線は、2次方程式$ ( y = ax^2 + bx + c ) $のグラフです。
  この曲線は放物線を表し、その形状は係数$ ( a, b, c ) $の値によって決まります。
• 横軸$ ( x ) $は変数$ ( x ) $の値を示し、縦軸$ ( y ) $はそれに対応する2次方程式の値$ ( y ) $を表します。

2. 解の計算:

• 2次方程式$ ( ax^2 + bx + c = 0 ) $の解は、$ ( x ) $の値で曲線が$ ( y = 0 ) $と交わる点です。
  これらの点が解を表します。
• 解は関数 solve_quadratic(a, b, c) を用いて計算されます。
  解が実数解を持つ場合は$ ( (x_1, x_2) ) $の形で表示され、重解の場合は同じ$ ( x ) $の値が二つの解として示されます。

3. 解のプロット:

• 計算された解は赤い点でグラフ上にプロットされます。
  これらの点は、2次方程式の解を示しています。
• 解の点が曲線上にある場合は、それが方程式の解であることを示します。

4. 曲線の特徴:

• 2次方程式の放物線は、係数$ ( a ) $によって開き方が異なります。
  $ ( a > 0 ) $の場合は上に凸の形状（頂点が最小値）、$ ( a < 0 ) $の場合は下に凸の形状（頂点が最大値）を示します。
• 係数$ ( b ) $は曲線の位置を左右にシフトさせ、係数$ ( c ) $は曲線の位置を上下にシフトさせます。

5. グラフの表示:

• グラフは横軸の範囲や解の位置によって、曲線の形状や解の配置が変わります。
• 横軸の範囲は np.linspace() 関数によって設定され、解が存在する範囲や曲線の描画範囲が決まります。

これらの要素を組み合わせることで、2次方程式の解とそのグラフを理解することができます。

グラフは曲線の形状や解の位置を直感的に示し、数学的な関係を視覚的に理解するのに役立ちます。

レナード-ジョーンズポテンシャル

レナード-ジョーンズポテンシャルは、原子間または分子間の相互作用を表現するための一般的なポテンシャルエネルギーモデルです。

このポテンシャルは、2つの粒子（原子や分子）間の距離に依存しており、分子動力学シミュレーションや物理化学で広く使用されています。

レナード-ジョーンズポテンシャル $ ( V(r) ) $は次の式で表されます：

$$
V(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right]
$$

ここで、各パラメータの意味は以下の通りです：

• $( r ) $: 分子間の距離（原子核間の最短距離からの距離）
• $( \sigma ) $: レナード-ジョーンズポテンシャルがゼロになる距離であり、分子間の最小距離を表す
• $( \epsilon ) $: ポテンシャルエネルギーの深さを決定する定数（吸引力または斥力の強さを示す）
• $( V(r) ) $: $( r ) $の関数として定義されるポテンシャルエネルギー

このポテンシャルエネルギー関数は、$( r ) $に対して典型的な二重極を持ちます。
距離$ ( r ) $が$ ( \sigma ) $以下の範囲では、ポテンシャルエネルギーは正の値を持ち、原子や分子が斥力を受けます。
一方、距離$ ( r ) $が$ ( \sigma ) $を超えると、ポテンシャルエネルギーは負の値を取り、原子や分子は引力を感じます。

レナード-ジョーンズポテンシャルは、分子間の相互作用を表す一般的な近似式であり、原子や分子の結合エネルギー、結合距離、相互作用の性質を理解するために広く用いられます。

物理学や化学の分野で、分子の構造や性質を研究する際に重要な役割を果たします。

ソースコード

以下は、レナード-ジョーンズポテンシャルを解いてグラフ化する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

import numpy as np
import matplotlib.pyplot as plt

def lennard_jones_potential(r, epsilon, sigma):
    """
    Calculate the Lennard-Jones potential energy for a given distance r.
    
    Parameters:
        r (float or numpy.ndarray): Distance between particles.
        epsilon (float): Depth of the potential energy well.
        sigma (float): Distance at which the potential energy is zero.
        
    Returns:
        float or numpy.ndarray: Potential energy at distance r.
    """
    # Calculate the terms in the LJ potential
    term1 = (sigma / r)**12
    term2 = (sigma / r)**6
    
    # LJ potential energy
    potential_energy = 4 * epsilon * (term1 - term2)
    
    return potential_energy

# Define parameters
epsilon = 1.0  # Depth of the potential energy well
sigma = 1.0    # Distance at which potential energy is zero

# Generate distance values
r_values = np.linspace(0.9, 3.0, 100)  # Range of distances to evaluate

# Calculate LJ potential energies
lj_potentials = lennard_jones_potential(r_values, epsilon, sigma)

# Plot the LJ potential
plt.figure(figsize=(8, 6))
plt.plot(r_values, lj_potentials, label='Lennard-Jones Potential')
plt.xlabel('Distance (r)')
plt.ylabel('Potential Energy (V)')
plt.title('Lennard-Jones Potential Energy')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)  # Horizontal line at V=0
plt.axvline(sigma, color='red', linestyle='--', label='Sigma (σ)')  # Vertical line at sigma
plt.legend()
plt.grid(True)
plt.show()

このコードでは、lennard_jones_potential 関数で与えられた距離$ ( r ) $に対するレナード-ジョーンズポテンシャルを計算し、numpy.linspace を使用して一定の範囲内で距離$ ( r ) $を生成します。
その後、計算されたポテンシャルエネルギーをグラフ化しています。

グラフでは、距離$ ( r ) $がレナード-ジョーンズポテンシャルの影響を受ける範囲を示すために、ポテンシャルがゼロになる距離$ ( \sigma ) $が赤色の破線で示されています。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算のためのライブラリであり、配列や数学関数を効率的に扱うことができます。
• matplotlib.pyplotはグラフ描画ライブラリであり、このプログラムではレナード-ジョーンズポテンシャルのグラフ化に使用されます。

2. レナード-ジョーンズポテンシャル関数の定義

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

def lennard_jones_potential(r, epsilon, sigma):
    """
    Calculate the Lennard-Jones potential energy for a given distance r.
    
    Parameters:
        r (float or numpy.ndarray): Distance between particles.
        epsilon (float): Depth of the potential energy well.
        sigma (float): Distance at which the potential energy is zero.
        
    Returns:
        float or numpy.ndarray: Potential energy at distance r.
    """
    # Calculate the terms in the LJ potential
    term1 = (sigma / r)**12
    term2 = (sigma / r)**6
    
    # LJ potential energy
    potential_energy = 4 * epsilon * (term1 - term2)
    
    return potential_energy

• lennard_jones_potential関数は、与えられた距離 r に対するレナード-ジョーンズポテンシャルエネルギーを計算します。
• r は粒子間の距離で、epsilonはポテンシャルエネルギーの深さを表し、sigmaはポテンシャルエネルギーがゼロになる距離を表します。
• term1とterm2はレナード-ジョーンズポテンシャルの計算に使用される項です。
• potential_energyはレナード-ジョーンズポテンシャルの値を計算し、その値を返します。

3. パラメータの定義

1
2

epsilon = 1.0  # Depth of the potential energy well
sigma = 1.0    # Distance at which potential energy is zero

• epsilonとsigmaはレナード-ジョーンズポテンシャルのパラメータであり、それぞれポテンシャルエネルギーの深さとゼロになる距離を表します。

4. 距離の範囲の生成

1

r_values = np.linspace(0.9, 3.0, 100)  # Range of distances to evaluate

• np.linspace関数を使用して、$0.9$から$3.0$までの範囲を等間隔で$100$点に分割した距離の配列 r_values を生成します。

5. レナード-ジョーンズポテンシャルの計算

1

lj_potentials = lennard_jones_potential(r_values, epsilon, sigma)

• 先ほど定義したlennard_jones_potential関数を使用して、各距離 r_values に対するレナード-ジョーンズポテンシャルエネルギーを計算し、lj_potentialsに格納します。

6. プロットの作成と表示

1
2
3
4
5
6
7
8
9
10

plt.figure(figsize=(8, 6))
plt.plot(r_values, lj_potentials, label='Lennard-Jones Potential')
plt.xlabel('Distance (r)')
plt.ylabel('Potential Energy (V)')
plt.title('Lennard-Jones Potential Energy')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)  # Horizontal line at V=0
plt.axvline(sigma, color='red', linestyle='--', label='Sigma (σ)')  # Vertical line at sigma
plt.legend()
plt.grid(True)
plt.show()

• matplotlib.pyplotを使用してグラフを描画します。
• plt.plotで距離 r_values に対するレナード-ジョーンズポテンシャルエネルギーをプロットします。
• 軸ラベルやタイトルを設定し、水平線と垂直線を追加しています。
• plt.show()でグラフを表示します。

このプログラムは、レナード-ジョーンズポテンシャルを計算し、距離に対するポテンシャルエネルギーの振る舞いを視覚化するためのものです。

グラフ解説

[実行結果]

このグラフは、レナード-ジョーンズポテンシャル（Lennard-Jones potential）を示しています。
レナード-ジョーンズポテンシャルは、原子間または分子間のポテンシャルエネルギーを表す一般的なモデルであり、分子動力学シミュレーションや物理化学で広く使用されています。

以下は、グラフの詳細説明です:

X軸（Distance (r)):

分子間の距離 $( r ) $を表します。
この距離は、原子間の最小距離から離れるにつれて増加します。
単位は任意の長さの単位（例えば、アングストロームなど）です。

Y軸（Potential Energy (V)):

レナード-ジョーンズポテンシャル $( V(r) ) $を表します。
このポテンシャルエネルギーは、分子間の距離 $( r )$ に対して計算され、分子間相互作用のエネルギーを示します。
単位はエネルギーの単位（例えば、エレクトロンボルトなど）です。

Lennard-Jones Potential Curve:

グラフ上の曲線は、レナード-ジョーンズポテンシャル $( V(r) ) $の値を$ ( r ) $の範囲で表します。
曲線は$ ( r ) $が大きくなるにつれてゼロに収束し、$( r ) $が小さいと急激に増加します。

Sigma (σ)の破線:

赤色の破線は、レナード-ジョーンズポテンシャルがゼロになる分子間の最小距離 $( \sigma ) $を示します。
この距離において、ポテンシャルエネルギーは最小値を持ちます。

グラフの特徴:

レナード-ジョーンズポテンシャル曲線は、分子間の距離に応じてポテンシャルエネルギーがどのように変化するかを示します。
$( r ) $が$ ( \sigma ) $以下の範囲では斥力が支配的で$、( r ) $が大きくなると引力が支配的になります。
曲線の形状は、分子間の相互作用に関する基本的な特性を示しています。

このグラフは、物理化学や分子動力学シミュレーションなどで、原子や分子間の相互作用を理解するための重要なツールとして使用されます。

プランクの放射則（Planck's radiation law）

プランクの放射則（Planck’s radiation law）は、物体の放射スペクトルを記述するための物理学の法則です。

この法則は、黒体放射の強度が波長と温度にどのように依存するかを示しています。

以下では、プランクの放射則を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 scipy.constants import h, c, k

# プランクの放射則を計算する関数
def planck_radiation(wavelength, temperature):
    numerator = 2 * h * c**2 / wavelength**5
    denominator = np.exp(h * c / (wavelength * k * temperature)) - 1
    return numerator / denominator

# 波長の範囲を定義 (単位: メートル)
wavelengths = np.linspace(1e-9, 3e-6, 1000)  # 1 nm から 3 μm の範囲

# 温度のリストを定義 (単位: ケルビン)
temperatures = [300, 600, 1200]  # 300 K, 600 K, 1200 K の場合を考える

# プランクの放射則を各温度で計算してグラフ化
plt.figure(figsize=(10, 6))
for temp in temperatures:
    intensity = planck_radiation(wavelengths, temp)
    plt.plot(wavelengths * 1e9, intensity, label=f'T = {temp} K')

plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Spectral Intensity (W / m^2 / nm / sr)', fontsize=12)
plt.title('Planck\'s Radiation Law', fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

このPythonコードでは、プランクの放射則を計算するための関数 planck_radiation を定義しました。

この関数は、指定された波長と温度に対する放射強度を計算します。

次に、計算した放射強度を異なる温度についてグラフ化します。

wavelengths は波長の範囲を定義し、temperatures は考慮する温度のリストです。
planck_radiation 関数を使用して各温度における放射スペクトルを計算し、結果をグラフにプロットしています。

[実行結果]

グラフは波長に対する放射強度を示しており、異なる温度に対する放射スペクトルの変化を視覚化しています。
この例は、プランクの放射則をPythonで計算してグラフ化する基本的な方法を示しています。

ソースコード解説

以下は、ソースコードの説明をします。

インポートとライブラリの準備

1
2
3

import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import h, c, k

• numpy パッケージは数値計算に使用されます。
• matplotlib.pyplot パッケージはグラフの作成や可視化に使用されます。
• scipy.constants からはプランク定数 (h)、光速度 (c)、ボルツマン定数 (k) を使用します。

プランクの放射則の関数定義

1
2
3
4

def planck_radiation(wavelength, temperature):
    numerator = 2 * h * c**2 / wavelength**5
    denominator = np.exp(h * c / (wavelength * k * temperature)) - 1
    return numerator / denominator

• planck_radiation 関数は、プランクの放射則に基づいて放射強度を計算します。
• wavelength は波長の配列、temperature は温度です。
• numerator は放射則の分子部分を計算し、denominator は分母部分を計算します。
• 関数は計算結果を返します。

波長の範囲の定義

1

wavelengths = np.linspace(1e-9, 3e-6, 1000)

• np.linspace を使用して、$1 nm $から$ 3 μm $の範囲を$ 1000 $等分した波長の配列を作成します。

温度のリストの定義

1

temperatures = [300, 600, 1200]

• グラフ化する各温度のリストを定義します。
  ここでは$ 300 K$、$600 K$、$1200 K $を考慮します。

グラフの作成とプロット

1
2
3
4
5
6
7
8
9
10
11

plt.figure(figsize=(10, 6))
for temp in temperatures:
    intensity = planck_radiation(wavelengths, temp)
    plt.plot(wavelengths * 1e9, intensity, label=f'T = {temp} K')

plt.xlabel('Wavelength (nm)', fontsize=12)
plt.ylabel('Spectral Intensity (W / m^2 / nm / sr)', fontsize=12)
plt.title('Planck\'s Radiation Law', fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

• plt.figure で新しい図を作成し、温度ごとにプランクの放射則を計算してプロットします。
• plot 関数を使用して波長に対する放射強度をプロットします。波長は$ nm $単位で表示されます。
• 軸ラベルやタイトルを設定し、凡例を表示します。
• 最後に show 関数でグラフを表示します。

このプログラムは、プランクの放射則に基づいて異なる温度での放射スペクトルを視覚化し、温度が高いほど短い波長での放射強度が増加することを示しています。

結果解説

[実行結果]

プランクの放射則を計算し、グラフ化することで表示される内容を詳しく説明します。

1. 横軸 (Wavelength - 波長):

• グラフの横軸は波長を表します。
• 単位はナノメートル ($nm$) で、範囲は$1 nm $から$ 3000 nm$ ($3 μm$) までです。
• 波長が短いほど青紫色に近く、長いほど赤色に近い光を示します。

2. 縦軸 (Spectral Intensity - スペクトル強度):

• グラフの縦軸は放射強度を表します。
• 単位はワット毎平方メートル毎ナノメートル毎立体ラジアン ($W/m^2/nm/sr$) です。
• 各波長における単位立体ラジアンあたりの放射強度が示されます。

3. 曲線 (Intensity vs. Wavelength - 強度対波長):

• 各曲線は異なる温度に対する放射スペクトルを示します。
• 温度が高いほど波長による放射強度のピークが高くなり、波長の変化に対する分布が異なります。
• 低温 ($300 K$) の場合はピークが波長の長い側にあり、高温 ($1200 K$) の場合はピークが波長の短い側にあります。

4. タイトルと凡例:

• グラフのタイトルは “Planck’s Radiation Law” で、プランクの放射則を示しています。
• 各曲線に対する凡例は、対応する温度 ($T$) を示しています。

このグラフは、温度が異なる条件下で物体が放射する光の波長ごとの強度分布を示しています。

温度が上昇すると、光のピークの位置が短い波長側に移動し、放射強度が増加することが観察されます。

これにより、物体の温度が光の放射特性にどのように影響を与えるかが可視化されます。

ベルヌーイ方程式

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

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

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

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

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

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

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

ベルヌーイ方程式を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

# 定数と流体の特性パラメータ
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. 定数と流体の特性パラメータの定義

1
2
3
4

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. グラフのプロット

1
2
3
4
5
6
7
8

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」として表示され、プロットされたデータの内容を説明しています。

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

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