有名な数学的方程式の１つとして、シュレーディンガー方程式があります。

これは量子力学における基本的な方程式であり、粒子の波動関数の時間変化を記述します。

一次元の自由粒子のシュレーディンガー方程式は以下のように表されます：

$$
i\hbar \frac{\partial}{\partial t} \Psi(x, t) = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} \Psi(x, t)
$$

ここで、$ ( \Psi(x, t) ) $は波動関数、$ ( x ) $は空間座標、$ ( t ) $は時間、$ ( \hbar ) $はディラック定数、$ ( m ) $は粒子の質量です。

この方程式を解いて、時間と空間に依存する波動関数をプロットします。

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

# 定数の設定
h_bar = 1.0  # ディラック定数
m = 1.0      # 粒子の質量

# 空間および時間の設定
x = np.linspace(-5, 5, 400)
t = np.linspace(0, 2*np.pi, 400)

# 波動関数の計算
X, T = np.meshgrid(x, t)
Psi = np.exp(-(X**2)/2) * np.exp(1j * T)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.imshow(np.abs(Psi)**2, extent=(-5, 5, 0, 2*np.pi), aspect='auto', cmap='hot')
plt.colorbar(label='Probability Density')
plt.title('Probability Density of a Free Particle')
plt.xlabel('Position')
plt.ylabel('Time')
plt.show()

このコードは、一次元の自由粒子のシュレーディンガー方程式を解き、時間と空間に依存する波動関数の確率密度をプロットします。

確率密度は波動関数の絶対値の二乗として与えられ、波動関数が存在する確率を表します。

[実行結果]

ソースコード解説

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

1. import numpy as np:

• NumPy ライブラリを np としてインポートします。
  NumPy は、数値計算を行うための強力なライブラリであり、配列操作や数学関数の計算に便利です。

2. import matplotlib.pyplot as plt:

• Matplotlib ライブラリの pyplot モジュールを plt としてインポートします。
  Matplotlib は、グラフを描画するためのライブラリであり、pyplot モジュールはその中でも一般的に使用されます。

3. h_bar = 1.0 および m = 1.0:

• シュレーディンガー方程式に登場する定数$ ( \hbar )$（ディラック定数）と粒子の質量$ ( m ) $の値を設定します。
  ここでは両方とも値を$ 1.0 $に設定しています。

4. x = np.linspace(-5, 5, 400) および t = np.linspace(0, 2*np.pi, 400):

• 空間座標$ ( x ) $および時間$ ( t ) $の範囲をそれぞれ -5 から 5 までの区間と 0 から 2π までの区間に分割し、それぞれ$ 400 $個の等間隔な点で表現します。
  これにより、グラフ上で表示される座標の範囲が決まります。

5. X, T = np.meshgrid(x, t):

• NumPy の meshgrid() 関数を使って、空間座標 x と時間 t を格子状の座標系に変換します。
  これにより、各点での波動関数の値を計算するための格子が作成されます。

6. Psi = np.exp(-(X**2)/2) * np.exp(1j * T):

• シュレーディンガー方程式の解である波動関数$ ( \Psi(x, t) ) $を計算します。
  ここで、np.exp() は指数関数を計算し、1j は虚数単位を表します。
  この波動関数は、時間と空間に依存する複素数で表現されます。

7. plt.figure(figsize=(10, 6)):

• Matplotlib で図を作成します。
  図のサイズを幅$10$インチ、高さ$6$インチに設定しています。

8. plt.imshow(np.abs(Psi)**2, extent=(-5, 5, 0, 2*np.pi), aspect='auto', cmap='hot'):

• imshow() 関数を使って、波動関数の絶対値の二乗である確率密度をイメージとして表示します。
  extent 引数は$ x $軸と$ y $軸の範囲を指定し、aspect 引数はアスペクト比を自動調整するように指定します。
  cmap 引数はカラーマップを指定します。

9. plt.colorbar(label='Probability Density'):

• カラーバーをグラフに追加し、確率密度の色に対応する値を表示します。

10. plt.title('Probability Density of a Free Particle'):

• グラフのタイトルを設定します。

11. plt.xlabel('Position') および plt.ylabel('Time'):

• $ x $軸と$ y $軸のラベルを設定します。

12. plt.show():

• グラフを表示します。

結果解説

下記のグラフは、一次元自由粒子のシュレーディンガー方程式の解に基づいて、時間と空間に依存する波動関数の確率密度を表しています。

[実行結果]

横軸は空間座標$ ( x ) $を表し、縦軸は時間$ ( t ) $を表しています。

グラフの各点における色の濃さは、その空間座標と時間における波動関数の確率密度を示しています。

色が濃い領域ほど確率密度が高く、粒子が存在する可能性が高いことを意味します。

このグラフは、波動関数の時間変化を視覚的に理解するのに役立ちます。

時間が経過するにつれて、波動関数の形状が変化し、特定の空間座標で粒子の存在確率が変動します。
特に、波動関数の振動が時間によってどのように進行し、確率密度がどのように変化するかが示されています。

宇宙の膨張 フリードマン方程式

宇宙の膨張を表すフリードマン方程式を考えます。

フリードマン方程式は、宇宙のスケールファクター $ ( a(t) ) $とエネルギー密度 $ ( \rho(t) ) $に関連します。

方程式は以下の通りです：

$$
H^2 = \left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G}{3}\rho - \frac{k}{a^2}
$$

ここで、$ ( H ) $はハッブル定数、$ ( G ) $は重力定数、$ ( k ) $は空間の曲率を表します。

また、ドットは時間に関する微分を表しています。

この方程式を解くには、初期条件やエネルギー密度などが必要です。

一般的には、標準的な宇宙論の初期条件を使用します。

以下は、Pythonを使用してこの方程式を解き、グラフ化する簡単な例です。

