コピペで学ぶポケモンPython微分積分学

ポケモンを題材にした微分積分学の学習カリキュラム案を作成しました。こちらは初期案（たたき台）としてご確認ください。

参考

# -*- coding: utf-8 -*-
# プログラム名 / Program Name: pokemon_shiny_probability_calcus.py
# 概要 / Overview:
# ポケモンの色違い確率を題材に、関数・極限・連続性・微分の概念をPythonで学ぶ
# Learn about functions, limits, continuity, and differentiation using shiny Pokémon appearance probability.

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, Limit, diff, lambdify

# --- 色違い確率の定義 / Define shiny probability ---
p = 1 / 4096  # 色違いが出る確率 / Probability of shiny Pokémon

# -------------------------------------
# 1. 関数 / Function:
# 捕まえた回数nに対して、少なくとも1匹色違いが出る確率
# Probability of getting at least one shiny after n encounters
# P(n) = 1 - (1 - p)^n
# -------------------------------------
def shiny_probability_function():
    n_vals = np.arange(1, 1001)
    P_n = 1 - (1 - p) ** n_vals

    # グラフ表示 / Plot graph
    plt.plot(n_vals, P_n)
    plt.title("Probability of At Least One Shiny Pokémon")
    plt.xlabel("Number of Encounters (n)")
    plt.ylabel("Probability")
    plt.grid()
    plt.show()

# -------------------------------------
# 2. 極限 / Limit:
# n→∞ で色違いが出る確率はどうなるか？
# What is the limit of shiny probability as n approaches infinity?
# -------------------------------------
def shiny_limit():
    n = symbols('n')
    f = 1 - (1 - p) ** n
    lim = Limit(f, n, np.inf).doit()
    print(f"[2] 極限 / Limit as n → ∞: {lim} → 必ず出る / Guaranteed shiny (Prob = 1)")

# -------------------------------------
# 3. 連続性 / Continuity:
# 関数 f(x) = 1 - (1 - p)^x は定義域内で連続
# The function is continuous because it's a composition of continuous functions
# -------------------------------------
def check_continuity():
    print("[3] f(x) = 1 - (1 - p)^x は連続関数（指数関数と定数の合成）")
    print("    f(x) = 1 - (1 - p)^x is continuous (composed of continuous functions)")

# -------------------------------------
# 4. 微分 / Differentiation:
# 関数 f(x) の導関数を求める / Find the derivative of f(x)
# f(x) = 1 - (1 - p)^x → f'(x) = ?
# -------------------------------------
def shiny_derivative():
    x = symbols('x')
    f = 1 - (1 - p) ** x
    df = diff(f, x)  # 微分 / Differentiate
    df_lamb = lambdify(x, df, 'numpy')  # 数値化 / Lambdify for plotting

    x_vals = np.linspace(1, 1000, 500)
    y_vals = df_lamb(x_vals)

    # グラフ表示 / Plot derivative
    plt.plot(x_vals, y_vals)
    plt.title("Rate of Change: Derivative of Shiny Probability")
    plt.xlabel("Number of Encounters (x)")
    plt.ylabel("f'(x): Derivative")
    plt.grid()
    plt.show()

    print(f"[4] 導関数 / Derivative f'(x) = {df}")

# -------------------------------------
# メイン関数 / Main function
# -------------------------------------
def main():
    print("=== Pokémon Shiny Probability × 関数・極限・微分 ===\n")
    shiny_probability_function()
    shiny_limit()
    check_continuity()
    shiny_derivative()

# --- 実行 / Run ---
if __name__ == "__main__":
    main()

# -*- coding: utf-8 -*-
# プログラム名 / Program Name: nn_backprop_pokemon_regression.py
# 概要 / Overview:
# 合成関数の微分、バックプロパゲーションの構造を学び、最後にポケモンの種族値に基づく回帰タスクをニューラルネットワークで実施。

import numpy as np
import matplotlib.pyplot as plt

# -------------------------------
# 1. 合成関数とチェインルールの学習
# Learn chain rule of composite functions
# -------------------------------
def composite_chain_rule():
    # y = f(g(x)), where f(u) = u^2, g(x) = 3x + 1
    x = 2
    g = lambda x: 3 * x + 1
    f = lambda u: u ** 2

    u = g(x)
    y = f(u)

    # 手計算での偏微分
    dy_dx = 6 * x + 2  # = f'(g(x)) * g'(x)

    print("[1] Composite Function: f(g(x)) = (3x+1)^2")
    print(f"    At x = 2 → f(g(x)) = {y}, dy/dx = {dy_dx}")

