クレスパラボロイドの方程式

クレスパラボロイドの方程式 $ z^2 = x^2 + y^2 $は、以下のように説明できます。

• この方程式は、$z$座標の2乗が、$x$座標と$y$座標の2乗和に等しいことを表しています。

• 左辺の$z^2$は、$z$座標の値を2乗したものです。

• 右辺の$x^2 + y^2$は、$x$座標と$y$座標の値をそれぞれ2乗し、その和を取ったものです。

• つまり、この方程式は、ある点$(x,y,z)$において、$z$座標の2乗値が、$x$,$y$座標の2乗和に一致するような点の集合を表しています。

• $z^2 = x^2 + y^2 $を満たす点は、原点$(0,0,0)$を頂点とする上に凸な放物面上にあります。

• また、この方程式は$z$座標について偶関数なので、$z > 0 $と$ z < 0 $の両領域で成り立ちます。

• したがって、この方程式は$z = 0 $の平面を対称面とする上下2つの放物面の集合を表すことになります。

• この形状は、滑らかで連続的で自己交差のない開放曲面となります。

• このような数学的性質から、クレスパラボロイドは工学や物理学の分野で重要な役割を果たしています。

つまり、クレスパラボロイドの方程式は、単純な代数式ながら、興味深い幾何学的形状を表す方程式なのです。

ソースコード

クレスパラボロイドの方程式は以下の通りです。

$$
z^2 = x^2 + y^2
$$

これをPythonでグラフ化するには、matplotlib.pylabライブラリを使います。

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

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

# グリッドデータの作成
x = np.arange(-2, 2, 0.1)
y = np.arange(-2, 2, 0.1)
X, Y = np.meshgrid(x, y)
Z = np.sqrt(X**2 + Y**2)

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

# 表面プロット
ax.plot_surface(X, Y, Z)
ax.plot_surface(X, Y, -Z)

# 軸ラベルとタイトル
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Cressida Paraboloid')

plt.show()

このコードでは、最初にnumpy関数で座標軸のデータを作成しています。
その後、plot_surfaceでクレスパラボロイドの上半分と下半分の曲面をプロットしています。

実行すると、以下のような3D曲面のグラフが表示されるはずです。

[実行結果]

曲面の中心部は原点で、$X$軸、$Y$軸に沿って放物線の形をしています。
$Z$軸方向には上下に対称な形状になります。

ソースコード解説

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

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

1
2
3

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

• NumPyライブラリをnpという名前で読み込みます。
  これは数値計算に使用します。
• Matplotlibライブラリをpltという名前で読み込みます。
  これは2D/3Dプロットに使用します。
• Matplotlib の3D プロット用モジュールAxes3Dをインポートします。

2. グリッドデータの作成

1
2
3
4

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

• np.arange関数で、$-2$から$2$までの$0.1$刻みの1次元配列$x$と$y$を作成します。
• np.meshgrid関数で、$x$と$y$の全組み合わせからなる2次元グリッド$X$、$Y$を生成します。
• 式$ z^2 = x^2 + y^2 $を満たす$Z$の値を計算し、$Z$配列を作成します。

3. 3Dプロットの準備

1
2

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

• plt.figure()で新しい図(Figure)オブジェクトを作成します。
• fig.add_subplot()で3Dプロット用のAxes3Dオブジェクトaxを作成します。

4. 表面プロット

1
2

ax.plot_surface(X, Y, Z)
ax.plot_surface(X, Y, -Z)

• $X$,$Y$,$Z$データからなるパラメトリック曲面をプロットします。
• 2行目で$-Z$を指定することで、下側の曲面もプロットしています。

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

1
2
3
4

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

• 各軸に$X$,$Y$,$Z$というラベルを設定します。
• プロットにCressida Paraboloidというタイトルを設定します。

6. プロットの表示

1

plt.show()

• この行で作成した図をウィンドウに表示します。

以上の手順で、クレスパラボロイドの3D曲面をPythonでプロットしています。

結果解説

[実行結果]

このグラフは3次元の直交座標系 ($X$,$Y$,$Z$軸)で表されています。
曲面は2つの部分から構成されており、1つ目は上に凸な放物面、2つ目はその放物面を下向きに反転させた面です。

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