この例では、ハッブル定数 $ ( H_0 ) $を$1$、重力定数 $ ( G ) $を$1$、曲率 $ ( k ) $を$0$と仮定します。

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
from scipy.integrate import odeint

# フリードマン方程式
def friedmann_equation(y, t):
    H0 = 1.0  # ハッブル定数
    G = 1.0   # 重力定数
    k = 0.0   # 曲率
    rho = y[0]  # エネルギー密度

    # フリードマン方程式
    dydt = [y[1], (8 * np.pi * G / 3) * rho - k / y[0]**2]
    return dydt

# 初期条件
initial_condition = [1.0, 0.0]  # 初期エネルギー密度、初期速度

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

# 方程式を解く
solution = odeint(friedmann_equation, initial_condition, t)

# 結果をプロット
plt.plot(t, solution[:, 0], label='Scale Factor (a(t))')
plt.xlabel('Time')
plt.ylabel('Scale Factor')
plt.title('Friedmann Equation Solution')
plt.legend()
plt.show()

この例では、初期エネルギー密度を$1$、初期速度を$0$としています。

得られた結果は、時間に対する宇宙のスケールファクターの変化を示しています。

これは、宇宙が膨張していく様子を表しています。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用して宇宙の膨張を表すフリードマン方程式を解き、その解をプロットするためのものです。

以下はソースコードの章立てと詳細な説明です：

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

1
2
3

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

• numpy：数値計算用のライブラリ。
  数学的な演算や配列操作に使用されます。
• matplotlib.pyplot：データの可視化のためのライブラリ。
  グラフの作成や表示に利用されます。
• scipy.integrate.odeint：常微分方程式を解くための関数が含まれています。

2. フリードマン方程式の定義

1
2
3
4
5
6
7
8
9

def friedmann_equation(y, t):
    H0 = 1.0  # ハッブル定数
    G = 1.0   # 重力定数
    k = 0.0   # 曲率
    rho = y[0]  # エネルギー密度

    # フリードマン方程式
    dydt = [y[1], (8 * np.pi * G / 3) * rho - k / y[0]**2]
    return dydt

• friedmann_equation 関数：

フリードマン方程式を定義します。
この関数は、時間 t と状態ベクトル y を受け取り、その微分を返します。

• H0：ハッブル定数。
  宇宙の膨張を示す定数。
• G：重力定数。
  重力の強さを示す定数。
• k：曲率。
  宇宙の空間が平坦か、閉じているか、開いているかを示すパラメータ。
• rho：エネルギー密度。
  宇宙のエネルギー分布を表す。
• dydt：微分方程式の右辺。
  宇宙の膨張を示すフリードマン方程式の微分方程式。

3. 初期条件の設定

1

initial_condition = [1.0, 0.0]  # 初期エネルギー密度、初期速度

• initial_condition：フリードマン方程式の初期条件。
  初期エネルギー密度と初期速度が設定されています。

4. 時間の範囲の設定

1

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

• t：時間の範囲。$0$から$10$までの範囲を等間隔に区切った$1000$のデータポイントを生成します。

5. 方程式の解の計算

1

solution = odeint(friedmann_equation, initial_condition, t)

• odeint 関数：常微分方程式を解くための関数。
  friedmann_equation 関数、初期条件 initial_condition、時間 t を引数として渡し、方程式の解を計算します。

6. 結果のプロット

1
2
3
4
5
6

plt.plot(t, solution[:, 0], label='Scale Factor (a(t))')
plt.xlabel('Time')
plt.ylabel('Scale Factor')
plt.title('Friedmann Equation Solution')
plt.legend()
plt.show()

• plt.plot：解をプロットします。
  横軸に時間 t、縦軸にスケールファクター a(t) をプロットします。
• plt.xlabel：X軸のラベルを設定。
• plt.ylabel：Y軸のラベルを設定。
• plt.title：グラフのタイトルを設定。
• plt.legend：凡例を表示。
• plt.show：グラフを表示します。

このコードは、フリードマン方程式を解いて宇宙の膨張を表すグラフを生成します。

初期条件やパラメータの変更によって、グラフの形状が変わります。

結果解説

[実行結果]

上記のPythonコードで生成されるグラフは、フリードマン方程式を解いた結果、時間に対する宇宙のスケールファクター$ (a(t)) $の変化を表しています。

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

1. 横軸（X軸） - 時間：

グラフの横軸は時間を表しています。
この例では、時間の範囲を$ (0) $から$ (10) $までとし、$ (1000) $の等間隔なデータポイントで構成されています。

2. 縦軸（Y軸） - スケールファクター $ (a(t)) $：

グラフの縦軸は宇宙のスケールファクター $ (a(t)) $を表しています。
スケールファクターは宇宙の膨張を示し、大きくなるほど宇宙が膨張していることを示します。

3. グラフの形状：

グラフは時間に対する宇宙の膨張を示しており、スケールファクターが増加していく様子を観察できます。
これは、ビッグバン理論に基づく宇宙の進化を示しています。
初期条件によってグラフの形状が変わります。

4. 初期条件：

この例では初期エネルギー密度を$ (1) $、初期速度を$ (0) $と仮定しています。
これらの初期条件に基づいて、時間が経過するにつれて宇宙がどのように膨張しているかが反映されています。

このグラフは、ビッグバンモデルにおける宇宙の膨張をシミュレートしたものであり、初期条件やパラメータの変更によって異なる結果が得られます。

常微分方程式

微分方程式の例として、単純な1階の常微分方程式を取り上げましょう。

以下の例では、微分方程式 $ ( \frac{dy}{dt} = -ky ) $を考えます。

この方程式は指数減衰を表します。

これを解いて、その結果をグラフ化します。

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

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

# 微分方程式の定義
def model(y, t, k):
    dydt = -k * y
    return dydt

# 初期条件
y0 = 5.0

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