# -------------------------------
# 2. バックプロパゲーションのステップ確認
# Manual backpropagation steps
# -------------------------------
def manual_backprop():
    # 入力 x = [1.0, 0.5], 重みとバイアスの初期化
    x = np.array([1.0, 0.5])
    W1 = np.array([[0.1, 0.3], [0.2, 0.4]])
    b1 = np.array([0.1, 0.2])
    W2 = np.array([0.3, 0.5])
    b2 = 0.1

    # 順伝播（Forward pass）
    z1 = np.dot(x, W1) + b1
    a1 = 1 / (1 + np.exp(-z1))  # sigmoid
    z2 = np.dot(a1, W2) + b2
    y = z2  # 恒等関数

    # 損失関数（2乗誤差） L = (y - t)^2
    t = 0.5
    loss = 0.5 * (y - t) ** 2

    print(f"[2] Forward Output y = {y:.3f}, Loss = {loss:.3f}")

    # 逆伝播（Backpropagation）
    dy = y - t         # ∂L/∂y
    dW2 = dy * a1      # ∂L/∂W2
    db2 = dy
    da1 = dy * W2      # ∂L/∂a1
    dz1 = da1 * a1 * (1 - a1)  # sigmoid微分
    dW1 = np.outer(x, dz1)    # ∂L/∂W1
    db1 = dz1

    print(f"[2] Gradients (Backpropagation):\n dW1 = {dW1},\n dW2 = {dW2},\n db1 = {db1}, db2 = {db2}")

# -------------------------------
# 3. ポケモンの種族値回帰をNNで実装
# Predict Atk from [HP, Spd] using simple neural network
# -------------------------------
def pokemon_regression():
    # 種族値データ（HP, Spd） → 攻撃（Atk）を予測
    data = {
        "Pikachu": [35, 90, 55],
        "Charizard": [78, 100, 84],
        "Bulbasaur": [45, 45, 49],
        "Snorlax": [160, 30, 110],
        "Mewtwo": [106, 130, 110],
        "Dragonite": [91, 80, 134]
    }

    X = np.array([[hp, spd] for hp, spd, atk in data.values()])
    y = np.array([atk for hp, spd, atk in data.values()])

    # 正規化
    X = (X - X.mean(axis=0)) / X.std(axis=0)
    y = (y - y.mean()) / y.std()

    # パラメータ初期化
    np.random.seed(0)
    W1 = np.random.randn(2, 4) * 0.1
    b1 = np.zeros(4)
    W2 = np.random.randn(4) * 0.1
    b2 = 0.0

    # 学習設定
    lr = 0.1
    epochs = 300
    losses = []

    for epoch in range(epochs):
        # 順伝播
        z1 = np.dot(X, W1) + b1
        a1 = 1 / (1 + np.exp(-z1))  # sigmoid
        z2 = np.dot(a1, W2) + b2
        y_pred = z2

        # 損失
        loss = np.mean((y_pred - y) ** 2)
        losses.append(loss)

        # 逆伝播
        dy = 2 * (y_pred - y) / len(y)
        dW2 = np.dot(a1.T, dy)
        db2 = np.sum(dy)
        da1 = np.dot(dy[:, None], W2[None, :])
        dz1 = da1 * a1 * (1 - a1)
        dW1 = np.dot(X.T, dz1)
        db1 = np.sum(dz1, axis=0)

        # パラメータ更新
        W2 -= lr * dW2
        b2 -= lr * db2
        W1 -= lr * dW1
        b1 -= lr * db1

    # 結果可視化
    plt.plot(losses)
    plt.title("Loss over Epochs (Pokemon Atk Regression)")
    plt.xlabel("Epoch")
    plt.ylabel("Loss")
    plt.grid()
    plt.show()

    print(f"[3] 最終損失 / Final Loss: {loss:.4f}")

# -------------------------------
# メイン関数
# -------------------------------
def main():
    print("=== Neural Network × Chain Rule × Backprop × Pokemon Regression ===")
    composite_chain_rule()
    manual_backprop()
    pokemon_regression()

