ヴァンデル・ポルト方程式

2つの振動子が相互作用する場合のヴァンデル・ポルト方程式があります。

この方程式は、化学や物理学、生物学などのさまざまな分野で用いられます。

2つの振動子$ (x) $と$ (y) $の間の相互作用を表すヴァンデル・ポルト方程式は次のように表されます:

ここで、$ (a) $、$ (b) $、$ (\alpha) $、$ (\epsilon) $、$ (\gamma) $はパラメータです。

この方程式をPythonで3Dグラフ化するには、以下のようなコードを使用できます。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

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

# ヴァンデル・ポルト方程式の定義
def vanderpol_equations(x, y, alpha, epsilon, gamma, a, b):
    dx_dt = alpha * (x - (x**3)/3 - y + a)
    dy_dt = epsilon * (x - gamma * y + b)
    return dx_dt, dy_dt

# パラメータの設定
alpha = 1.0
epsilon = 1.0
gamma = 2.0
a = 1.0
b = 1.0

# 時間のステップと初期値
dt = 0.01
steps = 10000

x_vals = np.zeros(steps + 1)
y_vals = np.zeros(steps + 1)

x_vals[0], y_vals[0] = 0.1, 0.1

# 数値積分（Euler法）
for i in range(steps):
    dx, dy = vanderpol_equations(x_vals[i], y_vals[i], alpha, epsilon, gamma, a, b)
    x_vals[i + 1] = x_vals[i] + dt * dx
    y_vals[i + 1] = y_vals[i] + dt * dy

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x_vals, y_vals, np.arange(steps + 1))

ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Time-axis')
ax.set_title('Van der Pol Oscillator')

plt.show()

このコードは、ヴァンデル・ポルト方程式を数値的に解き、結果を3Dプロットとして表示します。

これにより、2つの振動子の相互作用が示され、その振る舞いが視覚化されます。

[実行結果]

ソースコード解説

このコードは、ヴァンデル・ポルト方程式と呼ばれる非線形微分方程式を数値的に解き、その結果を3Dプロットとして表示しています。

以下がコードの詳細です。

1. ヴァンデル・ポルト方程式の定義:

• vanderpol_equations関数では、与えられた$ (x) $と$ (y) $の値に基づいて、ヴァンデル・ポルト方程式を計算しています。
  これらの方程式は2つの微分方程式で構成されており、$ (x) $および$ (y) $の時間変化を計算します。

2. パラメータの設定:

• ヴァンデル・ポルト方程式にはいくつかのパラメータがあり、このコードではそれらを設定しています。
  例えば、$ (\alpha) $、$ (\epsilon) $、$ (\gamma) $、$ (a) $、$ (b) $の値が与えられています。

3. 数値積分（Euler法）:

• ループを使用して、Euler法を使って微分方程式を数値的に解いています。
  初期値$ (x = 0.1) $および$ (y = 0.1) $から始め、ループを通じて微小時間$ (dt) $だけ変化を加えて新しい$ (x) $と$ (y) $の値を計算しています。

4. 3Dプロット:

• 最後に、Matplotlibを使用して、計算された$ (x) $と$ (y) $の値を3Dプロットとして可視化しています。
  このプロットは、$ (x) $軸と$ (y) $軸に対応する値が時間とともにどのように変化するかを示しています。

5. グラフの要素:

• X軸は振動子$ (x) $の値を表します。
• Y軸は振動子$ (y) $の値を表します。
• Z軸（時間軸）は時間を表し、振動子の値が時間とともにどのように変化するかを示します。

このコードは、ヴァンデル・ポルト方程式を数値的に解き、振動子の振る舞いを視覚化しています。

振動子の振る舞いはパラメータと初期値に依存するため、これらの値を変更すると異なる振る舞いが得られます。

グラフ解説

[実行結果]

このコードは、ヴァンデル・ポルト方程式と呼ばれる非線形微分方程式を解き、その結果を3Dプロットとして表示しています。

ヴァンデル・ポルト方程式は、2つの振動子$ (x) $と$ (y) $の相互作用を表しています。

ここではEuler法を用いて微分方程式を数値的に解き、結果を3Dグラフにプロットしています。

以下はグラフの要素についての詳細です。

1. 軸の意味:

• X軸は振動子$ (x) $の値を表します。
• Y軸は振動子$ (y) $の値を表します。
• Z軸（時間軸）は時間を表し、各点が時間の経過に伴って振動子の値がどのように変化するかを示しています。

2. グラフの特徴:

• ヴァンデル・ポルト方程式は非線形であり、2つの振動子$ (x) $と$ (y) $の間の相互作用が時間とともにどのように振る舞うかを表します。
• このグラフは、$ (x) $と$ (y) $の値が時間に応じてどのように変化するかを示しており、振動子の軌道が示されています。
  これらの軌道は、特定の初期値から始まり、振動子の相互作用によって時間の経過とともに変化します。

3. 振動の特性:

• ヴァンデル・ポルト方程式は、非線形な振動子系の一例です。
  このグラフは、振動子の非線形な相互作用により、複雑な振る舞いを示しています。
  振動子が周期的に振動する場合や、カオス的な振る舞いを示す場合があります。

このグラフは、2つの振動子が互いに相互作用するヴァンデル・ポルト方程式の振る舞いを視覚化しており、非線形ダイナミクスが示されています。

振動子がどのように相互作用し、時間とともにどのように変化するかを理解するのに役立ちます。