# パラメータ
k = 0.1

# 微分方程式を解く
solution = odeint(model, y0, t, args=(k,))

# 結果のグラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

このコードでは、SciPyのodeint関数を使用して微分方程式を解き、結果をグラフ化しています。

[実行結果]

指数減衰が確認できるでしょう。

パラメータ$ (k) $を変更することで解の挙動が変わります。

ソースコード解説

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

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

1
2
3

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

• numpy: 数値計算用ライブラリで、数値配列や数学関数などを提供します。
• matplotlib.pyplot: グラフの描画に使用されるライブラリ。
• scipy.integrate.odeint: 常微分方程式を数値的に解くための関数が含まれているSciPyライブラリの一部。

2. 微分方程式の定義:

1
2
3

def model(y, t, k):
   dydt = -k * y
   return dydt

• model関数は、微分方程式$ ( \frac{dy}{dt} = -ky ) $の右辺を定義しています。
  この関数は odeint 関数によって呼び出されます。

3. 初期条件の設定:

1

y0 = 5.0

• 微分方程式の初期条件$ (y(0)) $を設定します。

4. 時間の範囲の設定:

1

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

• 時間の範囲を$0$から$10$までの区間を$100$点で均等に区切って定義します。

5. パラメータの設定:

1

k = 0.1

• 微分方程式のパラメータ$ (k) $を設定します。

6. 微分方程式を解く:

1

solution = odeint(model, y0, t, args=(k,))

• odeint 関数を使用して微分方程式を解きます。
  初期条件 y0、時間の範囲 t、およびパラメータ k を指定しています。

7. 結果のグラフ化:

1
2
3
4
5
6
7
8

plt.figure(figsize=(8, 6))
plt.plot(t, solution, label='y(t)')
plt.title('指数減衰：$\\frac{dy}{dt} = -ky$')
plt.xlabel('時間 (t)')
plt.ylabel('y(t)')
plt.legend()
plt.grid(True)
plt.show()

• matplotlib を使用して、微分方程式の解 solution を時間に対してプロットし、グラフを表示します。
• グラフのタイトル、軸ラベル、凡例、グリッドが設定されています。

このプログラムは、指数減衰を示す微分方程式を数値的に解き、その解をグラフで視覚化する基本的な例です。

結果解説

[実行結果]

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

横軸 (X軸):時間$ (t) $

$0$から$10$の範囲を$100$点で等分しています。

縦軸 (Y軸):変数$ (y(t))$

微分方程式の解であり、指数減衰の挙動を示します。
初期条件$ (y(0)) $は$ 5.0 $としています。

グラフの線:

微分方程式を解いた結果の$ (y(t)) $の変化を表しています。
指数減衰の特徴的な形状で、初めは急速に減少し、後に緩やかな減少に移ります。

グラフのタイトル:指数減衰：$ \frac{dy}{dt} = -ky $

微分方程式の性質と解いたものが指数減衰を示すことを示しています。

X軸およびY軸のラベル:

時間（時間$ (t)$）および変数$ (y(t)) $の説明があります。

凡例 (Legend):$y(t)$

グラフ内の線が何を示しているかを示すラベルです。

グリッド:

グラフ内に格子状の線が描かれ、視認性を向上させています。

このグラフは、微分方程式の解を通して、指数減衰が時間とともにどのように進行するかを示しています。

指数減衰は、物理学や生態学などのさまざまな分野で重要な概念であり、この例ではPythonを用いて簡単に表現されています。

最適化問題

最適化問題の例として、特定の制約条件下での関数の最小化を考えてみましょう。

以下は、2変数の関数の最小化問題です。

この例では、scipy.optimize ライブラリを使用して最適化を行います。

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

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# 最小化する関数
def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

# 制約条件
constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

# 初期値
initial_guess = [0, 0]

# 最適化の実行
result = minimize(objective_function, initial_guess, constraints=constraints)

# 結果表示
print("最適解:", result.x)
print("最小値:", result.fun)

# グラフの作成
x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

# 制約条件を描画
plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

# 最小値を描画
plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

# 関数の等高線を描画
contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

# グラフの装飾
plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

# グラフの表示
plt.show()

このコードでは、関数$ x^2 + y^2 $を最小化する問題を考え、制約条件として$ x + y = 5 $を設定しています。

最適化結果や最小化される関数のグラフが表示されます。

[実行結果]

異なる関数や制約条件を試すことで、様々な最適化問題に対応できます。

ソースコード解説

以下にソースコードの説明を示します。

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

1
2
3

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

• numpy: 数値計算のためのライブラリ。
• scipy.optimize.minimize: 最適化問題を解くための関数。
• matplotlib.pyplot: グラフ描画のためのライブラリ。

2. 最小化する関数の定義:

1
2

def objective_function(x):
    return x[0]**2 + x[1]**2  # 例として x^2 + y^2 を最小化

• x は2要素の配列で、x[0] と x[1] が変数となります。
• $ x[0]^2 + x[1]^2 $を最小化する関数を定義しています。

3. 制約条件の設定:

1

constraints = ({'type': 'eq', 'fun': lambda x: x[0] + x[1] - 5})  # 制約条件: x + y = 5

• 制約条件として、$ x + y = 5 $を設定しています。

4. 初期値の設定:

1

initial_guess = [0, 0]

• 最適化の初期値を [0, 0] に設定します。

5. 最適化の実行:

1

result = minimize(objective_function, initial_guess, constraints=constraints)

• minimize 関数を使用して最適化を実行します。
• objective_function を最小化し、制約条件 constraints の下で最適解を見つけます。

6. 結果の表示:

1
2

print("最適解:", result.x)
print("最小値:", result.fun)

• 最適解と最小値を表示します。

7. グラフの作成:

1
2
3
4