if __name__ == "__main__":
    main()

# プログラム名: pokemon_function_types_plot.py
# 概要: 関数の種類をポケモンの設定に例えて可視化（攻撃力・経験値・命中率など）

import numpy as np
import matplotlib.pyplot as plt

# --- 描画範囲 / Define domain ---
x = np.linspace(-10, 10, 400)
x_positive = np.linspace(0.1, 10, 400)  # 対数・指数用（正の範囲）

# --- 各関数の定義 / Define functions ---
y_linear = 5 * x + 10                                # 一次関数: 攻撃力 = 5 * レベル + 10
y_quadratic = -1 * (x - 5)**2 + 100                  # 二次関数: 命中率 = 最適スピード(5)で最大
y_polynomial = 0.1*x**3 - 1.5*x**2 + 4*x + 30        # 多項式: 複雑な経験値カーブ
y_exponential = 2 ** x_positive                     # 指数関数: 回復量が倍増
y_log = np.log(x_positive) * 10 + 20                # 対数関数: 親密度の飽和
y_sin = 50 * np.sin(x) + 100                        # 三角関数: 鳴き声の波 or ダメージ変動

# --- グラフ描画設定 / Plot settings ---
plt.figure(figsize=(14, 10))

# ① 一次関数（攻撃力 vs レベル）
plt.subplot(2, 3, 1)
plt.plot(x, y_linear, label="Attack vs Level")
plt.title("Linear Function: Attack Increases per Level")
plt.xlabel("Level")
plt.ylabel("Attack")
plt.grid(True)

# ② 二次関数（命中率 vs 素早さ）
plt.subplot(2, 3, 2)
plt.plot(x, y_quadratic, color="orange", label="Accuracy vs Speed")
plt.title("Quadratic Function: Optimal Speed for Accuracy")
plt.xlabel("Speed")
plt.ylabel("Accuracy")
plt.grid(True)

# ③ 多項式関数（経験値カーブ）
plt.subplot(2, 3, 3)
plt.plot(x, y_polynomial, color="green", label="Exp vs Level")
plt.title("Polynomial Function: Complex EXP Growth")
plt.xlabel("Level")
plt.ylabel("Experience")
plt.grid(True)

# ④ 指数関数（回復アイテム使用）
plt.subplot(2, 3, 4)
plt.plot(x_positive, y_exponential, color="red", label="Recovery Effect")
plt.title("Exponential Function: Healing Grows Fast")
plt.xlabel("Number of Potions")
plt.ylabel("HP Recovery")
plt.grid(True)

# ⑤ 対数関数（親密度・なつき度）
plt.subplot(2, 3, 5)
plt.plot(x_positive, y_log, color="purple", label="Friendship vs Time")
plt.title("Logarithmic Function: Diminishing Growth")
plt.xlabel("Time")
plt.ylabel("Friendship Level")
plt.grid(True)

# ⑥ 三角関数（鳴き声・攻撃リズム）
plt.subplot(2, 3, 6)
plt.plot(x, y_sin, color="teal", label="Sound Wave / Damage")
plt.title("Trigonometric Function: Rhythm of Actions")
plt.xlabel("Time")
plt.ylabel("Volume or Damage")
plt.grid(True)

plt.tight_layout()
plt.show()

# -*- coding: utf-8 -*-
# プログラム名 / Program Name: pokemon_logistic_taylor_classification.py
# 概要 / Overview:
# ロジスティック関数のテイラー展開（10次）と、
# 素早さ(SPD)による勝敗予測をロジスティック回帰で分類する
# Visualizes 10th-order Taylor expansion of the logistic function,
# and models win probability based on Pokémon Speed (SPD) using logistic regression.

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, series, simplify, lambdify, exp
from sklearn.linear_model import LogisticRegression