サドル点

一般的な3次元プロットでよく使用される関数としては、$z = x^2 - y^2$（サドル点）があります。

以下にPythonのmatplotlibとnumpyを用いた実装例を示します。

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

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

#x, yの値を定義
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
x, y = np.meshgrid(x, y)

#z = x^2 - y^2という関数
z = x**2 - y**2

#3D描画
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z)

# グラフ表示
plt.show()

このプログラムを実行すると、3Dグラフが表示されます。

numpyのlinspace関数を使用して$x$と$y$の範囲を定義し、meshgrid関数で$x$と$y$の座標を生成しています。

その座標の組に対して関数$ (z = x^2 - y^2) $を適用して$z$値を求め、plot_surface関数で3Dグラフに描画しています。

[実行結果]

ソースコード解説

このコードは、PythonのライブラリであるNumPyとMatplotlibを使用して、$ (z = x^2 - y^2) $という関数の3Dプロットを作成しています。

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

1. import文:

• numpyとmatplotlib.pyplotから必要なモジュールをインポートしています。
  また、3Dプロットを作成するためのAxes3Dもインポートされています。

2. x, yの値の定義:

• np.linspaceを使って、$-5$から$5$までの範囲を$100$個の等間隔な点で$x$と$y$軸に割り当てています。
• meshgridを使って、$x$と$y$の格子点を作成しています。
  これにより、$x$と$y$の組み合わせのすべての点が作成されます。

3. zの計算:

• $ (z = x^2 - y^2) $の関数を表すため、各点における$ (x^2 - y^2) $の値を計算して$z$に割り当てています。

4. 3D描画の設定:

• plt.figure()を使用して新しい図を作成し、Axes3D関数を使用して3Dの軸を定義します。

5. ax.plot_surface:

• plot_surfaceを使って、$x$、$y$、および対応する$z$の値を用いて表面プロットを行っています。

6. グラフの表示:

• plt.show()を使用して、作成した3Dグラフを表示します。

このコード全体の目的は、$ (z = x^2 - y^2) $の関数を3Dプロットして、その曲面の形状を可視化することです。

$ x$と$y$が範囲内で変化すると、関数の値$z$がそれに応じて変化し、それが3次元のグラフとして表示されます。

シュレーディンガー方程式

理論物理学で使われる複雑な数式として、シュレーディンガー方程式を挙げてみましょう。

シュレーディンガー方程式は量子力学における基本的な方程式の1つです。

$$
[
i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + V(\mathbf{r}, t) \psi
]
$$

こちらの方程式を使って、2次元のハーモニック・オシレーターの確率密度をプロットするPythonの例を示します。

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

import numpy as np
import matplotlib.pyplot as plt

# 2Dハーモニック・オシレーターの確率密度関数
def harmonic_oscillator(x, y):
    return np.exp(-x**2 - y**2) * np.cos(2 * x) * np.sin(2 * y)**2

# メッシュグリッドの生成
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x, y)
Z = harmonic_oscillator(X, Y)

# 2Dプロット
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, Z, cmap='viridis')
plt.colorbar()
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.title('2D Harmonic Oscillator Probability Density')
plt.show()

このコードは、2次元のハーモニック・オシレーターの確率密度関数を定義し、それをプロットしています。

得られる2Dプロットは、確率密度が異なる位置においてどのように分布するかを示しています。

[実行結果]

ソースコード解説

このプログラムは、2次元のハーモニック・オシレーターの確率密度関数を計算し、その結果を2次元プロットで視覚化しています。

1. import文:

• numpyとmatplotlib.pyplotをインポートし、数値計算とプロット機能を使用できるようにします。

2. harmonic_oscillator関数:

• 2次元ハーモニック・オシレーターの確率密度関数を定義しています。
• この関数は、引数として$x座標$と$y座標$を受け取り、その位置における確率密度を計算します。
• 式には指数関数、三角関数が含まれており、座標に基づいて確率密度が計算されます。

3. メッシュグリッドの生成:

• np.linspace()を使用して、$-3$から$3$までの$x軸$と$y軸$の値を生成します。
• np.meshgrid()を使用して、これらの値からメッシュグリッドを生成します。
• 生成されたメッシュグリッドは、確率密度関数の各座標に対応します。

4. 2Dプロット:

• plt.figure()を使用して図のサイズを指定します。
• plt.contourf()を使用して、2次元の等高線プロットを作成します。
  等高線は確率密度を表し、カラーマップ（viridis）で表現されます。
• plt.colorbar()を使用して、カラーバーを表示し、色と確率密度の関連付けを提供します。
• plt.xlabel()とplt.ylabel()を使用して、$x軸$と$y軸$のラベルを指定します。
• plt.title()を使用して、プロットのタイトルを設定します。
• plt.show()でプロットを表示します。

このプログラムは、2次元のハーモニック・オシレーターの確率密度を計算し、等高線プロットとして視覚化しています。

等高線の濃淡が異なる領域での確率密度を示しており、それによって確率密度の分布を視覚的に理解できます。

グラフ解説

[実行結果]

このプロットは、2次元のハーモニック・オシレーターの確率密度関数を示しています。
ハーモニック・オシレーターは、量子力学や物理学で頻繁に使われるモデルの1つで、振動する粒子や原子の運動を表現します。

プロットされたグラフは、ハーモニック・オシレーターの確率密度を示しており、色の濃淡で異なる確率密度を表現しています。
色が濃い領域ほど高い確率密度を持ち、色が薄い領域ほど低い確率密度を持ちます。

