複素数関数 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

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

# グリッドを作成
real = np.linspace(-5, 5, 100)
imag = np.linspace(-5, 5, 100)
real, imag = np.meshgrid(real, imag)
z = real + 1j * imag

# 方程式
eqn = z**3 - 1  # 方程式の例

# 3Dグラフの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(np.real(z), np.imag(z), np.real(eqn), c=np.real(eqn), cmap='viridis')

# グラフのラベル設定
ax.set_xlabel('Real axis')
ax.set_ylabel('Imaginary axis')
ax.set_zlabel('Real part of the equation')
ax.set_title('Graph of a complex equation')

plt.show()

この例では、$ (z^3 - 1) $という複素数の3次方程式を使用しています。

これは複素数平面上の方程式です。

3Dプロットは、実部と虚部の関係を示し、方程式の解の振る舞いを可視化します。

[実行結果]

ソースコード解説

このコードは、複素平面上で$ (z^3 - 1) $のグラフを3Dで描画しています。

ここでの$ (z) $は、複素数$ (z = x + yi) $で表され、実部が$ (x) $軸、虚部が$ (y) $軸にマッピングされています。

1. グリッドの作成:

np.linspace を用いて、実数と虚数の範囲を指定し、それぞれ$ (x) $軸および$ (y) $軸の値を生成し、meshgrid を使って複素平面上の格子点を作成しています。

2. 方程式の定義:

z**3 - 1 という方程式を定義し、eqn に格納しています。これは$ (z^3 - 1 = 0) $という方程式を表しています。

3. 3Dグラフの描画:

matplotlib を用いて、3Dの散布図を描画しています。
scatter 関数を使い、複素数$ (z) $の実部と虚部に対応する点をプロットし、色は$ (z^3 - 1) $の実部によって変わります。

4. 軸ラベルとタイトルの設定:

set_xlabel、set_ylabel、set_zlabel を用いてそれぞれ$ (x) $軸、$ (y) $軸、$ (z) $軸のラベルを設定し、set_title でグラフのタイトルを設定しています。

5. グラフの表示:

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

このグラフは、複素平面上で$ (z^3 - 1) $の方程式を表し、その実部がゼロになる$ (z) $の値を表しています。

それにより得られる図形は、3次元曲線になります。

結果解説

[実行結果]

このグラフは、複素平面上で方程式$ ( z^3 - 1 ) $を表現しています。

グラフは複素数 $ ( z = x + yi ) $の実部、虚部、および方程式 $ ( z^3 - 1 ) $の実部を3Dで表現しています。

具体的には、複素平面内で実部が$ ( x ) $軸、虚部が$ ( y ) $軸、そして$ ( z^3 - 1 ) $の実部が$ ( z ) $軸にマッピングされています。

グラフ上の点は複素数 $ ( z ) $の実部と虚部に基づいて配置され、色は$ ( z^3 - 1 ) $の実部を示しています。

このグラフは$ ( z^3 - 1 ) $の実部が$ ( 0 ) $になるような$ ( z ) $の値を示しており、その結果得られる図形は特定の形を持つ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

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

# ロジスティック方程式（捕食者-被食者モデル）
def predator_prey_system(y, t, alpha, beta, delta, gamma):
    x, y = y
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# パラメータ設定
alpha = 0.1  # 捕食者の増加率
beta = 0.02  # 捕食者と被食者の相互作用率
delta = 0.02  # 被食者の増加率
gamma = 0.1  # 捕食者の減少率

# 初期値
initial_population = [40, 9]  # 捕食者と被食者の初期数

# 時間
t = np.linspace(0, 200, 1000)

# モデルの解を計算
sol = odeint(predator_prey_system, initial_population, t, args=(alpha, beta, delta, gamma))
predator, prey = sol[:, 0], sol[:, 1]

# グラフ化
plt.figure(figsize=(8, 6))
plt.plot(t, predator, label='Predator')
plt.plot(t, prey, label='Prey')
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Predator-Prey Dynamics')
plt.legend()
plt.grid(True)
plt.show()

