トロイダル座標系

トロイダル座標系は、3次元空間の特殊な曲面を記述するための座標系です。

その主な特徴は以下の通りです。

1. 渦巻き状の形状をしています。
   つまり、円柱状の中空の管が渦を巻いているような形状をしています。

2. 3つの座標変数 $(u, v, r)$で表されます。

   • $u$: 円周方向の角度($0$から$2π$まで)
   • $v$: 渦巻き方向の角度($0$から$2π$まで)
   • $r$: 渦巻き半径

3. 中心軸から一定の距離$r$にある円環状の曲線が、$v$角度分ずつ回転しながら螺旋状に巻きついた形をしています。

4. 電磁気学や流体力学などで、渦巻き状の磁場やフロー場を記述するのに使われます。

5. 座標変数$(u, v, r)$から直交座標$(x, y, z)$への変換式が定義されています。

6. トロイダル座標系はトーラス(ドーナツ型)の形状の物体を自然に記述できる座標系です。

トロイダル座標系は渦巻き状の曲面や中空の管状の物体の形状を簡潔に表現するのに適した特殊な曲線座標系なのです。

電磁気学や流体力学分野で重要な役割を果たしています。

トロイダル座標系を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
from mpl_toolkits.mplot3d import Axes3D

# トロイダル座標系のパラメータ範囲
c = 2  # トロイダル半径
a = 1  # 渦巻き半径
u_range = np.linspace(0, 2 * np.pi, 100)
v_range = np.linspace(0, 2 * np.pi, 100)

# トロイダル座標系の変換式
U, V = np.meshgrid(u_range, v_range)
X = (c + a * np.cos(V)) * np.cos(U)
Y = (c + a * np.cos(V)) * np.sin(U)
Z = a * np.sin(V)

# 3D プロットの設定
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, cmap='viridis', edgecolor='none')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Toroidal Coordinate System')

plt.show()

このコードでは、以下の手順で作業を行っています。

1. numpy と matplotlib をインポートし、3D プロットに必要な Axes3D も読み込みます。
2. トロイダル座標系のパラメータ c (トロイダル半径)、a (渦巻き半径)、および u、v の範囲を設定します。
3. np.meshgrid を使って、u、v の値の組み合わせを生成します。
4. トロイダル座標系の変換式を使って、各 (u, v) の組み合わせに対応する (x, y, z) の値を計算します。
5. matplotlib の plot_surface 関数を使って、計算した (x, y, z) の値から3D 曲面をプロットします。
6. 軸ラベルとタイトルを設定します。
7. plt.show() で最終的なプロットを表示します。

実行すると、次のようなトロイダル座標系の3D プロットが表示されるはずです。

[実行結果]

このプロットでは、トロイダル半径 c = 2、渦巻き半径 a = 1 の値を使用しています。

パラメータ c と a の値を変更すれば、異なる形状のトロイダルを描くことができます。

結果解説

[実行結果]

このグラフは、トロイダル座標系を3次元空間上に可視化したものです。
トロイダル座標系は、渦巻き状の曲面を表すのに適した座標系です。

グラフ上に表示されている曲面は、次の特徴を持っています:

1. 全体的な形状は、ドーナツ状のトーラス(環状)になっています。
   つまり、中心に空洞を持つ管状の構造をしています。

2. 曲面は渦を巻いた形状をしており、螺旋状に巻きついています。
   この螺旋の巻き具合は、トロイダル座標の$v$の値によって決まります。

3. 曲面の半径は一定ではありません。
   曲面の内側半径と外側半径が異なり、これがトロイダル座標の$r$の値に対応します。
   内側半径をトロイダル半径、外側半径を渦巻き半径と呼びます。

4. 曲面は滑らかに連続しており、自己交差することはありません。

5. 曲面の色は、viridisというカラーマップを使って割り当てられています。
   色の濃淡は、曲面の高さ($z$座標値)を表しています。

6. 軸は直交する$x$、$y$、$z$座標で構成され、それぞれにラベルが付けられています。

このようにグラフは、トロイダル座標系の本質的な性質である渦巻き状の管状構造をリアルに可視化しています。

電磁気学や流体力学の分野で重要な役割を果たすこのような特異な曲面を、直感的に理解することができます。

ボーアの周期律

ボーアの周期律は、ボーアの原子模型に基づいて、原子の電子配置と元素の性質の周期的変化を説明する法則です。

主な内容は以下の通りです。

1. 電子は特定の軌道上を回っており、その軌道半径と運動エネルギーは主量子数$n$によって決まる。

2. 電子は各軌道に$2n^2$個まで収容できる。
   $n=1$のK殻には$2$電子、$n=2$のL殻には$8$電子が入る。

3. 原子内の電子は最小のエネルギー準位から順に殻や軌道に入っていく。

4. 同じ主量子数の軌道では、エネルギー準位が低い順に$ s$, $p$, $d$, $f$ と入っていく。

5. 原子番号が増えるにつれ、電子配置が規則的に変化する。
   同じ電子配置を持つ元素は、化学的性質が似通っている。

6. 各周期では、アルカリ金属からハロゲンまで電子配置が変化し、性質も周期的に変化する。

7. 同じ族の元素は外側電子の数が同じなので、化学的性質が似ている。

このように、ボーアの周期律は電子配置に基づいて元素の性質を周期的に整理し、周期表の作成につながりました。

しかし完全ではなく、後にさらに発展した量子力学的な考え方に置き換えられています。

ソースコード

ボーアの原子模型に基づく原子のエネルギー準位は、次の式で与えられます：

$$
E_n = -\frac{k \cdot Z^2 \cdot R_H}{n^2}
$$

ここで、$(E_n) $は第$ (n) $準位のエネルギー、$(k) $はクーロン定数、$(Z) $は原子番号、$(R_H) $はリュードベリ定数、$(n) $は主量子数です。

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 bohr_energy(n, Z=1):
    k = 8.9875517923e9  # クーロン定数 [N m^2 / C^2]
    R_H = 2.1798741e-18  # リュードベリ定数 [J]
    return -(k * Z**2 * R_H) / n**2

# 原子番号
Z = 1  # 水素原子の場合

# 主量子数の範囲を定義
n_values = np.arange(1, 11)

# エネルギー準位を計算
energy_levels = [bohr_energy(n, Z) for n in n_values]

# グラフをプロット
plt.figure(figsize=(8, 6))
plt.plot(n_values, energy_levels, marker='o', linestyle='-')
plt.title('Bohr Model Energy Levels')
plt.xlabel('Principal Quantum Number (n)')
plt.ylabel('Energy (J)')
plt.grid(True)
plt.show()