ここで使用された確率密度関数は、指数関数や三角関数を組み合わせたもので、粒子の位置の特定の座標における確率を表現しています。

このプロットから、確率密度が異なる位置でどのように変化するかを視覚的に理解することができます。

$x軸$と$y軸$は、ハーモニック・オシレーターの位置座標を示しており、それぞれの座標での確率密度をカラーマップを用いて表現しています。

モビウス帯

モビウス帯は、帯状の特殊な幾何学的な面です。
これは通常の帯とは異なり、一周させると表面が反転してしまう特性を持っています。
この特性は、帯の表と裏の区別がなく、外側と内側が連続していることを示しています。

モビウス帯を作るには、帯の幅を持つストリップ（帯状の紙など）を取り、一端を半回転させて反対側とつなげます。

すると、通常の帯では表と裏がありますが、モビウス帯では一周しても片側だけになります。
これは3次元の空間における非常に興味深い形状であり、数学的にも重要です。

モビウス帯はトポロジー学や数学の教育において非常に興味深い対象であり、幾何学や物理学、コンピュータグラフィックスなど様々な分野で応用されています。

また、この特殊な形状は数学の美しさを示すための象徴的なものとしても知られています。

数式

一般的なモビウス帯は、次のようなパラメトリック方程式で表されます：

\begin{align*}
x(u, v) &= \left( r + v \cdot \cos\left(\frac{u}{2}\right) \right) \cdot \cos(u) \
\end{align*}
\begin{align*}
y(u, v) &= \left( r + v \cdot \cos\left(\frac{u}{2}\right) \right) \cdot \sin(u) \
\end{align*}
\begin{align*}
z(u, v) &= v \cdot \sin\left(\frac{u}{2}\right)
\end{align*}

ここで、$ (u) $は$ 0 $から$ (2\pi) $までのパラメータで、モビウス帯を一周させるための角度を表しています。

$ (v) $は帯の幅を決定するパラメータで、通常は$ (-w) $から$ (w) $の範囲で変化します。
そして$ (r) $は帯の中心からの距離を表します。

このパラメトリック方程式では、$ (u) $と$ (v) $の値の範囲によってモビウス帯の形状が変化します。

$ (u) $を$ 0 $から$ (2\pi) $まで変化させると、モビウス帯が一周しますが、一周した際に表面が反転することが特徴です。

また、$ (v) $を変化させることで帯の幅を調整することができます。

グラフ化

モビウス帯を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

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

# モビウス帯のパラメータ
r, w = 1, 0.2  # 半径と帯の幅
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(-w, w, 50)
U, V = np.meshgrid(u, v)
X = (r + V * np.cos(U / 2)) * np.cos(U)
Y = (r + V * np.cos(U / 2)) * np.sin(U)
Z = V * np.sin(U / 2)

# 3Dプロット
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

# グラフの設定
ax.set_title('Mobius Strip')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

このコードは、モビウス帯のパラメータを使用して3Dプロットを行い、Matplotlibで可視化します。

実行すると、モビウス帯の形状を表示するグラフが生成されます。

[実行結果]

ソースコード解説

このコードは、PythonのNumPyとMatplotlibを使用してモビウス帯をプロットしています。

各行の役割や手順を見ていきましょう。

モビウス帯のパラメータ設定

1
2
3
4
5
6
7

r, w = 1, 0.2  # 半径と帯の幅
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(-w, w, 50)
U, V = np.meshgrid(u, v)
X = (r + V * np.cos(U / 2)) * np.cos(U)
Y = (r + V * np.cos(U / 2)) * np.sin(U)
Z = V * np.sin(U / 2)

• r はモビウス帯の半径、w は帯の幅を指定します。
• u は$0$から$2π$までを$100$等分した配列を作り、v は -w から w までを$50$等分した配列を作ります。
• np.meshgrid を使用して、u と v のメッシュグリッドを作成します。
  これにより、U と V の2つの行列が得られます。
• X, Y, Z は、パラメータ式を使ってモビウス帯の各点の座標を計算します。

3Dプロット

1
2
3

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

• fig は図を表し、ax はその図上に配置される3Dサブプロットを表します。
• ax.plot_surface で X, Y, Z のデータからモビウス帯の表面をプロットします。

グラフの設定と表示

1
2
3
4
5
6

ax.set_title('Mobius Strip')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

• ax.set_title, ax.set_xlabel などでグラフのタイトルや軸ラベルを設定します。
• plt.show() でプロットを表示します。

このコードは、NumPyを使ってモビウス帯のパラメータを設定し、Matplotlibを使ってそれを3Dプロットしています。

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
24
25
26
27
28
29
30
31
32
33
34

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigs

# 定数
N = 1000  # 格子点の数
mass = 1.0  # 粒子の質量
hbar = 1.0  # 換算プランク定数
L = 1.0  # 井戸の幅
x = np.linspace(0, L, N)  # 位置

# 差分化
dx = x[1] - x[0]
T = hbar**2 / (2.0 * mass * dx**2) * (diags([1, -2, 1], [-1, 0, 1], shape=(N, N)).toarray())

# エネルギー固有値と波動関数の計算
values, vecs = eigs(T, k=30, which='SM')
idx = values.real.argsort()
values = values.real[idx]
vecs = vecs.real[:, idx]