このコードは、捕食者（Predator）と被食者（Prey）の数の時間変化を示します。

[実行結果]

捕食者と被食者の相互作用を表すパラメータを調整することで、弱肉強食の関係を模倣することができます。

ソースコード解説

このコードは、捕食者-被食者（Predator-Prey）モデルを用いて、捕食者と被食者の数の時間変化をシミュレーションしています。

これを行うために、次の手順が含まれています：

1. 関数の定義:

predator_prey_system 関数では、捕食者と被食者の数の時間変化を表す微分方程式（ロジスティック方程式）が定義されています。
この方程式は、捕食者と被食者の増減率、相互作用率に基づいて、それぞれの数の変化を記述します。

2. パラメータ設定:

モデル内の各パラメータ（捕食者の増加率、相互作用率など）が定義され、モデルに適用されます。
これらのパラメータは生態系内での捕食者と被食者の相互作用を表現します。

3. 初期値の設定:

捕食者と被食者の初期値が設定されます。
この値は、シミュレーションの開始時点での捕食者と被食者の数を表します。

4. 時間の定義:

シミュレーションする時間軸が定義されます。
この場合、$0$から$200$までの時間を$1000$ステップで区切っています。

5. 微分方程式の数値解析:

odeint を使用して微分方程式を解きます。
初期値、時間、および定義した微分方程式の関数が渡され、解析された結果が sol に格納されます。

6. グラフ化:

得られた結果をグラフに描画します。
時間に対する捕食者と被食者の数がそれぞれ線グラフとして表示され、どのように増減するかが視覚化されます。

このコードを実行することで、捕食者と被食者の数の時間変化を表すグラフが表示されます。

時間が経過するにつれて、捕食者と被食者の数がどのように変化するかが可視化されます。

結果解説

[実行結果]

このグラフは、捕食者（Predator）と被食者（Prey）の数の時間変化を示しています。
捕食者と被食者の関係を表すロジスティック方程式を用いて、彼らの相互作用をモデル化しました。

• 初期値に基づき、捕食者の数は増加し、被食者の数もそれに応じて増加します。
• しかし、被食者の増加に伴って捕食者の数が増え、それによって被食者の数が減少します。
• その後、捕食者の数が減少すると被食者の数が再び増加します。

このように、捕食者と被食者の数は時間とともに振動する関係を持ちます。

捕食者と被食者の相互作用率や増加率などのパラメータを調整すると、この関係性がどのように変化するかを観察できます。
このようなモデルは、生態系における捕食者と被食者の関係性やそのダイナミクスを理解するために使用されます。

栄養学 最適化問題

栄養学の最適化問題の一例として、食事計画を最適化する問題を考えてみましょう。

例えば、特定の栄養素を最大化または最小化するという目標を持った食事プランを作成することができます。

以下のコードは、特定の栄養素（タンパク質、脂質、炭水化物など）を最大化する食事プランを最適化する問題を扱います。

SciPyのminimize関数を使用します。

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

from scipy.optimize import minimize

# 目的関数：最大化する栄養素
def objective_function(x):
    # 栄養素の吸収率（例としてランダムな値を使用）
    protein_absorption = 0.8
    fat_absorption = 0.7
    carb_absorption = 0.9
    
    # 栄養素の含有量（仮の値を使用）
    protein_content = x[0]
    fat_content = x[1]
    carb_content = x[2]
    
    # 含有する栄養素の量に対して吸収率を乗算して目的関数を設定
    return -(protein_content * protein_absorption + fat_content * fat_absorption + carb_content * carb_absorption)

# 初期推定値
initial_guess = [50, 80, 200]  # タンパク質、脂質、炭水化物の初期量

# 制約条件
constraints = ({'type': 'ineq', 'fun': lambda x: 200 - sum(x)},  # 合計カロリー制約
               {'type': 'ineq', 'fun': lambda x: x[0] - 50},     # タンパク質最低量制約
               {'type': 'ineq', 'fun': lambda x: x[1] - 60},     # 脂質最低量制約
               {'type': 'ineq', 'fun': lambda x: x[2] - 150})    # 炭水化物最低量制約

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