1. 原点$(0, 0, 0)$を通る曲面
2. $X$軸、$Y$軸に関して点対称の形状
3. $X=0$平面と$Y=0$平面に沿った交線は放物線
4. $Z>0$の領域は上に凸な放物面
5. $Z<0$の領域は下に凸な放物面
6. $|X|$、$|Y|$が大きくなるほど$Z$の絶対値も大きくなる
7. $Z$軸方向に上下に無限に広がる開曲面

この形状はクレスパラボロイド(Cressida Paraboloid)と呼ばれ、数学的には$z^2 = x^2 + y^2 $の方程式で表されます。

グラフの範囲は$-2 \leqq x \leqq 2$、$-2 \leqq y \leqq 2$、$-2 \leqq z \leqq 2$となっていますが、実際の曲面はこの範囲を超えて無限に広がっています。

曲面は滑らかな放物面なので、光源を適切に設定すれば美しい光沢のある3D描画が可能です。
この形状はさまざまな工学分野で利用されています。

ロリツ環面

パラメトリック方程式を使って複雑な3次元曲線やフラクタル構造を描画することができます。

ここでは、ロリツ環面という美しい曲面をパラメトリック方程式で表現する例を示します。

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

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

def lortiz_surface(u, v, a=1, b=1, c=1):
    x = a * np.cos(u) * np.cosh(v)
    y = b * np.sin(u) * np.cosh(v)
    z = c * np.sinh(v) * np.cos(u/2)
    return x, y, z

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

X, Y, Z = lortiz_surface(U, V)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='plasma')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

このコードの説明:

1. lortiz_surface関数は、ロリツ環面のパラメトリック方程式を実装しています。
   引数uとvはパラメータ、a、b、cはスケーリング係数です。

2. np.linspaceを使って、パラメータuとvの値を生成しています。

3. np.meshgridを使って、uとvの値から格子点を生成しています。

4. lortiz_surface関数を呼び出して、格子点における曲面上の点の座標(X, Y, Z)を計算しています。

5. plot_surface関数を使って、計算した座標から曲面をプロットしています。
   rstrideとcstrideは、プロットするポリゴンの密度を制御するパラメータです。
   cmap='plasma'でカラーマップを設定しています。

実行結果は以下のようになります。

[実行結果]

ロリツ環面は、非常に複雑で美しい形状をしています。
パラメトリック方程式を変更したり、異なる方程式を用いることで、さまざまな3次元曲線やフラクタル構造を生成できます。

パラメトリック方程式は、3次元グラフィックスやコンピューターグラフィックス、数学的モデリングなど、さまざまな分野で利用されています。

結果解説

[実行結果]

ロリツ環面のグラフを詳しく説明します。

このグラフは、ロリツ環面と呼ばれる3次元曲面を表しています。
ロリツ環面は、パラメトリック方程式によって定義される曲面で、以下の3つの方程式で表されます。

$$
x = a * cos(u) * cosh(v)
$$
$$
y = b * sin(u) * cosh(v)
$$
$$
z = c * sinh(v) * cos(u/2)
$$
ここで、$u$ と$v$ はパラメータ、$a$、$b$、$c$ はスケーリング係数です。

グラフ中の曲面は、$u$ と $v$ のパラメータ値を$-π$から$+π$まで変化させながら、上記の3つの方程式から計算された$(x, y, z)$の座標値をプロットしたものです。

曲面の形状は、パラメータ$u$ が変化すると環状に回転し、$v $が変化するとくびれた部分ができる非常に複雑な構造になっています。
曲面全体は、ねじれた環状の形をしています。

カラーマップには”plasma”が使われており、曲面の高さ($z$座標値)に応じて青から赤へと連続的に色が変化するよう設定されています。
青い部分が低く、赤い部分が高くなっています。

背景は黒で、曲面の輪郭が浮き出るようになっており、3次元構造をよりわかりやすく表現できます。

$X$、$Y$、$Z$軸のラベルも付けられており、それぞれの座標軸の向きと大きさを把握しやすくなっています。

このようにパラメトリック方程式を用いることで、複雑で美しい3次元曲面をコンピューターグラフィックスで表現することができます。

シェルピンスキーのガスケット

シェルピンスキーのガスケットは、数学的なフラクタルの一種であり、無限に続く小さな三角形から構成される図形です。

このガスケットはポーランドの数学者カジミェシュ・シェルピンスキーによって考案されました。