# 波動関数の正規化とプロット
for n in range(3):  # 最初の3つのエネルギー状態をプロット
    psi = vecs[:, n] / np.sqrt(dx)
    plt.plot(x, psi, label='E = %.2f' % values[n])

plt.title('Wavefunctions in a 1D Infinite Square Well')
plt.xlabel('x')
plt.ylabel('psi')
plt.legend(loc='upper right')
plt.grid(True)

plt.show()

このコードでは最初の3つのエネルギー状態の波動関数をプロットします。

エネルギーが増えるにつれて波のノード（波の振幅がゼロとなる点）の数が増えることが確認できます。

これは量子力学における特從性の一つであり、粒子のエネルギー状態を視覚化することができます。

[実行結果]

ソースコード解説

このコードは、1次元の無限長方形井戸内の粒子の波動関数を計算し、その結果をグラフ化するものです。

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

1
2
3
4

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import eigs

• numpyは数値計算のためのライブラリで、配列操作や数学関数を提供します。
• matplotlib.pyplotはグラフ描画のためのライブラリで、グラフの作成やカスタマイズを行います。
• scipy.sparse.diagsは疎行列を扱うためのライブラリで、特定のパターンの疎行列を生成します。
• scipy.sparse.linalg.eigsは疎行列の固有値問題を解くためのライブラリで、固有値と固有ベクトルを計算します。

2. 定数の設定

1
2
3
4
5

N = 1000  # 格子点の数
mass = 1.0  # 粒子の質量
hbar = 1.0  # 換算プランク定数
L = 1.0  # 井戸の幅
x = np.linspace(0, L, N)  # 位置

• Nは格子点の数、massは粒子の質量、hbarは換算プランク定数、Lは井戸の幅を表します。
• np.linspaceは$0$から L までの範囲を N 分割した値の配列を生成します。
  これにより、位置 x の配列が作成されます。

3. 差分化

1
2

dx = x[1] - x[0]
T = hbar**2 / (2.0 * mass * dx**2) * (diags([1, -2, 1], [-1, 0, 1], shape=(N, N)).toarray())

• dxは隣接する位置の差を表します。
• diags関数は特定のパターンの疎行列を生成します。
  ここでは、2階微分の差分演算子を表す疎行列$ (T) $を生成しています。

4. エネルギー固有値と波動関数の計算

1
2
3
4

values, vecs = eigs(T, k=30, which='SM')
idx = values.real.argsort()
values = values.real[idx]
vecs = vecs.real[:, idx]

• eigs関数を使用して疎行列$ (T) $の固有値と対応する固有ベクトル（波動関数）を計算します。
• k=30は最初の$30個$の固有値と固有ベクトルを取得することを意味します。
• which='SM'は最小の固有値から計算を行うことを指定します。

5. 波動関数の正規化とプロット

1
2
3
4
5
6
7
8
9
10
11

for n in range(3):  # 最初の3つのエネルギー状態をプロット
    psi = vecs[:, n] / np.sqrt(dx)
    plt.plot(x, psi, label='E = %.2f' % values[n])

plt.title('Wavefunctions in a 1D Infinite Square Well')
plt.xlabel('x')
plt.ylabel('psi')
plt.legend(loc='upper right')
plt.grid(True)

plt.show()

• 取得した固有ベクトル（波動関数）を正規化し、最初の3つのエネルギー状態に対応する波動関数をグラフ化しています。
• plt.plotで各波動関数を位置$ (x) $に対する波動関数の値$ (\psi) $としてプロットしています。
• plt.title、plt.xlabel、plt.ylabelはそれぞれグラフのタイトルと軸ラベルを設定し、plt.legendで凡例を表示しています。
• plt.grid(True)はグリッド線を表示し、plt.show()でグラフを表示しています。

グラフ解説

[実行結果]

このグラフは、1次元の無限長方形井戸内の粒子の波動関数を計算し、その結果を表示したものです。

まず、N個の格子点で井戸の幅が L である領域を考えています。
その後、この領域を離散化するために、差分化を行っています。
具体的には、2階微分を表す行列$ (T) $を作成しています。
この行列は、離散的な位置$ (x)$ における運動エネルギーの演算子（運動量の二乗に比例する部分）を表します。

次に、eigs関数を使用して$ (T) $の固有値と対応する固有ベクトル（波動関数）を計算しています。
ここで k=30 は、最初の$30$個の固有値と固有ベクトルを取得することを意味します。
そして which='SM' は、最小の固有値から計算を行うことを指定しています。

最後に、最初の3つのエネルギー状態に対応する波動関数を取り出し、それぞれを正規化してグラフ化しています。
波動関数は各エネルギー状態における粒子の振る舞いを示し、横軸が位置$ (x) $、縦軸が波動関数$ (\psi) $の値を表しています。

タイトルは「Wavefunctions in a 1D Infinite Square Well」であり、各波動関数のラベルには対応するエネルギー値 $ (E) $が表示されており、それぞれのエネルギー状態が異なるラインで示されています。

ポートフォリオ最適化問題

ポートフォリオ最適化問題の例を、3つの資産（ビットコイン、イーサリアム、米ドル）で考え、3Dグラフ化します。

以下のPythonコードでは、それぞれの資産への投資比率を調整し、ポートフォリオのリターンを最大化するように最適化します。

各投資比率は$0$から$1$の間の値をとり、全体の和は$1$となる制約があります。

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

1

!pip install cvxpy matplotlib numpy

それでは、3Dでのグラフ化を含む以下のコードをご覧ください。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

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