# 最適化結果
optimal_meal_plan = result.x
maximized_nutrient = -result.fun

print(f"Optimal Meal Plan (Protein, Fat, Carbs): {optimal_meal_plan}")
print(f"Maximized Nutrient: {maximized_nutrient}")

このコードは、最大化したい栄養素を設定し、各栄養素の含有量を変数として最適化します。

制約条件を設定し、合計カロリーと各栄養素の最低摂取量を制限しています。

結果として、最適な食事プランと最大化された栄養素量が出力されます。

[実行結果]

Optimal Meal Plan (Protein, Fat, Carbs): [ 35.  45. 135.]
Maximized Nutrient: 181.0

ソースコード解説

このコードは、栄養学の観点から最適な食事プランを見つける最適化問題を扱います。

1. 目的関数の定義:

• objective_function関数は、最大化したい栄養素（タンパク質、脂質、炭水化物など）を表しています。
• 各栄養素の含有量とそれぞれの吸収率を考慮し、それらの積の総和を最大化することがこの関数の目的です。

2. 初期推定値:

• initial_guessは、最適化アルゴリズムによって使用される栄養素の初期量を示しています。

3. 制約条件:

• constraintsでは、合計カロリーと各栄養素の最低摂取量に関する制約条件を設定しています。
• この場合、合計カロリーは200未満でなければならず、タンパク質、脂質、炭水化物はそれぞれ50g、60g、150g以上である必要があります。

4. 最適化の実行:

• minimize関数を使用して、定義した目的関数と制約条件を考慮して最適な栄養素を見つけます。
• ここではCOBYLA法を使用して最適化を行っています。

5. 最適化結果の表示:

• 最適な食事プラン（タンパク質、脂質、炭水化物の含有量）と、それによって最大化された栄養素の量が出力されます。

結果解説

この結果は、特定の栄養素（タンパク質、脂質、炭水化物など）を最大化する食事プランを見つけたことを示しています。

[実行結果]

Optimal Meal Plan (Protein, Fat, Carbs): [ 35.  45. 135.]
Maximized Nutrient: 181.0

• 最適な食事プラン:
  最適な食事プランでは、タンパク質が35g、脂質が45g、炭水化物が135g含まれています。
  これは、目標としていた栄養素を最大化するために、それぞれの栄養素が含まれる量を示しています。
  これらの値は最適化アルゴリズムが見つけた最良のバランスを表しています。

• 最大化された栄養素量:
  最大化された栄養素の値は181.0です。
  この値は、目的関数が最大化した栄養素の総量を示しています。
  最適化の目標は、この数値をできるだけ高くすることでした。

宇宙工学 最適化問題

宇宙工学の最適化問題の一例として、宇宙船の軌道設計を考えてみましょう。

例えば、地球周回軌道への衛星の打ち上げを考えます。

具体的には、衛星の打ち上げコストと到達する軌道の高度、または到達するまでの時間などの複数の要素を最適化する問題があります。

ここでは、簡単な例として、地球周回軌道への衛星の打ち上げコストと軌道の高度の関係を考えてみます。

軌道の高度が増すにつれて打ち上げコストが増えると仮定します。

以下のコードは、この問題を最適化するためにScipyライブラリを使用しています。

具体的には、打ち上げコストを最小化するための最適な軌道高度を探します。

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.optimize import minimize

# 目的関数：打ち上げコストの最小化
def launch_cost(height):
    # 仮想的な打ち上げコスト関数（ここでは単純な2次関数）
    return (height - 500) ** 2  # 高度と500kmの差の二乗をコストと仮定

# 初期推定値
initial_guess = 1000  # 初期の軌道高度の推定値

# 最適化の実行
result = minimize(launch_cost, initial_guess, method='BFGS')

# 最適化結果
optimal_height = result.x
min_launch_cost = result.fun

# グラフ化
heights = np.linspace(400, 600, 100)  # 400kmから600kmまでの高度を考慮
costs = (heights - 500) ** 2  # 打ち上げコスト関数