x_vals = np.linspace(-2, 7, 100)
y_vals = np.linspace(-2, 7, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = X**2 + Y**2

• グラフ描画用の x, y 座標を生成します。
• Z には$ X^2 + Y^2 $の関数値が格納されます。

8. 制約条件の描画:

1

plt.plot(x_vals, 5 - x_vals, label='制約条件: x + y = 5', color='red')

• 制約条件を赤い線で描画します。

9. 最小値の描画:

1

plt.scatter(result.x[0], result.x[1], color='green', marker='x', label='最小値')

• 最小値の点を緑の “x” マークで描画します。

10. 等高線プロットの描画:

1
2

contour = plt.contour(X, Y, Z, levels=15, cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)

• 関数の等高線プロットを描画します。

11. グラフの装飾:

1
2
3
4
5

plt.title('2変数関数の最小化と制約条件')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)

• グラフのタイトルや軸ラベル、凡例、グリッドを設定します。

12. グラフの表示:

1

plt.show()

• これまでの設定で作成したグラフを表示します。

このコードは、2変数の関数を最小化し、その結果をグラフに可視化しています。

特に、制約条件下での最小値を探索している点がポイントです。

結果解説

[実行結果]

この最適化問題は、関数$ ( f(x, y) = x^2 + y^2 ) $を最小化するというもので、制約条件として$ ( x + y = 5 ) $が課されています。

最適化アルゴリズムは、この関数を最小化する$ ( x ) と ( y ) $の値を見つけることを試みます。

結果は以下の通りです：

• 最適解: [2.5, 2.5]
• 最小値: 12.5

最適解 $ [2.5, 2.5] $は、制約条件 $ ( x + y = 5 ) $を満たす中で$ ( f(x, y) ) $を最小化する点です。

最小値 $ 12.5 $は、この最小化された関数の値です。

グラフには、等高線プロットが表示されています。

等高線は、関数の高さを示し、最小値がグリーンの “x” マークで示されています。

また、制約条件 $ ( x + y = 5 ) $を赤い線で示しています。

最適解がこの赤い線上にあり、かつ関数を最小化する点であることが確認できます。

この例は非常にシンプルですが、最適化問題の基本的な理解を提供しています。

実際の応用では、より複雑な関数や複数の制約条件がありますが、同様のアプローチが使われます。

2変数の定積分

例題として、2変数の定積分を考えてみましょう。

以下の関数を考えます：

$$
f(x, y) = x^2 + y^2
$$

そして、この関数を$ [x \in [-2, 2], y \in [-2, 2]] $の範囲で積分してみます。

Pythonで積分を計算し、結果をグラフ化するためには、SciPyと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

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

# 積分する関数
def integrand(x, y):
    return x**2 + y**2

# 積分範囲
x_range = (-2, 2)
y_range = (-2, 2)

# 積分計算
result, _ = integrate.dblquad(integrand, x_range[0], x_range[1], lambda x: y_range[0], lambda x: y_range[1])

print("積分の結果:", result)

# グラフ化
x_vals = np.linspace(x_range[0], x_range[1], 100)
y_vals = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(x_vals, y_vals)
Z = integrand(X, Y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', rstride=5, cstride=5, alpha=0.7, linewidth=0.5, edgecolors='w')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('2D Integral Example')

plt.show()

このコードでは、dblquadを使用して2変数の積分を計算しています。

そして、plot_surfaceメソッドを使用して関数の3Dプロットを作成しています。

プログラムを実行すると、積分の結果が表示され、関数がどのように分布しているかがグラフで視覚化されます。

[実行結果]

ソースコード解説

このソースコードは、Pythonを使用して2変数の関数を定積分し、その結果を3Dグラフ化するサンプルコードです。

以下に、コードの各部分の説明を示します。

1
2
3
4

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

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

numpyは数値計算を効率的に行うためのライブラリ、scipyは科学技術計算のためのライブラリ、matplotlibはグラフの描画のためのライブラリ、Axes3Dは3Dグラフを描画するためのサブモジュールです。

1
2
3

# 積分する関数
def integrand(x, y):
    return x**2 + y**2

integrandという関数を定義しています。

この関数は、2変数$ [x, y] $に対して$ [x^2 + y^2] $を返します。

これが積分対象の関数です。

1
2
3

# 積分範囲
x_range = (-2, 2)
y_range = (-2, 2)

x_rangeおよびy_rangeで積分の範囲を指定しています。

この例では、$ [x \in [-2, 2]] $および$ [y \in [-2, 2]] $の範囲で積分を行います。

1
2

# 積分計算
result, _ = integrate.dblquad(integrand, x_range[0], x_range[1], lambda x: y_range[0], lambda x: y_range[1])

integrate.dblquad関数を使用して2重積分を計算しています。
dblquadは2変数関数の定積分を計算する関数で、第一引数に積分対象の関数integrand、次に$ [x] $の範囲、そしてその後に$ [y] $の範囲を指定します。
積分結果はresultに格納されます。

1

print("積分の結果:", result)

計算結果をコンソールに出力しています。

1
2
3
4
5
6

# グラフ化
x_vals = np.linspace(x_range[0], x_range[1], 100)
y_vals = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(x_vals, y_vals)
Z = integrand(X, Y)

グラフを描画するために、$ [x] $と$ [y] $の値を生成しています。
linspace関数を使用して、指定された範囲内に等間隔の点を生成しています。
meshgrid関数を使用して、これらの点から座標行列$ [X, Y] $を作成し、それを用いて関数の値$ [Z] $を計算しています。

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, Z, cmap='viridis', rstride=5, cstride=5, alpha=0.7, linewidth=0.5, edgecolors='w')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('2D Integral Example')

plt.show()

最後に、matplotlibを使用して3Dグラフを描画しています。
plot_surfaceメソッドを使用して、$ [X, Y, Z] $の3次元座標で関数をプロットし、その他のオプションでグラフのスタイルを調整しています。