# ビットコイン、イーサリアム、米ドルのリターン
r_btc = 7.0
r_eth = 5.0
r_usd = 2.0

# 投資比率の変数 (0 ≦ x, y ≦ 1)
x = cp.Variable() # ビットコインへの投資比率
y = cp.Variable() # イーサリアムへの投資比率

# 目的関数 - ポートフォリオリターンの最大化
objective = cp.Maximize(r_btc * x + r_eth * y + r_usd * (1 - x - y))

# 制約条件
constraints = [x >= 0, y >= 0, x + y <= 1]

# 最適化問題の定義
prob = cp.Problem(objective, constraints)

# 3Dグラフ描画のための範囲とグリッドを作成
btc_range = np.linspace(0, 1, 100)
eth_range = np.linspace(0, 1, 100)
b, e = np.meshgrid(btc_range, eth_range)
usd_range = 1 - b - e

# 最適化問題を解く
result = prob.solve()

print("最大ポートフォリオリターン: ", result)
print("最適なビットコインへの投資比率: ", x.value)
print("最適なイーサリアムへの投資比率: ", y.value)
print("最適な米ドルへの投資比率: ", 1 - x.value - y.value)

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

# ポートフォリオのリターンをZ軸としてプロット
Z = r_btc * b + r_eth * e + r_usd * usd_range
mask = np.where(usd_range < 0, np.nan, Z)

ax.plot_surface(b, e, mask, rstride=1, cstride=1, color='b', alpha=0.2, zorder=-1)

# 最大化問題の点をプロット
ax.scatter(x.value, y.value, result, color="r", s=100)

ax.set_xlabel('Bitcoin')
ax.set_ylabel('Ethereum')
ax.set_zlabel('Return')

plt.show()

ちなみになりますが、多資産のポートフォリオ最適化は非常に複雑であり、多くの場合、相関関係、リスク等の他の要因を考慮する必要があります。

上記のコードは理論上の模範であり、実際の投資判断には利用しないでください。

[実行結果]

ソースコード解説

このPythonスクリプトは、cvxpyを使用して投資ポートフォリオの最適化問題を解き、3Dグラフを作成しています。

1. モジュールのインポート:

cvxpy、numpy、matplotlibの必要なモジュールをインポートしています。
cvxpyは最適化問題を解くためのライブラリです。

2. リターンの設定:

ビットコイン、イーサリアム、米ドルのリターンをそれぞれ変数r_btc、r_eth、r_usdに設定しています。

3. 変数の定義:

投資比率を示す変数x（ビットコインへの投資比率）とy（イーサリアムへの投資比率）をcvxpyの変数として定義しています。

4. 目的関数と制約条件の定義:

最大化したいポートフォリオリターンを目的関数として定義し、投資比率が$0$以上で、合計が$1$以下であるという制約条件を設定しています。

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

cvxpyを使用して最適化問題を定義し、.solve()メソッドを使って問題を解いています。

6. 最適化結果の出力:

解をresultに格納し、最適な投資比率を表示しています。

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

ビットコインとイーサリアムへの投資比率をX軸とY軸に、ポートフォリオのリターンをZ軸にした3Dグラフを作成しています。
赤い点は最大化されたポートフォリオの点を示しています。

このコードは、ビットコインとイーサリアムの投資比率を最適化し、ポートフォリオのリターンを視覚化しています。

結果解説

[実行結果]

この3Dグラフは、ビットコイン、イーサリアム、米ドルという3つの異なる資産に対する投資比率（X軸はビットコイン、Y軸はイーサリアム）と、それらの投資比率から導かれるポートフォリオのリターン（Z軸）の関係を示しています。

頂上の赤い点は、投資比率を最適化することによって達成可能な最大ポートフォリオリターンを示しています。

これは、ビットコイン、イーサリアム、米ドルへの最適な投資比率を用いた場合のリターンです。

このグラフを見ると、一部の領域が存在しないことに気づくかもしれません。

これは、ビットコイン、イーサリアム、米ドルへの投資比率の合計が100%を超えることはできないためです。

これら3つの資産への投資合計が1（または100%）を超える投資比率の組み合わせは不可能であり、そのためグラフから除外されています。

この3Dグラフはあくまで理論的なもので、実際の市場環境や投資戦略、リスク許容度によって最適なポートフォリオは変動します。

だからと言って、このグラフが無意味であるわけではありません。

次元の異なる複数の資産を持つポートフォリオのリターンを視覚化する一助となり、リスクとリターンのトレードオフを理解する上で有用なツールです。

最適化問題（3Dグラフ）

以下は、Pythonで3Dグラフを用いて最適化問題を解いて可視化する例です。

この例では、CVXPYを使用して、変数 $x$, $y$ の2つを持つ最適化問題を解きます。

目的関数は$ x^2 + y^2 $で、制約条件は$ x + y <= 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

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

# 変数の定義
x = cp.Variable()
y = cp.Variable()

# 目的関数の定義
objective = cp.Minimize(x**2 + y**2)

# 制約条件の定義
constraints = [x + y <= 1]

# 最適化問題の定義
prob = cp.Problem(objective, constraints)

# 最適化問題の解決
result = prob.solve()

# 解を出力
print("Optimal value:", result)
print("Optimal x:", x.value)
print("Optimal y:", y.value)

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