このコードでは、主量子数$ (n) $ の値に対するエネルギー準位を計算し、それをグラフ化しています。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• NumPyはPythonで数値計算を行うためのライブラリです。
• MatplotlibはPythonでグラフ描画を行うためのライブラリです。
  このスクリプトでは、Matplotlibのpyplotモジュールを使用しています。

2. ボーアのエネルギー計算関数:

1
2
3
4

def bohr_energy(n, Z=1):
   k = 8.9875517923e9  # クーロン定数 [N m^2 / C^2]
   R_H = 2.1798741e-18  # リュードベリ定数 [J]
   return -(k * Z**2 * R_H) / n**2

• bohr_energy 関数は、ボーアの原子模型に基づいて与えられた主量子数$ (n) $に対する原子のエネルギー準位を計算します。
  デフォルトでは水素原子の場合を考慮しており、原子番号$ (Z) $が$1$となっています。

3. 原子番号の設定:

1

Z = 1  # 水素原子の場合

• Z は原子の原子番号を表し、ここでは水素原子の場合を考慮しています。

4. 主量子数の範囲の定義:

1

n_values = np.arange(1, 11)

• np.arange 関数を使って、主量子数$ (n) $の値を$1$から$10$までの範囲で定義しています。

5. エネルギー準位の計算:

1

energy_levels = [bohr_energy(n, Z) for n in n_values]

• bohr_energy 関数を使って、与えられた主量子数$ (n) $に対するエネルギー準位を計算し、リスト energy_levels に格納しています。

6. グラフのプロット:

1
2
3
4
5
6
7

plt.figure(figsize=(8, 6))
plt.plot(n_values, energy_levels, marker='o', linestyle='-')
plt.title('Bohr Model Energy Levels')
plt.xlabel('Principal Quantum Number (n)')
plt.ylabel('Energy (J)')
plt.grid(True)
plt.show()

• Matplotlibを使ってグラフを作成し、エネルギー準位を表示しています。
• plt.figure(figsize=(8, 6)) でグラフのサイズを設定します。
• plt.plot で n_values に対する energy_levels をプロットし、マーカー ‘o’ を使用してデータ点を示し、直線でデータ点をつなぎます。
• グラフのタイトル、軸ラベルを設定し、グリッドを表示しています。
• plt.show() でグラフを表示します。

結果解説

[実行結果]

このグラフは、ボーアの原子模型に基づいて計算された原子のエネルギー準位を示しています。

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

1. x軸:

• x軸は主量子数$ (n) $を表しています。
  主量子数は、電子が原子の中心から離れた軌道またはエネルギー準位を示します。
  主量子数が大きくなるにつれて、電子のエネルギー準位が上昇します。

2. y軸:

• y軸はエネルギーを表しています。
  エネルギーはボーアの原子模型によって予測された値であり、主量子数$ (n) $に依存します。
  負の値を持つことに注意してください。
  これは、電子が原子核に束縛されているため、系全体のエネルギーが負の値となるためです。

3. グラフのプロット:

• グラフには、主量子数$ (n) $に対するエネルギー準位がプロットされています。
  各点は、特定の主量子数$ (n) $に対応するエネルギーを示しています。
• マーカー ‘o’ が使用されていますが、これはデータ点を示すための円を意味します。
• 線はデータ点をつなぐための直線を表し、データの間の関係を視覚化します。

4. タイトル:

• グラフのタイトルは “Bohr Model Energy Levels” となっており、ボーアの原子模型に基づくエネルギー準位を示していることを示しています。

5. 軸ラベル:

• x軸は “Principal Quantum Number (n)” とラベル付けされており、主量子数$ (n) $を示します。
• y軸は “Energy (J)” とラベル付けされており、エネルギーをジュール単位で示しています。

6. グリッド:

• グラフには背景にグリッドが表示されており、データ点の位置を見やすくし、目盛り線に沿ってデータを読み取るのを助けます。

このグラフは、主量子数$ (n) $とエネルギーの関係を視覚化し、ボーアの原子模型に基づく原子のエネルギー準位の理解を支援します。

クラーク-ラムダ分布（Clark-Lambda Distribution）

クラーク-ラムダ分布（Clark-Lambda Distribution）は、確率分布の一種であり、リスク分析や信号処理などの分野で使用されます。

この分布の確率密度関数（PDF）は次のように定義されます:

$$
f(x; \lambda) = \frac{2}{\pi} \frac{\lambda}{\lambda^2 + x^2}
$$

ここで、$(x)$は確率変数であり、$(\lambda)$は分布のパラメータです。

以下は、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

import numpy as np
import matplotlib.pyplot as plt

def clark_lambda_distribution(x, lam):
    return (2 / np.pi) * (lam / (lam**2 + x**2))

# パラメータの設定
lam = 1.0  # λの値

# xの範囲を定義
x_values = np.linspace(-10, 10, 1000)

# 確率密度関数（PDF）を計算
pdf_values = clark_lambda_distribution(x_values, lam)

# グラフをプロット
plt.figure(figsize=(8, 6))
plt.plot(x_values, pdf_values, label=f'λ={lam}')
plt.title('Clark-Lambda Distribution')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.show()

このコードは、$-10$から$10$までの範囲でClark-Lambda分布の確率密度関数を計算し、グラフ化します。

[実行結果]

Lambdaパラメータ（λ）を変更することで、異なる分布を観察できます。

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• NumPyはPythonで数値計算を行うためのライブラリです。
• MatplotlibはPythonでグラフ描画を行うためのライブラリです。
  このスクリプトでは、Matplotlibのpyplotモジュールを使用しています。

2. Clark-Lambda分布の確率密度関数を定義する関数:

1
2

def clark_lambda_distribution(x, lam):
   return (2 / np.pi) * (lam / (lam**2 + x**2))

• clark_lambda_distribution 関数は、Clark-Lambda分布の確率密度関数を計算します。
  引数として x（確率変数）と lam（分布のパラメータ）を受け取ります。

3. パラメータの設定:

1

lam = 1.0  # λの値

• lam はClark-Lambda分布のパラメータ$ ( \lambda ) $の値を示しています。

4. xの範囲を定義:

1

x_values = np.linspace(-10, 10, 1000)

• np.linspace 関数を使って、$-10$から$10$の範囲を等間隔で$1000$点に分割した x の値を生成します。
  これにより、確率密度関数のグラフが滑らかになります。

5. 確率密度関数（PDF）を計算:

1

pdf_values = clark_lambda_distribution(x_values, lam)

• 先ほど定義した clark_lambda_distribution 関数を使って、与えられた x の値とパラメータ$ ( \lambda ) $に基づいてClark-Lambda分布の確率密度関数の値を計算します。