最終的にplt.show()でグラフを表示します。

結果解説

[実行結果]

積分の結果が$ 42.666666666666664 $と表示されました。

これは、関数$ [ f(x, y) = x^2 + y^2 ] $を範囲$ [x \in [-2, 2], y \in [-2, 2]] $で積分した値です。
この積分は、与えられた領域で関数の値を足し合わせることに相当します。

この場合、積分結果が数値で表示されるため、具体的な意味が直感的には理解しにくいかもしれません。
しかし、この数値は範囲内で関数が取る面積や体積に相当します。

積分が正の値であることから、関数$ [ f(x, y) = x^2 + y^2 ] $が指定された範囲で正の値を取り、範囲内の領域の面積が正であることを示しています。

グラフは、積分対象の関数$ [ f(x, y) = x^2 + y^2 ] $の3Dプロットを示しています。
この関数は放物面であり、中心が原点で、$ [ x ] $と$ [ y ] $の値が大きくなるほど関数の値が増加します。
したがって、グラフは中心から外側に向かって放物線が立体的に広がる形状をしています。

この例での積分とグラフは、数学的な計算と視覚的な理解を結びつけるものであり、積分が関数の形状や領域に対する面積や体積の計算にどのように寄与するかを示しています。

人口動態

人口動態の例題として、シンプルな人口成長モデルを考えてみましょう。

ここでは、離散的な時間ステップでの人口の変化を表すシンプルな差分方程式を使用します。

例として、次のような差分方程式を考えます：

$$
[ P_{t+1} = P_t + r \cdot P_t \cdot (1 - \frac{P_t}{K}) ]
$$

ここで、$ (P_t) $は時刻$ (t) $での人口、$ (r) $は増加率、$ (K) $は環境収容力（持続可能な最大人口）です。

以下に、この人口成長モデルを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

import matplotlib.pyplot as plt

def population_growth_model(initial_population, growth_rate, carrying_capacity, generations):
    population = [initial_population]

    for t in range(1, generations):
        next_population = population[-1] + growth_rate * population[-1] * (1 - population[-1] / carrying_capacity)
        population.append(next_population)

    return population

# パラメータ設定
initial_population = 100  # 初期人口
growth_rate = 0.1        # 増加率
carrying_capacity = 500   # 環境収容力
generations = 50          # 世代数

# 人口成長モデルの実行
population_over_time = population_growth_model(initial_population, growth_rate, carrying_capacity, generations)

# 結果のグラフ化
plt.plot(range(generations), population_over_time, marker='o', linestyle='-', color='b')
plt.title('人口成長モデル')
plt.xlabel('世代')
plt.ylabel('人口')
plt.grid(True)
plt.show()

このコードでは、初期人口が$100$人で、増加率が$0.1$、環境収容力が$500$人の場合の人口成長モデルをシミュレートしています。

結果を折れ線グラフで表示しています。

[実行結果]

グラフを通じて、初期急激な成長があり、次第に環境収容力に向かって人口が収束していく様子がわかります。

ソースコード解説

このソースコードは、人口成長モデルを定義し、それを用いて人口の時間変化をシミュレートし、最終的には結果をMatplotlibを使用してグラフ化するものです。

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

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

1

import matplotlib.pyplot as plt

Matplotlibのpyplotモジュールをインポートしています。

これはグラフの描画に使用されます。

2. 人口成長モデルの定義

1
2
3
4
5
6
7
8

def population_growth_model(initial_population, growth_rate, carrying_capacity, generations):
    population = [initial_population]

    for t in range(1, generations):
        next_population = population[-1] + growth_rate * population[-1] * (1 - population[-1] / carrying_capacity)
        population.append(next_population)

    return population

population_growth_model関数は、差分方程式に基づくシンプルな人口成長モデルを実装しています。
初期人口から始まり、各世代ごとに人口を計算し、結果をリストとして返します。

3. パラメータ設定

1
2
3
4

initial_population = 100  # 初期人口
growth_rate = 0.1        # 増加率
carrying_capacity = 500   # 環境収容力
generations = 50          # 世代数

モデルのパラメータを設定しています。
初期人口、増加率、環境収容力、およびシミュレーションの世代数が含まれます。

4. 人口成長モデルの実行

1

population_over_time = population_growth_model(initial_population, growth_rate, carrying_capacity, generations)

先に定義した人口成長モデルを実行し、時間経過に伴う人口の変化をリストとして取得します。

5. 結果のグラフ化

1
2
3
4
5
6

plt.plot(range(generations), population_over_time, marker='o', linestyle='-', color='b')
plt.title('人口成長モデル')
plt.xlabel('世代')
plt.ylabel('人口')
plt.grid(True)
plt.show()

Matplotlibを使用して人口の時間変化をグラフに描画します。
横軸は世代数、縦軸は人口数です。
グリッドが表示され、グラフのタイトルと軸ラベルが設定されています。
グラフは点線（'o'）で表示され、色は青（'b'）です。

このコード全体では、シンプルな人口成長モデルを定義し、それを用いて特定のパラメータでの人口変化を可視化しています。

簡単な数理解析

方程式の解を求める簡単な例題を取り上げます。

以下の2次方程式を解く問題を考えましょう。

$$
[ax^2 + bx + c = 0]
$$

ここで、$ (a) $,$ (b) $,$ (c) $は任意の実数です。

解の公式を用いて、この方程式の解を求める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 cmath  # 複素数の計算が必要な場合に使います

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

    # 判別式が非負なら実数解が存在
    if discriminant >= 0:
        root1 = (-b + discriminant**0.5) / (2*a)
        root2 = (-b - discriminant**0.5) / (2*a)
        return root1, root2
    else:
        # 判別式が負の場合、複素数解が存在
        complex_root1 = (-b + cmath.sqrt(discriminant)) / (2*a)
        complex_root2 = (-b - cmath.sqrt(discriminant)) / (2*a)
        return complex_root1, complex_root2