シェルピンスキーのガスケットは、次のような特徴を持っています。

1. 階層的な構造:

シェルピンスキーのガスケットは、再帰的なプロセスによって構築されます。
最初に正三角形を描き、その三辺の中点を結んで中央の小さな三角形を作ります。
そして、再帰的にこのプロセスを繰り返し、それぞれの小さな三角形に対して同じ操作を行います。
これにより、ガスケット全体が階層的な構造を持ちます。

2. 自己相似性:

シェルピンスキーのガスケットは、自己相似性を示します。
つまり、ガスケット全体が小さな部分にも似た形状やパターンを持っています。
これは、再帰的なプロセスが同じ操作を繰り返し適用することによって生じます。

3. フラクタル性:

シェルピンスキーのガスケットはフラクタルの一種であり、無限に細かい構造が繰り返されるため、無限の解像度で見るとその複雑さが保たれます。
これにより、ガスケットの外形の長さや面積は有限ですが、その中の構造は無限に続きます。

シェルピンスキーのガスケットは、数学的な美しさと興味深さを持ち、コンピューターグラフィックスやフラクタル幾何学などの分野で広く使用されています。

ソースコード

シェルピンスキーのガスケットをPythonで解くためには、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

import numpy as np
import matplotlib.pyplot as plt

def midpoint(p1, p2):
    return (0.5 * (p1[0] + p2[0]), 0.5 * (p1[1] + p2[1]))

def sierpinski(points, depth):
    colors = ['r', 'g', 'b']
    for i in range(3):
        plt.plot([points[i][0], points[(i+1)%3][0]], [points[i][1], points[(i+1)%3][1]], color=colors[i])
    if depth > 0:
        sierpinski([points[0], midpoint(points[0], points[1]), midpoint(points[0], points[2])], depth-1)
        sierpinski([points[1], midpoint(points[0], points[1]), midpoint(points[1], points[2])], depth-1)
        sierpinski([points[2], midpoint(points[2], points[1]), midpoint(points[0], points[2])], depth-1)

def main():
    points = [(0, 0), (1, 0), (0.5, np.sqrt(3)/2)]
    plt.figure(figsize=(6, 6))
    sierpinski(points, 5)  # 深さ5のシェルピンスキーのガスケットを生成
    plt.axis('equal')
    plt.axis('off')
    plt.show()

if __name__ == "__main__":
    main()

このコードでは、まず基本となる三角形の頂点を定義し、その三角形の各辺の中点を計算して新たな三角形を作り出す再帰関数を定義しています。

そして、再帰関数を使用して指定した深さまでシェルピンスキーのガスケットを描画します。

[実行結果]

ソースコード解説

このソースコードは、Pythonとmatplotlibを使用してシェルピンスキーのガスケットを生成し、グラフ化するものです。
以下は、各部分の詳細な説明です。

1. midpoint 関数:

1
2

def midpoint(p1, p2):
   return (0.5 * (p1[0] + p2[0]), 0.5 * (p1[1] + p2[1]))

この関数は、二つの点 p1 と p2 の中点を計算します。
それぞれの点は (x, y) の座標を持つタプルで表されます。
中点の座標は、各次元の座標の平均値を取ることで計算されます。

2. sierpinski 関数:

1
2
3
4
5
6
7
8

def sierpinski(points, depth):
   colors = ['r', 'g', 'b']
   for i in range(3):
       plt.plot([points[i][0], points[(i+1)%3][0]], [points[i][1], points[(i+1)%3][1]], color=colors[i])
   if depth > 0:
       sierpinski([points[0], midpoint(points[0], points[1]), midpoint(points[0], points[2])], depth-1)
       sierpinski([points[1], midpoint(points[0], points[1]), midpoint(points[1], points[2])], depth-1)
       sierpinski([points[2], midpoint(points[2], points[1]), midpoint(points[0], points[2])], depth-1)

この関数は、与えられた三角形の頂点を受け取り、指定された深さまでシェルピンスキーのガスケットを描画します。
各三角形の辺は異なる色で描かれます。
再帰的なアプローチが使用され、各三角形はその頂点の中点を使って三つの小さな三角形に分割されます。

3. main 関数:

1
2
3
4
5
6
7