6. グラフをプロット:

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(x_values, pdf_values, label=f'λ={lam}')
plt.title('Clark-Lambda Distribution')
plt.xlabel('x')
plt.ylabel('Probability Density')
plt.legend()
plt.grid(True)
plt.show()

• Matplotlibを使ってグラフを作成し、Clark-Lambda分布の確率密度関数を表示します。
• plt.figure(figsize=(8, 6)) でグラフのサイズを設定します。
• plt.plot で x_values と pdf_values をプロットし、ラベルにはパラメータ$ ( \lambda ) $の値が表示されます。
• plt.title, plt.xlabel, plt.ylabel を使ってタイトルや軸のラベルを設定します。
• plt.legend() で凡例を表示します。
• plt.grid(True) でグリッドを有効にし、視覚的な参照を改善します。
• plt.show() でグラフを表示します。

結果解説

[実行結果]

グラフに表示される内容の詳細な説明です。

1. x軸:

• x軸は確率変数 $(x)$を表します。
  この確率変数はClark-Lambda分布におけるランダム変数であり、確率密度関数が評価される値の範囲が表示されます。
  ここでは、$-10$から$10$までの範囲が選択されています。

2. y軸:

• y軸は確率密度関数（PDF）の値を表します。

Clark-Lambda分布の確率密度関数は、$(x)$に関する確率密度を示し、y軸の値はその密度を表します。

3. ラベル:

• グラフにはClark-Lambda分布の確率密度関数が、指定されたLambdaパラメータ（λ）でプロットされています。
  ラベルにはこのパラメータの値が表示され、グラフ上で異なるパラメータの効果を直感的に理解するのに役立ちます。

4. 凡例:

• 複数のLambdaパラメータを使用する場合、それぞれの分布が区別できるように凡例が表示されます。
  この例では、唯一のLambdaパラメータが使用されていますが、複数の分布を比較する場合には凡例が重要です。

5. タイトル:

• グラフのタイトルには、プロットされている分布の種類が示されます。
  この場合、「Clark-Lambda Distribution」と表示されます。

6. グリッド:

• グラフには背景にグリッドが表示されます。
  これにより、目盛り線に沿ってデータを読み取るのが容易になります。

このグラフは、Clark-Lambda分布の確率密度関数を直感的に理解しやすくするためのものであり、異なるLambdaパラメータを使用することで、分布の形状がどのように変化するかを視覚的に把握できます。

オイラー-ラグランジュ方程式

オイラー-ラグランジュ方程式は、力学系の運動を記述する非線形常微分方程式です。

以下は、Pythonコードの例です。

このコードでは、scipy.integrateモジュールを使って常微分方程式を数値的に解き、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
from scipy.integrate import odeint

# オイラー-ラグランジュ方程式
def euler_lagrange(y, t, m, k, g, l):
    theta, theta_dot = y
    dtheta_dt = theta_dot
    dtheta_dot_dt = -(g/l)*np.sin(theta) - (k/(m*l**2))*theta
    return [dtheta_dt, dtheta_dot_dt]

# 初期条件
theta0 = np.pi/4
theta_dot0 = 0
y0 = [theta0, theta_dot0]

# パラメータ
m = 1.0  # 質量
k = 1.0  # ばね定数
g = 9.81 # 重力加速度
l = 1.0  # 長さ

# 時間範囲
t = np.linspace(0, 20, 1000)

# 常微分方程式を解く
sol = odeint(euler_lagrange, y0, t, args=(m, k, g, l))

# グラフ描画
plt.figure(figsize=(10, 6))
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.plot(t, sol[:, 1], 'g', label='theta_dot(t)')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.title('Euler-Lagrange Equation')
plt.show()

このコードを実行すると、時間tに対する変数thetaとtheta_dotのグラフが描画されます。
初期条件やパラメータを変更すると、異なる振る舞いが観察できます。

オイラー-ラグランジュ方程式は、非線形性があるため、解の振る舞いが複雑になる場合があります。
このコードでは、単振り子の運動を例として示しましたが、さまざまな力学系への応用が可能です。

必要に応じて、グラフの範囲やラベル、タイトルなどを調整してください。
また、より詳細な解析や可視化が必要な場合は、適宜コードを拡張できます。

オイラー-ラグランジュ方程式は、力学系の運動を記述する基本的な方程式であり、物理学や工学分野で広く使われています。

[実行結果]

ソースコード解説

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

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

1
2
3

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

• NumPyライブラリをnpという名前で、数値計算に使用します。
• matplotlibライブラリをpltという名前で、グラフ描画に使用します。
• scipy.integrateモジュールからodeinitを、常微分方程式の数値解を求めるために使用します。

2. オイラー-ラグランジュ方程式の定義

1
2
3
4
5

def euler_lagrange(y, t, m, k, g, l):
    theta, theta_dot = y
    dtheta_dt = theta_dot
    dtheta_dot_dt = -(g/l)*np.sin(theta) - (k/(m*l**2))*theta
    return [dtheta_dt, dtheta_dot_dt]

• この関数では、オイラー-ラグランジュ方程式を定義しています。
• 入力引数yには、角度thetaと角速度theta_dotの値が渡されます。
• tは時間、m、k、g、lはそれぞれ質量、バネ定数、重力加速度、長さです。
• 方程式の右辺は、dtheta_dtとdtheta_dot_dtで表されています。
• 関数の出力は、dtheta_dtとdtheta_dot_dtのリストです。

3. 初期条件と系のパラメータの設定

1
2
3
4
5
6
7
8

theta0 = np.pi/4
theta_dot0 = 0
y0 = [theta0, theta_dot0]

m = 1.0  # 質量
k = 1.0  # ばね定数
g = 9.81 # 重力加速度
l = 1.0  # 長さ

• 初期条件として、theta0とtheta_dot0をそれぞれπ/4(約45度)と0(静止状態)に設定しています。
• y0は初期値のリストです。
• m、k、g、lには、質量、バネ定数、重力加速度、長さの値を代入しています。

4. 時間範囲の設定と常微分方程式の解法

1
2

t = np.linspace(0, 20, 1000)
sol = odeint(euler_lagrange, y0, t, args=(m, k, g, l))

• tには、$0$から$20$秒までの$1000$個の時間点が等間隔で設定されます。
• odeiintを使って、オイラー-ラグランジュ方程式の数値解sol が計算されます。
• 関数euler_lagrangeと初期値y0、時間範囲tが引数として渡されます。
• さらに、m、k、g、lの値もargsで渡されます。

5. グラフの描画

1
2
3
4
5
6
7
8