# 例として、a=1, b=-3, c=2の場合を考える
a = 1
b = -3
c = 2

roots = solve_quadratic_equation(a, b, c)
print("方程式の解:", roots)

このコードでは、解の公式を用いて方程式の解を計算しています。

判別式が非負なら実数解が存在し、判別式が負の場合には複素数解が存在します。

上記の例では、$ (a=1) $,$ (b=-3) $,$ (c=2) $の2次方程式の解を求めています。

[実行結果]

方程式の解: (2.0, 1.0)

ソースコード解説

このソースコードは、2次方程式$ (ax^2 + bx + c = 0) $の解を求めるための関数を定義し、その関数を具体的な例で呼び出しています。

以下に、コードの章立てと詳細な説明を示します。

1. 複素数計算のためのライブラリの導入:

1

import cmath  # 複素数の計算が必要な場合に使います

コードの冒頭で cmath ライブラリを導入しています。
これは、複素数を扱うためのライブラリで、実数解だけでなく複素数解も考慮できるようになります。

2. 2次方程式の解を求める関数の定義:

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

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

   # 判別式が非負なら実数解が存在
   if discriminant >= 0:
       root1 = (-b + discriminant**0.5) / (2*a)
       root2 = (-b - discriminant**0.5) / (2*a)
       return root1, root2
   else:
       # 判別式が負の場合、複素数解が存在
       complex_root1 = (-b + cmath.sqrt(discriminant)) / (2*a)
       complex_root2 = (-b - cmath.sqrt(discriminant)) / (2*a)
       return complex_root1, complex_root2

この関数は、2次方程式の係数 (a), (b), (c) を引数として受け取り、解を返すものです。
判別式を計算し、その値に基づいて解の計算方法を分岐させています。

• 判別式が非負の場合（実数解が存在する場合）：
  実数解を計算して root1 および root2 に代入します。
• 判別式が負の場合（複素数解が存在する場合）：
  複素数解を計算して complex_root1 および complex_root2 に代入します。

最終的に、計算された解が返されます。

3. 例として方程式を解く:

1
2
3
4
5
6
7

# 例として、a=1, b=-3, c=2の場合を考える
a = 1
b = -3
c = 2

roots = solve_quadratic_equation(a, b, c)
print("方程式の解:", roots)

具体的な係数$ (a) $,$ (b) $,$ (c) $を与えて、solve_quadratic_equation 関数を呼び出し、得られた解を出力します。

これにより、与えられた2次方程式の解が実数か複素数かに応じて正しく計算され、結果が出力されます。

仮想通貨価格のモデリング (Geometric Brownian Motion)

一般的な仮想通貨価格のモデリングにおいて、Geometric Brownian Motion (GBM) がよく用いられます。

以下はその例です。

$$
dS = \mu S dt + \sigma S dW
$$

• $ (S) $は仮想通貨の価格
• $ (\mu) $はドリフト（平均収益率）
• $ (\sigma) $はボラティリティ（価格変動の度合い）
• $ (W) $はウィーナープロセス（ブラウン運動）

以下はこのモデルを用いた 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

import numpy as np
import matplotlib.pyplot as plt

def generate_price(mu, sigma, S0, days, dt=0.01):
    np.random.seed(42)
    t = np.arange(0, days, dt)
    n = len(t)
    
    S = np.zeros(n)
    S[0] = S0
    
    for i in range(1, n):
        dW = np.random.normal(0, np.sqrt(dt))
        dS = mu * S[i-1] * dt + sigma * S[i-1] * dW
        S[i] = S[i-1] + dS
    
    return t, S

# パラメータの設定
mu = 0.001  # ドリフト
sigma = 0.02  # ボラティリティ
S0 = 100  # 初期価格
days = 365  # シミュレーション期間

# 価格の生成
t, price_simulation = generate_price(mu, sigma, S0, days)

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(t, price_simulation, label='仮想通貨価格')
plt.title('仮想通貨価格のGeometric Brownian Motionモデル')
plt.xlabel('時間')
plt.ylabel('価格')
plt.legend()
plt.show()

このコードでは、GBMモデルに基づいて仮想通貨の価格を生成し、Matplotlibを使用して時間と価格のグラフを描画しています。

パラメータ$ (\mu) $、$ (\sigma) $、$ (S0) $を変更することで異なる結果が得られます。

[実行結果]

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算のためのライブラリで、行列計算や乱数生成などができます。
• matplotlib.pyplotはグラフの描画に使用されます。

2. 価格生成関数の定義

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

def generate_price(mu, sigma, S0, days, dt=0.01):
    np.random.seed(42)
    t = np.arange(0, days, dt)
    n = len(t)
    
    S = np.zeros(n)
    S[0] = S0
    
    for i in range(1, n):
        dW = np.random.normal(0, np.sqrt(dt))
        dS = mu * S[i-1] * dt + sigma * S[i-1] * dW
        S[i] = S[i-1] + dS
    
    return t, S

• generate_price関数はGBMモデルに基づいて仮想通貨の価格を生成する役割を担います。
• muはドリフト（平均収益率）、sigmaはボラティリティ（価格変動の度合い）、S0は初期価格、daysはシミュレーション期間です。
• dtは時間ステップを表し、デフォルトでは$ 0.01 $に設定されています。
• 生成された価格データと対応する時間軸を返します。

3. パラメータの設定

1
2
3
4

mu = 0.001  # ドリフト
sigma = 0.02  # ボラティリティ
S0 = 100  # 初期価格
days = 365  # シミュレーション期間

• 仮想通貨価格のシミュレーションに使用されるパラメータを設定します。

4. 価格の生成

1

t, price_simulation = generate_price(mu, sigma, S0, days)

