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

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

ナーシム–シュレーダー方程式(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 について、ジュリア集合をプロットします。
まず、xy の範囲を設定し、それらの組み合わせについて 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$が使われています。
  • ジュリア集合自体は非常に複雑な構造を持ち、フラクタル的な自己相似性が見られます。
    つまり、拡大するとほぼ同じ模様が現れる性質があります。
  • 白い領域(ジュリア集合)の周りに同心円状の渦巻き模様があり、そこから外に向かうにつれて発散速度が速くなっていきます。
  • グラフ右側のカラーバーは、発散速度(反復回数)と色の対応関係を示しています。

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

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

シンプソンのパラドックス

シンプソンのパラドックス

シンプソンのパラドックスは、データを解釈する際に生じる混乱を示す統計学上のパラドックスです

具体的には、集団全体部分集団の間で、同じデータを用いて異なる結論が導かれることを示します。
このパラドックスは、条件付き確率の誤解やサンプルサイズの影響に関連しています。

例えば、ある特定の属性を持つ集団を部分集団として考えた場合、その属性に関する統計的傾向が逆転することがあります。
つまり、部分集団での割合が高いのに、全体集団での割合が低くなる場合があります。

このパラドックスは、統計的推論データ解釈において、慎重な考慮が必要であることを強調します。

サンプルソース

シンプソンのパラドックスは、3つの異なる群のデータを結合した場合に、それぞれのグループの平均値が異なるにも関わらず、結合した全体の平均値が逆にそれらの平均値よりも高い場合を指します。

以下に、シンプソンのパラドックスを示す簡単な例を用いてPythonでグラフ化する方法を示します。

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

# グループ1: 平均値=5, 標準偏差=1, サンプル数=100
group1_mean = 5
group1_std = 1
group1_samples = np.random.normal(group1_mean, group1_std, 100)

# グループ2: 平均値=7, 標準偏差=1, サンプル数=100
group2_mean = 7
group2_std = 1
group2_samples = np.random.normal(group2_mean, group2_std, 100)

# グループ3: 平均値=9, 標準偏差=1, サンプル数=100
group3_mean = 9
group3_std = 1
group3_samples = np.random.normal(group3_mean, group3_std, 100)

# 結合したデータ
combined_data = np.concatenate((group1_samples, group2_samples, group3_samples))

# プロット
plt.figure(figsize=(10, 6))

# グループ1のヒストグラム
plt.hist(group1_samples, bins=20, alpha=0.5, label='Group 1')

# グループ2のヒストグラム
plt.hist(group2_samples, bins=20, alpha=0.5, label='Group 2')

# グループ3のヒストグラム
plt.hist(group3_samples, bins=20, alpha=0.5, label='Group 3')

# 結合したデータのヒストグラム
plt.hist(combined_data, bins=20, alpha=0.5, label='Combined', color='black')

plt.xlabel('Value')
plt.ylabel('Frequency')
plt.title('Simpson\'s Paradox Example')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、3つの異なる群のデータを生成し、それぞれをヒストグラムとしてプロットします。

結合した全体のデータもヒストグラムとしてプロットされ、シンプソンのパラドックスを可視化します。

[実行結果]

ソースコード解説

ソースコードの部分について詳細に説明します。

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

1
2
import numpy as np
import matplotlib.pyplot as plt
  • numpy は数値計算を行うためのPythonライブラリで、np としてインポートされています。
  • matplotlib.pyplot はグラフを描画するためのPythonライブラリで、plt としてインポートされています。

2. グループのデータ生成:

1
2
3
group1_mean = 5
group1_std = 1
group1_samples = np.random.normal(group1_mean, group1_std, 100)
  • group1_meangroup1_std はグループ1のデータの平均値標準偏差を表します。
  • np.random.normal() 関数を使用して、平均値group1_mean標準偏差group1_std正規分布に従うデータを生成します。
  • これにより、グループ1のデータが group1_samples として格納されます。
    $100$個のサンプルが生成されます。

3. 同様に、グループ2とグループ3のデータも生成されます。

4. データの結合:

1
combined_data = np.concatenate((group1_samples, group2_samples, group3_samples))
  • np.concatenate() 関数を使用して、3つのグループのデータを結合して1つの配列にします。
  • これにより、全体のデータが combined_data として格納されます。

5. プロット:

1
2
3
4
5
6
7
8
9
10
11
plt.figure(figsize=(10, 6))
plt.hist(group1_samples, bins=20, alpha=0.5, label='Group 1')
plt.hist(group2_samples, bins=20, alpha=0.5, label='Group 2')
plt.hist(group3_samples, bins=20, alpha=0.5, label='Group 3')
plt.hist(combined_data, bins=20, alpha=0.5, label='Combined', color='black')
plt.xlabel('Value')
plt.ylabel('Frequency')
plt.title('Simpson\'s Paradox Example')
plt.legend()
plt.grid(True)
plt.show()
  • plt.figure(figsize=(10, 6)) は、図のサイズを設定します。
  • plt.hist() 関数は、ヒストグラムを描画します。
    alpha パラメータは透明度を設定し、label パラメータは凡例のラベルを設定します。
  • plt.xlabel()plt.ylabel() は、x軸とy軸のラベルを設定します。
  • plt.title() は図のタイトルを設定します。
  • plt.legend() は、凡例を表示します。
  • plt.grid(True) は、グリッド線を表示します。
  • plt.show() は、図を表示します。

このコードは、3つの異なるグループのデータを生成し、それらをヒストグラムとしてプロットしています。

さらに、結合したデータのヒストグラムもプロットされ、シンプソンのパラドックスが可視化されています。

結果解説

[実行結果]

このグラフでは、シンプソンのパラドックスの例を示しています。

各グループのヒストグラム:

  • “Group 1”, “Group 2”, “Group 3” というラベルが付けられた3つのヒストグラムが表示されています。
  • 各グループは、それぞれ異なる平均値標準偏差を持ちます。
  • 各グループのヒストグラムは、それぞれのグループが持つデータの分布を示しています。
    ヒストグラムは、データの値の範囲をいくつかの階級(ビン)に分割し、各階級に含まれるデータの頻度(個数)を表示します。

結合したデータのヒストグラム:

  • “Combined” というラベルが付けられた黒いヒストグラムが表示されています。
  • このヒストグラムは、3つの異なるグループからのデータを結合した全体のデータの分布を示しています。
  • 結合したデータのヒストグラムは、各グループのデータの組み合わせから得られる全体のデータの分布を示しています。

シンプソンのパラドックスの本質は、個々のグループの平均値が異なる方向に移動する一方で、結合したデータの平均値がそれらの平均値よりも大きくなることです。

これは、データが持つ特定の相関構造関係性によって引き起こされることがあります。

グロース・ピタエフスキー(Gross-Pitaevskii)方程式

グロース・ピタエフスキー(Gross-Pitaevskii)方程式

グロース・ピタエフスキー(Gross-Pitaevskii)方程式は、ボース・アインシュタイン凝縮体の基底状態を記述する非線形シュレディンガー方程式です。

ここでは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
from scipy.fftpack import fft2, ifft2

# パラメータの設定
N = 128 # グリッドの次元
L = 10.0 # 領域の大きさ
dx = L / N # 空間グリッドの間隔
dt = 0.001 # 時間ステップ
g = 1.0 # 非線形係数
Nt = 1000 # 時間ステップ数

# 初期条件
x = np.arange(-N / 2, N / 2) * dx
X, Y = np.meshgrid(x, x)
psi0 = np.exp(-(X ** 2 + Y ** 2) / 2) / np.sqrt(2 * np.pi)

# 運動方程式の計算
psi = psi0.copy()
for n in range(Nt):
psi_fft = fft2(psi)
psi_fft *= np.exp(-0.5j * dt * (4 * np.pi ** 2 / L ** 2) * (np.fft.fftfreq(N, d=dx) ** 2 + np.fft.fftfreq(N, d=dx).reshape(-1, 1) ** 2))
psi_fft -= 1j * dt * g * np.abs(psi) ** 2 * psi_fft
psi = ifft2(psi_fft)

# グラフ化
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(np.abs(psi) ** 2, cmap='coolwarm', extent=[-L / 2, L / 2, -L / 2, L / 2])
ax.set_xlabel('x', fontsize=14)
ax.set_ylabel('y', fontsize=14)
ax.set_title('Gross-Pitaevskii Equation', fontsize=16)
fig.colorbar(im, ax=ax)
plt.show()

このコードでは、まずパラメータと初期条件を設定します。
初期条件としてガウス分布を与えています。

次に、時間発展の計算ループに入ります。
各ステップで、まず波動関数のフーリエ変換を計算し、運動エネルギー項をフーリエ空間で適用します。
続いて、非線形項を実空間で適用します。
最後に逆フーリエ変換を行い、新しい波動関数を得ます。

計算が終了したら、波動関数の確率密度を2次元プロットで可視化しています。
初期のガウス分布が、非線形項の影響で変形していく様子が分かります。

グロース・ピタエフスキー方程式は、ボース・アインシュタイン凝縮体の基底状態を記述する重要な方程式です。
超流動性量子渦などの興味深い現象を記述できます。

[実行結果]

ソースコード解説

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

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

1
2
3
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft2, ifft2
  • NumPy: 数値計算のための基本ライブラリ
  • Matplotlib: グラフィックスライブラリ
  • scipy.fftpack: フーリエ変換のためのサブライブラリ

2. パラメータの設定

1
2
3
4
5
6
N = 128  # グリッドの次元
L = 10.0 # 領域の大きさ
dx = L / N # 空間グリッドの間隔
dt = 0.001 # 時間ステップ
g = 1.0 # 非線形係数
Nt = 1000 # 時間ステップ数
  • $N$: 空間グリッドの次元数
  • $L$: 計算領域の大きさ
  • $dx$: 空間グリッドの間隔
  • $dt$: 時間ステップ幅
  • $g$: グロース・ピタエフスキー方程式の非線形係数
  • $Nt$: 時間発展の計算ステップ数

3. 初期条件の設定

1
2
3
x = np.arange(-N / 2, N / 2) * dx
X, Y = np.meshgrid(x, x)
psi0 = np.exp(-(X ** 2 + Y ** 2) / 2) / np.sqrt(2 * np.pi)
  • $x$: 1次元の空間グリッド
  • $X, Y$: 2次元の空間グリッド
  • $psi0$: 初期波動関数(ガウス分布)

4. 運動方程式の時間発展計算

1
2
3
4
5
6
psi = psi0.copy()
for n in range(Nt):
psi_fft = fft2(psi)
psi_fft *= np.exp(-0.5j * dt * (4 * np.pi ** 2 / L ** 2) * (np.fft.fftfreq(N, d=dx) ** 2 + np.fft.fftfreq(N, d=dx).reshape(-1, 1) ** 2))
psi_fft -= 1j * dt * g * np.abs(psi) ** 2 * psi_fft
psi = ifft2(psi_fft)
  • $psi$0から$psi$にコピー
  • 時間ループ
    • 波動関数のフーリエ変換
    • 運動エネルギー項の適用(フーリエ空間)
    • 非線形ポテンシャル項の適用(実空間)
    • 逆フーリエ変換で新しい波動関数を得る

5. 結果の可視化

1
2
3
4
5
6
7
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(np.abs(psi) ** 2, cmap='coolwarm', extent=[-L / 2, L / 2, -L / 2, L / 2])
ax.set_xlabel('x', fontsize=14)
ax.set_ylabel('y', fontsize=14)
ax.set_title('Gross-Pitaevskii Equation', fontsize=16)
fig.colorbar(im, ax=ax)
plt.show()
  • 図とプロット領域の設定
  • 確率密度の2次元プロット
  • 軸ラベル、タイトル設定
  • カラーバーの追加
  • グラフ表示

このコードは、グロース・ピタエフスキー方程式の数値解の時間発展を計算し、確率密度分布を可視化しています。

主な手順は、(1)初期条件と計算パラメータの設定、(2)時間発展の計算ループ、(3)結果の可視化、となっています。

計算では、フーリエ変換と実空間での演算を組み合わせて行われています。

結果解説

[実行結果]

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

このグラフは、グロース・ピタエフスキー方程式の数値解の確率密度分布を可視化したものです。

  • $x軸$と$y軸$はそれぞれ空間座標を表しており、単位は長さの次元を持っています。
    軸のスケールは$-5$から$5$までの範囲を示しています。

  • 色の濃淡は確率密度の大きさを表しています。
    赤が最も高密度、青が最も低密度の領域を示します。

  • 初期条件としてガウス分布を与えているため、時刻$t=0$では中心に赤い高密度領域があり、周りに低密度の青い領域が広がっています。

  • 時間が経過するにつれて、非線形項の影響で確率密度分布が変形しています。
    中心付近の高密度領域が縮んだり、複数の高密度領域(ピーク)が現れたりしています。

  • 非線形シュレディンガー方程式の解は一般に複雑な構造を示します。
    この解の時間発展では、確率密度が連続的に変形しながら複雑な模様を描いています。

  • グラフの右側にカラーバーがあり、確率密度の大きさと色の対応関係が示されています。

このように、グロース・ピタエフスキー方程式の数値解を可視化することで、ボース・アインシュタイン凝縮体の基底状態がどのように時間発展するかを確認できます。

特に非線形項の効果が顕著に現れている様子が分かります。

ブラック・ショールズ方程式

ブラック・ショールズ方程式

ブラック・ショールズ方程式は、オプション価格を計算する際に使用される重要な式です。

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

# パラメータの設定
S0 = 100 # 現在の株価
K = 105 # 権利行使価格
T = 1.0 # 残存期間(年)
r = 0.05 # 無リスク利子率
sigma = 0.2 # ボラティリティ

# ブラック・ショールズ方程式
def black_scholes(S0, K, T, r, sigma):
from scipy.stats import norm
d1 = (np.log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
call = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
return call

# グラフ化するための株価の範囲
S = np.arange(60, 140, 1)

# コールオプション価格を計算
call_values = [black_scholes(s, K, T, r, sigma) for s in S]

# グラフの描画
plt.figure(figsize=(8, 6))
plt.plot(S, call_values)
plt.xlabel('Stock Price')
plt.ylabel('Call Option Price')
plt.title('Black-Scholes Option Pricing')
plt.grid(True)
plt.show()

このコードでは、以下の手順を実行しています。

  1. 必要なライブラリ(numpy、matplotlib、scipy)をインポートします。
  2. オプションのパラメータ(株価、権利行使価格、残存期間、無リスク利子率、ボラティリティ)を設定します。
  3. ブラック・ショールズ方程式を実装した関数black_scholesを定義します。
  4. グラフ化する株価の範囲を設定します。
  5. その株価範囲に対してコールオプション価格を計算します。
  6. matplotlibを使ってグラフを描画します。

実行すると、横軸が株価、縦軸がコールオプション価格のグラフが表示されます。

[実行結果]

このグラフを見ることで、株価オプション価格の関係を視覚的に確認できます。

ソースコード解説

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

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

1
2
import numpy as np
import matplotlib.pyplot as plt

このコードでは、NumPyMatplotlibをインポートしています。
NumPyは数値計算に使用され、Matplotlibはグラフィックスを描画するために使用されます。

2. パラメータの設定

1
2
3
4
5
S0 = 100  # 現在の株価
K = 105 # 権利行使価格
T = 1.0 # 残存期間(年)
r = 0.05 # 無リスク利子率
sigma = 0.2 # ボラティリティ

ここでは、ブラック・ショールズ方程式に必要なパラメータを設定しています。

  • S0現在の株価
  • K権利行使価格
  • Tオプションの残存期間(年数)
  • r無リスク利子率
  • sigmaボラティリティ(株価の変動率)

3. ブラック・ショールズ方程式の実装

1
2
3
4
5
6
def black_scholes(S0, K, T, r, sigma):
from scipy.stats import norm
d1 = (np.log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
call = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
return call

この関数は、ブラック・ショールズ方程式を実装しています。

引数として、株価、権利行使価格、残存期間、無リスク利子率、ボラティリティを受け取ります。

関数内部では、d1d2を計算し、それらを使ってコールオプション価格を計算しています。

4. グラフ化するための株価の範囲の設定

1
S = np.arange(60, 140, 1)

このコードは、グラフ化する株価の範囲を設定しています。

np.arange(60, 140, 1)は、$60$から$139$までの整数を$1$つ飛ばしで生成します。

5. コールオプション価格の計算

1
call_values = [black_scholes(s, K, T, r, sigma) for s in S]

この行では、リストコンプリヘンションを使って、Sのそれぞれの値に対してコールオプション価格を計算しています。

black_scholes関数を呼び出し、その結果をリストに格納しています。

6. グラフの描画

1
2
3
4
5
6
7
plt.figure(figsize=(8, 6))
plt.plot(S, call_values)
plt.xlabel('Stock Price')
plt.ylabel('Call Option Price')
plt.title('Black-Scholes Option Pricing')
plt.grid(True)
plt.show()

最後に、Matplotlibを使ってグラフを描画しています。

  • plt.figure(figsize=(8, 6))は、グラフのウィンドウサイズを設定しています。
  • plt.plot(S, call_values)は、株価(S)とコールオプション価格(call_values)をプロットしています。
  • plt.xlabelplt.ylabelは、$x軸$と$y軸$のラベルを設定しています。
  • plt.titleはグラフのタイトルを設定しています。
  • plt.grid(True)は、グリッド線を表示します。
  • plt.show()は、グラフを表示します。

このコードを実行すると、株価コールオプション価格の関係を示すグラフが表示されます。

結果解説

[実行結果]

このグラフは、株価(横軸)コールオプション価格(縦軸)の関係を示しています。

横軸は、株価の範囲を表しています。
このコードでは、株価が$60$から$139$までの値をプロットしています。

縦軸は、その株価に対応するコールオプション価格を表しています。
コールオプション価格は、ブラック・ショールズ方程式を用いて計算されています。

グラフ上の曲線は、株価コールオプション価格の関係を示しています。
曲線の形状から、以下の特徴が読み取れます。

1. 原点近くでは緩やかな上昇 :

株価が権利行使価格を下回る場合、オプションは価値がほとんどありません。
そのため、曲線の始まりは緩やかです。

2. 権利行使価格付近で急な上昇 :

株価が権利行使価格に近づくにつれ、オプションの価値が急激に上昇します。
これは、オプションの値が内在価値(権利行使価格と株価の差)に依存するためです。

3. 権利行使価格を超えた後は緩やかな上昇 :

株価が権利行使価格を上回ると、オプションの価値は時間的価値のみに依存するため、曲線の傾きは緩やかになります。

このようにグラフを見ることで、株価の変化に対するオプション価格の挙動を視覚的に把握できます。

オプション取引においては、この関係を理解することが重要です。

アトキンソン・シフト方程式

アトキンソン・シフト方程式

アトキンソン・シフト方程式は、以下のように表されます。

$$
y’ = y(1 - y)
$$

ここで、$y$はある生物の個体数の割合を表します。

この方程式は、初期値$y(0)=y0$を与えることで解くことができます。

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 atkinson(y0, t):
'''
アトキンソン・シフト方程式を解く関数
y0: 初期値
t: 時間
'''
return y0 * np.exp(t) / (y0 * (np.exp(t) - 1) + 1)

# 時間の範囲を設定
t = np.linspace(0, 10, 1000)

# 様々な初期値でプロット
y0_values = [0.1, 0.3, 0.5, 0.7, 0.9]
for y0 in y0_values:
y = atkinson(y0, t)
plt.plot(t, y, label=f'y0 = {y0}')

plt.xlabel('Time')
plt.ylabel('Population Fraction')
plt.title('Atkinson Shift Equation')
plt.legend()
plt.show()

このコードでは、まずアトキンソン・シフト方程式を解析的に解く関数atkinsonを定義しています。

次に、時間の範囲を設定し、様々な初期値$y0$で方程式を解いて可視化しています。

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

[実行結果]

横軸が時間、縦軸が個体数の割合を表しています。

初期値が異なると、個体数の変化の様子が大きく変わることがわかります。

このグラフから、アトキンソン・シフト方程式は、一定の割合で飽和する振る舞いを示すことがわかります。

ソースコード解説

このソースコードは以下の構成になっています。

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

1
2
import numpy as np
import matplotlib.pyplot as plt

NumPyMatplotlibライブラリをインポートしています。
NumPyは数値計算、Matplotlibはグラフ描画に使用されます。

2. 関数の定義

1
2
3
4
5
6
7
def atkinson(y0, t):
'''
アトキンソン・シフト方程式を解く関数
y0: 初期値
t: 時間
'''
return y0 * np.exp(t) / (y0 * (np.exp(t) - 1) + 1)

アトキンソン・シフト方程式 $y’ = y(1 - y)$ を解析的に解く関数 atkinson を定義しています。

引数は初期値 $y0$ と時間 $t$ です。
返り値は、その時間$t$における方程式の解の値になります。

3. 時間範囲の設定

1
t = np.linspace(0, 10, 1000)

時間の範囲を$0$から$10$まで$1000$点に分けた等間隔の値を作成しています。

この値を使って解の挙動をプロットします。

4. 初期値のリストの作成

1
y0_values = [0.1, 0.3, 0.5, 0.7, 0.9]

プロットする初期値のリストを作成しています。

この例では$0.1$, $0.3$, $0.5$, $0.7$, $0.9$の$5$つの初期値で解の挙動をプロットします。

5. プロットのループ

1
2
3
for y0 in y0_values:
y = atkinson(y0, t)
plt.plot(t, y, label=f'y0 = {y0}')

初期値のリストから$1$つずつ初期値を取り出し、対応する解の値を計算して、Matplotlibでプロットしています。

引数labelは凡例の値です。

6. 軸ラベル、タイトル、凡例の設定

1
2
3
4
plt.xlabel('Time')
plt.ylabel('Population Fraction')
plt.title('Atkinson Shift Equation')
plt.legend()

グラフの横軸ラベル、縦軸ラベル、タイトル、凡例を設定しています。

7. グラフの表示

1
plt.show()

作成したグラフを表示する命令です。

このコードではまず、アトキンソン・シフト方程式を解析的に解く関数を定義し、次にその関数を使って様々な初期値における解の挙動をプロットしています。

最後にグラフの体裁を整え、表示しています。

グラフ解説

[実行結果]

このグラフは、アトキンソン・シフト方程式の解の挙動を示しています。

  • 横軸は時間 $t$を表しています。
    時間範囲は$0$から$10$まで等間隔で$1000$点プロットされています。
  • 縦軸は個体数の割合 $y$を表しています。
    個体数の割合は$0$から$1$の範囲にあります。
  • グラフには5本の曲線が描かれています。
    それぞれ初期値$y0$が$0.1$、$0.3$、$0.5$、$0.7$、$0.9$のときの解の挙動を表しています。
  • 曲線の色は異なり、凡例によって各曲線に対応する初期値が示されています。

これらの曲線の形状から、以下の特徴が読み取れます:

  1. すべての解曲線は滑らかな単調増加関数となっています。
  2. 初期値$y0$が大きいほど、曲線の立ち上がりが早くなっています。
  3. すべての曲線は時間がたつにつれて$y=1$に漸近していきます。
    つまり、十分時間が経てば個体数の割合は1 (100%) に収束します。
  4. 初期値が$0.5$のときは、曲線は$y=0.5$を通過しています。
    初期値が$0.5$より大きい(小さい)ときは、曲線は常に$y>0.5 (y<0.5)$となっています。

このようにアトキンソン・シフト方程式は、初期値に応じて個体数の割合が徐々に増加し、最終的には100%に収束する振る舞いを記述しています。

生物の個体群の成長をモデル化する際に使われる方程式です。

グラフを見ることで、異なる初期値での解の挙動の違いを視覚的に捉えることができ、この方程式の性質を理解しやすくなっています。

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

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

レナード・ジョーンズポテンシャルは、原子間の相互作用をモデル化するために使用されるポテンシャル関数です。

一般的な形式は次のようになります:

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

ここで、$ ( V(r) ) $はポテンシャルエネルギー、$ ( r ) $は原子間の距離、$ ( \epsilon ) $はポテンシャルの深さを表すパラメータ、$ ( \sigma ) $は原子間の最小距離(ポテンシャルがゼロとなる距離)を表すパラメータです。

Pythonでこのポテンシャルをグラフ化するために、適切なパラメータを選び、原子間距離の範囲でポテンシャルエネルギーを計算し、プロットします。

以下にサンプルコードを示します。

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

# レナード・ジョーンズポテンシャルの定義
def lennard_jones_potential(r, epsilon, sigma):
return 4 * epsilon * ((sigma / r)**12 - (sigma / r)**6)

# パラメータの設定
epsilon = 1.0 # ポテンシャルの深さ
sigma = 1.0 # 最小距離

# 距離の範囲
r_values = np.linspace(0.9 * sigma, 3 * sigma, 1000)

# ポテンシャルの計算
potential_values = lennard_jones_potential(r_values, epsilon, sigma)

# グラフ化
plt.plot(r_values, potential_values)
plt.xlabel('Distance (σ)')
plt.ylabel('Potential Energy (ε)')
plt.title('Lennard-Jones Potential')
plt.grid(True)
plt.show()

このコードは、レナード・ジョーンズポテンシャルを計算し、原子間距離の範囲でポテンシャルエネルギーをプロットします。
パラメータ$ ( \epsilon ) $と$ ( \sigma ) $を変更することでポテンシャルの形状を調整することができます。

[実行結果]

ソースコード解説

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

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

1
2
import numpy as np
import matplotlib.pyplot as plt
  • NumPy: 数値計算を行うためのライブラリ。
    数値的な操作や配列演算を効率的に行うために使用されます。
  • Matplotlib: グラフやプロットを作成するためのライブラリ。
    様々な種類の図やグラフを描画する機能を提供します。

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

1
2
def lennard_jones_potential(r, epsilon, sigma):
return 4 * epsilon * ((sigma / r)**12 - (sigma / r)**6)
  • lennard_jones_potential: レナード・ジョーンズポテンシャルを計算する関数。
    r原子間距離epsilonポテンシャルの深さを表すパラメータ、sigma最小距離を表すパラメータです。

3. パラメータの設定:

1
2
epsilon = 1.0  # ポテンシャルの深さ
sigma = 1.0 # 最小距離
  • epsilon: ポテンシャルの深さを示すパラメータ。
    ポテンシャルの最小値となるポテンシャルエネルギーを表します。
  • sigma: 原子間の最小距離を示すパラメータ。
    ポテンシャルがゼロとなる距離を表します。

4. 距離の範囲の設定:

1
r_values = np.linspace(0.9 * sigma, 3 * sigma, 1000)
  • np.linspace: 指定した範囲内で等間隔の数値を生成するNumPyの関数。
    ここでは、0.9 * sigmaから3 * sigmaの範囲で$1000$個の点を生成しています。

5. ポテンシャルの計算:

1
potential_values = lennard_jones_potential(r_values, epsilon, sigma)
  • lennard_jones_potential関数を使用して、距離の範囲に対するレナード・ジョーンズポテンシャルの値を計算します。

6. グラフの作成:

1
2
3
4
5
6
plt.plot(r_values, potential_values)
plt.xlabel('Distance (σ)')
plt.ylabel('Potential Energy (ε)')
plt.title('Lennard-Jones Potential')
plt.grid(True)
plt.show()
  • plt.plot: 指定された$x座標$と$y座標$のデータから折れ線グラフを作成します。
  • plt.xlabel: $x軸$のラベルを設定します。
  • plt.ylabel: $y軸$のラベルを設定します。
  • plt.title: グラフのタイトルを設定します。
  • plt.grid: グリッド線を表示するかどうかを設定します。
  • plt.show(): グラフを表示します。

これにより、レナード・ジョーンズポテンシャルの距離に対するポテンシャルエネルギーの関係を示すグラフが生成されます。

グラフ解説

[実行結果]

このグラフは、レナード・ジョーンズポテンシャルを可視化したものです。
ポテンシャルエネルギー $(( V ))$は、原子間の距離 $(( r ))$の関数として表されます。
以下に、グラフに表示される内容を詳しく説明します。

x軸(Distance (σ)):

原子間の距離 $(( r ))$を表します。
単位はレナード・ジョーンズポテンシャルのパラメータ$ ( \sigma ) $と同じです。
これは最小距離を$1$としたときの相対的な距離を表します。

y軸(Potential Energy (ε)):

ポテンシャルエネルギー$(( V ))$を表します。
単位はレナード・ジョーンズポテンシャルのパラメータ$ ( \epsilon ) $と同じです。
これは相互作用する原子間で生じるポテンシャルエネルギーの深さを示します。

レナード・ジョーンズポテンシャルの形状:

グラフは、原子間距離に対するポテンシャルエネルギーの関係を表しています。
この関数は、原子間距離が大きくなるにつれて急速に減少し、原子間距離が最小距離に近づくにつれてポテンシャルエネルギーが極小値を持ちます。

その後、原子間距離が大きくなるにつれてポテンシャルエネルギーがゆっくりと増加します。
これは、原子が非常に近いときと非常に遠いときには引力が働き、一定の距離で反発力が働くというレナード・ジョーンズポテンシャルの特徴を示しています。

このグラフを用いることで、レナード・ジョーンズポテンシャルによって定義される原子間の相互作用のエネルギーを直感的に理解することができます。

ヒル方程式

ヒル方程式

ヒル方程式は、生化学や医学などの分野で、酵素反応速度受容体のリガンド結合などの現象を記述するのに使われる非線形方程式です。

ヒル方程式は通常、次のように表されます:

$$
f(x) = \frac{x^n}{K^n + x^n}
$$

ここで、$ ( f(x) ) $は応答関数(応答の強さを示す)、$ ( x ) $は入力、$ ( n ) $はヒル係数(曲線の傾きを示す)、$ ( K ) $は半数飽和定数(入力が半数の最大効果をもたらす点)です。

以下は、Pythonでヒル方程式を解いてグラフ化するコードです。

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

# ヒル方程式
def hill_equation(x, n, K):
return x**n / (K**n + x**n)

# パラメータ設定
x_values = np.linspace(0, 10, 100) # xの範囲
n = 2 # ヒル係数
K = 5 # 半数飽和定数

# ヒル方程式の計算
y_values = hill_equation(x_values, n, K)

# グラフの描画
plt.plot(x_values, y_values)
plt.title('Hill Equation')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()

このコードでは、指定された範囲内の$ ( x ) $についてヒル方程式を計算し、結果をグラフ化します。

ヒル係数 $ ( n ) $と半数飽和定数 $ ( K ) $を変えることで、曲線の形状を調整することができます。

[実行結果]

ソースコード解説

以下にソースコードの詳細な説明を示します:

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

1
2
import numpy as np
import matplotlib.pyplot as plt
  • numpy:数値計算を行うためのPythonライブラリです。
  • matplotlib.pyplot:グラフの描画に使用されるMatplotlibのサブモジュールです。

2. ヒル方程式の定義:

1
2
def hill_equation(x, n, K):
return x**n / (K**n + x**n)
  • hill_equation 関数は、ヒル方程式を定義しています。
    この関数は、入力 $ (x) $、ヒル係数 $ (n)$、半数飽和定数 $ (K) $を引数として受け取り、ヒル方程式 $ (f(x) = \frac{x^n}{K^n + x^n}) $を計算します。

3. パラメータの設定:

1
2
3
x_values = np.linspace(0, 10, 100)  # xの範囲
n = 2 # ヒル係数
K = 5 # 半数飽和定数
  • x_values:x軸の値を設定します。
    ここでは、$0$から$10$までの範囲を等間隔で$100$個の点に区切った値を設定しています。
  • nヒル係数を設定します。
    この値は、ヒル方程式の曲線の傾きを制御します。
  • K半数飽和定数を設定します。
    この値は、入力が半数の最大効果をもたらす点を制御します。

4. ヒル方程式の計算:

1
y_values = hill_equation(x_values, n, K)
  • hill_equation 関数を使用して、ヒル方程式の値を計算します。
    入力$ (x) $の範囲に対応する$ (f(x)) $の値を計算し、y_values に格納します。

5. グラフの描画:

1
2
3
4
5
6
plt.plot(x_values, y_values)
plt.title('Hill Equation')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(True)
plt.show()
  • plt.plot(x_values, y_values)ヒル方程式の曲線をプロットします。
  • plt.title('Hill Equation'):グラフのタイトルを設定します。
  • plt.xlabel('x'):$x軸$のラベルを設定します。
  • plt.ylabel('f(x)'):$y軸$のラベルを設定します。
  • plt.grid(True):グリッドを表示します。
  • plt.show():グラフを表示します。

このプログラムは、指定されたパラメータでヒル方程式を計算し、その結果をグラフ化しています。

結果解説

[実行結果]

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

x軸:

入力$ (x) $の値を示します。
この範囲は、$0$から$10$までの$100$個の等間隔の点で設定されています。
この範囲内の$ (x) $の値に対して、ヒル方程式が計算されます。

y軸:

応答関数$ (f(x)) $の値を示します。
これは、ヒル方程式によって計算された応答の強さを表します。

タイトル:

グラフのタイトルは「Hill Equation」となっています。
これは、このグラフがヒル方程式の曲線を表していることを示しています。

x軸ラベル:

$x軸$のラベルは「x」となっています。
これは、グラフの$x軸$が入力$ (x) $を表していることを示しています。

y軸ラベル:

$y軸$のラベルは「f(x)」となっています。
これは、グラフの$y軸$がヒル方程式の応答関数 $ (f(x)) $を表していることを示しています。

グリッド:

グラフにはグリッドが表示されています。
これは、目盛りの補助として、グラフのデータポイントの位置を視覚的に把握するのに役立ちます。

このグラフは、ヒル方程式の形状や曲線の特性を視覚的に理解するのに役立ちます。

エルミート多項式

エルミート多項式

エルミート多項式は、次の微分方程式を満たす多項式です:

$$
\frac{d^2 H_n(x)}{dx^2} - 2x \frac{dH_n(x)}{dx} + 2n H_n(x) = 0
$$

この微分方程式をPythonで解き、グラフ化するために、SciPyodeint関数を使用します。

以下は、このタスクを実行する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
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# エルミート多項式の微分方程式
def hermite_eq(y, x, n):
dydx = [y[1], 2*x*y[1] - 2*n*y[0]]
return dydx

# 初期条件と解く範囲の設定
x = np.linspace(-5, 5, 1000) # xの範囲
y0 = [1, 0] # 初期条件 H_0(x) = 1, H_0'(x) = 0

# エルミート多項式の計算
n_values = [0, 1, 2, 3] # nの値
for n in n_values:
y = odeint(hermite_eq, y0, x, args=(n,))
plt.plot(x, y[:, 0], label=f'n={n}')

plt.title("Hermite Polynomials")
plt.xlabel("x")
plt.ylabel("H_n(x)")
plt.legend()
plt.grid(True)
plt.show()

このコードは、エルミート多項式を計算し、$n=0,1,2,3$に対する結果をグラフ化します。

グラフの$x軸$は$-5$から$5$の範囲で、$y軸$はエルミート多項式 $(H_n(x))$の値です。

[実行結果]

ソースコード解説

以下にコードの詳細な説明を示します。

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

1
2
3
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
  • numpyは数値計算用のライブラリであり、配列操作や数値計算機能を提供します。
  • scipy.integrate.odeintは常微分方程式の数値積分を行う関数です。
  • matplotlib.pyplotはグラフの描画を行うためのライブラリです。

2. エルミート多項式の微分方程式の定義:

1
2
3
def hermite_eq(y, x, n):
dydx = [y[1], 2*x*y[1] - 2*n*y[0]]
return dydx
  • hermite_eq関数はエルミート多項式微分方程式を定義しています。
    この微分方程式は2階の常微分方程式です。

3. 初期条件と解く範囲の設定:

1
2
x = np.linspace(-5, 5, 1000)  # xの範囲
y0 = [1, 0] # 初期条件 H_0(x) = 1, H_0'(x) = 0
  • xはグラフの描画範囲であり、$-5$から$5$の範囲を$1000$等分して設定されています。
  • y0は微分方程式の初期条件です。
    ここでは、$ (H_0(x) = 1) $と$ (H_0’(x) = 0) $を設定しています。

4. エルミート多項式の計算とグラフ化:

1
2
3
4
n_values = [0, 1, 2, 3]  # nの値
for n in n_values:
y = odeint(hermite_eq, y0, x, args=(n,))
plt.plot(x, y[:, 0], label=f'n={n}')
  • n_valuesエルミート多項式の次数$ (n) $の値のリストです。
    ここでは、$0$から$3$までの整数が設定されています。
  • odeint関数を用いて微分方程式を解き、エルミート多項式の値を計算しています。
    結果は変数yに格納されます。
  • plt.plot()関数を用いて、計算されたエルミート多項式のグラフを描画します。
    ラベルにはnの値が表示されます。

5. グラフの装飾:

1
2
3
4
5
6
plt.title("Hermite Polynomials")
plt.xlabel("x")
plt.ylabel("H_n(x)")
plt.legend()
plt.grid(True)
plt.show()
  • グラフのタイトル、$x軸$ラベル、$y軸$ラベルを設定します。
  • plt.legend()で凡例を表示し、plt.grid(True)でグリッドを表示します。
  • 最後にplt.show()でグラフを表示します。

これにより、エルミート多項式微分方程式を解いて計算し、それらをグラフ化する完全なプログラムが得られます。

結果解説

[実行結果]

このプログラムは、エルミート多項式$ (H_n(x)) $を計算し、その結果をグラフ化しています。
以下はグラフに表示される内容の詳細です:

1. x軸の範囲:

$-5$から$5$の間で$1000$個の等間隔の点が生成されます。
これはグラフの横軸に対応します。

2. y軸の値:

エルミート多項式 $ (H_n(x)) $の値が計算され、それぞれの$n$に対して異なる線がプロットされます。
各線は異なる色で表されます。

3. nの値:

$n$はエルミート多項式の次数を表します。
プログラムでは、$n$が$0$から$3$までの$4$つの値を取ります。
これにより、それぞれのエルミート多項式が表示されます。

4. グラフのタイトル:

グラフのタイトルには「Hermite Polynomials」と表示されます。

5. x軸とy軸のラベル:

x軸とy軸にはそれぞれ「x」と「H_n(x)」というラベルが付けられます。

6. 凡例:

各線の凡例には、対応する$n$の値が表示されます。
これにより、各線がどのエルミート多項式に対応しているかがわかります。

7. グリッド線:

グラフにはグリッド線が表示され、視覚的な参照を提供します。