plt.plot(heights, costs, label='Launch Cost')
plt.scatter(optimal_height, min_launch_cost, color='red', label='Optimal Height')
plt.xlabel('Orbital Height (km)')
plt.ylabel('Launch Cost')
plt.title('Optimization of Satellite Launch Height')
plt.legend()
plt.grid(True)
plt.show()

print(f"Optimal Height for Minimum Launch Cost: {optimal_height} km")
print(f"Minimum Launch Cost: {min_launch_cost}")

このコードは、打ち上げコストを最小化するための最適な軌道高度を見つけ、その結果をグラフ化しています。

グラフでは、軌道高度に対する打ち上げコストの関係が示され、最適な軌道高度が赤い点として示されます。

[実行結果]

ソースコード解説

このコードは、宇宙工学における最適化問題の一例を示しています。

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

1. 目的関数の定義:

• launch_cost関数は、仮想的な打ち上げコスト関数を定義しています。
  ここでは、高度と500kmの差の二乗をコストとしています。

2. 初期推定値の設定:

• initial_guessは、最適化アルゴリズムの初期の推定値として1000kmを設定しています。

3. 最適化の実行:

• minimize関数を使って、launch_cost関数を最小化することで最適な軌道高度を探します。
  ここではBFGS法を用いて最適化を行っています。

4. 最適化結果の取得:

• 最適な軌道高度とその時の打ち上げコストを取得しています。

5. グラフ化:

• heightsは400kmから600kmまでの範囲の高度を考慮するための配列です。
• costsはそれぞれの高度における打ち上げコストを示します。
• plt.plot()で高度に対する打ち上げコストの関係を表すグラフをプロットしています。
• plt.scatter()で最小の打ち上げコストを達成する最適な高度を赤い点として表示しています。
• その他、グラフのタイトルや軸ラベル、凡例、グリッドを設定しています。

6. 結果の表示:

• 最適な軌道高度とその時の最小打ち上げコストを出力しています。

結果解説

[実行結果]

この結果は、打ち上げコストを最小化するための最適な軌道高度を見つけたことを示しています。
結果として得られた最適な軌道高度は約500kmです。
しかし、数値計算の誤差のため、厳密には499.99999999kmという値が得られています。

さらに、この最適な高度での打ち上げコストは非常に小さく、ほぼゼロに近い値となっています（6.48426847315948e-17）。
これは浮動小数点の誤差などによるもので、ほぼゼロに等しい非常に小さな値です。

グラフでは、軌道高度が高くなるにつれて打ち上げコストが増加することが示されています。
そして、その関係を表す2次関数的なカーブが描かれており、最小値が500km付近であることが視覚的に示されています。
赤い点は最適な高度での最小コストを示しています。

リカードの比較優位

貿易理論には、複数の方程式やモデルが存在しますが、最も基本的なものの一つにリカードの比較優位に関連する方程式があります。

リカードの比較優位理論は、異なる国が特定の商品を生産する際の優位性を分析します。

これには特に次の数式が関連しています。

ここで、$ (P_x) $と$ (P_y) $はそれぞれ商品$ (X) $と商品$ (Y) $の価格を示し、$ (A_x) $と$ (A_y) $は生産コストまたは生産に要する労働時間などを示します。

この不等式が成り立つ場合、つまり商品$ (X) $の国内価格が国際価格よりも低いとき、その国は商品$ (X) $の生産に比較的適しており、他国との間で特化と貿易を行うことで両国とも利益を最大化できる可能性が高くなります。

この方程式はリカードの比較優位を示すものであり、貿易理論における基本的なモデルの一部です。

例題

リカードの比較優位に関連した最適化問題の一例として、2つの国（国Aと国B）が2つの商品（商品Xと商品Y）を生産する場合を考えてみましょう。

以下はその問題設定です。

国Aで商品Xの生産に1時間かかり、商品Yの生産に2時間かかるとします。
国Bでは商品Xの生産に3時間かかり、商品Yの生産に1時間かかるとします。
また、国際価格における商品Xと商品Yの価格比率が2:1であるとします。