• 先ほど定義したgenerate_price関数を呼び出し、仮想通貨価格のデータを生成します。

5. グラフの描画

1
2
3
4
5
6
7

plt.figure(figsize=(10, 6))
plt.plot(t, price_simulation, label='仮想通貨価格')
plt.title('仮想通貨価格のGeometric Brownian Motionモデル')
plt.xlabel('時間')
plt.ylabel('価格')
plt.legend()
plt.show()

• matplotlib.pyplotを使用して、生成された仮想通貨価格データをグラフに描画します。
• グラフには時間に対する価格の変動が可視化されます。
• グラフのタイトルや軸ラベル、凡例なども設定されています。

グラフ解説

[実行結果]

このソースコードによって生成されるグラフは、Geometric Brownian Motion（幾何学ブラウン運動）モデルに基づいて仮想通貨の価格変動を模擬したものです。

以下に、グラフの内容を詳しく説明します。

グラフの構成要素:

1. X軸（時間）:

• グラフの横軸は時間を表しています。
• daysパラメータで指定されたシミュレーション期間内の時間軸がプロットされます。

2. Y軸（価格）:

• グラフの縦軸は仮想通貨の価格を表しています。
• 価格の単位は、モデル内での相対的な値です。

3. 曲線（仮想通貨価格）:

• プロットされた曲線は、Geometric Brownian Motionモデルに基づいて生成された仮想通貨の価格変動を示しています。
• ブラウン運動の性質に基づき、確率的な変動が取り入れられており、連続的かつランダムな変動が観察されます。

4. タイトル:

• グラフのタイトルは「仮想通貨価格のGeometric Brownian Motionモデル」です。
• 使用されている価格変動モデルがGBMであることが示唆されています。

5. 軸ラベル:

• X軸とY軸には、「時間」と「価格」という軸ラベルがそれぞれ設定されています。

6. 凡例:

• グラフには「仮想通貨価格」という凡例が設定されており、プロットされた曲線が何を表しているかを示しています。

内容の解釈:

• グラフは時間の経過に伴う仮想通貨価格の変動を視覚化しています。
• ブラウン運動のランダムな性質により、価格は一定のトレンドを持ちつつも、確率的に変動しています。
• 初期価格$（(S0)）$から始まり、時間が進むにつれて価格が変動していることが観察できます。

このグラフは、GBMモデルに基づく単純な価格変動のシミュレーションであり、実際の市場の複雑な動きを模倣しているわけではありません。

ただし、モデルの理解やパラメータの調整を通じて、価格変動の基本的な性質を観察するのに役立ちます。

ランダムウォークモデル

ランダムウォークモデルは、価格や変数がランダムに変動するモデルで、特に金融市場の価格変動などをモデル化するのによく使用されます。

ランダムウォークの基本的な方程式は以下のように表されます。

$$
[ P_t = P_{t-1} + \epsilon_t ]
$$

ここで、$ ( P_t ) $は時刻$ ( t ) $における価格、$ ( \epsilon_t ) $は平均がゼロのランダムな変数（ノイズ）です。

以下は、このランダムウォークモデルを用いて価格の変動をシミュレートし、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

import numpy as np
import matplotlib.pyplot as plt

# ランダムウォークに基づく価格変動のモデル
def random_walk(initial_price, num_steps, volatility):
    prices = [initial_price]
    for _ in range(num_steps):
        price = prices[-1] + np.random.normal(0, volatility)
        prices.append(price)
    return prices

# 初期価格
initial_price = 100

# ステップ数
num_steps = 100

# ボラティリティ（価格変動の程度）
volatility = 1.0

# 価格変動のシミュレーション
price_simulation = random_walk(initial_price, num_steps, volatility)

# グラフを描画
plt.plot(range(num_steps + 1), price_simulation)
plt.title('価格変動のランダムウォーク')
plt.xlabel('ステップ')
plt.ylabel('価格')
plt.grid(True)
plt.show()

このコードでは、ランダムウォークモデルに基づき、指定したステップ数分の価格変動をシミュレーションしています。

ボラティリティの値を変更することで、価格変動の程度を制御できます。

[実行結果]

ソースコード解説

このPythonソースコードは、ランダムウォークモデルを使用して価格変動をシミュレーションし、それをグラフで視覚化するものです。

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算をサポートするライブラリで、主に配列や行列などの操作に利用されます。
• matplotlib.pyplotはグラフ描画のためのライブラリです。

2. ランダムウォーク関数の定義:

1
2
3
4
5
6

def random_walk(initial_price, num_steps, volatility):
   prices = [initial_price]
   for _ in range(num_steps):
       price = prices[-1] + np.random.normal(0, volatility)
       prices.append(price)
   return prices

• random_walk関数は、初期価格、ステップ数、ボラティリティを受け取り、ランダムウォークモデルに基づいて価格変動をシミュレーションします。
• np.random.normal(0, volatility)は平均が0で標準偏差がvolatilityの正規分布からランダムな変数を生成しています。

3. パラメータの設定:

1
2
3

initial_price = 100
num_steps = 100
volatility = 1.0

• initial_priceは初期価格、num_stepsはステップ数、volatilityは価格変動の程度を表します。

4. 価格変動のシミュレーション:

1

price_simulation = random_walk(initial_price, num_steps, volatility)

• random_walk関数を使用して価格変動をシミュレーションし、その結果をprice_simulationに格納します。

5. グラフの描画:

1
2
3
4
5
6

plt.plot(range(num_steps + 1), price_simulation)
plt.title('価格変動のランダムウォーク')
plt.xlabel('ステップ')
plt.ylabel('価格')
plt.grid(True)
plt.show()

• plt.plotで価格変動のランダムウォークを描画します。
• plt.titleでグラフにタイトルを追加し、plt.xlabelとplt.ylabelで軸のラベルを指定します。
• plt.grid(True)でグリッドを表示し、plt.show()でグラフを表示します。