# -------------------------------
# 1. ロジスティック関数と10次テイラー展開の比較 / Logistic function vs 10th-order Taylor expansion
# -------------------------------
def logistic_taylor():
    x = symbols('x')
    f = 1 / (1 + exp(-x))  # ロジスティック関数 / Logistic function
    taylor_10 = series(f, x, 0, 11)  # 10次までのテイラー展開 / 10th-order Taylor expansion
    taylor_10 = simplify(taylor_10.removeO())

    f_np = lambdify(x, f, 'numpy')
    taylor_np = lambdify(x, taylor_10, 'numpy')
    x_vals = np.linspace(-6, 6, 400)

    plt.plot(x_vals, f_np(x_vals), label="Sigmoid", color='blue')
    plt.plot(x_vals, taylor_np(x_vals), label="Taylor (10th)", color='red', linestyle='--')
    plt.title("Logistic Function vs 10th Taylor Approximation")
    plt.xlabel("x")
    plt.ylabel("σ(x)")
    plt.legend()
    plt.grid()
    plt.show()

    return str(taylor_10)

# -------------------------------
# 2. ロジスティック回帰：素早さ（SPD）に基づく勝敗予測 / Logistic regression based on SPD
# -------------------------------
def logistic_regression_pokemon():
    # 仮想データ: SPD と WIN/LOSE / Synthetic data: SPD and WIN/LOSE
    data = {
        "Pikachu":     [90, 1],
        "Charizard":   [100, 1],
        "Bulbasaur":   [45, 0],
        "Snorlax":     [30, 0],
        "Mewtwo":      [130, 1],
        "Dragonite":   [80, 1],
        "Geodude":     [20, 0],
        "Machamp":     [55, 0],
        "Gengar":      [110, 1],
        "Jigglypuff":  [20, 0]
    }

    X = np.array([[v[0]] for v in data.values()])  # SPD
    y = np.array([v[1] for v in data.values()])    # WIN(1)/LOSE(0)

    model = LogisticRegression()
    model.fit(X, y)

    x_plot = np.linspace(0, 150, 300).reshape(-1, 1)
    y_proba = model.predict_proba(x_plot)[:, 1]

    plt.scatter(X, y, label="Battle Outcome", color='black')
    plt.plot(x_plot, y_proba, label="Logistic Regression", color='green')
    plt.title("Logistic Regression: Win Probability vs Speed")
    plt.xlabel("Speed (SPD)")
    plt.ylabel("Probability of Winning")
    plt.grid()
    plt.legend()
    plt.show()

    coef = model.coef_[0][0]
    intercept = model.intercept_[0]
    return coef, intercept

# 実行 / Execute
taylor_expression = logistic_taylor()
coef, intercept = logistic_regression_pokemon()
taylor_expression, coef, intercept

# -*- coding: utf-8 -*-
# プログラム名: pokemon_multivariable_extrema.py
# 概要:
# ポケモンの種族値(HP, SPD)と攻撃力を使って、偏微分・全微分・極値・条件付き極値を学ぶ

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, diff, solve, Eq, lambdify, Matrix, hessian

# --- シンボル定義 / Define symbols ---
x, y, λ = symbols('x y λ')  # x: HP, y: SPD

# --- 関数定義 / Define function ---
# ホントは正しくないけど、ダミーで作った攻撃力をモデル化した関数 z = f(x, y) / Attack power based on HP and SPD
f = -0.01*(x - 100)**2 - 0.02*(y - 80)**2 + 120  # 凸関数で最大値あり

# 1. 偏微分 / Partial derivatives
df_dx = diff(f, x)
df_dy = diff(f, y)

# 2. 極値（臨界点）/ Find critical point
critical = solve([df_dx, df_dy], [x, y])
x0, y0 = critical[0]
z0 = f.subs({x: x0, y: y0})

# 3. ヘッセ行列と正負判定（極値判定）/ Hessian matrix to classify extrema
H = hessian(f, (x, y)).subs({x: x0, y: y0})
H_np = np.array(H).astype(np.float64)
eigenvals = np.linalg.eigvals(H_np)
extrema_type = "最大" if all(val < 0 for val in eigenvals) else "最小 or 鞍点"

# 4. 全微分 / Total differential
dz = df_dx * (x - x0) + df_dy * (y - y0)

# 5. 合成関数の微分（u = x+y, v = x−y）/ Composite function: f(u(x,y), v(x,y))
u = x + y
v = x - y
g = f.subs({x: u/2 + v/2, y: u/2 - v/2})  # 合成関数
dg_du = diff(g, u)
dg_dv = diff(g, v)

# 6. 陰関数定理の利用 / Implicit function theorem
# g(x, y) = f(x, y) - C = 0 の形で y を x の関数として微分
C = z0
g_implicit = f - C
dy_dx_implicit = -diff(g_implicit, x) / diff(g_implicit, y)