x_vals = np.linspace(-1, 1, 100)
y_vals = np.linspace(-1, 1, 100)
x_mesh, y_mesh = np.meshgrid(x_vals, y_vals)
z = x_mesh**2 + y_mesh**2

ax.plot_surface(x_mesh, y_mesh, z, cmap='viridis', alpha=0.8)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Objective function')

# 最適解をプロット
opt_x = x.value
opt_y = y.value
opt_z = opt_x**2 + opt_y**2
ax.scatter(opt_x, opt_y, opt_z, color='red', s=100, label='Optimal point')

plt.title('Optimization with CVXPY and 3D Plotting')
plt.legend()
plt.show()

このコードでは、$x$と$y$の値を変化させた際の2変数関数$ x^2 + y^2 $を3Dプロットし、最適化された点を赤い点で示しています。

[実行結果]

ソースコード解説

このコードは、CVXPYを使用して最適化問題を解き、その結果を3Dプロットで可視化するものです。

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

ライブラリのインポート

• cvxpy as cp: CVXPYをcpとしてインポートします。
• numpy as np: NumPyをnpとしてインポートします。
• matplotlib.pyplot as plt: Matplotlibのpyplotモジュールをpltとしてインポートします。
• mpl_toolkits.mplot3d: Matplotlibの3Dプロットを利用するためのモジュールです。

変数の定義

• xとyはCVXPYの変数として定義されます。
  これらの変数は最適化問題の変数として使用されます。

目的関数の定義

• objectiveは、最小化する目的関数を表しています。
  ここでは、$ (x^2 + y^2) $を最小化するように設定されています。

制約条件の定義

• constraintsは、制約条件のリストです。
  この場合、$ (x + y \leq 1) $という制約条件が定義されています。

最適化問題の定義

• probは、最適化問題を表します。
  objectiveとconstraintsを引数として渡して最適化問題を定義します。

最適化問題の解決

• prob.solve()は、定義した最適化問題を解きます。
  最適な変数の値を見つけ、最小値を求めます。

解の出力

• 求めた最適な変数の値と最適値を表示します。

3Dプロット

• plt.figure(): 新しい図を作成します。
• fig.add_subplot(111, projection='3d'): 3Dのサブプロットを追加します。
• x_valsとy_valsは、それぞれ$-1$から$1$までの$100個$の等間隔の値を持つ配列です。
  これを使ってメッシュグリッドを生成します。
• x_meshとy_meshはx_valsとy_valsを使って作成された格子状のデータです。
  それらに基づいて目的関数の値zを計算します。
• ax.plot_surface(): メッシュグリッドを使用して3Dの曲面をプロットします。
  cmapはカラーマップ、alphaは透明度を指定します。
• ax.set_xlabel()、ax.set_ylabel()、ax.set_zlabel(): 各軸のラベルを設定します。

最適解のプロット

• ax.scatter(): 最適解を赤い点でプロットします。
  これは、最適解の$ (x) $と$ (y) $の値に対応する目的関数の値を表します。

グラフの表示

• plt.title(): グラフのタイトルを設定します。
• plt.legend(): 凡例を追加します。
• plt.show(): グラフを表示します。

結果解説

[実行結果]

このコードは、CVXPYを使用して2変数の最適化問題を解き、その結果を3Dグラフで可視化しています。

最適化問題は、目的関数 $ (x^2 + y^2) $を最小化するというもので、制約条件は$ (x + y \leq 1) $です。
これは、2変数$ (x) $と$ (y) $の値を持つ平面上で$ (x^2 + y^2) $を最小化する問題です。

プログラムでは、変数$ (x) $と$ (y) $を定義し、目的関数および制約条件を設定しました。
その後、CVXPYを使って最適化問題を解決し、最適解の$ (x) $と$ (y) $の値を得ています。

3Dグラフでは、平面上の$ (x) $と$ (y) $の値を範囲指定してメッシュグリッドを生成し、それに対応する$ (x^2 + y^2) $の値を計算してプロットしています。
これにより、2変数関数$ (x^2 + y^2) $の3D表面が描かれます。

最適解は、求めた$ (x) $と$ (y) $の値を赤い点で表現しています。
これは、目的関数を最小化する最適な $ (x) $と$ (y) $の値を示しています。

つまり、このグラフは、2変数関数を3Dで可視化し、最適解を視覚的に示しています。

ベッセル関数

ベッセル関数は物理や数学で幅広く利用される特殊関数の一つです。

ベッセル関数の数式として、第一種ベッセル関数 $ ( J_n(x) ) $の例を示します：

$$
J_n(x) = \frac{1}{\pi} \int_{0}^{\pi} \cos(n\theta - x \sin \theta) , d\theta
$$

このベッセル関数をPythonでグラフ化するために、SciPyライブラリのベッセル関数を利用します。

以下がコードになります：

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn

x = np.linspace(0, 20, 1000)  # xの範囲を設定

# 第一種ベッセル関数の計算（n=0, 1, 2）
for n in range(3):
    y = jn(n, x)  # ベッセル関数の計算
    plt.plot(x, y, label=f'J{n}(x)')  # グラフのプロット

plt.title('Bessel Functions')
plt.xlabel('x')
plt.ylabel('Jn(x)')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、scipy.special.jn を使用してベッセル関数を計算し、それをMatplotlibを使ってグラフ化しています。

[実行結果]

ソースコード解説

このコードは、Pythonでベッセル関数（Bessel Functions）を計算し、Matplotlibを使用してグラフ化するものです。

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