plt.figure(figsize=(10, 6))
plt.plot(t, sol[:, 0], 'b', label='theta(t)')
plt.plot(t, sol[:, 1], 'g', label='theta_dot(t)')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.title('Euler-Lagrange Equation')
plt.show()

• 新しい図ウィンドウが、サイズ$ (10, 6) $インチで作成されます。
• sol[:, 0]は、thetaの解の時間変化を表す配列です。
• sol[:, 1]は、theta_dotの解の時間変化を表す配列です。
• それぞれ青と緑の線で描画されます。
• 凡例が配置され、ラベルが付けられます。
• 横軸のラベルは’Time’、縦軸のラベルは’Solution’です。
• グラフのタイトルは’Euler-Lagrange Equation’に設定されています。
• plt.show()で、グラフが表示されます。

このコードでは、オイラー-ラグランジュ方程式を使って単振り子の運動を記述し、その解の時間変化をグラフで可視化しています。

初期条件やパラメータを変更することで、異なる振る舞いを観察できます。

結果解説

[実行結果]

グラフには、2つの曲線が表示されています。

1. 青色の曲線: theta(t)

この曲線は、時間tに対する角度theta(単位はラジアン) の変化を表しています。
theta(t)は、オイラー-ラグランジュ方程式の第1変数です。
この例では、初期値theta0 = π/4 = 0.785ラジアン(約45度)から始まっています。
曲線の挙動は、単振り子の振動を表しています。

2. 緑色の曲線: theta_dot(t)

この曲線は、時間tに対する角速度theta_dot(単位はラジアン/秒) の変化を表しています。
theta_dot(t)は、オイラー-ラグランジュ方程式の第2変数です。
この例では、初期値theta_dot0 = 0 (静止状態)から始まっています。
曲線の挙動は、単振り子の振動の速さの変化を表しています。

グラフの横軸は時間tを表し、単位は秒です。
縦軸は、青曲線ではthetaの値(ラジアン)、緑曲線ではtheta_dotの値(ラジアン/秒)を表します。

曲線の挙動は、以下のパラメータによって決まります。

• m = 1.0 (質量)
• k = 1.0 (バネ定数)
• g = 9.81 (重力加速度)
• l = 1.0 (単振り子の長さ)

これらのパラメータを変更すると、曲線の振る舞いが変化します。

つまり、このグラフは単振り子の運動をオイラー-ラグランジュ方程式で記述し、その解の時間変化を可視化したものです。

力学の分野で重要な役割を果たすこの方程式の振る舞いを、視覚的に確認できます。

ケーリー方程式

ケーリー方程式は以下の方程式で表されます。

$$
y^2 = x^3 + ax + b
$$

ここで、$a$、$b$は定数です。

この方程式は以下の特徴を持っています。

• 左辺の$y^2$は2次式、右辺の$x^3 + ax + b$は3次式なので、この方程式は代数曲線を表しています。
• 特に、この曲線は楕円曲線と呼ばれる特殊な形状をしています。
  楕円は円に似た閉じた曲線です。
• $a$、$b$の値を変えると、楕円の形状が変わります。
• 楕円曲線上の点の集合は、加算やスカラー倍といった代数的な演算を定義できる代数構造になっています。
• この性質から、楕円曲線は暗号理論などの応用分野で重要な役割を果たしています。

つまり、ケーリー方程式は単に楕円曲線の方程式を表しているだけでなく、楕円曲線の持つ豊富な代数的構造を含んでいるのです。

簡単な方程式ながら、非常に深遠な数学的性質を持っているということができます。

ソースコード

Pythonでケーリー方程式を描画するコードは以下のようになります。

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

import matplotlib.pyplot as plt
import numpy as np

# 定数a、bを設定
a = 1
b = 1

# x、yの範囲を設定
x_range = np.linspace(-10, 10, 1000)
y_range = np.linspace(-10, 10, 1000)

# グリッド上の各点でy^2 = x^3 + ax + bを計算
X, Y = np.meshgrid(x_range, y_range)
Z = Y**2 - (X**3 + a*X + b)

# 等高線を描画
plt.contour(X, Y, Z, [0])
plt.title(f"Weierstrass Elliptic Curve: y^2 = x^3 + {a}x + {b}")
plt.xlabel("x")
plt.ylabel("y")
plt.axhline(y=0, color='k')
plt.axvline(x=0, color='k')
plt.show()

このコードでは、最初に定数$a$、$b$を設定しています。

次に、$x$、$y$の範囲を設定し、それぞれの点でケーリー方程式の値を計算しています。

その後、plt.contour関数を使って等高線を描画しています。
等高線の値が$0$の場合がちょうどケーリー方程式を満たす点になります。

最後に、グラフのタイトルと軸のラベルを設定し、x軸とy軸を描画しています。

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

[実行結果]

このグラフは、$a=1$、$b=1$のときのケーリー方程式を表しています。
$a$や$b$の値を変更すると、楕円曲線の形状が変わります。

グラフ解説

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

1. 曲線の形状

• グラフには曲線が描かれています。
  これがケーリー方程式 y^2 = x^3 + ax + b (ここでは a=1, b=1) を満たす点の集合です。
• 曲線の形状は、楕円に似た閉じた曲線になっています。
  これは楕円曲線の典型的な形状です。
• 曲線は原点$ (0, 0) $を通っています。これは、
  原点がケーリー方程式を自動的に満たすためです。

2. 対称性

• 曲線は x 軸と y 軸に対して点対称になっています。
  これは、ケーリー方程式が x と y に関して対称であることに由来します。
• また、曲線は原点に関しても点対称になっています。

3. 交点

• 曲線は x 軸上の 3 点$ (-2.49, 0) $、$ (0, 0) $、$ (2.49, 0) $で交わっています。
  これらの点は、方程式 y^2 = x^3 + x + 1 を満たす y=0 の解になっています。

4. 漸近線

• 曲線は x 軸や y 軸に漸近線を持ちません。
  つまり、無限遠で軸に近づくことはありません。

5. グリッド線

• グラフには x 軸と y 軸が黒色の太い線で描かれています。
  これにより、曲線の位置を視覚的に捉えやすくなっています。

6. 軸のラベル

• x 軸のラベルは “x”、y 軸のラベルは “y” と付けられています。

7. タイトル

• グラフのタイトルは “Weierstrass Elliptic Curve: y^2 = x^3 + x + 1” と表示されています。
  これにより、このグラフがケーリー方程式の特定のケース (a=1, b=1) を表していることがわかります。

このように、この単純な曲線グラフからも、ケーリー方程式や楕円曲線の様々な性質を読み取ることができます。

リンデロフ問題

