栄養学最適化問題 SciPy

January 8, 2024

栄養学最適化問題

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

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

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

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

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

January 7, 2024

宇宙工学最適化問題

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

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

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

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

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

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

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

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付近であることが視覚的に示されています。
赤い点は最適な高度での最小コストを示しています。

リカードの比較優位 SciPy

January 6, 2024

リカードの比較優位

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

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

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

ここで、$ (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ライブラリを使って最適化を行い、結果をグラフ化します。

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を生産することで生じる貿易利益を表しています。

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

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

January 5, 2024

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

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

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

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

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

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

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つの振動子が互いに相互作用するヴァンデル・ポルト方程式の振る舞いを視覚化しており、非線形ダイナミクスが示されています。

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

サドル点

January 4, 2024

サドル点

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

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

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次元のグラフとして表示されます。

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

January 3, 2024

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

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

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

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

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

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

モビウス帯

January 2, 2024

モビウス帯

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

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

すると、通常の帯では表と裏がありますが、モビウス帯では一周しても片側だけになります。
これは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を使用して三次元プロットを行います。

以下はモビウス帯を描画するサンプルコードです。

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を使用してモビウス帯をプロットしています。

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

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

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 のデータからモビウス帯の表面をプロットします。

グラフの設定と表示

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 の値を変更することで、異なる半径や帯の幅のモビウス帯を描画することができます。

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

January 1, 2024

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

シュレーディンガー方程式は量子力学における基本的な方程式の一つであり、粒子の波動関数の時間変化を記述します。

しかし、シュレーディンガー方程式を正確に解くことは通常、解析的には不可能であるため、数値的な近似解を求めることが一般的です。

以下に、一次元の無限井戸ポテンシャル下でのシュレーディンガー方程式の時間依存性を無視した定常状態の解を求めるPythonのコードを示します。

このコードは、一次元の無限井戸の中を動く粒子の波動関数を計算し、それをプロットします。

無限井戸は物理学における理想化されたモデルで、井戸の中を動く粒子は井戸の壁を超えて脱出することはありません。

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. ライブラリのインポート

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. 定数の設定

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. エネルギー固有値と波動関数の計算

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. 波動関数の正規化とプロット

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

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

December 31, 2023

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

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

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

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

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

1	!pip install cvxpy matplotlib numpy

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

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グラフ） CVXPY

December 30, 2023

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

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

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

目的関数は$ x^2 + y^2 $で、制約条件は$ x + y <= 1 $です。

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で可視化し、最適解を視覚的に示しています。