このコードは、単純なランダムウォークモデルを使用して価格変動をシミュレーションし、それをグラフで視覚化するサンプルです。

グラフ解説

[実行結果]

上記のグラフは、ランダムウォークモデルに基づく価格変動のシミュレーションを表しています。

以下はグラフの詳細な説明です。

横軸 (X軸):

ステップ数を示しています。
各ステップは時間の経過を表しており、時刻$ ( t ) $が$ 0 $から$ 100 $まで変化します。

縦軸 (Y軸):

価格を示しています。初期価格は$ 100 $からスタートしており、各ステップでランダムな変動が加わり、価格が変動します。

グラフの形状:

各ステップでの価格変動がランダムであるため、グラフはランダムに上下に動きます。
ノイズ（ランダムな変数）が価格に加わり、その結果として価格がランダムウォークする様子が観察されます。

ボラティリティの影響:

ボラティリティが$ 1.0 $と設定されているため、価格変動が比較的大きくなっています。
ボラティリティが増加すると、価格の変動が激しくなります。

このグラフは、ランダムウォークモデルを用いて生成された価格変動のシミュレーションであり、実際の金融市場の価格変動などをモデル化するために使用される手法の一つです。

天体の動きは、通常、楕円軌道を描くケプラーの法則によって記述されます。

以下に、楕円軌道上を動く天体の運動を表す方程式と、その方程式をグラフ化するPythonコードの例を示します。

$$
r(\theta) = \frac{a(1 - e^2)}{1 + e \cdot \cos(\theta)}
$$

• $ ( r(\theta) ) $は楕円軌道上の距離を表します。
• $ ( a ) $は長半径（セミメジャーアクシス）を表します。
• $ ( e ) $は離心率を表します。
• $ ( \theta ) $は楕円軌道上の角度を表します。

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 elliptical_orbit(theta, a, e):
    return (a * (1 - e**2)) / (1 + e * np.cos(theta))

# パラメータの設定
a = 5  # 長半径（セミメジャーアクシス）
e = 0.5  # 離心率

# 角度の範囲を指定
theta_values = np.linspace(0, 2 * np.pi, 1000)

# 楕円軌道上の距離を計算
r_values = elliptical_orbit(theta_values, a, e)

# 楕円軌道をグラフ化
plt.polar(theta_values, r_values, label='楕円軌道 (a={}, e={})'.format(a, e))
plt.title('楕円軌道のグラフ')
plt.legend()
plt.show()

このコードでは、ケプラーの法則に基づいた楕円軌道の方程式を用いて、指定された長半径$ (a) $と離心率$ (e) $の楕円軌道を極座標系でグラフ化しています。

[実行結果]

パラメータ$ (a) $と$ (e) $を変更することで、異なる楕円軌道を描画できます。

ソースコード解説

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

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

1
2

import numpy as np
import matplotlib.pyplot as plt

• numpyは数値計算をサポートするライブラリで、数学的な操作に使用されます。
• matplotlib.pyplotはグラフの描画に使用されます。

2. 楕円軌道の方程式の定義:

1
2

def elliptical_orbit(theta, a, e):
   return (a * (1 - e**2)) / (1 + e * np.cos(theta))

• elliptical_orbit関数は、楕円軌道の方程式を表しています。
  引数として角度 theta、長半径 a、離心率 e を取ります。

3. パラメータの設定:

1
2

a = 5  # 長半径（セミメジャーアクシス）
e = 0.5  # 離心率

• aは楕円軌道の長半径（セミメジャーアクシス）を指定します。
• eは楕円軌道の離心率を指定します。

4. 角度の範囲を指定:

1

theta_values = np.linspace(0, 2 * np.pi, 1000)

• linspace関数は、指定された範囲内で等間隔の数値を生成します。
  ここでは、$0$から$ (2\pi) $までの$1000$個の角度を生成しています。

5. 楕円軌道上の距離を計算:

1

r_values = elliptical_orbit(theta_values, a, e)

• 先ほど定義した elliptical_orbit 関数を使用して、各角度における楕円軌道上の距離を計算します。

6. 楕円軌道をグラフ化:

1
2
3
4

plt.polar(theta_values, r_values, label='楕円軌道 (a={}, e={})'.format(a, e))
plt.title('楕円軌道のグラフ')
plt.legend()
plt.show()

• plt.polar関数は極座標系でのプロットを行います。
  楕円軌道の角度と距離を指定して描画します。
• plt.titleでグラフにタイトルを追加し、plt.legendで凡例を表示します。
• plt.showでグラフを表示します。

このコードは、指定されたパラメータに基づいて楕円軌道を描画し、その形状と特性を可視化します。

グラフ解説

[実行結果]

上記のグラフは、ケプラーの法則に基づく楕円軌道を表しています。

このグラフでは、パラメータとして$ (a = 5) $（長半径が5単位）および$ (e = 0.5) $（離心率が0.5）を使用しています。

グラフの特徴:

1. 楕円の形状:

グラフは楕円の形状を示しており、楕円軌道が中心から離れたり近づいたりすることを表しています。

2. 焦点とセンター:

楕円の中心に焦点があり、この焦点からの距離が離心率によって変化します。
離心率が$0$に近いほど、楕円は円に近くなります。

3. 長半径と短半径:

長半径は楕円の中心から最も遠い点までの距離で、この例では$ (a = 5) $です。
楕円の中心から最も近い点までの距離は$ (a(1 - e)) $になります。

4. 楕円の回転:

グラフは楕円が楕円軌道上を回転していることを示しています。
楕円がどの方向に回転しているかは、楕円軌道上の角度$ ( \theta ) $によって定まります。

これらの要素により、楕円軌道がどのように形成され、どのような特性を持っているかがグラフ上で可視化されています。