リンデロフ問題は、混沌理論における代表的な例題です。

以下のような非線形2階常微分方程式で表されます。

$$
dx/dt = -y^2 - x
$$
$$
dy/dt = x
$$

この方程式は非常にシンプルですが、その解の振る舞いは複雑で予測不可能な様相を呈します。

初期値が僅かに異なっても、その後の軌道は指数関数的に発散してしまいます。
つまり、初期値のわずかな違いが、長時間後には大きな違いとなって現れるのです。
これが混沌の典型的な性質です。

位相空間軌跡を見ると、特徴的な蝶々形状の引力領域(アトラクター)が現れます。
この引力領域(アトラクター)は非常に複雑な構造を持ち、無限に細かい自己相似性を持っています。

また時間発展を見ると、振動が徐々に発散し、最終的に振幅が無限大になっていく様子がわかります。

このようにリンデロフ問題は、非線形力学系が持つ混沌的振る舞いの本質を簡単な方程式で表現しているため、混沌理論の基礎的な例題として広く知られています。

ソースコード

リンデロフ問題は、非線形の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
29
30
31
32
33
34
35
36
37
38

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

# リンデロフ方程式
def dU_dt(U, t):
    x, y = U
    return [-y**2 - x, x]

# 初期条件
x0 = 1.0
y0 = 1.0
U0 = [x0, y0]

# 時間範囲
t = np.linspace(0, 10, 10000)

# 数値解法
sol = odeint(dU_dt, U0, t)

# グラフ化
fig, ax = plt.subplots(1, 2, figsize=(12, 5))

# 位相空間軌跡
ax[0].plot(sol[:, 0], sol[:, 1])
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
ax[0].set_title('Phase Space Trajectory')

# 時間発展
ax[1].plot(t, sol[:, 0], label='x')
ax[1].plot(t, sol[:, 1], label='y')
ax[1].set_xlabel('t')
ax[1].set_title('Time Evolution')
ax[1].legend()

plt.tight_layout()
plt.show()

このコードでは、scipy.integrate.odeintを使ってリンデロフ方程式を数値的に解いています。

[実行結果]

最初のグラフは、位相空間軌跡($x$, $y$平面上の軌跡)を示しています。
リンデロフ問題の解は、特徴的な「蝶々」の形状を持つ軌道を描きます。

2番目のグラフは、$x(t)$と$y(t)$の時間発展を示しています。
振動が徐々に発散し、最終的に振幅が無限大になることがわかります。
これは、リンデロフ問題が非線形不安定であり、初期値の微小な摂動が時間とともに指数関数的に増幅されることを意味しています。

このように、リンデロフの問題は混沌系の典型的な振る舞いを示す簡単な例となっています。

ソースコード解説

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

1. インポート

まず必要なモジュールをインポートしています。

• numpy - 数値計算に使用
• matplotlib.pyplot - グラフの作成に使用
• scipy.integrate - 常微分方程式の数値解法に使用

2. リンデロフ方程式の定義

リンデロフ方程式は非線形の2次元常微分方程式で、次のように定義されています。

1
2
3

def dU_dt(U, t):
    x, y = U
    return [-y**2 - x, x]

この関数は、状態ベクトル U = [x, y] と時間 t を受け取り、x と y の微分値 dx/dt、dy/dt をタプルで返します。

3. 初期条件の設定

リンデロフ方程式を解くための初期条件を設定しています。

1
2
3

x0 = 1.0
y0 = 1.0
U0 = [x0, y0]

4. 時間範囲の設定

方程式を解く時間範囲を設定しています。
この例では t = 0 から t = 10 までの範囲を 10000 個の点で分割しています。

1

t = np.linspace(0, 10, 10000)

5. 数値解法

scipy.integrate.odeint 関数を使って、リンデロフ方程式を数値的に解いています。

1

sol = odeint(dU_dt, U0, t)

この関数は、方程式 dU_dt、初期値 U0、時間範囲 t を受け取り、各時間ステップにおける x と y の値を含む行列 sol を返します。

6. グラフ化

matplotlib.pyplot を使って、2つのグラフを作成しています。

• 1つ目のグラフは、位相空間軌跡 (x, y 平面上の軌跡) を表しています。
• 2つ目のグラフは、時間発展 x(t) と y(t) を表しています。

グラフのタイトル、軸ラベル、凡例などの設定を行っています。

7. グラフの表示

最後に plt.show() を実行して、作成したグラフを表示しています。

このコードでは、リンデロフ方程式を数値的に解き、その解の振る舞いを2種類のグラフで可視化しています。

位相空間軌跡からは混沌的な振る舞いが見て取れ、時間発展からは振動が徐々に発散していく様子がわかります。

このようにして、リンデロフ問題が混沌系の典型例であることを確認できます。

ヤング–ミルズ方程式

ヤング–ミルズ方程式は、場の理論で重要な方程式の一つですが、これをPythonで解くのは一般的に困難です。

なぜなら、この方程式は非線形偏微分方程式であり、厳密な解析解が存在しない場合がほとんどだからです。

ただし、近似的な数値解法を用いて方程式を解くことは可能です。
ただ、この方法でも一般的には複雑な解析の要素が残ります。

以下に、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.integrate import odeint

# ヤング–ミルズ方程式の定義
def young_mills_equations(y, x):
    # 例として単純な非線形微分方程式を使用
    dydx = [y[1], -np.sin(y[0])]
    return dydx

# 初期条件
y0 = [0, 1]  # y[0] が θ(x), y[1] が dθ/dx

# 積分範囲
x = np.linspace(0, 10, 1000)

# ヤング–ミルズ方程式を解く
solution = odeint(young_mills_equations, y0, x)

# 解をプロット
plt.plot(x, solution[:, 0], label='θ(x)')
plt.plot(x, solution[:, 1], label='dθ/dx')
plt.xlabel('x')
plt.ylabel('Solution')
plt.title('Solution to Young-Mills Equations')
plt.legend()
plt.grid(True)
plt.show()

このコードは、ヤング–ミルズ方程式の代わりに単純な非線形微分方程式を使用していますが、このアプローチをヤング–ミルズ方程式に拡張することができます。

結果をプロットすることで、数値解の挙動を視覚化することができます。

[実行結果]

ソースコード解説

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

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

1
2
3

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

• numpy は数値計算を行うためのライブラリであり、npとしてインポートされます。
• matplotlib.pyplot はグラフの描画に使用されるライブラリです。
• scipy.integrate モジュールから odeint 関数が使用され、常微分方程式を数値的に解きます。

2. ヤング–ミルズ方程式の定義:

1
2
3

def young_mills_equations(y, x):
    dydx = [y[1], -np.sin(y[0])]
    return dydx