• numpy を np としてインポートし、行列計算や数値演算を行うための機能を利用します。
• matplotlib.pyplot を plt としてインポートし、グラフ描画のための機能を利用します。

2. SciPyからベッセル関数のインポート:

• scipy.special モジュールから jn 関数をインポートしています。
  これはベッセル関数を計算するための関数です。

3. データの準備:

• np.linspace(0, 20, 1000) によって、$0$から$20$の範囲を等間隔に$1000$点取得し、x に格納します。

4. ベッセル関数の計算とプロット:

• for ループで range(3) を使用し、n が0から2まで変化します。
• jn(n, x) を使って、x に対する第一種ベッセル関数（n=0, 1, 2）を計算します。
• plt.plot(x, y, label=f'J{n}(x)') により、各 n のベッセル関数をプロットします。label を使って凡例を設定しています。

5. グラフの装飾:

• plt.title、plt.xlabel、plt.ylabel を使って、タイトルや軸ラベルを設定します。
• plt.legend() で凡例を表示します。
• plt.grid(True) でグリッドを表示します。

6. グラフの表示:

• plt.show() により、グラフを表示します。

このコードは、ベッセル関数を計算し、それをグラフ化しています。

それぞれのベッセル関数の振る舞いを可視化し、それが x の値に対してどのように変化するかを示しています。

グラフ解説

[実行結果]

このグラフは、第一種ベッセル関数（Bessel Function of the First Kind）の値を示しています。
第一種ベッセル関数は微分方程式理論や物理学、工学など幅広い分野で重要な役割を果たす特殊関数です。

解析的な式で表されるこの関数は、x（横軸）と Jn(x)（縦軸）の関係を示しており、このグラフでは n=0、n=1、n=2 の3つのベッセル関数が描かれています。

• n=0のベッセル関数（青の線）は、xが$0$のときに$1$に近づき、その後減少していきます
  xが大きくなるにつれて振動が小さくなり、$0$に収束します。
• n=1のベッセル関数（オレンジの線）は、xが$0$のときに$0$になり、xが増加するにつれて振動して$±1$の間を移動します。
• n=2のベッセル関数（緑の線）は、xが0のときに$0$になり、xが増加するにつれて振動が増幅し、2回の極小値と1回の極大値を持ちます。

xが大きくなるほど、ベッセル関数の振る舞いが規則的になり、特定のパターンに収束していく様子が観察されます。

このグラフは、数学的な特殊関数の性質を視覚的に理解するのに役立ちます。

乳がんの診断（ランダムフォレストモデル）

乳がんの診断（良性、悪性）を予測するためのランダムフォレストモデルを作ります。

scikit-learnの乳がんデータセットを利用します。

データセットの分割とランダムフォレストモデルの訓練を行い、最後にテストデータを使って予測の精度を計測・表示します。

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

from sklearn.datasets import load_breast_cancer
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

# データの読み込み
data = load_breast_cancer()
X = data.data
y = data.target

# データの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# モデルの訓練
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)

# テストデータでの精度
accuracy = clf.score(X_test, y_test)
print(f'Model accuracy: {accuracy*100:.2f}%')

[実行結果]

Model accuracy: 96.49%

次に、診断の重要度を評価し、それを可視化します。

この可視化は、どの特徴が乳がんの診断に最も重要であるかを理解するのに役立ちます。

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

import numpy as np

# 特徴量の重要度
importances = clf.feature_importances_

# 重要度が高い順にソート
indices = np.argsort(importances)[::-1]

# 特徴量の名前をソートした順に並べ替え
names = [data.feature_names[i] for i in indices]

# プロットの作成
plt.figure(figsize=(10, 5))
plt.title("Feature Importance")
plt.bar(range(X.shape[1]), importances[indices])
plt.xticks(range(X.shape[1]), names, rotation=90)

plt.show()

[実行結果]

最後に、テストデータの正解ラベルと予測ラベルを比較します。

これによりモデルがどの程度のパフォーマンスを発揮しているかを視覚的に理解することができます。

1
2
3
4
5
6
7
8
9

# テストデータの予測
y_pred = clf.predict(X_test)

# 結果の表示
plt.figure(figsize=(15, 5))
plt.plot(y_test, 'o', markersize=10)
plt.plot(y_pred, 'v', markersize=10)
plt.legend(['True', 'Predicted'])
plt.show()

[実行結果]

以上のステップで、乳がん診断のためのランダムフォレストモデルの学習と評価を行い、結果を視覚化しました。

ソースコード解説

このコードは乳がんデータセットを使用し、ランダムフォレストを使って腫瘍が良性か悪性かを分類する機械学習モデルを訓練しています。

1. データの読み込みと分割:

• load_breast_cancerで乳がんデータセットを読み込み、特徴量（X）とターゲット（y）を取得します。
• train_test_splitを使ってデータをトレーニング用とテスト用に分割します。

2. モデルの訓練:

• ランダムフォレスト分類器（RandomForestClassifier）を使って、トレーニングデータでモデルを訓練します。

3. モデルの評価:

• テストデータで訓練されたモデルの性能を評価し、精度（accuracy）を表示します。

4. 特徴量の重要度のプロット:

• feature_importances_を使用して、各特徴量の重要度を取得します。
• 重要度の高い特徴量を降順にソートし、特徴量の名前とともにバーのプロットを作成します。

5. 予測結果の可視化:

• テストデータでの実際のターゲット値（y_test）と予測値（y_pred）をプロットして比較します。