この場合、どの国がどの商品を生産すべきか、そして両国間での貿易が有益かを最適化する問題を考えることができます。

以下はPythonでこの問題を解くコードです。
ここではScipyライブラリを使って最適化を行い、結果をグラフ化します。

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
from scipy.optimize import minimize

# 各国の生産にかかる時間（単位：時間）
# [国Aの商品X生産時間, 国Aの商品Y生産時間]
# [国Bの商品X生産時間, 国Bの商品Y生産時間]
production_times = np.array([[1, 2], [3, 1]])

# 国際価格比率（商品X価格 / 商品Y価格）
international_price_ratio = 2

# 目的関数：貿易利益の最大化（国Aと国Bの生産量）
def objective_function(production_quantities):
    trade_profit = -np.sum(production_quantities * production_times)
    return trade_profit

# 制約条件：国際価格比率と生産量の関係
def constraint_function(production_quantities):
    return international_price_ratio - production_quantities[0] / production_quantities[1]

# 初期推定値
initial_guess = [1, 1]

# 最適化の実行
constraint = {'type': 'eq', 'fun': constraint_function}
result = minimize(objective_function, initial_guess, constraints=constraint)

# 最適化結果
optimal_production_quantities = result.x
trade_profit = -result.fun

# グラフ化
labels = ['国Aの商品X', '国Aの商品Y']
fig, ax = plt.subplots()
ax.bar(labels, optimal_production_quantities)
ax.set_ylabel('生産量')
ax.set_title('最適生産量')
plt.show()

print(f"最適な生産量: {optimal_production_quantities}")
print(f"貿易利益: {trade_profit}")

このコードは、Scipyのminimizeを使って貿易利益を最大化する生産量を見つけ、結果を棒グラフで表示します。

最適な生産量と貿易利益が表示されます。

ただし、このコードは簡略化されたモデルですので、実際の貿易問題には多くの要因が含まれる場合があります。

[実行結果]

ソースコード解説

このコードは、2つの国（国Aと国B）が2つの商品（商品Xと商品Y）を生産する場合の貿易最適化問題を解くためのPythonスクリプトです。

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

1. 各国の生産時間と国際価格比率の設定:

• production_timesは2つの国での商品の生産時間を表す行列です。
  例えば、国Aでは商品Xの生産に1時間かかり、商品Yの生産に2時間かかります。
• international_price_ratioは商品Xと商品Yの価格比率を示します。

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

• objective_functionは貿易利益を最大化するための目的関数です。
  生産量と生産時間を掛け合わせ、その合計をマイナスにしています。
• constraint_functionは国際価格比率と生産量の関係を示す制約条件です。

3. 最適化の実行:

• minimize関数を使って、目的関数を最大化するための最適化を行います。
  初期推定値はinitial_guessで設定されており、制約条件はconstraintに定義されています。

4. 最適化結果の取得:

• result.xは最適な生産量を示しています。
• result.funは最適な目的関数の値で、貿易利益を示しています。

5. グラフ化:

• 最適な生産量を棒グラフとして表示しています。

最終的に、このコードは貿易利益を最大化するための最適な生産量を計算し、それをグラフ化して表示しています。

また、最適な生産量と貿易利益をコンソールに出力しています。

結果解説

[実行結果]

最適な生産量は、国Aが商品Xをおよそ1.37兆単位、商品Yをおよそ0.69兆単位生産することを示しています。
一方、国Bは商品Xをほぼ0単位、商品Yをおよそ0.69兆単位生産することを示しています。

この生産量の結果、貿易利益は755兆USDに相当します。
この利益は、国際価格比率と生産コストの違いを利用して、国Aが商品Xを生産し、国Bが商品Yを生産することで生じる貿易利益を表しています。

しかし、結果が非常に大きな数字になっており、生産量のバランスや制約条件の設定に誤差がある可能性があります。
このような場合、モデルや数値計算の調整が必要になるかもしれません。

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

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) $が表示されており、それぞれのエネルギー状態が異なるラインで示されています。