• young_mills_equations 関数は、ヤング–ミルズ方程式（あるいは代替的な非線形微分方程式）を定義します。
  この例では、単純な非線形微分方程式が使用されています。
  yは未知関数のベクトル、xは独立変数です。
• ヤング–ミルズ方程式は実際には別のものであり、ここではその代替として非線形微分方程式が使われています。

3. 初期条件の設定:

1

y0 = [0, 1]  # y[0] が θ(x), y[1] が dθ/dx

• 初期条件 y0 は、未知関数とその導関数の初期値を指定します。
  ここでは、θ(x) の初期値が$0$、dθ/dx の初期値が$1$と設定されています。

4. 積分範囲の設定:

1

x = np.linspace(0, 10, 1000)

• np.linspace 関数を使って、積分の範囲を設定します。
  この場合、$0$から$10$の範囲を$1000$点で区切ります。

5. ヤング–ミルズ方程式の解を計算:

1

solution = odeint(young_mills_equations, y0, x)

• odeint 関数を使用して、ヤング–ミルズ方程式（または代替的な非線形微分方程式）の解を計算します。

初期条件 y0 と積分範囲 x が与えられます。

6. 解のプロット:

1
2
3
4
5
6
7
8

plt.plot(x, solution[:, 0], label='θ(x)')
plt.plot(x, solution[:, 1], label='dθ/dx')
plt.xlabel('x')
plt.ylabel('Solution')
plt.title('Solution to Young-Mills Equations')
plt.legend()
plt.grid(True)
plt.show()

• 解 solution をプロットします。
  solution[:, 0] はθ(x)、solution[:, 1] はdθ/dxに対応します。
• plt.xlabel と plt.ylabel はそれぞれ$x軸$と$y軸$のラベルを設定します。
• plt.title でグラフのタイトルを設定します。
• plt.legend で凡例を表示し、plt.grid でグリッド線を表示します。
• plt.show でグラフを表示します。

このプログラムは、ヤング–ミルズ方程式やその他の非線形微分方程式を解いて、その解をグラフ化する基本的な手法を示しています。

結果解説

[実行結果]

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

1. x軸：

変数$ x $の範囲です。
この範囲は np.linspace(0, 10, 1000) を用いて、$0$から$10$までの区間を$1000$等分した値です。
この範囲内で微分方程式を解いています。

2. y軸（左）：

θ(x)の値です。
ヤング–ミルズ方程式では、θ(x)は場の理論の中で重要な量の一つであり、その解釈は理論の背後にある対称性や場の構造を理解するのに役立ちます。

3. y軸（右）：

dθ/dxの値です。
これはθ(x)の導関数であり、θ(x)がどのように変化しているかを示します。
ヤング–ミルズ方程式において、場の構造がどのように変化するかを理解するのに重要な情報を提供します。

4. タイトル：

グラフのタイトルは「Solution to Young-Mills Equations（ヤング–ミルズ方程式の解）」です。
これは、解がヤング–ミルズ方程式から得られたものであることを示しています。

5. 凡例：

プロットされている曲線の説明です。
“θ(x)”はθ(x)の値を表し、”dθ/dx”はdθ/dxの値を表します。

6. グリッド：

グラフ内にはグリッド線が表示されており、データの視覚的な解釈を補助します。

ナーシム–シュレーダー方程式

ナーシム–シュレーダー方程式（Nagumo–Schaffer Equation）は、神経科学や興奮伝播のモデルとして広く使用される非線形微分方程式の一つです。

ナーシム–シュレーダー方程式は、以下の形式の非線形微分方程式で表します：

$$
\frac{du}{dt} = u - u^3 - v + a
$$

$$
\frac{dv}{dt} = b(u - cv)
$$

ここで、$ (u) $と$ (v) $は時間$ (t) $の関数であり、$ (a) $、$(b)$、$(c) $は定数です。

これらの方程式は興奮伝播のモデルとして広く使用されます。

第1式は$ (u) $の時間変化を、第2式は$ (v) $の時間変化を表します。
$(a) $は$ (u) $のリセット項、$(b) $は興奮の増幅率、$(c) $は相互抑制の効果を示します。

ソースコード

以下では、この方程式を Python を使用して解き、結果をグラフ化します。

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

1
2
3

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

次に、ナーシム–シュレーダー方程式を定義します。

1
2
3

def nagumo_schaffer(t, y, a, b, c):
    return [y[0] - y[0]**3 - y[1] + a,
            b * (y[0] - c * y[1])]

この方程式には3つのパラメータ$ (a)$、$(b)$、$(c) $があります。

これらは、方程式の特性を制御します。

次に、初期条件と時間の範囲を設定します。

1
2

y0 = [0.1, 0.1]  # 初期条件 [u, v]
t_span = (0, 50)  # シミュレーション時間の範囲

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

1
2
3
4

a = 0.2
b = 0.2
c = 3.0
sol = solve_ivp(nagumo_schaffer, t_span, y0, args=(a, b, c), t_eval=np.linspace(0, 50, 1000))

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

1
2
3
4
5
6
7
8

plt.plot(sol.t, sol.y[0], label='u(t)')
plt.plot(sol.t, sol.y[1], label='v(t)')
plt.xlabel('Time')
plt.ylabel('Value')
plt.title('Nagumo-Schaffer Equation Solution')
plt.legend()
plt.grid(True)
plt.show()

これにより、ナーシム–シュレーダー方程式の解が時間とともにどのように振る舞うかを示すグラフが生成されます。

結果解説

[実行結果]

このグラフは、ナーシム–シュレーダー方程式の解を示しています。

方程式は2つの変数$ (u(t)) $と$ (v(t)) $で構成されており、時間$ (t) $に関する解をプロットしています。

具体的には、以下の内容がグラフに表示されます：

1. 時間(Time):

x軸は時間を表しています。この例では、$0$から$50$の範囲で、$1000$の等間隔の時間ステップがプロットされます。

2. 値(Value):

y軸は$ (u(t)) $と$ (v(t)) $の値を表します。
それぞれの変数の時間変化がプロットされます。

3. (u(t)) 曲線:

ナーシム–シュレーダー方程式の解$ (u(t)) $の時間変化がプロットされます。
この曲線は興奮性の変化を示しています。

4. (v(t)) 曲線:

ナーシム–シュレーダー方程式の解$ (v(t)) $の時間変化がプロットされます。
この曲線は抑制性の変化を示しています。

5. タイトル(Title):

グラフのタイトルは “Nagumo-Schaffer Equation Solution” となっています。

6. 凡例(Legend):