このコードは、ランダムフォレストを使った分類タスクの全体的な流れを示しています。

データの読み込みから特徴量の重要度の可視化まで、モデル構築から評価までの一連のステップを示しています。

結果解説

[実行結果]

このグラフは、ランダムフォレストクラシフィケーションモデルを用いて特徴量の重要度を視覚化したものです。

グラフの各バーは一つひとつの特徴量を表し、その長さ（高さ）はその特徴量のランダムフォレスト内での相対的な重要度を示しています。

特徴量の重要度は、その特徴が何回スプリットに利用され、それらの分割がどれだけ不純物（通常はジニ不純度またはエントロピー）を減らすのに寄与したかを基づいて計算されます。

x軸上には特徴量の名前が表示され、これらはその重要度に基づいて左から右へと並べられています。

最も重要な特徴量が左側に表示され、重要度が低いものが右側へと続きます。

ばらつきの大きい特徴量は、対象とする問題において重要な役割を果たすことが多いです。

これらの特徴量を特定することで、データの理解が深まり、より効果的な特徴エンジニアリング戦略を考える上でも有効です。

なお、このソースコードは乳がんの診断に関するデータセットを対象としているため、特徴量は乳がんの診断に関連する一連の医学的測定値を指します。

例えば、細胞の形やサイズ、細胞核の形状などがこれに該当します。

[実行結果]

このグラフは、予測結果と実際の結果を視覚化したものです。

テストデータに対して予測がどのように行われ、その結果がどの程度実際の目標値に一致しているかを判断するためのものです。

“True”ライン：

これはy_testの値をプロットしたもので、テストデータの真のラベル（クラス）がどのように分布しているかを示しています。
これは「正解」の分布を示す基準線となります。

“Predicted”ライン：

これはy_predの値をプロットしたもので、モデルがテストデータをどのように予測したかを示します。
これはモデルの「予測」の分布を示しており、こちらが“True”ラインとどの程度一致しているかが、モデルの性能を評価するための主要な視覚的指標となります。

点のy軸の位置はクラスラベルを示し、x軸の位置は各データポイントを示しています。

ここでは二項分類問題を扱っており、各クラス（良性、悪性）はそれぞれ$0$と$1$でラベル付けされています。

したがって、予測が真のラベルと一致する場合、“True”点と”Predicted”の点は同じy座標に配置されます。

逆に、予測が真のラベルと一致しない場合、“True”点と”Predicted”点は異なるy座標に存在します。

このタイプのグラフはモデルの具体的なミスを視覚的に把握するのに役立ちます。

例えば、一部のデータポイントで予測が一貫して外れている場合、それはそのデータポイントが何らかの特定のパターンを有しており、それがモデルにとって難しいか、または学習データ内で不足している可能性を示しています。

温度による酵素の活性

生物学的な例題として、温度による酵素の活性 を挙げます。

酵素は、生物の体内で化学反応を促進するタンパク質です。

酵素の活性は、温度によって変化します。

一般的に、酵素の活性は、温度が上昇するにつれて増加し、ある温度でピークを迎え、その後は減少します。

この現象をPythonで解いて、グラフ化するには、以下のコードを使用します。

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

import numpy as np
import matplotlib.pyplot as plt

# 酵素の活性の定義
def enzyme_activity(temperature):
    if temperature < 20:
        return 0
    elif temperature <= 30:
        return temperature - 20
    else:
        return 100 - (temperature - 30)**2

# グラフの表示
temperatures = np.linspace(10, 40, 100)
activities = [enzyme_activity(temperature) for temperature in temperatures]
plt.plot(temperatures, activities)
plt.xlabel('温度(℃)')
plt.ylabel('酵素の活性')
plt.show()

このコードを実行すると、以下のグラフが表示されます。

[実行結果]

このグラフは、酵素の活性を温度の関数として示しています。

温度が20℃から30℃の範囲では、酵素の活性は増加していますが、温度が30℃を超えると、酵素の活性は減少しています。

このように、Pythonを使えば、生物学的な現象を解いて、グラフ化することができます。

この例では、酵素の活性のモデルとして、単純な線形モデルを使用しています。

実際の酵素の活性は、温度だけでなく、pHやイオン濃度などの条件によっても変化します。

ソースコード解説

このコードは、酵素の活性を温度に基づいてモデル化し、その関係をグラフで可視化しています。

1. import文では、NumPyとMatplotlibライブラリをインポートしています。
   これらは数値計算やグラフ描画のために使用されます。

2. enzyme_activity関数は、温度を引数として受け取り、その温度に基づいて酵素の活性を計算します。

• 温度が20℃未満の場合、活性は$0$となります。
• 温度が20℃以上30℃以下の場合、活性は温度から20を引いた値となります。
• 温度が30℃より高い場合、活性は100から(温度から30を引いた値)の二乗を引いた値となります。

3. temperaturesは、10℃から40℃までの温度を100分割した配列を作成しています。

4. activitiesは、enzyme_activity関数を使って各温度に対する活性を計算し、リストに格納しています。

5. plt.plot()は、温度に対する酵素活性のグラフを作成します。
   x軸には温度(℃)、y軸には酵素の活性が表示されます。

6. plt.xlabel()とplt.ylabel()は、x軸とy軸のラベルを設定しています。

7. plt.show()は、グラフを表示します。

このコード全体は、温度によって酵素の活性がどのように変化するかを表したシンプルなグラフを生成しています。