def main():
   points = [(0, 0), (1, 0), (0.5, np.sqrt(3)/2)]
   plt.figure(figsize=(6, 6))
   sierpinski(points, 5)  # 深さ5のシェルピンスキーのガスケットを生成
   plt.axis('equal')
   plt.axis('off')
   plt.show()

この関数は、プログラムのエントリーポイントとして機能します。
まず、基本となる三角形の頂点が定義され、sierpinski 関数が呼び出されて指定された深さまでのシェルピンスキーのガスケットを生成します。
最後に、グラフが表示されます。

4. if __name__ == "__main__"::

これは、Pythonのスクリプトが直接実行された場合に main 関数が呼び出されることを確認するための条件です。
これにより、このスクリプトが他のスクリプトからインポートされた場合には main 関数が呼び出されないようになります。

このソースコードは、シェルピンスキーのガスケットを再帰的な手法を使用して生成し、matplotlibを使用してグラフ化することで、シェルピンスキーの美しい幾何学的な図形を視覚化しています。

グラフ解説

[実行結果]

このプログラムは、シェルピンスキーのガスケットを生成し、matplotlibを使用してグラフ化します。
以下は、プログラムが生成する内容の詳細な説明です。

1. 三角形の定義:
   最初に、シェルピンスキーのガスケットを構成する基本となる三角形の頂点が定義されます。
   この例では、座標が (0, 0)、(1, 0)、(0.5, √3/2) の三つの頂点を持つ正三角形が使用されます。

2. 中点の計算:
   次に、与えられた三角形の各辺の中点を計算する midpoint 関数が定義されます。
   この関数は、二つの頂点の座標を受け取り、それらの中点の座標を返します。

3. 再帰的な描画:
   sierpinski 関数は、与えられた三角形と深さに基づいて、シェルピンスキーのガスケットを再帰的に描画します。
   具体的には、与えられた三角形の各辺の中点を使用して新しい三角形を作成し、それぞれの新しい三角形に対して再帰的に同じ処理を行います。
   深さが$0$に達するまでこの処理を繰り返します。

4. グラフの描画:
   main 関数では、matplotlibを使用してグラフを描画します。
   sierpinski 関数を呼び出してシェルピンスキーのガスケットを描画し、plt.axis('equal')とplt.axis('off')を使用して、図形のアスペクト比を保持し、軸を非表示にします。

このプログラムは、再帰的なアプローチを使用して、シェルピンスキーのガスケットを生成し、美しい幾何学的な図形を作成します。

コッホ雪片（Koch snowflake）

コッホ雪片（Koch snowflake）は、フラクタル幾何学の代表的な図形の1つです。
スウェーデンの数学者Helge von Kochによって提案されました。
コッホ雪片は、次のような手順で構築されます。

1. 始点の設定:

最初に、正三角形（すなわち、$3$つの等しい長さの辺と$3$つの$60$度の内角を持つ図形）を描きます。

2. 分割:

各辺を$3$等分し、中央の$1/3$地点に新しい頂点を追加します。
これにより、正三角形の各辺が4つの小さな正三角形に分割されます。

3. 内側の辺の追加:

新しい頂点を中心として、外側の辺の中心に向かって、正三角形の$1/3$の長さの辺を追加します。
これにより、外側の辺が内側に折り返されます。

4. 再帰:

生成された小さな三角形に対して、同じ手順を再帰的に適用します。
各ステップで、図形はより細かい詳細を持つより複雑な形状に進化します。

5. 極限に向かって:

無限に続くこのプロセスにより、図形の周囲は収束し、コッホ雪片と呼ばれる図形が形成されます。

コッホ雪片は、境界の長さが無限になるにつれて、形状が複雑になり、自己相似性（自分自身の一部が全体と似ている性質）を持つフラクタルとして知られています。

このようにして、コッホ雪片は独特の幾何学的美しさと複雑さを持ち、数学や美術など様々な分野で興味深い研究対象として扱われています。

ソースコード

コッホ雪片（Koch snowflake）を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
38

import numpy as np
import matplotlib.pyplot as plt