各プロットに対応する凡例が表示されます。
u(t)は$ u(t)$、v(t) は$ v(t) $を表しています。

7. グリッド(Grid):

グリッド線が表示され、目盛りを読み取りやすくします。

このグラフは、ナーシム–シュレーダー方程式の解が時間とともにどのように振る舞うかを可視化しており、その結果が興奮性と抑制性の両方の側面を捉えています。

フェルマーの最終定理

フェルマーの最終定理は「$n>2$のとき、$x^n + y^n = z^n$を満たす自然数$x,y,z$は存在しない」というものです。

これをPythonで表現し、グラフ化するには以下のようなコードを使うことができます。

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

import matplotlib.pyplot as plt
import numpy as np

# パラメータ設定
n = 3  # 検証する指数
max_val = 100  # x,y,zの最大値

# グリッド上の全ての点をチェック
x = np.arange(max_val)
y = np.arange(max_val)
X, Y = np.meshgrid(x, y)
Z = X**n + Y**n

# 方程式を満たす点をプロット
plt.figure(figsize=(8, 6))
plt.scatter(X[Z <= max_val**n], Y[Z <= max_val**n], s=1, c='b')
plt.scatter(X[Z > max_val**n], Y[Z > max_val**n], s=1, c='r')
plt.xlim(0, max_val)
plt.ylim(0, max_val)
plt.title(f'Fermat\'s Last Theorem (n={n})')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

このコードでは、指数$n$と$x$,$y$,$z$の最大値を設定できます。

そして、グリッド上の全ての$(x,y)$の組み合わせについて、$x^n + y^n$が$max＿val^n$を超えるかどうかをチェックしています。

方程式を満たす点$(x^n + y^n <= max＿val^n)$は青色で、満たさない点は赤色でプロットされます。

グラフの例(n=3, max_val=100)

このグラフから、$n=3$の場合、$x^3 + y^3 = z^3$を満たす自然数の組$(x,y,z)$が存在しないことが分かります。
$max＿val$を大きくすれば、より大きな値についてもフェルマーの最終定理が成り立つことを確認できます。

ソースコード解説

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

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

1
2

import matplotlib.pyplot as plt
import numpy as np

• matplotlib.pyplot: グラフを描画するためのライブラリ
• numpy: 数値計算ライブラリ

2. パラメータ設定

1
2

n = 3  # 検証する指数
max_val = 100  # x,y,zの最大値

• $n$には検証したい指数(この例では$n=3$)を設定します
• $max_val$には、$x$,$y$,$z$の最大値を設定します

3. グリッド上の全ての点をチェック

1
2
3
4

x = np.arange(max_val)
y = np.arange(max_val)
X, Y = np.meshgrid(x, y)
Z = X**n + Y**n

• $x, y$にそれぞれ$0$から$max＿val-1$までの値を格納した1次元NumPy配列を作成
• np.meshgrid(x, y)で、$x$,$y$の全ての組み合わせからなる2次元グリッドX,Yを生成
• X**n + Y**nで、各グリッド点$(x,y)$における$x^n + y^n$の値を計算し、$Z$にNumPy配列として格納

4. 方程式を満たす点をプロット

1
2
3
4
5
6
7
8
9

plt.figure(figsize=(8, 6))
plt.scatter(X[Z <= max_val**n], Y[Z <= max_val**n], s=1, c='b')
plt.scatter(X[Z > max_val**n], Y[Z > max_val**n], s=1, c='r')
plt.xlim(0, max_val)
plt.ylim(0, max_val)
plt.title(f'Fermat\'s Last Theorem (n={n})')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

• plt.figure(figsize=(8, 6))で新しい図を作成し、サイズを指定
• 1つ目のplt.scatter():
  • X[Z <= max_val**n]は、$x^n + y^n <= max＿val^n$を満たす$x$の値
  • Y[Z <= max_val**n]は、$x^n + y^n <= max＿val^n$を満たす$y$の値
  • s=1で点のサイズを1に設定
  • c=’b’で点の色を青に設定
• 2つ目のplt.scatter():
  • X[Z > max_val**n]は、$x^n + y^n > max_val^n$を満たす$x$の値
  • Y[Z > max_val**n]は、$x^n + y^n > max_val^n$を満たす$y$の値
  • s=1で点のサイズを1に設定
  • c=’r’で点の色を赤に設定