# 7. 条件付き極値（Lagrange）/ Lagrange multiplier method
# 制約: x + y = 180（総合種族値制限）
L = f + λ * (180 - x - y)
Lx = diff(L, x)
Ly = diff(L, y)
Lλ = diff(L, λ)
lagrange_sol = solve([Lx, Ly, Lλ], [x, y, λ])
x_l, y_l = lagrange_sol[0][0:2]
f_l = f.subs({x: x_l, y: y_l})

# --- 結果出力 / Print results ---
print("=== Pokémon × Multivariable Calculus ===\n")
print(f"[1] 偏微分 df/dx: {df_dx}, df/dy: {df_dy}")
print(f"[2] 極値候補: HP = {x0}, SPD = {y0}, Atk = {z0:.2f}（{extrema_type}）")
print(f"[3] 全微分（変化量近似）: dz = {dz}")
print(f"[4] 合成関数の偏微分: ∂f/∂u = {dg_du}, ∂f/∂v = {dg_dv}")
print(f"[5] 陰関数定理による dy/dx: {dy_dx_implicit}")
print(f"[6] 条件付き極値（x + y = 180）→ HP = {x_l}, SPD = {y_l}, Atk = {f_l:.2f}")

# -*- coding: utf-8 -*-
# プログラム名: pokemon_double_integral_stats.py
# 概要:
# ポケモンのHPとAtkの仮想的な連続確率分布に基づき、2重積分と累次積分を使って確率を求める

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, exp, integrate

# --- 変数定義 / Define variables ---
x, y = symbols('x y')  # x: HP, y: Atk

# --- 2変数の正規化された確率密度関数（仮定） ---
# HP (x): 0~200, Atk (y): 0~200 を対象とした分布モデル
λx = 1/100
λy = 1/80
f_xy = λx * λy * exp(-λx * x - λy * y)

# --- 累次積分：領域 x=50~150, y=60~140 の確率 ---
inner = integrate(f_xy, (y, 60, 140))            # yについて先に積分
prob = integrate(inner, (x, 50, 150))            # xについて積分（累次積分完了）