def koch_snowflake(order, length):
    if order == 0:
        # 基本の三角形の頂点座標
        return np.array([[0, 0], [length/2, length*np.sqrt(3)/2], [length, 0]])
    else:
        # 前のオーダーの座標
        points = koch_snowflake(order-1, length)
        new_points = [points[0]]
        for i in range(len(points)-1):
            # 現在の辺を3等分する点の座標を計算
            p1 = points[i]
            p2 = points[i+1]
            dx = p2[0] - p1[0]
            dy = p2[1] - p1[1]
            new_points.append(p1 + np.array([dx/3, dy/3]))
            new_points.append(p1 + np.array([dx/2 + np.sqrt(3)*dy/6, dy/2 - np.sqrt(3)*dx/6]))
            new_points.append(p1 + np.array([2*dx/3, 2*dy/3]))
            new_points.append(p2)
        return np.vstack(new_points)

def plot_koch_snowflake(order, length):
    snowflake = koch_snowflake(order, length)
    plt.figure(figsize=(8, 8))
    plt.title("Koch Snowflake (order = {})".format(order))
    plt.plot(snowflake[:, 0], snowflake[:, 1], 'b-')
    plt.axis('equal')
    plt.axis('off')
    plt.show()

# コッホ雪片のオーダーと辺の長さ
order = 4
length = 300

# コッホ雪片のプロット
plot_koch_snowflake(order, length)

このコードは、指定されたオーダーと辺の長さに基づいてコッホ雪片を生成し、matplotlibを使用してグラフ化します。

[実行結果]

ソースコード解説

以下は、コードの詳細な説明です。

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは、数値計算を行うためのPythonライブラリです。
• matplotlib.pyplotは、グラフ描画機能を提供するMatplotlibライブラリのサブモジュールです。

2. コッホ雪片を生成する関数の定義:

1
2

def koch_snowflake(order, length):
   ...

• orderとlengthを引数として受け取り、指定されたオーダーと辺の長さに基づいてコッホ雪片の頂点座標を計算します。
• 再帰的なアプローチを使用して、コッホ雪片の各オーダーの頂点座標を計算します。

3. コッホ雪片をプロットする関数の定義:

1
2

def plot_koch_snowflake(order, length):
   ...

• orderとlengthを引数として受け取り、指定されたオーダーと辺の長さで生成されたコッホ雪片をプロットします。
• koch_snowflake関数を呼び出してコッホ雪片の頂点座標を取得し、それをmatplotlibを使用して描画します。

4. コッホ雪片のオーダーと辺の長さの指定:

1
2

order = 4
length = 300

• order変数には、コッホ雪片のオーダーを指定します。
  オーダーが大きいほど、より詳細なコッホ雪片が生成されます。
• length変数には、コッホ雪片の辺の長さを指定します。

5. コッホ雪片のプロット:

1

plot_koch_snowflake(order, length)

• plot_koch_snowflake関数を呼び出して、指定されたオーダーと辺の長さでコッホ雪片をプロットします。

これらの手順によって、指定されたオーダーと辺の長さでコッホ雪片が生成され、matplotlibを使用してグラフ化されます。

グラフ解説

[実行結果]

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

1. コッホ雪片の生成:

koch_snowflake関数は、指定されたオーダー（order）と辺の長さ（length）に基づいて、コッホ雪片の頂点の座標を計算します。
コッホ雪片は再帰的な方法で生成され、与えられたオーダーに応じて詳細な構造が形成されます。

2. グラフの描画:

plot_koch_snowflake関数は、生成されたコッホ雪片をmatplotlibを使用してプロットします。
プロットされたコッホ雪片は青色の線で表され、その形状はコッホ曲線の特徴的な三角形の形状を持ちます。
グラフのタイトルは、「Koch Snowflake (order = 指定されたオーダー)」となります。

3. グラフの設定:

plt.figure関数によって、図のサイズが設定されます。
plt.title関数によって、グラフにタイトルが追加されます。
plt.plot関数は、コッホ雪片の頂点座標を線でつなぎ、コッホ雪片を描画します。
plt.axis('equal')によって、$x$軸と$y$軸の比率が等しくなり、コッホ雪片の歪みがなくなります。
plt.axis('off')は、$x$軸と$y$軸の目盛りや枠線を非表示にします。

4. 表示される内容:

グラフに表示される内容は、生成されたコッホ雪片の形状です。
コッホ雪片は、フラクタル図形の一種であり、無限に複雑な形状を持つ三角形です。
このグラフは、コッホ雪片のオーダーに応じた細かい構造を視覚的に示しています。

バーニー龍（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$軸はシステムの状態を表します。
• 真の状態と推定された状態の値が示されています。

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

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