• plt.xlim(0, max_val)、plt.ylim(0, max_val)でx軸、y軸の範囲を設定
• plt.title(f’Fermat's Last Theorem (n={n})’)でグラフのタイトルを設定
• plt.xlabel(‘x’)、plt.ylabel(‘y’)で軸ラベルを設定
• plt.show()でグラフを表示

このコードでは、指定した範囲の全ての$(x,y)$の組み合わせについて、$x^n + y^n$の値が$max＿val^n$を超えるかどうかをチェックし、超えない点(方程式を満たす点)を青色で、超える点を赤色でプロットしています。

また、適切な軸範囲、タイトル、ラベルを設定してグラフを描画しています。

結果解説

[実行結果]

このグラフには以下の情報が表示されています。

1. 点の色

青い点: $x^n + y^n $の値が$max＿val^n$以下の$(x,y)$の組
赤い点: $x^n + y^n $の値が$max＿val^n$を超える$(x,y)$の組

2. 軸

x軸: x値 ($0$から$max＿val$まで)
y軸: y値 ($0$から$max＿val$まで)

3. タイトル

グラフのタイトルには “Fermat’s Last Theorem (n=3)” と表示されています。
つまり、このグラフはフェルマーの最終定理を$n=3$の場合について可視化したものです。

4. 点の分布

原点$(0, 0)$を中心に、青い点が放射状に広がっています。
赤い点は原点から離れた領域に存在しています。
これは、小さな値の$x$,$y$に対しては$x^3 + y^3$が$max＿val^3$以下になりやすいが、$x$,$y$が大きくなるにつれて$max＿val^3$を超えやすくなることを示しています。

5. グラフの形状

点の分布から、青い領域は滑らかなカーブ状の境界線を持っていることが分かります。
この境界線は、$x^3 + y^3 = max＿val^3 $という3次曲面を表しています。

6. 点の密度

原点付近では青い点が密集していますが、原点から離れるにつれて点の密度が低下しています。
これは、小さな値の$x$,$y$の組み合わせの方が多いためです。

このグラフを見ることで、$n=3$の場合にはフェルマーの最終定理が視覚的に裏付けられていることが分かります。

青い点は$x^n + y^n = z^n$を満たす可能性のある$(x,y)$の組を表し、赤い点はその組み合わせが存在しないことを示しています。

ジュリア集合

ジュリア集合は、複素関数の反復による軌道の振る舞いを示す集合です。

ジュリア集合の収束性の条件を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

import numpy as np
import matplotlib.pyplot as plt

def julia(c, x, y, max_iter=100):
    z = x + 1j * y
    for i in range(max_iter):
        z = z ** 2 + c
        if np.abs(z) > 2:
            return i
    return max_iter

def plot_julia(c, zoom=1, x_range=(-2, 2), y_range=(-1.5, 1.5), max_iter=100):
    x = np.linspace(x_range[0] / zoom, x_range[1] / zoom, 1000)
    y = np.linspace(y_range[0] / zoom, y_range[1] / zoom, 1000)
    X, Y = np.meshgrid(x, y)

    Z = np.zeros_like(X, dtype=int)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = julia(c, X[i, j], Y[i, j], max_iter)

    plt.figure(figsize=(10, 8))
    plt.imshow(Z, cmap='hot', extent=(x_range[0], x_range[1], y_range[0], y_range[1]))
    plt.colorbar()
    plt.xlabel('Re(c)', fontsize=14)
    plt.ylabel('Im(c)', fontsize=14)
    plt.title(f'Julia Set for c={c:.2f}', fontsize=16)
    plt.show()

# 例: c = -0.8 + 0.156j のジュリア集合
plot_julia(-0.8 + 0.156j, zoom=1, max_iter=100)

このコードでは、まず julia 関数を定義しています。
この関数は、与えられた複素数 c と初期値 x, y について、反復 z = z**2 + c を最大 max_iter 回行い、発散するまでの反復回数を返します。

次に、plot_julia 関数を定義しています。
この関数は、与えられた c について、ジュリア集合をプロットします。
まず、x と y の範囲を設定し、それらの組み合わせについて julia 関数を呼び出して発散までの反復回数を求めます。
その結果を行列 Z に格納し、plt.imshow でプロットしています。

最後に、plot_julia(-0.8 + 0.156j, zoom=1, max_iter=100) と指定して、ジュリア集合をプロットしています。

実行結果は、複素平面上のジュリア集合を可視化した画像になります。

[実行結果]

色が濃いほど発散までの反復回数が多く、白い部分が発散しない領域(ジュリア集合自体)になります。
このジュリア集合は、フラクタル構造を持っていることがわかります。

ジュリア集合は、カオス理論や複素動力学系の研究で重要な役割を果たしています。
また、その美しいフラクタル構造から、コンピュータグラフィックスなどでも活用されています。

ソースコード解説

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

1 : ジュリア集合の計算関数

1
2
3
4
5
6
7

def julia(c, x, y, max_iter=100):
    z = x + 1j * y
    for i in range(max_iter):
        z = z ** 2 + c
        if np.abs(z) > 2:
            return i
    return max_iter

• julia関数は、与えられた複素数cと初期値(x, y)について、反復z = z**2 + cを最大max_iter回行います。
• 反復の過程で$z$の絶対値が$2$を超えれば(発散した場合)、その反復回数を返します。
• max_iter回まで発散しなければmax_iterを返します。

2 : ジュリア集合の可視化関数

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

def plot_julia(c, zoom=1, x_range=(-2, 2), y_range=(-1.5, 1.5), max_iter=100):
    x = np.linspace(x_range[0] / zoom, x_range[1] / zoom, 1000)
    y = np.linspace(y_range[0] / zoom, y_range[1] / zoom, 1000)
    X, Y = np.meshgrid(x, y)

    Z = np.zeros_like(X, dtype=int)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = julia(c, X[i, j], Y[i, j], max_iter)

    plt.figure(figsize=(10, 8))
    plt.imshow(Z, cmap='hot', extent=(x_range[0], x_range[1], y_range[0], y_range[1]))
    plt.colorbar()
    plt.xlabel('Re(c)', fontsize=14)
    plt.ylabel('Im(c)', fontsize=14)
    plt.title(f'Julia Set for c={c:.2f}', fontsize=16)
    plt.show()

• plot_julia関数は、与えられたcについてジュリア集合をプロットします。
• まず、$x$、$y$の範囲を設定し、それらの組み合わせについてjulia関数を呼び出します。
• 結果をZ行列に格納し、plt.imshowでプロットします。
• 軸ラベル、タイトル、カラーバーを設定し、グラフを表示します。

3 : プログラムの実行

1
2

# 例: c = -0.8 + 0.156j のジュリア集合
plot_julia(-0.8 + 0.156j, zoom=1, max_iter=100)

• plot_julia関数を呼び出し、c = -0.8 + 0.156jのジュリア集合をプロットしています。
• zoomは$1$(デフォルト)、max_iterは100回に設定されています。

このコードは、ジュリア集合の計算と可視化を行っています。主な手順は以下の通りです。

1. julia関数で、与えられた$c$と初期値$ (x, y) $について、反復z = z**2 + cを行い、発散するまでの反復回数を求めます。

2. plot_julia関数で、複素平面上の$(x, y)$の範囲を設定し、それぞれの点についてjulia関数を呼び出します。

3. 得られた発散までの反復回数を行列Zに格納します。

4. 行列Zをプロットし、発散速度(反復回数)を色で可視化します。

5. 軸ラベル、タイトル、カラーバーを設定してグラフを表示します。

最終的に、ジュリア集合の複雑で自己相似的な構造を可視化したグラフが得られます。

結果解説

[実行結果]

このグラフはジュリア集合の様子を可視化したものです。

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

• グラフの横軸は複素数$c$の実部を、縦軸は$c$の虚部を表しています。
  つまり、複素平面上の点がプロットされています。
• 背景の色は、ある初期値$(x, y)$からジュリア集合の反復$z = z**2 + c$を行ったときの発散する速さを表しています。
  赤や黄色は早く発散し、青や紫は遅く発散することを意味します。
• 一方で、白い領域は発散しない点の集まりです。
  これがジュリア集合自体を表しています。
• ジュリア集合の形状は、与えられた$c$の値によって大きく異なります。
  このグラフでは$c = -0.8 + 0.156j$が使われています。
• ジュリア集合自体は非常に複雑な構造を持ち、フラクタル的な自己相似性が見られます。
  つまり、拡大するとほぼ同じ模様が現れる性質があります。
• 白い領域(ジュリア集合)の周りに同心円状の渦巻き模様があり、そこから外に向かうにつれて発散速度が速くなっていきます。
• グラフ右側のカラーバーは、発散速度(反復回数)と色の対応関係を示しています。

このように、ジュリア集合は複素動力学系における重要な概念で、それ自体が不規則でフラクタル的な美しい構造を持っています。

この可視化により、ジュリア集合の様々な性質を確認できます。