# --- 可視化（領域のハッチ表示） ---
x_vals = np.linspace(0, 200, 100)
y_vals = np.linspace(0, 200, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = λx * λy * np.exp(-λx * X - λy * Y)

plt.contourf(X, Y, Z, levels=50, cmap='plasma')
plt.colorbar(label="f(x, y)")
plt.title("Joint PDF of HP and Atk (Laplace-like)")
plt.xlabel("HP (x)")
plt.ylabel("Attack (y)")
plt.axvline(50, color='white', linestyle='--')
plt.axvline(150, color='white', linestyle='--')
plt.axhline(60, color='white', linestyle='--')
plt.axhline(140, color='white', linestyle='--')
plt.grid(alpha=0.4)
plt.tight_layout()
plt.show()

# --- 結果出力 / Print result ---
print(f"[確率] HP ∈ [50, 150] かつ Atk ∈ [60, 140] に入る確率: {prob.evalf():.4f}")

# -*- coding: utf-8 -*-
# プログラム名: pokemon_shiny_probability_grad.py
# 内容: 色違いが出る確率のPyTorchによる微分計算

import torch

# --- 色違い確率 / Shiny rate ---
p = 1 / 4096

# --- 捕獲回数（連続変数として）/ Number of encounters ---
n = torch.tensor(100.0, requires_grad=True)  # requires_grad=True で微分対象

# --- 色違いが出る確率 P(n) = 1 - (1 - p)^n ---
P = 1 - (1 - p) ** n

# --- 自動微分 / Compute gradient ---
P.backward()

# --- 出力 / Display ---
print("=== PyTorch × Shiny Probability Differentiation ===")
print(f"Number of attempts (n): {n.item()}")
print(f"Shiny probability P(n): {P.item():.6f}")
print(f"dP/dn (Gradient): {n.grad.item():.6f}")

# プログラム名: gaussian_integral_pokemon_stats.py
# 概要: ポケモンの種族値を正規分布とガウス積分の視点から解析・可視化

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad

# --- 平均と標準偏差を仮定 / Assume mean and std of Attack stat ---
mu = 80      # 平均攻撃種族値 / Mean attack base stat
sigma = 15   # 標準偏差 / Standard deviation

# --- x軸の定義 / X-axis range ---
x = np.linspace(mu - 4*sigma, mu + 4*sigma, 400)
pdf = norm.pdf(x, mu, sigma)  # 正規分布の確率密度関数 / PDF of normal distribution

# --- プロット / Plot PDF ---
plt.figure(figsize=(10, 6))
plt.plot(x, pdf, label="Normal Distribution of Attack Stat", color="purple")
plt.title("Gaussian Distribution of Pokemon Attack Stats")
plt.xlabel("Attack Base Stat")
plt.ylabel("Probability Density")
plt.grid(True)

# --- 指定区間の確率（例：攻撃力100〜120）/ Probability for Attack in [100,120] ---
a, b = 100, 120
area, _ = quad(lambda x: norm.pdf(x, mu, sigma), a, b)

# --- 塗りつぶし / Shade area under curve ---
x_fill = np.linspace(a, b, 200)
plt.fill_between(x_fill, norm.pdf(x_fill, mu, sigma), alpha=0.4, color="orange", label=f"P({a} ≤ Attack ≤ {b}) ≈ {area:.3f}")

plt.legend()
plt.show()

# --- 全体積確認（ガウス積分）/ Total integral check over all real numbers ---
total, _ = quad(lambda x: norm.pdf(x, mu, sigma), -np.inf, np.inf)
print(f" ガウス積分による全体確率: {total:.5f}（理論上は1.0）")
print(f" 攻撃力が {a}〜{b} の範囲にいる確率 ≈ {area:.3f}")

#pikachu_thunderbolt_with_noise_power.py
# --- ライブラリのインポート / Import libraries ---
import numpy as np
import matplotlib.pyplot as plt

# --- 時間軸（0〜10ms）/ Time axis ---
t = np.linspace(0, 0.01, 1000)  # 秒単位（1000サンプル）

# --- 電圧定義 / Voltage definition ---
pulse = 100000 * np.exp(-500 * t)                # 主パルス（10万ボルト減衰）
noise = 300 * np.sin(2 * np.pi * 8000 * t)      # 高周波ノイズ
voltage = pulse + noise                          # 合成波形

# --- 抵抗を仮定（例: 100Ω）/ Assume resistive load ---
R = 100  # Ω（電気を受ける相手の抵抗）

# --- 瞬時電力 P(t) = V(t)^2 / R ---
power = voltage**2 / R

# --- 累積エネルギー（電力量）E = ∫ P(t) dt ---
energy = np.trapz(power, t)  # [J] ジュール（ワット秒）

# --- グラフ1: 電圧波形 / Voltage waveform ---
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(t * 1000, voltage, label='Thunderbolt Signal (V)', color='orange')
plt.title("Pikachu's Thunderbolt (with Noise)")
plt.xlabel("Time (ms)")
plt.ylabel("Voltage (V)")
plt.grid(True)
plt.legend()

# --- グラフ2: 電力波形 / Instantaneous Power ---
plt.subplot(2, 1, 2)
plt.plot(t * 1000, power, label='Instantaneous Power (W)', color='red')
plt.xlabel("Time (ms)")
plt.ylabel("Power (W)")
plt.title(f"Instantaneous Power – Total Energy: {energy:.2f} J")
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

# プログラム名: voltorb_moment_of_inertia.py
# 内容: ビリリダマを球として z軸周りの慣性モーメントを三重積分で求める

import sympy as sp

# --- 記号定義 ---
r, theta, phi, R, rho = sp.symbols('r theta phi R rho', real=True, positive=True)

# --- 極座標変換 x² + y² = r² sin²θ ---
integrand = (r**2 * sp.sin(theta)**2) * rho * r**2 * sp.sin(theta)

# --- 積分範囲 ---
# φ: 0 to 2π, θ: 0 to π, r: 0 to R
I = sp.integrate(
    sp.integrate(
        sp.integrate(integrand, (r, 0, R)),
        (theta, 0, sp.pi)
    ),
    (phi, 0, 2 * sp.pi)
)

I = sp.simplify(I)
print("[Voltorb 球モデルによる z軸慣性モーメント]")
print("I_z = ", I)