More than 1 year has passed since last update.

Pythonとコピペと生成AIで学ぶ名問の森

Last updated at 2025-07-25Posted at 2025-07-17

はじめに

「名門の森」は、物理体系的理解と演習を重視した問題集です。本記事では、この「名門の森」をPythonと生成AI（ChatGPTなど）を活用しながら、より実践的・再現可能な方法で学習・再構築していくことを目的としています。

◆ 学習の進め方

学習の基本方針は次の通りです：

Google Colab 上で Python コードを使いながら、数式やグラフ、計算結果をインタラクティブに確認する。ブラウザとGoogleアカウントがあれば利用でき、ノートブック形式で記録と再現が可能です。
わからないところや式展開の途中などは、生成AI（例：ChatGPT）に式や質問を入力して補完する。途中式の導出、変数定義、定義の再確認などに活用します。
理解した内容や確認できた数式・コード・考察は、自分の言葉でプレーンテキストまたはMarkdown＋コードセル付きのノート形式で整理する。これにより、検索・再利用・共有が容易になります。

問題1：放物運動

# ============================================================
# プログラム名：monkey_hunting_simulation.py
# Program Title: Monkey Hunting Simulation – Does the projectile hit?
# ============================================================

import numpy as np
import matplotlib.pyplot as plt

# ---------------------------------------------------------------
# ■ 問題（日本語）
# 地上の射手が、木の上にいるサルを狙って物体を角度θで発射する。
# 同時にサルは枝から落下（自由落下）を始める。
# 射手がサルの初期位置を正確に狙えば、サルに命中するかを確認せよ。
#
# ■ Problem (English)
# A projectile is launched at a monkey hanging from a tree.
# The monkey lets go and begins to fall at the same moment the projectile is fired.
# If the initial velocity vector is aimed exactly at the monkey’s initial position,
# will the projectile hit the monkey?
# ---------------------------------------------------------------

# -------------------------------
# PARAMETERS / パラメータ設定
# -------------------------------

g = 9.8             # gravitational acceleration [m/s^2] / 重力加速度
v0 = 20.0           # initial speed [m/s] / 発射初速度
a = 10.0            # horizontal distance to monkey [m] / 水平方向距離
b = 5.0             # vertical height of monkey [m] / 鉛直方向高さ

# Compute angle theta so that the initial velocity vector points to (a, b)
# θ = arctan(b / a)
theta = np.arctan2(b, a)

# Velocity components
v0x = v0 * np.cos(theta)
v0y = v0 * np.sin(theta)

# Time to impact (based on straight-line motion)
t_hit = np.sqrt(a**2 + b**2) / v0

# Time samples from 0 to t_hit
t = np.linspace(0, t_hit, 100)

# -------------------------------
# TRAJECTORIES / 軌道計算
# -------------------------------

# P's trajectory: projectile motion
x_p = v0x * t
y_p = v0y * t - 0.5 * g * t**2

# Q's trajectory: free fall from (a, b)
x_q = np.full_like(t, a)
y_q = b - 0.5 * g * t**2

# -------------------------------
# PLOTTING / プロット描画
# -------------------------------

plt.figure(figsize=(8, 6))
plt.plot(x_p, y_p, label='Projectile P', color='blue')
plt.plot(x_q, y_q, label='Falling monkey Q', color='red', linestyle='--')
plt.plot(a, 0, 'ko', label='Impact Point')  # 衝突点
plt.plot(0, 0, 'bo', label='Launch Point')
plt.plot(a, b, 'ro', label='Monkey Start Position')
plt.xlabel('Horizontal distance x [m]')
plt.ylabel('Vertical height y [m]')
plt.title('Monkey Hunting Simulation')
plt.grid(True)
plt.axis('equal')
plt.legend()
plt.tight_layout()
plt.show()

# -------------------------------
# RESULT / 結果表示
# -------------------------------

# Compute vertical difference at x = a
# t_hit is the time both reach same horizontal point
y_p_at_hit = v0y * t_hit - 0.5 * g * t_hit**2
y_q_at_hit = b - 0.5 * g * t_hit**2
diff = abs(y_p_at_hit - y_q_at_hit)

print("---- Simulation Result ----")
print(f"Time to impact: {t_hit:.3f} s")
print(f"Projectile y at impact: {y_p_at_hit:.3f} m")
print(f"Monkey y at impact:     {y_q_at_hit:.3f} m")
print(f"Vertical difference:     {diff:.5f} m")

if diff < 0.01:
    print("✅ HIT: The projectile hits the falling monkey.")
else:
    print("❌ MISS: The projectile misses the monkey.")

問題2：放物運動

import numpy as np
import pandas as pd

# Program Name: Projectile Motion Calculation Based on Meimon no Mori p.12
print("Program Name: Projectile Motion Calculation Based on Meimon no Mori p.12\n")

# 定数設定
g = 9.8               # 重力加速度 [m/s^2]
e = 0.8               # 反発係数
theta_deg = 45        # 発射角度 [度]
theta = np.radians(theta_deg)  # ラジアンに変換
v0 = 10.0             # 初速度 [m/s]
l = 5.0               # 壁までの距離 [m]

# 計算
# 最高点の高さ h1
h1 = (v0**2) * (np.sin(theta)**2) / (2 * g)

# 反射後の最高点の高さ h2
h2 = (e * v0 * np.sin(theta))**2 / (2 * g)

# AB間の時間 t2
t1 = l / (v0 * np.cos(theta))
t2 = 2 * v0 * np.sin(theta) / g

# BC = e * v0 * cosθ * (t2 - t1)
vx = v0 * np.cos(theta)
BC = e * vx * (t2 - t1)

# OB = e^2 * v0^2 * sin(2θ) / g
OB = (e**2) * v0**2 * np.sin(2 * theta) / g

# 合計距離 l = OB + BC
l_total = OB + BC

# 結果のデータフレーム作成
results = pd.DataFrame({
    "項目": [
        "初速度 v0 [m/s]",
        "発射角度 θ [deg]",
        "最高点の高さ h1 [m]",
        "反射後の高さ h2 [m]",
        "AB間の時間差 t2 - t1 [s]",
        "水平方向距離 BC [m]",
        "着地点 OB [m]",
        "合計距離 l = OB + BC [m]"
    ],
    "値": [
        v0,
        theta_deg,
        round(h1, 3),
        round(h2, 3),
        round(t2 - t1, 3),
        round(BC, 3),
        round(OB, 3),
        round(l_total, 3)
    ]
})

# Google Colabで表示
results

問題3：放物運動

# ============================================================
# Program: reflection_projectile_theoretical_vs_simulation.py
# Purpose: 名門の森に基づく「斜方投射＋多段反射」運動の可視化と理論比較
# ============================================================

import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# PARAMETERS / パラメータ設定
# -----------------------------
g = 9.8                   # gravity [m/s^2]
h = 2.0                   # initial height [m]
e = 0.8                   # coefficient of restitution
theta_deg = 45            # launch angle in degrees
theta = np.radians(theta_deg)

# -----------------------------
# Step 1: 初速度 v0 = sqrt(2gh)
# -----------------------------
v0 = np.sqrt(2 * g * h)

# -----------------------------
# Step 2: 速度成分
# -----------------------------
v0x = v0 * np.cos(theta)
v0y = v0 * np.sin(theta)

# -----------------------------
# Step 3: 時間と高さを各バウンドについて計算
# -----------------------------
n_bounce = 3  # バウンド回数
colors = ['blue', 'green', 'orange', 'purple']
x_all = []
y_all = []

x_last = 0
for i in range(n_bounce):
    # 各バウンドの速度
    vx_i = (e ** i) * v0x
    vy_i = (e ** i) * v0y
    t_flight = 2 * vy_i / g  # 滞空時間

    # 時間分割
    t = np.linspace(0, t_flight, 100)
    x = x_last + vx_i * t
    y = vy_i * t - 0.5 * g * t**2

    x_all.append(x)
    y_all.append(y)
    x_last = x[-1]  # 次の初期位置に

# -----------------------------
# Step 4: 理論値の計算
# -----------------------------
# T_inf = 2 * e * v0 * sinθ / (g * (1 - e))
T_inf = (2 * e * v0 * np.sin(theta)) / (g * (1 - e))

# OB_theory = (4 * e * h * sinθ) / (1 - e^2)
OB_theory = (4 * e * h * np.sin(theta)) / (1 - e**2)

# -----------------------------
# Step 5: プロット描画
# -----------------------------
plt.figure(figsize=(10, 5))
for i in range(n_bounce):
    plt.plot(x_all[i], y_all[i], label=f'Trajectory {i+1}', color=colors[i])

plt.axhline(0, color='gray', linestyle='--')
plt.title("Multiple Reflection Projectile ")
plt.xlabel("Horizontal distance x [m]")
plt.ylabel("Vertical height y [m]")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# -----------------------------
# Step 6: 結果出力
# -----------------------------
print("===== Simulation Summary =====")
print(f"Initial speed v₀           = {v0:.3f} m/s")
print(f"Total flight time (T_inf)  = {T_inf:.3f} s")
print(f"Theoretical OB distance     = {OB_theory:.3f} m")
print(f"Simulated x_final           = {x_last:.3f} m")
print("===============================")

問題4：剛体の釣り合い


# ============================================================
# Program: moment_and_force_balance.py
# Purpose: Calculation and Visualization of Moment and Force Balance for a Rigid Body on an Inclined Plane Based on Meimon no Mori
# ============================================================

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# -----------------------------
# Parameters / パラメータ設定
# -----------------------------
g = 9.8           # gravity acceleration [m/s^2]
M = 1.0           # mass of the rod [kg], assumed for calculation
l = 1.0           # length unit [m], assumed for scaling
theta_deg = 53.13 # angle θ in degrees (approx. for sinθ=4/5, cosθ=3/5)
theta = np.radians(theta_deg)

# Trigonometric values from triangle ABC (3:4:5)
cos_theta = 3/5
sin_theta = 4/5

# -----------------------------
# Step 1: Moment about point B
# -----------------------------
# Mg * 5l * cosθ = R * 10l * sinθ
R = (M * g * 5 * l * cos_theta) / (10 * l * sin_theta)
R = (3/8) * M * g  # Simplified as per derivation

# -----------------------------
# Step 2: Vertical force balance
# -----------------------------
N = M * g

# -----------------------------
# Step 3: Horizontal force balance
# -----------------------------
F = R

# -----------------------------
# Step 4: Net reaction force at point B
# -----------------------------
net_force = np.sqrt(N**2 + F**2)
net_force_simplified = (np.sqrt(73) / 8) * M * g

# -----------------------------
# Step 5: Alternative moment about point A (verification)
# -----------------------------
# N * 10l * cosθ = Mg * 5l * cosθ + F * 10l * sinθ
moment_A_left = N * 10 * l * cos_theta
moment_A_right = M * g * 5 * l * cos_theta + F * 10 * l * sin_theta
moment_A_check = np.isclose(moment_A_left, moment_A_right)

# -----------------------------
# Step 6: DataFrame for results
# -----------------------------
results = pd.DataFrame({
    "項目 / Item": [
        "支持反力 R [N] / Support force R",
        "垂直抗力 N [N] / Normal force N",
        "静止摩擦力 F [N] / Friction force F",
        "合力の大きさ [N] / Net reaction force",
        "A点まわりのモーメント平衡 / Moment balance at A"
    ],
    "値 / Value": [
        round(R, 3),
        round(N, 3),
        round(F, 3),
        round(net_force, 3),
        "Satisfied" if moment_A_check else "Not satisfied"
    ]
})

# -----------------------------
# Step 7: Vector visualization with adjusted scale
# -----------------------------
plt.figure(figsize=(10, 6))

# Origin at point B (0,0)
origin = np.array([0, 0])

# Vectors for forces at point B
N_vector = np.array([0, N])  # Normal force (vertical)
F_vector = np.array([F, 0])  # Friction force (horizontal)
R_vector = np.array([-R * sin_theta, R * cos_theta])  # Support force at angle
net_force_vector = N_vector + F_vector  # Net reaction force at B

# Plot vectors with adjusted scale
scale = 20.0  # Adjusted scale (larger value = shorter arrows)
plt.quiver(*origin, *N_vector, color='blue', scale=scale, label='Normal force N', width=0.005, alpha=0.8)
plt.quiver(*origin, *F_vector, color='red', scale=scale, label='Friction force F', width=0.005, alpha=0.8)
plt.quiver(*origin, *R_vector, color='green', scale=scale, label='Support force R', width=0.005, alpha=0.8)
plt.quiver(*origin, *net_force_vector, color='purple', scale=scale, label='Net reaction force', width=0.005, alpha=0.8)

# Set axis limits based on maximum force magnitude
max_force = max(N, F, net_force, np.linalg.norm(R_vector)) * 1.2  # 20% margin
plt.xlim(-max_force, max_force)
plt.ylim(-max_force, max_force)

# Plot settings
plt.title("Force Vectors at Point B")
plt.xlabel("Horizontal Force [N]")
plt.ylabel("Vertical Force [N]")
plt.grid(True)
plt.legend()
plt.axis('equal')  # Equal scaling for x and y axes
plt.tight_layout()

# Show plot
plt.show()

# -----------------------------
# Step 8: Display results
# -----------------------------
print("===== Calculation Summary =====")
print(results.to_string(index=False))
print("===============================")

問題5：剛体の釣り合い

# ============================================================
# Program: sliding_vs_tipping_analysis.py
# Purpose: Sliding vs. Tipping condition analysis for rigid body on incline
# 名門の森の式に基づき、滑り条件・転倒条件を数式で比較
# ============================================================

import numpy as np
import pandas as pd

# -----------------------------
# Parameters / パラメータ設定
# -----------------------------
g = 9.8           # gravity [m/s^2]
M = 1.0           # mass [kg]
a = 0.5           # center of mass to base (horizontal) [m]
b = 1.0           # center of mass to base (vertical) [m]
mu = 0.6          # friction coefficient
theta_deg = 30    # incline angle [deg]
theta = np.radians(theta_deg)  # convert to radians

# -----------------------------
# Step 1: Sliding condition (摩擦による滑り)
# F1 = Mg (μ cosθ - sinθ)
F1 = M * g * (mu * np.cos(theta) - np.sin(theta))

# -----------------------------
# Step 2: Tipping condition (モーメントによる転倒)
# F2 = (Mg / 2b) * (a cosθ - b sinθ)
F2 = (M * g / (2 * b)) * (a * np.cos(theta) - b * np.sin(theta))

# -----------------------------
# Step 3: Comparison
# → F1 > F2 → tipping occurs first
# → F2 > F1 → sliding occurs first
if F1 > F2:
    mode = "Tipping occurs first (転倒が先)"
elif F2 > F1:
    mode = "Sliding occurs first (滑りが先)"
else:
    mode = "Simultaneous (同時に発生)"

# -----------------------------
# Output as table
# -----------------------------
result_data = {
    "Incline θ (deg)": theta_deg,
    "F1 (Sliding condition)": round(F1, 3),
    "F2 (Tipping condition)": round(F2, 3),
    "Result": mode
}

df = pd.DataFrame([result_data])
print("=== Sliding vs. Tipping Analysis ===")
print(df.to_string(index=False))

問題7：運動方程式

import math

def motion_times_with_quadratic(h, g=9.8):
    """
    Computes t1 and t2 using the quadratic formula for vertical motion,
    then returns total time t = t1 + t2.
    """
    # Time to reach height h
    t1 = math.sqrt(2 * h / g)

    # Velocity at top
    v0 = math.sqrt(2 * g * h)

    # Quadratic formula:
    # (1/2)g * t² - v₀ * t - h = 0
    a = 0.5 * g
    b = -v0
    c = -h

    discriminant = b**2 - 4*a*c
    if discriminant < 0:
        raise ValueError("No real solution for t2 (check your inputs)")

    sqrt_discriminant = math.sqrt(discriminant)

    # Only positive root is physically meaningful
    t2 = (-b + sqrt_discriminant) / (2 * a)

    total_time = t1 + t2
    return t1, t2, total_time

# --- Example ---
h = 2.0  # Height in meters
t1, t2, total = motion_times_with_quadratic(h)

print(f"Height h = {h} m")
print(f"t1 (rise time): {t1:.4f} s")
print(f"t2 (fall time): {t2:.4f} s")
print(f"Total time t1 + t2: {total:.4f} s")

問題13：運動量保存則

import numpy as np

# ==== パラメータ定義 / Parameter definitions ====
m = 1.0     # 質点の質量 / Mass of particle
M = 3.0     # リングの質量 / Mass of ring
e = 0.8     # 反発係数 / Coefficient of restitution
v0 = 1.0    # 初期速度（質点） / Initial velocity of mass m
V0 = 0.0    # 初期速度（リング） / Initial velocity of M

# ==== 行列の構築 / Matrix setup ====
# 行列Aは速度の更新を表す2x2行列
A = np.array([
    [(m - e*M) / (m + M), (1 + e)*M / (m + M)],
    [(1 + e)*m / (m + M), (M - e*m) / (m + M)]
])

# 初期状態ベクトル v = [v; V]
v = np.array([v0, V0])

# ==== シミュレーション回数と記録 / Iteration ====
N = 10  # 衝突回数 / Number of collisions
print("n\tv_n\t\tV_n")
for n in range(N + 1):
    print(f"{n}\t{v[0]:.5f}\t{v[1]:.5f}")
    v = A @ v  # 行列による速度の更新

【説明】

質点とリングの運動量保存＋反発係数の2式を使って、次の状態への線形変換として行列 $A$ を定義。
各衝突ごとに $v_n$, $V_n$ を更新して出力。
$n \to \infty$ で両者の速度が等しくなり、一体運動になる（$u = \frac{m}{m+M}v_0$ に収束）。

問題14：運動量保存則

# Program Name: projectile_vector_analysis.py
# Creation Date: 20250722
# Overview: Analyze projectile motion vector and 2D impact coordinates
# Usage: Run to compute vector difference and x/y positions

import numpy as np

# ==== パラメータ定義 / Parameters ====
v0 = 10.0       # 初速度 [m/s] / Initial speed
theta_deg = 45  # 発射角 [度] / Launch angle in degrees
g = 9.8         # 重力加速度 [m/s^2]

# ==== 単位変換 / Unit conversion ====
theta = np.deg2rad(theta_deg)

# ==== 相対速度の大きさ / Relative velocity magnitude ====
# v1 = (v0*cosθ, 0)
# v2 = (-v0*cosθ, 2*v0*cosθ)
ux = -2 * v0 * np.cos(theta)
uy = 2 * v0 * np.cos(theta)
u = np.sqrt(ux**2 + uy**2)

print(f"Relative speed u = {u:.3f} m/s")

# ==== 位置の式 / Position equations ====
x = -(v0**2 / (2 * g)) * np.sin(2 * theta)
y = (3 * v0**2 / (2 * g)) * np.sin(2 * theta)

print(f"x = {x:.3f} m")
print(f"y = {y:.3f} m")

問題15：保存則

# Program Name: momentum_energy_quadratic_solution.py
# Creation Date: 20250722
# Overview: Solve quadratic equation from momentum and energy conservation for elastic collision
# Usage: Run to compute post-collision velocity using two-body elastic collision formula

import sympy as sp

# ==== 記号定義 / Symbol definitions ====
m, M, v0, v = sp.symbols('m M v0 v')

# ==== 運動量保存とエネルギー保存の式から導出される2次方程式 / Quadratic from momentum and energy ====
# (m + M) * v**2 - 2 * m * v0 * v + (m - M) * v0**2 = 0
eq = (m + M) * v**2 - 2 * m * v0 * v + (m - M) * v0**2

# ==== 解く / Solve quadratic ====
sol = sp.solve(eq, v)
print("Solutions for v:")
for s in sol:
    print(sp.simplify(s))

# ==== 速度（向きあり）と速さ（絶対値） / Velocity and speed ====
v_expr = (m - M) / (m + M) * v0
v_abs_expr = abs(m - M) / (m + M) * v0

print("\nVelocity v =")
sp.pprint(sp.simplify(v_expr))

print("\nSpeed |v| =")
sp.pprint(sp.simplify(v_abs_expr))

問題17：保存則

# Program Name: elastic_collision_matrix_form.py
# Creation Date: 20250722
# Overview: Solve 2-body elastic collision using matrix formulation from momentum and restitution
# Usage: Run to compute post-collision velocities using matrix method

import sympy as sp

# ==== 記号定義 / Symbol definitions ====
m, M, g, h = sp.symbols('m M g h')
v, V, v0 = sp.symbols('v V v0')

# ==== エネルギー保存から初速度 v0 を導出 / Initial velocity from potential energy ====
v0_expr = sp.sqrt(2 * g * h)

# ==== 運動量保存則 / Momentum conservation ====
eq1 = sp.Eq(m * v + M * V, m * v0)

# ==== 反発係数 e = 1 の関係式 / Restitution relation ====
eq2 = sp.Eq(v - V, -(v0))

# ==== 行列形式で解く / Solve as matrix ====
A = sp.Matrix([[m, M], [1, -1]])
b = sp.Matrix([m * v0, -v0])
sol = A.LUsolve(b)

v_sol = sol[0].subs(v0, v0_expr)
V_sol = sol[1].subs(v0, v0_expr)

print("Post-collision velocity of small mass v:")
sp.pprint(sp.simplify(v_sol))

print("\nPost-collision velocity of large mass V:")
sp.pprint(sp.simplify(V_sol))

問題18：保存則

# Program Name: pendulum_cart_conservation.py
# Creation Date: 20250722
# Overview: Solve pendulum inside a box using momentum and energy conservation
# Usage: Run to compute velocity v and V at lowest point of swing

import sympy as sp

# ==== 記号定義 / Symbol definitions ====
m, M, g, l, theta = sp.symbols('m M g l theta', positive=True)
v, V = sp.symbols('v V')  # それぞれの速度 / Velocities

# ==== 力学的エネルギー保存 / Energy conservation ====
# m*g*(l - l*cosθ) = 1/2*m*v² + 1/2*M*V²
energy_eq = sp.Eq(m * g * (l - l * sp.cos(theta)), (1/2) * m * v**2 + (1/2) * M * V**2)

# ==== 運動量保存 / Momentum conservation ====
# 水平方向の全運動量は保存 / m*v = M*V
momentum_eq = sp.Eq(m * v, M * V)

# ==== V を消去して v を解く / Eliminate V to solve for v ====
V_expr = sp.solve(momentum_eq, V)[0]
energy_sub = energy_eq.subs(V, V_expr)
v_sol = sp.simplify(sp.solve(energy_sub, v)[0])

# ==== V も求める / Back-substitute for V ====
V_sol = V_expr.subs(v, v_sol)

print("v (mass P at bottom):")
sp.pprint(v_sol)

print("\nV (box speed at same moment):")
sp.pprint(V_sol)

問題27：等速円運動

import math

# パラメータ設定（例）
m = 1.0       # 質量 [kg]
g = 9.81      # 重力加速度 [m/s^2]
r = 2.0       # 曲率半径 [m]
h = 0.5       # 重心の高さ [m]
d = 0.3       # 幅 [m]
mu = 0.4      # 摩擦係数

# 滑り出す速度の境界 v2
v2 = math.sqrt(g * d * r / h)

# 滑り出す速度 v1 の条件からの不等式チェック
v1_limit = math.sqrt(mu * g * r)

# μ < d/h の条件チェック
condition = mu < d / h

print(f"v2 (滑り出す速度の境界) = {v2:.3f} m/s")
print(f"v1の上限 = {v1_limit:.3f} m/s")
print(f"μ < d/h の条件は {condition}")

# v1 < v2 であれば滑り出さない
if v1_limit < v2:
    print("車輪は滑り出さない条件を満たしています。")
else:
    print("車輪は滑り出す可能性があります。")

import math

def banking_angle(v, r, g=9.8):
    """
    Calculate the banking angle θ [radians, degrees] using:
    tan(θ) = v² / (g * r)
    """
    tan_theta = v**2 / (g * r)
    theta_rad = math.atan(tan_theta)         # radians
    theta_deg = math.degrees(theta_rad)      # degrees
    return theta_rad, theta_deg

# Example usage
v = 20.0  # velocity in m/s
r = 50.0  # radius in m
theta_rad, theta_deg = banking_angle(v, r)

print(f"θ (radians) = {theta_rad:.4f}")
print(f"θ (degrees) = {theta_deg:.2f}")

問題30：円運動

import math

def motion_on_incline(u, r, g=9.8, alpha_deg=30):
    """
    Calculate time t and horizontal distance x for a particle on an incline
    after the string is cut.
    
    Parameters:
    u          : initial speed along incline [m/s]
    r          : radius of circular motion before release [m]
    g          : gravitational acceleration [m/s^2]
    alpha_deg  : incline angle in degrees
    
    Returns:
    t : time to reach point A [s]
    x : horizontal displacement [m]
    """
    alpha = math.radians(alpha_deg)
    
    # Compute discriminant inside sqrt
    term_under_sqrt = u**2 + 4 * math.sqrt(3) * g * r * math.sin(alpha)
    
    # Time of flight (positive root of quadratic equation)
    t = (u + math.sqrt(term_under_sqrt)) / (2 * g * math.sin(alpha))
    
    # Horizontal motion
    x = (r / 2) * (math.sqrt(3) / (4 * g * math.sin(alpha))) * (u + math.sqrt(term_under_sqrt))
    
    return t, x

# --- Example usage ---
u = 5.0     # Initial velocity along incline [m/s]
r = 2.0     # Radius [m]

t, x = motion_on_incline(u, r)

print(f"Time to reach point A: {t:.4f} s")
print(f"Horizontal distance x: {x:.4f} m")

問題37：単振動

$$
\omega = \sqrt{\frac{2k}{m}}
$$

$$
u = \frac{v_0}{2} \cos\left(\sqrt{\frac{2k}{m}} t\right)
$$

$$
v_A = \frac{v_0}{2} + u = \frac{v_0}{2}\left(1 + \cos\left(\sqrt{\frac{2k}{m}} t\right)\right)
$$

import numpy as np
import matplotlib.pyplot as plt

def relative_velocity(v0, k, m, t_array):
    """
    Calculates v_A(t) = v0/2 * (1 + cos(sqrt(2k/m) * t))
    
    Parameters:
    v0       : Initial velocity [m/s]
    k        : Spring constant [N/m]
    m        : Mass [kg]
    t_array  : numpy array of time values [s]
    
    Returns:
    v_A      : numpy array of velocity values [m/s]
    """
    omega = np.sqrt(2 * k / m)
    v_A = (v0 / 2) * (1 + np.cos(omega * t_array))
    return v_A

# --- Parameters ---
v0 = 1.0      # Initial velocity [m/s]
k = 2.0       # Spring constant [N/m]
m = 1.0       # Mass [kg]
t = np.linspace(0, 10, 500)  # Time from 0 to 10 seconds

# --- Calculate and Plot ---
vA = relative_velocity(v0, k, m, t)

plt.figure(figsize=(8, 4))
plt.plot(t, vA, label=r'$v_A(t)$')
plt.title("Velocity of A over Time")
plt.xlabel("Time [s]")
plt.ylabel("Velocity [m/s]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題38

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# --- Constants ---
M = 1.0           # Mass [kg]
k = 20.0          # Spring constant [N/m]
omega = 5.0       # Angular velocity [rad/s]

# Effective spring constant: L = k - M * omega^2
L = k - M * omega**2
omega_0 = np.sqrt(L / M)  # Natural frequency of oscillation

# --- Define the differential equation ---
# Let x' = v, v' = -omega_0^2 * x
def harmonic_oscillator(t, y):
    x, v = y
    dxdt = v
    dvdt = -omega_0**2 * x
    return [dxdt, dvdt]

# --- Initial conditions and time array ---
x0 = 0.1   # Initial displacement [m]
v0 = 0.0   # Initial velocity [m/s]
t_span = (0, 10)  # Time interval [s]
t_eval = np.linspace(*t_span, 1000)

# --- Solve the differential equation ---
sol = solve_ivp(harmonic_oscillator, t_span, [x0, v0], t_eval=t_eval)

# --- Plot the results ---
plt.figure(figsize=(10, 4))
plt.plot(sol.t, sol.y[0], label='x(t): Displacement')
plt.plot(sol.t, sol.y[1], label='v(t): Velocity', linestyle='--')
plt.title("Harmonic Oscillation of a Spring in Rotational System")
plt.xlabel("Time [s]")
plt.ylabel("Displacement / Velocity")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# --- Print additional information ---
print(f"Angular frequency ω₀ = {omega_0:.3f} [rad/s]")
print(f"Period T = {2 * np.pi / omega_0:.3f} [s]")

問題39：単振動

(1) 最大摩擦力に達したときの静止摩擦係数の導出

ばねが自然長から変位 $x = d$ に伸びたとき、最大摩擦力に達して静止状態にある：

$$
kd = \mu_0 mg \Rightarrow \mu_0 = \frac{kd}{mg}
$$

(2) 力の式と単振動の中心

物体 M が左に滑るとき、摩擦力は右向き（ばねは左へ引く）。合力：

$$
F = -kx + \mu mg = -k \left( x - \frac{\mu mg}{k} \right)
$$

つまり、振動の中心は

$$
x_c = \frac{\mu mg}{k}
$$

(3) 両端 $x_0, x_1$ の平均が中心

$$
x_c = \frac{x_0 + x_1}{2}
\Rightarrow \mu = \frac{k(x_0 + x_1)}{2mg}
$$

(4) 摩擦熱を考慮したエネルギー保存

$$
\frac{1}{2}kx_0^2 = \frac{1}{2}kx_1^2 + \mu mg(x_0 - x_1)
$$

これを整理すると：

$$
\mu = \frac{k(x_0^2 - x_1^2)}{2mg(x_0 - x_1)} = \frac{k(x_0 + x_1)}{2mg}
$$

(5) 左右に滑るときの振動中心

左に滑る：摩擦力右向き、振動中心は

$$
x_c = \frac{x_0 + x_1}{2} = 0.5d
$$

右に滑る：摩擦力左向き、合力は

$$
F = -kx - \mu mg = -k \left( x + \frac{\mu mg}{k} \right)
\Rightarrow x_c' = -\frac{\mu mg}{k}
$$

(6) 減衰の終わりの位置まで

例：$x_1 = -2.5d$ から振動中心 $x_c' = -0.5d$ を中心に振動 → 振幅 = 2d

右端まで $x_2 = 1.5d$、左へ $x_3 = -0.5d$

そして、赤線部（|x| ≦ d）を超えられなくなると摩擦により静止。

(7) グラフ

横軸：時間 $t / t_1$
縦軸：位置 $x / d$
減衰しながら最終的に停止
振動の中心が周期ごとに移動

各点の時間と位置のまとめ

時刻 $t$	振動中心 $x_c$	始点 $x_{\text{start}}$	終点 $x_{\text{end}}$	振幅	備考
$t = 0$	–	$x_0 = 3.5d$	–	–	初期位置（静止）
$t = t_1$	$x_c = 0.5d$	$x_0 = 3.5d$	$x_1 = -2.5d$	$3d$	左向き滑り
$t = 2t_1$	$x_c' = -0.5d$	$x_1 = -2.5d$	$x_2 = 1.5d$	$2d$	右向き滑り
$t = 3t_1$	$x_c = 0.5d$	$x_2 = 1.5d$	$x_3 = -0.5d$	$1d$	左向き滑り
$t = 4t_1$	静止（赤線範囲）	$x_3 = -0.5d$	–	–	停止

# Program Name: damped_oscillation_generalized.py
# Creation Date: 20250722
# Overview: Plot generalized damped spring-mass motion with friction, all positions/time/labels defined via variables
# Usage: Adjust constants m, k, d, and displacement coefficients. Run to visualize annotated oscillation diagram.

!pip install matplotlib

import matplotlib.pyplot as plt
import numpy as np

# --- System Constants ---
mass = 1.0                     # 質量 [kg]
spring_const = 4.0             # ばね定数 [N/m]
unit_disp = 1.0                # 基準変位 d [m]
gravity = 9.8                  # 重力加速度 [m/s^2]
mu_static = spring_const * unit_disp / (mass * gravity)  # 静止摩擦係数 μ₀

# --- Period & Time ---
T = 2 * np.pi * np.sqrt(mass / spring_const)    # 周期 T
half_T = T / 2                                   # 半周期 t₁

# --- Displacement Coefficients ---
x_coeffs = {
    "x0": 3.5,
    "x1": -2.5,
    "x2": 1.5,
    "x3": -0.5,
    "stop": -0.5
}

# --- Time Coefficients ---
t_coeffs = {
    "t0": 0,
    "t1": 1,
    "t2": 2,
    "t3": 3,
    "t4": 4
}

# --- Data Points List ---
points = [
    {"label": "x0", "t_coeff": t_coeffs["t0"], "x_coeff": x_coeffs["x0"]},
    {"label": "x1", "t_coeff": t_coeffs["t1"], "x_coeff": x_coeffs["x1"]},
    {"label": "x2", "t_coeff": t_coeffs["t2"], "x_coeff": x_coeffs["x2"]},
    {"label": "x3", "t_coeff": t_coeffs["t3"], "x_coeff": x_coeffs["x3"]},
    {"label": "Stop", "t_coeff": t_coeffs["t4"], "x_coeff": x_coeffs["stop"]}
]

# --- Extract Times, Positions, Labels ---
times = [p["t_coeff"] * half_T for p in points]
positions = [p["x_coeff"] * unit_disp for p in points]
labels = [f"${p['label']} = {p['x_coeff']}d$\nt = {p['t_coeff'] * half_T:.2f}s" for p in points]

# --- Center Line Coefficients ---
center_coeffs = [0.5, -0.5]   # ±μmg/k の中心位置係数

# --- Plot ---
plt.figure(figsize=(10, 5))
plt.plot(times, positions, 'o-', color='tab:red', linewidth=2, markersize=8)

# Annotate points
for t, x, label in zip(times, positions, labels):
    plt.text(t, x + 0.25, label, ha='center', fontsize=10)

# Center lines
for c in center_coeffs:
    y = c * unit_disp
    plt.axhline(y, linestyle='--', color='blue', linewidth=1)
    plt.text(times[-1] + 0.3, y, f"center = {c:.1f}d", va='center', color='blue')

# Stopping region boundaries ±d
for lim, name in zip([unit_disp, -unit_disp], ["+d limit", "−d limit"]):
    plt.axhline(lim, linestyle=':', color='green', linewidth=1)
    plt.text(times[-1] + 0.3, lim, name, va='center', color='green')

# Axis
plt.title("Damped Oscillation with Generalized Parameters", fontsize=14)
plt.xlabel("Time [s]", fontsize=12)
plt.ylabel("Position [m]", fontsize=12)
plt.grid(True)
plt.xticks(times)
plt.yticks(np.arange(-3.5 * unit_disp, 4 * unit_disp, unit_disp))
plt.tight_layout()
plt.show()

問題41：単振動

問題の概要

力 $F = -kx$ （ここでは $k = \frac{\mu M g}{l}$）の復元力が働くため、質点（板）は単振動をする。
この場合の運動方程式は

$$
M \frac{d^2 x}{dt^2} = -k x
$$

となり、ここで

$$
k = \frac{\mu M g}{l}
$$

微分方程式

$$
M \frac{d^2 x}{dt^2} + k x = 0
$$

または

$$
\frac{d^2 x}{dt^2} + \frac{k}{M} x = 0
$$

解の形

これは単振動の微分方程式で、一般解は

$$
x(t) = A \cos(\omega t) + B \sin(\omega t)
$$

ここで角振動数 $\omega$ は

$$
\omega = \sqrt{\frac{k}{M}} = \sqrt{\frac{\mu g}{l}}
$$

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
mu = 0.3    # 動摩擦係数の例
M = 1.0     # 質量（kg）
g = 9.81    # 重力加速度（m/s^2）
l = 0.5     # 長さ（m）

k = mu * M * g / l
omega = np.sqrt(k / M)

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

# 振幅d
d = 0.1  # 振幅の例（m）

# 変位x(t)
x = d * np.cos(omega * t)

# プロット
plt.plot(t, x)
plt.xlabel("Time (s)")
plt.ylabel("Displacement (m)")
plt.title("Simple Harmonic Motion")
plt.grid(True)
plt.show()

バネ振り子などに見られる、以下のような単振動を考えます：

$$
x(t) = A \cos(\omega t + \phi)
$$

または、初期位置 $x(0) = x_0$、初期速度 $v(0) = v_0$ を与えると、解は：

$$
x(t) = x_0 \cos(\omega t) + \frac{v_0}{\omega} \sin(\omega t)
$$

import math
import numpy as np
import matplotlib.pyplot as plt

def shm_with_initial_conditions(x0, v0, k, m, t_max=10, num_points=500):
    """
    Simulates simple harmonic motion given initial conditions:
    x(0) = x0, v(0) = v0
    
    Parameters:
    x0         : initial position [m]
    v0         : initial velocity [m/s]
    k          : spring constant [N/m]
    m          : mass [kg]
    t_max      : time range to simulate [s]
    num_points : number of time points
    
    Returns:
    t, x       : time array and displacement array
    """
    omega = math.sqrt(k / m)
    t = np.linspace(0, t_max, num_points)
    
    x = x0 * np.cos(omega * t) + (v0 / omega) * np.sin(omega * t)
    
    return t, x

# --- Parameters ---
x0 = 0.1    # initial position [m]
v0 = 0.0    # initial velocity [m/s]
k = 2.0     # spring constant [N/m]
m = 1.0     # mass [kg]

# --- Run simulation ---
t, x = shm_with_initial_conditions(x0, v0, k, m)

# --- Plot ---
plt.figure(figsize=(8, 4))
plt.plot(t, x, label=f'x0={x0}, v0={v0}')
plt.title('Simple Harmonic Motion (with Initial Conditions)')
plt.xlabel('Time [s]')
plt.ylabel('Displacement [m]')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題43

import numpy as np

# --- 定数 ---
R = 6.4e6        # 地球の半径 [m]
g = 9.8          # 地表での重力加速度 [m/s^2]
h = 1.0e6        # 深さ [m]（例として 1000 km の深さ）

# --- 振動周期の半分（BA間） ---
T_half = np.pi * np.sqrt(R / g)    # [s]
T_min = T_half / 60                # 分に換算

# --- 振幅 A = sqrt(R^2 - h^2) ---
A = np.sqrt(R**2 - h**2)

# --- 最大速度 v_max = A * 2π / T = A * sqrt(g / R) ---
v_max = A * np.sqrt(g / R)

# --- 結果表示 ---
print(f"半周期 T/2 = {T_half:.2f} [s] ≈ {T_min:.2f} [min]")
print(f"振幅 A = {A:.2f} [m]")
print(f"最大速度 v_max = {v_max:.2f} [m/s]")

問題44

# --- 基本定数 ---
m_ice = 200               # 氷の質量 [g]
Lf = 336                  # 融解熱 [J/g]
c_water = 4.2             # 水の比熱 [J/g・K]
delta_T_ice = 50          # 0℃ → 50℃ の温度変化
Q1 = 600 * (124 - 12)     # 熱量1 BC間 [J]
Q2 = 600 * (199 - 124)    # 熱量2 CD間 [J]
Q3 = 600 * (12 - 0)       # 熱量3 AB間 [J]

# --- 容器の比熱計算（CD間） ---
# Q2 = m_ice*c_water*ΔT + C*ΔT → C = ?
C = (Q2 - m_ice * c_water * delta_T_ice) / delta_T_ice
print(f"(1) 容器の熱容量 C = {C:.1f} J/K")

# --- 水の比熱計算（AB間） ---
# Q3 = m*c1*ΔT + C*ΔT → c1 = ?
delta_T_AB = 15
c1 = (Q3 - C * delta_T_AB) / (m_ice * delta_T_AB)
print(f"(2) 水の比熱 c1 = {c1:.1f} J/g・K")

# --- 銅の比熱計算（(3)の熱授受） ---
# 水と容器が放出した熱 = 銅が得た熱
T_hot = 50
T_final = 47.7
T_cold = -10
delta_T_copper = T_final - T_cold

Q_lost = m_ice * c_water * (T_hot - T_final) + C * (T_hot - T_final)
m_copper = 90

c_cu = Q_lost / (m_copper * delta_T_copper)
print(f"(3) 銅の比熱 c_cu = {c_cu:.3f} J/g・K")

問題45

# --- 必要な定数 ---
rho0 = 1.20        # 外気の密度 [kg/m^3]
V = 500            # 気球の体積 [m^3]
W = 180            # 吊り下げ物体の質量 [kg]
T0 = 280           # 外気の温度 [K]
T1 = 400           # 内部の温度 [K]
g = 9.8            # 重力加速度 [m/s^2]
P0 = 1.0e5         # 地上での気圧 [Pa]

# --- (1) 内部の空気密度の算出 ---
rho = (rho0 * V * g - W * g) / (V * g)
print(f"(1) 内部の空気密度 rho = {rho:.3f} kg/m³")

# --- (1b) 温度T1の確認（すでに与えられている） ---
T1_calc = (rho0 / rho) * T0
print(f"(1b) T1の確認: T1 = {T1_calc:.1f} K")

# --- (2) 浮力条件からρ₁（新しい密度）を計算 ---
W_inner = 18     # 気球の重量 [kg]
w_inner = 180    # 中身の吊り下げ物 [kg]

rho1 = (T1 * (W_inner - w_inner)) / (V * (T1 - T0))
print(f"(2) 新しい内部密度 rho₁ = {rho1:.2f} kg/m³")

# --- (3) 内部の気圧 P1 を計算 ---
P1 = (rho1 / rho0) * P0
print(f"(3) 内部気圧 P1 = {P1:.2e} Pa")

# --- (4) 上昇高度 h を求める ---
h = 2 * (P0 - P1) / ((rho0 + rho1) * g)
print(f"(4) 上昇高度 h ≒ {round(h)} m")

問題46：分子運動

(1) 分子が壁と衝突する際の運動量の変化

分子の速度ベクトル $v$ のうち、壁面に垂直な成分：$v \cos \theta$
衝突によって反転するので、運動量の変化は：

$$
\Delta p = mv \cos \theta - (-mv \cos \theta) = 2mv \cos \theta
$$

(2) 衝突回数（時間あたり）

分子は $vt$ の距離を進む
壁との往復距離 $2r \cos \theta$ ごとに1回衝突
よって、時間 $t$ あたりの衝突回数：

$$
\frac{vt}{2r \cos \theta}
$$

(3) 分子1個あたりの力積と全体の力

単位時間に壁に与える力積（1個）：

$$
\Delta p \cdot \text{衝突回数} = 2mv \cos \theta \cdot \frac{vt}{2r \cos \theta} = \frac{mv^2 t}{r}
$$

分子 $N$ 個全体の力積：

$$
Ft = N \cdot \frac{mv^2 t}{r} \Rightarrow F = \frac{Nmv^2}{r}
$$

(4) 圧力 $P$ の導出

球の表面積 $A = 4\pi r^2$
圧力は $P = \frac{F}{A}$ より：

$$
P = \frac{F}{4\pi r^2} = \frac{Nmv^2}{4\pi r^3}
$$

球の体積 $V = \frac{4}{3}\pi r^3$ を使うと：

$$
P = \frac{Nmv^2}{3V}
$$

(5) 理想気体の状態方程式と温度の関係

理想気体の状態方程式：

$$
PV = nRT = \frac{N}{N_A} RT
$$

圧力の式より：

$$
PV = \frac{1}{3} Nmv^2
$$

よって：

$$
\frac{1}{2} mv^2 = \frac{3}{2} \cdot \frac{R}{N_A} T
$$

(6) 内部エネルギー $U$

1分子あたりの平均運動エネルギー：

$$
\frac{1}{2} mv^2 = \frac{3}{2} k_B T \quad (k_B = \frac{R}{N_A})
$$

全体の内部エネルギー：

$$
U = N \cdot \frac{1}{2} mv^2 = \frac{3}{2} nRT
$$

# Program Name: kinetic_theory_pressure_energy.py
# Creation Date: 20250722
# Overview: Calculates pressure and internal energy of an ideal gas based on kinetic theory using input constants.
# Usage: Set constants N, m, v, r, T, n, then run to compute P and U with print output.

!pip install numpy

import numpy as np

# --- Constants ---
N = 6.0e23              # Number of molecules
m = 4.65e-26            # Mass of one molecule [kg] (e.g., Nitrogen N2)
v = 500.0               # Average speed of molecule [m/s]
r = 0.1                 # Radius of spherical container [m]
T = 300.0               # Temperature [K]
n = 1.0                 # Amount of substance [mol]
R = 8.314               # Gas constant [J/(mol·K)]
NA = 6.022e23           # Avogadro constant [1/mol]
pi = np.pi

# --- Derived Constants ---
V = (4/3) * pi * r**3             # Volume of the container [m^3]
A = 4 * pi * r**2                 # Surface area of the container [m^2]
kB = R / NA                       # Boltzmann constant [J/K]

# --- (1) Momentum change per collision ---
delta_p = 2 * m * v              # Momentum change [kg·m/s]

# --- (2) Collision count per time ---
collision_count = v / (2 * r)    # Collisions per second

# --- (3) Impulse per second (Force per particle) ---
Ft_one = delta_p * collision_count
F_total = N * Ft_one

# --- (4) Pressure ---
P = F_total / A
P_via_volume = N * m * v**2 / (3 * V)

# --- (5) Average kinetic energy per particle ---
KE_avg = (3 / 2) * kB * T

# --- (6) Internal energy of the gas ---
U = N * (1/2) * m * v**2
U_from_RT = (3 / 2) * n * R * T

# --- Print Results ---
print(f"Container Volume V       = {V:.4e} m^3")
print(f"Container Area A         = {A:.4e} m^2")
print(f"Momentum change Δp       = {delta_p:.4e} kg·m/s")
print(f"Collision frequency      = {collision_count:.2f} times/s")
print(f"Force on one particle Ft = {Ft_one:.4e} N·s/s")
print(f"Total Force F            = {F_total:.4e} N")
print(f"Pressure P (via force)   = {P:.4e} Pa")
print(f"Pressure P (via volume)  = {P_via_volume:.4e} Pa")
print(f"Kinetic Energy per mol.  = {KE_avg:.4e} J")
print(f"Internal Energy U        = {U:.4e} J")
print(f"Internal Energy U (nRT)  = {U_from_RT:.4e} J")

問題47：分子運動

(1) 分子がピストンに衝突する際の速度変化

ピストン速度：$u$
分子衝突前速度：$v_x$
弾性衝突 ⇒ 相対速度反転 ⇒

$$
v_x' = - (v_x - u) + u = v_x - 2u
$$

(2) 運動エネルギーの減少

$$
\Delta E = \frac{1}{2}mv_x^2 - \frac{1}{2}m(v_x - 2u)^2
= 2m v_x u - 2m u^2
= 2m v_x u \left(1 - \frac{u}{v_x} \right) \approx 2m v_x u
$$

（$u \ll v_x$ の近似）

(3) 単位時間あたりの衝突回数

距離 $2L$ ごとに衝突 ⇒ 時間 $\Delta t$ における回数：

$$
\frac{v_x \Delta t}{2L}
$$

(4) 単位時間あたりのエネルギー減少量

1回あたりの減少 × 衝突回数：

$$
\Delta e = 2m v_x u \cdot \frac{v_x \Delta t}{2L}
= \frac{m v_x^2 u}{L} \Delta t
$$

(5) 容器の体積変化 $\Delta V$

断面積 $S$, 高さ $L$ より体積 $V = SL$

$$
\Delta V = S \cdot u \cdot \Delta t \Rightarrow \frac{\Delta V}{V} = \frac{u}{L} \Delta t
$$

(6) 多数の分子の平均化（全体のエネルギー変化）

$$
\Delta E = - \Delta e = - m \bar{v}_x^2 \cdot \frac{u}{L} \Delta t
= - m \bar{v}_x^2 \cdot \frac{\Delta V}{V}
$$

(7) 三平方より：

$$
v^2 = v_x^2 + v_y^2 + v_z^2 \Rightarrow \bar{v}_x^2 = \frac{1}{3} \bar{v}^2
$$

(8) 内部エネルギーの変化

$$
\Delta E = - \frac{1}{3} m \bar{v}^2 \cdot \frac{\Delta V}{V}
$$

また、内部エネルギー $E = \frac{1}{2} m \bar{v}^2$ より：

$$
\Delta E = - \frac{2}{3} E \cdot \frac{\Delta V}{V} \tag{①}
$$

(9) エネルギーと温度の関係

比例関係：$E = aT$

$$
\Delta E = a \Delta T \Rightarrow \frac{\Delta E}{E} = \frac{\Delta T}{T} \tag{②}
$$

①と②より：

$$
\frac{\Delta T}{T} = - \frac{2}{3} \cdot \frac{\Delta V}{V} \Rightarrow \Delta T = - \frac{2T}{3V} \Delta V \tag{③}
$$

微分方程式と断熱変化の式導出（続き）

(10) 微分形式での温度変化：

先ほどの温度変化式：

$$
\Delta T = -\frac{2T}{3V} \Delta V
$$

を微小変化で書くと：

$$
\frac{dT}{dV} = -\frac{2T}{3V}
$$

(11) 両辺を積分：

両辺を変数分離して積分：

$$
\int \frac{1}{T} dT = -\frac{2}{3} \int \frac{1}{V} dV
$$

よって：

$$
\ln T = -\frac{2}{3} \ln V + \ln C
\Rightarrow \ln (T V^{2/3}) = \ln C
$$

両辺を指数化：

$$
T V^{2/3} = \text{const.}
$$

(12) よって、断熱変化の式：

$$
T V^{\gamma - 1} = \text{const.}
\quad \text{ただし } \gamma = \frac{C_P}{C_V}
$$

単原子分子の場合：

$$
\gamma = \frac{5}{3} \Rightarrow T V^{2/3} = \text{const.}
$$

# Program Name: adiabatic_cooling_due_to_piston_collision.py
# Creation Date: 20250722
# Overview: Simulate energy loss and temperature drop due to molecular collisions with a moving piston
# Usage: Set initial values for m, vx, u, L, S, T. Run the script to calculate ΔE, ΔT based on piston motion.

import math

# --- Constants (all variables defined explicitly) ---
m = 4.65e-26          # mass of a molecule (kg) - approx. mass of one nitrogen molecule
vx = 500              # average velocity of molecule in x-direction (m/s)
u = 0.01              # piston velocity (m/s)
L = 0.1               # container length (m)
S = 0.01              # container cross-sectional area (m^2)
T = 300               # initial temperature (K)
dt = 0.01             # small time interval (s)

# --- Derived values ---
V = S * L                         # container volume [m^3]
delta_V = S * u * dt             # change in volume
delta_V_over_V = delta_V / V     # relative volume change

# average v_x^2 from 3D thermal motion: vx² = (1/3) * v²
v2_avg = 3 * vx ** 2             # mean square speed
vx2_avg = v2_avg / 3             # ⟨vx²⟩

# (6) Energy decrease
delta_E = -m * vx2_avg * delta_V_over_V

# (8) Total energy
E = 0.5 * m * v2_avg

# (9) Temperature drop from ΔE/E = ΔT/T
delta_T = - (2 * T / 3) * (delta_V_over_V)


# --- Output ---
print(f"Initial Volume V: {V:.4e} m³")
print(f"Volume Change ΔV: {delta_V:.4e} m³")
print(f"Relative Volume Change ΔV/V: {delta_V_over_V:.4e}")
print(f"Energy Decrease ΔE: {delta_E:.4e} J")
print(f"Initial Internal Energy E: {E:.4e} J")
print(f"Temperature Change ΔT: {delta_T:.4f} K")
print(f"Final Temperature: {T + delta_T:.2f} K")

# Program Name: adiabatic_relation.py
# Creation Date: 20250722
# Overview: Derive and confirm the adiabatic relation T * V^(gamma - 1) = constant
# Usage: Define initial T and V, then compute and confirm constancy of T * V^(gamma - 1)

import numpy as np

# --- Constants ---
Cv = 3/2  # Heat capacity at constant volume for monoatomic gas [R units]
Cp = 5/2  # Heat capacity at constant pressure [R units]
gamma = Cp / Cv  # Adiabatic index γ = Cp / Cv

# Initial state
T0 = 300.0  # Initial temperature [K]
V0 = 1.0    # Initial volume [m^3]

# --- Function to compute T from V using adiabatic relation ---
def temperature(V, T_ref=T0, V_ref=V0, gamma=gamma):
    """Compute T from V using T*V^(γ-1)=const"""
    return T_ref * (V_ref / V)**(gamma - 1)

# --- Generate volume range and compute temperature ---
V_values = np.linspace(0.5 * V0, 2.0 * V0, 100)
T_values = temperature(V_values)

# --- Confirm constancy of T * V^(γ - 1) ---
TV_gamma_minus1 = T_values * V_values**(gamma - 1)

# --- Print result sample ---
print("Adiabatic index gamma =", gamma)
print("Sample results (T, V, T*V^(γ-1)):")
for i in range(0, len(V_values), 20):
    print(f"V = {V_values[i]:.3f} m³, T = {T_values[i]:.2f} K, T*V^(γ-1) = {TV_gamma_minus1[i]:.2f}")

# --- Optional plot ---
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5))
plt.plot(V_values, T_values, label="T(V)", color='red')
plt.title("Adiabatic Expansion/Compression: Temperature vs Volume")
plt.xlabel("Volume [m³]")
plt.ylabel("Temperature [K]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題49：熱力学

# 定数・パラメータの設定（仮の値）
P0 = 100    # 圧力 [kPa]
S = 0.01    # ピストン断面積 [m^2]
L = 0.1     # ピストン初期高さ [m]
n = 1       # モル数
R = 8.314   # 気体定数 [J/(mol·K)]
T0 = 300    # 初期温度 [K]
M = 0.5     # 重りの質量 [kg]
g = 9.81    # 重力加速度 [m/s^2]

# 温度T1の計算
T1 = T0 + (M * g * L) / (n * R)

# 温度T2の計算
T2 = (3/2) * T1

# W2の計算
W2 = 0.5 * n * R * T1

# Q2の計算
Q2 = (5/4) * (n * R * T0 + M * g * L)

# T3の計算（例）
factor = (3/2) ** (5/3)
T3 = factor * T1

# W3の計算
W3 = (3/2) ** 2 * ((3/2) ** (2/3) - 1) * (n * R * T0 + M * g * L)

print(f"T1 = {T1:.2f} K")
print(f"T2 = {T2:.2f} K")
print(f"W2 = {W2:.2f} J")
print(f"Q2 = {Q2:.2f} J")
print(f"T3 = {T3:.2f} K")
print(f"W3 = {W3:.2f} J")

問題50

# --- 定数定義 ---
R = 8.314    # [J/(mol·K)] 気体定数
T0 = 300.0   # [K] 初期温度
Q = (5/12) * R * T0   # 与えられた熱量 Q（オの3分の1の時間で）

# --- オ：定圧変化での T'（7/6 T0） ---
T1 = (7/6) * T0
print(f"(1) T' = T1 = {T1:.2f} K")

# 状態方程式より l = (7/9) L
L = 1.0   # シリンダーの全長（仮）[m]
l = (7/9) * L
print(f"(2) ピストンの位置 l = {l:.3f} m")

# --- カ：ピストンが底まで到達時の温度 T''（4/3 T0） ---
T2 = (4/3) * T0
print(f"(3) T'' = T2 = {T2:.2f} K")

# 加えた熱量 q1（定圧変化）
q1 = (5/2) * R * (T2 - T1)
print(f"(4) q1 = {q1:.2f} J/mol")

# 残りの熱量 q2（定積変化）
q2 = Q * (2/3) - q1
print(f"(5) q2 = {q2:.2f} J/mol")

# q2 = nCvΔT より T3 を求める
# q2 = (3/2) * R * (T3 - T2)
T3 = q2 / ((3/2) * R) + T2
print(f"(6) 最終温度 T3 = {T3:.2f} K = {(29/18)*T0:.2f} K (検算)")

問題52：熱力学

# Program Name: tv_curve_max_gradient.py
# Creation Date: 20250722
# Overview: Find maximum temperature of T(V) using gradient (derivative) method
# Usage: Run the script to compute and print volume V where T is maximum

import numpy as np

# ==== パラメータ定義 / Parameter definitions ====
P0 = 1.0       # 初期圧力 [Pa] / Initial pressure
V0 = 1.0       # 初期体積 [m^3] / Initial volume
n = 1.0        # モル数 / Amount of substance
R = 8.314      # 気体定数 [J/(mol·K)] / Ideal gas constant

# ==== T(V)の勾配法による最大値探索 / Maximize T(V) by gradient method ====
# T(V) = (P0 / (n R V0)) * V * (3V0 - V)
# dT/dV = (P0 / (n R V0)) * (3V0 - 2V) = 0 → V = 3V0 / 2

def T(V):
    return (P0 / (n * R * V0)) * V * (3 * V0 - V)

def dT_dV(V):
    return (P0 / (n * R * V0)) * (3 * V0 - 2 * V)

# 初期値と学習率 / Initial guess and learning rate
V = V0  # 初期体積から開始
alpha = 0.01
eps = 1e-8
max_iter = 1000

for _ in range(max_iter):
    grad = dT_dV(V)
    if abs(grad) < eps:
        break
    V += alpha * grad

T_max = T(V)

# 結果出力 / Print result
print(f"Maximum temperature T_max ≈ {T_max:.5f} [K] at V ≈ {V:.5f} [m³]")

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
P0 = 100
V0 = 1.0

# 体積Vの範囲
V = np.linspace(0.5 * V0, 3.0 * V0, 500)

# Qの計算
Q = (P0 / (2 * V0)) * (V - V0) * (11 * V0 - 4 * V)

# 最大値を探す
max_index = np.argmax(Q)
max_V = V[max_index]
max_Q = Q[max_index]

# プロット
plt.figure(figsize=(8, 5))
plt.plot(V, Q, label=r'$Q = \frac{P_0}{2V_0}(V - V_0)(11V_0 - 4V)$')
plt.axvline(x=V0, color='gray', linestyle='--', label=r'$V_0$')
plt.scatter(max_V, max_Q, color='red', label=f'Max Q at V={max_V:.2f}')
plt.xlabel('Volume $V$')
plt.ylabel('Heat $Q$')
plt.title('Heat $Q$ vs Volume $V$')
plt.legend()
plt.grid(True)
plt.show()

問題53

# --- 定数定義 ---
p0 = 1.0e5       # 外圧 [Pa]
S = 0.01         # ピストン断面積 [m^2]
M = 1.0          # ピストン質量 [kg]
g = 9.8          # 重力加速度 [m/s^2]
k = 100.0        # ばね定数 [N/m]
n = 1.0          # モル数
R = 8.314        # 気体定数 [J/mol·K]

l0 = 0.40        # 初期長さ [m]
l1 = 0.30        # 最終長さ [m]

# --- 仕事 W の計算（赤色面積の計算） ---
W = (p0 * S + 0.5 * M * g) * (l0 - l1)
print(f"仕事 W = {W:.2f} J")

# --- ΔU の計算（理想気体の内部エネルギー変化） ---
delta_U = (3/2) * (p0 * S * (l0 - l1)) + (3/2) * M * g * l0
print(f"内部エネルギー変化 ΔU = {delta_U:.2f} J")

# --- 加えた熱量 Q の計算（Q = ΔU + W） ---
Q = delta_U + W
print(f"加えられた熱量 Q = {Q:.2f} J")

# --- ばねのエネルギー項を含む別解の検証 ---
spring_energy = 0.5 * k * (l0 - l1)**2
rhs = M * g * (l0 - l1) + p0 * S * (l0 - l1)
lhs = W + spring_energy

print(f"\n【別解検証】")
print(f"W + 1/2 k(l0 - l1)^2 = {lhs:.2f} J")
print(f"= Mg(l0 - l1) + p0S(l0 - l1) = {rhs:.2f} J")

問題54

# --- 定数 ---
R = 8.314    # 気体定数 [J/mol·K]
T = 300      # 初期温度 [K]
V = 1.0      # 容器の体積 [m^3]
n = 1.0      # 総モル数 [mol]

# --- 各領域のモル数 ---
nA = (3/5) * n
nB = (2/5) * n

# --- (1) 混合後の温度 T' の計算（内部エネルギー保存より） ---
T_prime = ( (3/2)*nA*R*2*T + (3/2)*nB*R*T ) / ( (3/2)*n*R )
# または T_prime = (8/5)*T
print(f"(1) 混合後の温度 T' = {T_prime:.2f} K")

# --- (2) 混合後の圧力 P' の計算（状態方程式） ---
V_total = 2 * V + V  # 合計体積
P_prime = n * R * T_prime / V_total
print(f"(2) 混合後の圧力 P' = {P_prime:.2f} Pa")

# --- 検算：公式 P' = (8nRT)/(15V) でも計算 ---
P_check = (8 * n * R * T) / (15 * V)
print(f"(検算) P' = {P_check:.2f} Pa")

問題55：熱力学

(1) 真空への膨張（温度一定）

初期の状態方程式：

$$
P = \frac{nRT}{V}
$$

真空に体積 $2V$ まで膨張 ⇒ 等温条件での圧力：

$$
P' = \frac{nRT}{2V} = \frac{P}{2}
$$

(2) 温度変化の確認（断熱膨張）

断熱条件：$PV^\gamma = \text{const}$ より、

$$
P V^\gamma = P'(2V)^\gamma \Rightarrow P' = \frac{P}{2^\gamma}
$$

状態方程式より温度 $T'$：

$$
T' = \frac{P' \cdot 2V}{nR} = \frac{2P'}{nR} \cdot V = \frac{2}{2^\gamma} \cdot \frac{PV}{nR} = \frac{T}{2^{\gamma - 1}}
$$

単原子分子の場合 $\gamma = \frac{5}{3}$ よって：

$$
T' = \frac{T}{2^{2/3}} < T
$$

(3) 圧力一定の定圧変化による温度変化

断熱でなく、定圧 $P$ のままピストンで押された場合：

仕事：

$$
W = P \Delta V = \frac{nRT}{V} \Delta V
$$

状態方程式を使って：

$$
P (2V - \Delta V) = nR T'' \Rightarrow T'' = \frac{P}{nR} (2V - \Delta V)
= \frac{nRT}{nR V} (2V - \Delta V) = \frac{T}{V} (2V - \Delta V)
$$

一方、内部エネルギー変化は：

$$
\Delta U = \frac{3}{2} nR (T'' - T)
$$

よって、エネルギー保存：

$$
\frac{3}{2} nR (T'' - T) = \frac{nRT}{V} \Delta V
$$

解くと：

$$
T'' = \frac{7}{5}T, \quad \Delta V = \frac{3}{5}V
$$

# Program Name: gas_expansion_analysis.py
# Creation Date: 20250722
# Overview: Print symbolic results for gas expansion cases: isothermal, adiabatic, and isobaric
# Usage: Run the script to display derived expressions for pressure, temperature, and volume changes

from sympy import symbols, Eq, solve, Rational, pprint, init_printing

# 日本語表示が崩れないよう数式表示を設定 / Set pretty printing
init_printing()

# ==== 変数定義 / Define variables ====
P, V, T, n, R, gamma, Pp, Tp, Vp, Tpp, dV = symbols('P V T n R γ P\' T\' V\' T\'\' ΔV')

# ==== (1) 等温膨張 / Isothermal expansion ====
Pp_expr = P / 2
print("(1) 等温膨張後の圧力 P'：")
pprint(Eq(Pp, Pp_expr))
print()

# ==== (2) 断熱膨張 / Adiabatic expansion ====
Pp_ad_expr = P / (2 ** gamma)
Tp_expr = T / (2 ** (gamma - 1))
print("(2) 断熱膨張後の温度 T'：")
pprint(Eq(Tp, Tp_expr))
print()

# ==== γ = 5/3 の場合（単原子分子）/ For monoatomic gas γ = 5/3 ====
Tp_value = Tp_expr.subs(gamma, Rational(5, 3))
print("γ = 5/3 のときの T'：")
pprint(Eq(Tp, Tp_value))
print()

# ==== (3) 定圧変化 / Isobaric expansion ====
Tpp_expr = T / V * (2 * V - dV)
delta_U = Rational(3, 2) * n * R * (Tpp - T)
W_expr = P * dV
eq_energy = Eq(delta_U, W_expr.subs(P, n * R * T / V))

print("(3) 定圧変化のエネルギー保存式：")
pprint(eq_energy)
print()

# ==== T'' を解く / Solve for T'' ====
sol = solve(eq_energy, Tpp)[0]
print("定圧変化後の温度 T''：")
pprint(Eq(Tpp, sol))
print()

# ==== T'' = 7/5 T, ΔV = 3/5 V の確認 / Verify final values ====
Tpp_check = T * Rational(7, 5)
dV_check = Rational(3, 5) * V
print("最終結果の確認：")
pprint(Eq(Tpp, Tpp_check))
pprint(Eq(dV, dV_check))

問題57

状態方程式の微小変化：

$$
\frac{\Delta P}{P} + \frac{\Delta V}{V} = \frac{\Delta T}{T}
$$

エネルギー保存より（断熱変化, $Q = 0$）：

$$
\Delta U = W = -P \Delta V
$$

単原子分子気体の内部エネルギー：

$$
\Delta U = \frac{3}{2} nR \Delta T = \frac{3}{2} \frac{PV}{T} \Delta T
$$

これらを連立すると：

$$
\Delta T = -\frac{2T}{3V} \Delta V
$$

# --- 定数設定 ---
P = 1.0e5       # 初期圧力 [Pa]
V = 1.0         # 初期体積 [m^3]
T = 300         # 初期温度 [K]
delta_V = 0.01  # 微小な体積変化 [m^3]

# --- ΔT の理論値計算 ---
delta_T = - (2 * T / (3 * V)) * delta_V
print(f"ΔT = {delta_T:.2f} K")

# --- ΔP の理論値計算（状態方程式の微小変化から） ---
# ΔP/P + ΔV/V = ΔT/T から ΔP = P * (ΔT/T - ΔV/V)
delta_P = P * (delta_T / T - delta_V / V)
print(f"ΔP = {delta_P:.2f} Pa")

# --- ΔU と W（= -PΔV）のチェック ---
nR = P * V / T
delta_U = (3/2) * nR * delta_T
W = -P * delta_V

print(f"内部エネルギー変化 ΔU = {delta_U:.2f} J")
print(f"仕事 W = {W:.2f} J")
print(f"Q = ΔU - W = {delta_U - W:.2e} J（断熱なので ≈ 0）")

問題62

import numpy as np
import pandas as pd

# --- Constants ---
m = 6                  # mass [kg]
g = 10                 # gravity [m/s^2]
lambda_cm = 10        # wavelength [cm]
lambda_m = lambda_cm / 100  # wavelength [m]
mass_density = 1.8e-3 / (120 * 1e-2)  # linear density [kg/m]
rho = mass_density
L_string = 0.3         # string length [m]

# --- Wave speed ---
v = np.sqrt(m * g / rho)

# --- Frequency ---
f = v / lambda_m

# --- Case: 4x mass ---
m1 = 4 * m
lambda1 = 2 * lambda_m
f1 = np.sqrt(4 * m * g / rho) / lambda1

# --- Next resonant length ---
half_lambda1_cm = lambda1 * 100 / 2
next_resonant_length = 40  # [cm]
move_distance = next_resonant_length - 35  # [cm]

# --- Case: 4x linear density (due to 2x diameter) ---
rho2 = 4 * rho
lambda2 = lambda_m / 2
f2 = np.sqrt(m * g / rho2) / lambda2
belly_spacing_cm = lambda2 * 100 / 2
belly_count = 30 / belly_spacing_cm

# --- Results Table ---
results = {
    "項目": [
        "波の速さ [m/s]",
        "周波数 [Hz]",
        "質量4倍のときの周波数 [Hz]",
        "次の共振長さ [cm]",
        "移動距離 [cm]",
        "線密度4倍のときの周波数 [Hz]",
        "腹の間隔 [cm]",
        "腹の数"
    ],
    "値": [
        v,
        f,
        f1,
        next_resonant_length,
        move_distance,
        f2,
        belly_spacing_cm,
        belly_count
    ]
}

# --- Display DataFrame ---
df = pd.DataFrame(results)
df

問題66

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# Wave parameters
A = 1.0           # Amplitude
f = 1.0           # Frequency [Hz]
λ = 1.0           # Wavelength [m]
L = 2.0           # Wall position [m]
v = f * λ         # Wave speed [m/s]

# Space and time setup
x = np.linspace(0, L, 500)
t_values = np.linspace(0, 2, 100)  # 2 seconds total

# Create the figure and axes
fig, ax = plt.subplots(figsize=(10, 6))
line1, = ax.plot([], [], 'r--', label='Incident wave y1')
line2, = ax.plot([], [], 'b:', label='Reflected wave y2')
line3, = ax.plot([], [], 'k-', label='Standing wave y_total', linewidth=2)

ax.set_xlim(0, L)
ax.set_ylim(-2 * A, 2 * A)
ax.set_xlabel("Position x [m]")
ax.set_ylabel("Displacement y")
ax.set_title("Wave Motion Animation")
ax.grid(True)
ax.legend()

# Initialization function
def init():
    line1.set_data([], [])
    line2.set_data([], [])
    line3.set_data([], [])
    return line1, line2, line3

# Update function for animation
def update(t):
    y1 = A * np.sin(2 * np.pi * (f * t - x / λ))  # Incident wave
    y2 = -A * np.sin(2 * np.pi * (f * t + (x - 2 * L) / λ))  # Reflected wave (fixed-end reflection)
    y_total = y1 + y2  # Standing wave
    line1.set_data(x, y1)
    line2.set_data(x, y2)
    line3.set_data(x, y_total)
    return line1, line2, line3

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=t_values, init_func=init, blit=True, interval=50)

plt.show()

# Install required package
!pip install matplotlib

# Import libraries
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML

# Wave parameters
A = 1.0           # Amplitude
f = 1.0           # Frequency [Hz]
λ = 1.0           # Wavelength [m]
L = 2.0           # Wall position [m]
v = f * λ         # Wave speed [m/s]

# Space and time setup
x = np.linspace(0, L, 500)
t_values = np.linspace(0, 2, 100)  # 2 seconds

# Define the standing wave using trigonometric identity
def standing_wave(x, t):
    return 2 * A * np.sin(2 * np.pi * (L - x) / λ) * np.cos(2 * np.pi * (f * t - L / λ))

# Plot setup
fig, ax = plt.subplots(figsize=(10, 5))
line, = ax.plot([], [], 'k-', label='Standing wave', linewidth=2)
ax.set_xlim(0, L)
ax.set_ylim(-2*A, 2*A)
ax.set_xlabel("Position x [m]")
ax.set_ylabel("Displacement y")
ax.set_title("Standing Wave Animation")
ax.grid(True)
ax.legend()

# Initialization
def init():
    line.set_data([], [])
    return line,

# Update each frame
def update(t):
    y = standing_wave(x, t)
    line.set_data(x, y)
    return line,

# Create animation
ani = animation.FuncAnimation(fig, update, frames=t_values, init_func=init, blit=True)

# Display animation in Colab
HTML(ani.to_jshtml())

import numpy as np
import matplotlib.pyplot as plt

# Parameters
A = 1.0           # Base amplitude
λ = 1.0           # Wavelength
L_values = np.linspace(0, 2, 500)  # Vary L from 0 to 2λ

# Amplitude of the resulting wave
A_result = 2 * A * np.abs(np.sin(2 * np.pi * L_values / λ))

# Plot
plt.figure(figsize=(10, 5))
plt.plot(L_values, A_result, label=r"$A' = 2A |\sin(\frac{2\pi L}{\lambda})|$")
plt.axhline(2*A, color='gray', linestyle='--', label='Maximum Amplitude')
plt.title("Amplitude of Resulting Wave vs. Wall Position L")
plt.xlabel("Wall Position L [m]")
plt.ylabel("Resulting Amplitude A'")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題70

import numpy as np
import matplotlib.pyplot as plt

# Constants
f0 = 440         # Emitted frequency in Hz
V = 340          # Speed of sound in m/s
v = 50           # Speed of source in m/s
d = 100          # Distance to observer in meters

# Angle from 0 to 180 degrees
theta_deg = np.linspace(0, 180, 500)
theta_rad = np.radians(theta_deg)

# Full expression for arrival time difference T using cosine law and binomial approximation
# Exact P'A using cosine law
PA_prime = np.sqrt(d**2 + (v / f0)**2 - 2 * d * v * np.cos(theta_rad) / f0)
t2 = 1 / f0 + PA_prime / V
t1 = d / V  # first wave arrives at distance d
T_exact = t2 - t1

# Approximate using binomial expansion (first-order)
T_approx = (1 / f0) + (1 / V) * ((d * (1 - (v * np.cos(theta_rad)) / (f0 * d))) - d)

# Observed frequency f = 1 / T
f_exact = 1 / T_exact
f_approx = 1 / T_approx

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(theta_deg, f_exact, label="Exact Doppler Frequency", color='blue')
plt.plot(theta_deg, f_approx, label="Approximation", color='red', linestyle='dashed')
plt.axhline(f0, color='gray', linestyle='--', label='Emitted Frequency $f_0$')
plt.title("Observed Frequency vs. Angle θ (Moving Source)")
plt.xlabel("Angle θ [degrees]")
plt.ylabel("Observed Frequency [Hz]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題72

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

# ---------------------------
# Part 1: Plot lens behavior
# ---------------------------

# Lens parameters
f = 8  # focal length in cm

# Object distances: just beyond focal point to 100 cm
a_vals = np.linspace(8.1, 100, 500)
b_vals = 1 / (1 / f - 1 / a_vals)        # image distance using lens formula
magnifications = -b_vals / a_vals        # magnification

# Plotting image distance
plt.figure(figsize=(10, 5))
plt.plot(a_vals, b_vals, label="Image Distance b(a)", color="orange")
plt.axvline(f, color='gray', linestyle='--', label='Focal Length f')
plt.title("Image Distance vs. Object Distance (Convex Lens)")
plt.xlabel("Object Distance a [cm]")
plt.ylabel("Image Distance b [cm]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# Plotting magnification
plt.figure(figsize=(10, 5))
plt.plot(a_vals, magnifications, color='green', label="Magnification M = -b/a")
plt.axhline(0, color='gray', linestyle='--')
plt.title("Magnification vs. Object Distance")
plt.xlabel("Object Distance a [cm]")
plt.ylabel("Magnification M")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

# ---------------------------
# Part 2: Solve lens equation
# ---------------------------

# Define symbol
a = sp.Symbol('a')

# Given equation: 1/a + 1/(50 - a) = 1/8
equation = 1/a + 1/(50 - a) - 1/8

# Solve symbolically
solutions = sp.solve(equation, a)
solutions = [sp.N(sol) for sol in solutions]  # convert to numeric

# Compute required shift
required_shift = abs(solutions[1] - solutions[0])

solutions, required_shift

問題73

# Re-declaring constants to ensure full context
a = 4.0             # slope a
b = 1.6             # intercept b
a_prime = 5.0       # slope a'
b_prime = 0.6       # intercept b'
time_per_div = 5.0e-6  # seconds per division
c = 3.0e8           # speed of light in m/s

# Rotation measurement
a_div = 4.0                          # divisions
t = a_div * time_per_div            # total time

teeth = 200                         # number of teeth on the wheel
rotation_per_sec = 1 / (t / teeth)  # full rotations per second

# Water delay effect
d = 500  # one-way distance in meters
delay_water = (2 * d / (c / 1.3)) - (2 * d / c)
div_per_delay = delay_water / time_per_div

# Equation solving: am + b = a'm + b'
m = (b_prime - b) / (a - a_prime)
total_div = a * m + b
l = 0.5 * c * total_div * time_per_div

{
    "Time for 4 divisions [s]": t,
    "Rotation speed [rot/s]": rotation_per_sec,
    "Extra time due to water [s]": delay_water,
    "Extra divisions due to water": div_per_delay,
    "m (interference index)": m,
    "Total divisions (am + b)": total_div,
    "Distance l [m]": l
}

追加47

# Given percent changes
dV_over_V = 0.03
dT_over_T = -2/3 * dV_over_V  # from thermodynamic relation

# Solve for pressure change using: ΔP/P + ΔV/V = ΔT/T
dP_over_P = dT_over_T - dV_over_V

# Output
print("Volume change ΔV/V:", dV_over_V)
print("Temperature change ΔT/T:", dT_over_T)
print("Pressure change ΔP/P:", dP_over_P)
print("Pressure change (%):", dP_over_P * 100)

追加73

# Parameters
N = 720     # number of teeth (example value)
l = 5000    # one-way distance in meters (example)
f0 = 625    # measured rotation frequency in Hz

# Calculate speed of light
c = 4 * N * l * f0

# Output
print("Estimated speed of light c =", c, "m/s")

問題1：波の干渉

import numpy as np
import matplotlib.pyplot as plt

# --- Parameters ---
A = 1.0              # Amplitude
λ = 2.0              # Wavelength
v = 1.0              # Wave velocity
T = λ / v            # Period
ω = 2 * np.pi / T    # Angular frequency
k = 2 * np.pi / λ    # Wave number
δ1 = 0.0             # Initial phase of wave 1
δ2 = 0.0             # Initial phase of wave 2

# --- Space and Time ---
x = np.linspace(0, 4 * λ, 1000)
t = 0  # fixed time snapshot

# --- Component Waves ---
y1 = A * np.sin(ω * t - k * x + δ1)
y2 = A * np.sin(ω * t + k * x + δ2)

# --- Resulting Standing Wave ---
y = y1 + y2

# --- Plot ---
plt.figure(figsize=(10, 4))
plt.plot(x, y1, '--', label='y1 (right-going wave)')
plt.plot(x, y2, '--', label='y2 (left-going wave)')
plt.plot(x, y, label='Standing Wave y = y1 + y2', color='red')
plt.xlabel("x")
plt.ylabel("Displacement")
plt.title("Standing Wave Formed by Two Opposite Traveling Waves")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()


import numpy as np

# Define lambda symbolically for generality
λ = 1  # Wavelength (set as 1 unit for simplicity)

# Define function to compute x based on m
def x_value(m):
    numerator = 9 - (m + 0.5)**2
    denominator = 2 * m + 1
    return (numerator / denominator) * λ

# Evaluate for m = 0 to 5
results = {}
for m in range(6):
    x = x_value(m)
    results[m] = x

results

問題2：波動

import numpy as np
import matplotlib.pyplot as plt

# 定数
λ = 1.0     # 波長 [任意単位]
T = 1.0     # 周期 [s]

# 角度 θ を 0〜90度（ラジアン）で定義
theta_deg = np.linspace(1, 89, 500)
theta_rad = np.deg2rad(theta_deg)

# vx, vy を計算
vx = λ / (T * np.sin(theta_rad))
vy = λ / (T * np.cos(theta_rad))

# グラフ描画
plt.figure(figsize=(10, 5))
plt.plot(theta_deg, vx, label=r'$v_x = \frac{\lambda}{T \sin\theta}$', color='blue')
plt.plot(theta_deg, vy, label=r'$v_y = \frac{\lambda}{T \cos\theta}$', color='green')
plt.xlabel(r'$\theta$ [degrees]')
plt.ylabel('Speed component')
plt.title('Wave velocity components vs θ')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題3：光の干渉

AP:

$$
AP = \sqrt{l^2 + \left(x + \frac{d}{2}\right)^2} \approx l \left{1 + \frac{1}{2} \left(\frac{x + \frac{d}{2}}{l}\right)^2 \right}
$$
BP:

$$
BP = \sqrt{l^2 + \left(x - \frac{d}{2}\right)^2} \approx l \left{1 + \frac{1}{2} \left(\frac{x - \frac{d}{2}}{l}\right)^2 \right}
$$
差：

$$
AP - BP \approx \frac{1}{2l} \left[ \left(x + \frac{d}{2}\right)^2 - \left(x - \frac{d}{2}\right)^2 \right] = \frac{dx}{l}
$$

def approx_ap_bp_difference(x, d, l):
    """
    Approximate AP - BP ≈ dx / l using second-order expansion.
    
    Parameters:
    x : float - small displacement from center
    d : float - separation between paths
    l : float - original length
    
    Returns:
    AP, BP : approximate distances
    diff   : AP - BP ≈ dx / l
    """
    term_ap = (x + d/2) / l
    term_bp = (x - d/2) / l

    ap = l * (1 + 0.5 * term_ap**2)
    bp = l * (1 + 0.5 * term_bp**2)
    diff = ap - bp

    # Also return theoretical dx/l value
    theory = (d * x) / l

    return ap, bp, diff, theory

# --- Example values ---
x = 0.01  # [m]
d = 0.005 # [m]
l = 1.0   # [m]

ap, bp, diff, theory = approx_ap_bp_difference(x, d, l)

print(f"AP ≈ {ap:.6f} m")
print(f"BP ≈ {bp:.6f} m")
print(f"AP - BP ≈ {diff:.6e} m")
print(f"dx/l (theory) = {theory:.6e} m")

問題6：波動

平方根の近似展開：

$$
(1 + x)^{a} \approx 1 + ax \quad \text{（ただし } |x| \ll 1\text{ のとき）}
$$

import math

# --- 例1: √(4 + 0.31) を近似 ---
def approx_sqrt_near(n_base, delta):
    """
    Approximate sqrt(n_base + delta) using first-order binomial expansion:
    sqrt(a + x) ≈ sqrt(a) * (1 + (x / (2a)))
    """
    approx = math.sqrt(n_base) * (1 + delta / (2 * n_base))
    exact = math.sqrt(n_base + delta)
    return approx, exact

n_base = 4
delta = 0.31
approx1, exact1 = approx_sqrt_near(n_base, delta)

print(f"Approx sqrt(4 + 0.31): {approx1:.6f}")
print(f"Exact   sqrt(4 + 0.31): {exact1:.6f}")


# --- 例2: √59 ≈ 8 * sqrt(1 - 5/64) を近似 ---
def approx_sqrt_59():
    """
    Approximate sqrt(59) using:
    sqrt(59) = 8 * sqrt(1 - 5/64) ≈ 8 * (1 - 5/128)
    """
    x = -5 / 64
    approx = 8 * (1 + 0.5 * x)
    exact = math.sqrt(59)
    return approx, exact

approx2, exact2 = approx_sqrt_59()

print(f"\nApprox sqrt(59): {approx2:.6f}")
print(f"Exact   sqrt(59): {exact2:.6f}")


# --- 例3: 波長・干渉式の計算（d = 682λ / 59） ---
def calc_d(lambda_nm=682, m=30, n2=509, n1=118):
    """
    Calculate d using the formula:
    d = (lambda * m * sqrt(n2/n1 - 0.5)) [nm]
    """
    ratio = n2 / n1
    value = math.sqrt(ratio - 0.5)
    d_nm = lambda_nm * m * value / 59
    return d_nm / 1000  # convert to μm

d_micrometers = calc_d()
print(f"\nResulting d ≈ {d_micrometers:.3f} μm")

問題7：光の干渉

# 基本パラメータ
F = 100         # 焦点距離 [cm]
λ = 656e-9      # 波長 [m]
d = 1 / 400     # スリット間隔 [m]（格子定数：400本/m）

# 明線の間隔 Δx を計算
delta_x = F * λ / d  # [m]
delta_x_cm = delta_x * 100  # cm に変換

print(f"明線の間隔 Δx = {delta_x_cm:.2f} cm")

# --- ドップラー効果による波長変化 ---
# Δx の変化量（問題で与えられた）:
delta_lambda = 0.011 * d / F  # 式：Δx - Δx' = (F/d)(λ - λ') → λ - λ' = Δλ

# 新しい波長
lambda_prime = λ - delta_lambda

# 元の振動数 f と観測された f' の計算（ドップラー公式）
c = 3.0e8  # 光速 [m/s]
f = c / λ
f_prime = c / lambda_prime

# 星雲の接近速度 v を計算（ドップラー効果逆計算）
v = c * (f_prime - f) / f_prime

print(f"変化後の波長 λ' = {lambda_prime:.2e} m")
print(f"観測された周波数 f' = {f_prime:.2e} Hz")
print(f"星雲の接近速度 v = {v:.2f} m/s")

問題9

# Given constants
C = 1500 / (91.67 - 20.00)                   # Heat capacity [J/K]
R = 8.3                                      # Gas constant [J/(mol·K)]
delta_T = 90.74 - 20.00                      # Temperature change [K]

# Solve for mol of Argon gas
lhs_energy = 1500e3
x = lhs_energy / ((3/2) * R * delta_T + C * delta_T)

# Results for (2)
mol_argon = x

# For (3)
n = 1 + (178 * 6.3e-7) / 0.20                # Refractive index at current density
n_minus_1 = n - 1
rho = 2.2e-2 / 2.5e-4                        # Current density [kg/m^3]
M = 2.5e-2                                   # Molar mass of argon [kg/mol]
rho_0 = M / 2.24e-2                          # Standard density [kg/m^3]

# Use formula: n0 = 1 + (rho0 / rho) * (n - 1)
n0 = 1 + (rho_0 / rho) * n_minus_1

mol_argon, n, n0

問題13：静電気単振動

等速運動に移るまでの時間：

$$
t_2 = \frac{x_1 - x_4}{V} = \frac{2(x_1 - x_2)}{V} = \frac{2mg}{aqV}(\mu_1 - \mu_2)
$$

def calc_time(m, g, a, q, V, mu1, mu2):
    """
    単振動から等速運動に移るまでの時間 t2 を計算する。
    
    引数:
    m -- 質量 [kg]
    g -- 重力加速度 [m/s^2]
    a -- 加速度 [m/s^2]
    q -- 係数（バネ定数などの比例定数）
    V -- ベルトの速度 [m/s]
    mu1 -- 動摩擦係数
    mu2 -- 静止摩擦係数
    
    戻り値:
    t2 -- 時間 [s]
    """
    t2 = (2 * m * g / (a * q * V)) * (mu1 - mu2)
    return t2

# --- 使用例 ---
m = 0.5       # kg
g = 9.8       # m/s^2
a = 2.0       # m/s^2
q = 1.5       # 任意の定数（例: バネ定数に相当）
V = 0.3       # m/s
mu1 = 0.4     # 動摩擦係数
mu2 = 0.2     # 静止摩擦係数

t2 = calc_time(m, g, a, q, V, mu1, mu2)
print(f"ベルト上で等速運動に移るまでの時間 t2 = {t2:.4f} 秒")

# --- Constants ---
R = 8.3                # Gas constant [J/mol·K]
C = 1500 / 71.67       # Heat capacity [J/K]
delta_T = 90.74 - 20.00

# --- (2) Calculate mol number x ---
numerator = 1.5e3
denominator = (3/2) * R * delta_T + C * delta_T
x = numerator / denominator  # mol
print(f"Mol number x ≈ {x:.3e} mol")

# --- (3) Calculate refractive index n ---
n = 1 + (178 * 6.3e-7) / 0.20
print(f"Refractive index n ≈ {n:.6f}")

# --- Density calculation ---
M = 1.0  # Assume molar mass of argon = 1.0 (will cancel later)
V = 2.5e-4
rho = x * M / V
print(f"Density rho ≈ {rho:.3e} kg/m³")

# --- Reference density and final refractive index n0 ---
rho0 = 1.0  # arbitrary unit since it cancels with M
n_minus_1 = n - 1
n0 = 1 + (rho0 / rho) * n_minus_1
print(f"Final refractive index n0 ≈ {n0:.6f}")

問題15：静電気

電位差：

$$
V = \left( \frac{b - a}{ab} \right) kQ
$$

電気容量 $C$：

$$
C = \frac{Q}{V} = \frac{ab}{k(b - a)}
$$

$k = \frac{1}{4\pi\varepsilon}$ のとき：

$$
C = \frac{4\pi\varepsilon ab}{b - a}
$$

import math

def calc_capacitance(a, b, epsilon):
    """
    導体間の電気容量Cを計算する（球導体 or 平行板近似）。
    
    引数:
    a -- A の位置（内側導体）[m]
    b -- B の位置（外側導体）[m]
    epsilon -- 誘電率 [F/m]
    
    戻り値:
    C -- 電気容量 [F]
    """
    C = (4 * math.pi * epsilon * a * b) / (b - a)
    return C

# --- 使用例 ---
a = 0.01       # A の半径 [m]
b = 0.02       # B の半径 [m]
epsilon = 8.854e-12  # 真空の誘電率 [F/m]

cap = calc_capacitance(a, b, epsilon)
print(f"電気容量 C = {cap:.3e} F")

問題19：キャパシタ

以下は、画像に基づいたコンデンサーの電気容量と電圧、静電エネルギー変化を計算する Python コードです。

総容量：

$$
C = \left( 1 + (\varepsilon_r - 1) \frac{x}{l} \right) C_0
$$

電圧：

$$
V_1 = \frac{Q_0}{C} = \frac{l}{l + (\varepsilon_r - 1)x} V
$$

外力のした仕事（エネルギー変化）：

$$
W_1 = \frac{(\varepsilon_r - 1) C_0 V^2 x}{2[l + (\varepsilon_r - 1)x]}
$$

def calc_capacitance_voltage_energy(epsilon_r, x, l, C0, V):
    """
    誘電体挿入によるコンデンサーの変化を計算

    引数:
    epsilon_r : 誘電率の比 (εr)
    x         : 誘電体の長さ [m]
    l         : コンデンサー全体の長さ [m]
    C0        : 元の電気容量 [F]
    V         : 元の電圧 [V]

    戻り値:
    C         : 合成容量 [F]
    V1        : 誘電体挿入後の電圧 [V]
    W1        : 外力による仕事（エネルギー変化）[J]
    """
    # 合成容量
    C = (1 + (epsilon_r - 1) * x / l) * C0

    # 電荷Q0は変化しない
    Q0 = C0 * V

    # 誘電体挿入後の電圧
    V1 = Q0 / C

    # 外力のした仕事（静電エネルギー変化）
    W1 = ((epsilon_r - 1) * C0 * V**2 * x) / (2 * (l + (epsilon_r - 1) * x))

    return C, V1, W1


# --- 使用例 ---
epsilon_r = 5.0     # 誘電率
x = 0.01            # 誘電体の長さ [m]
l = 0.02            # 全長 [m]
C0 = 1e-12          # 元の電気容量 [F]
V = 100.0           # 元の電圧 [V]

C, V1, W1 = calc_capacitance_voltage_energy(epsilon_r, x, l, C0, V)

print(f"合成容量 C = {C:.3e} F")
print(f"電圧 V1 = {V1:.2f} V")
print(f"外力による仕事 W1 = {W1:.3e} J")

問題22：キャパシタ

与えられる漸化式：

$$
x_{n+1} = \frac{1}{2} x_n + V_0 \tag{1}
$$

または変形して：

$$
x_{n+1} - 2V_0 = \frac{1}{2}(x_n - 2V_0) \tag{2}
$$

これは公比 $\frac{1}{2}$ の等比数列として収束：

$$
x_n = 2V_0 - \frac{V_0}{2^{n-1}} \quad \Rightarrow \quad x_\infty = 2V_0
$$

import matplotlib.pyplot as plt

def capacitor_sequence(V0, n_terms):
    """
    Generate the sequence x_n defined by:
    x_{n+1} = 0.5 * x_n + V0
    starting with x1 = V0
    """
    x_vals = [V0]  # x1 = V0
    for n in range(1, n_terms):
        x_next = 0.5 * x_vals[-1] + V0
        x_vals.append(x_next)
    return x_vals

# --- Parameters ---
V0 = 5.0           # Initial voltage [V]
n_terms = 20       # Number of iterations

# --- Generate sequence ---
x_sequence = capacitor_sequence(V0, n_terms)

# --- Final value (limit) ---
x_limit = 2 * V0

# --- Plot ---
plt.figure(figsize=(8, 4))
plt.plot(range(1, n_terms + 1), x_sequence, 'o-', label='xₙ')
plt.axhline(x_limit, color='red', linestyle='--', label=f'Limit = {x_limit:.1f}')
plt.title("Convergence of Capacitor Voltage xₙ")
plt.xlabel("n (operation count)")
plt.ylabel("xₙ [V]")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題29

import numpy as np
import matplotlib.pyplot as plt

# Data points: Current I [mA], Voltage V [V]
current = np.array([10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150])
voltage = np.array([1.638, 1.634, 1.630, 1.625, 1.620, 1.615, 1.610, 1.605, 1.600, 
                    1.595, 1.590, 1.585, 1.580, 1.575, 1.570])

# Linear regression (least squares method)
coefficients = np.polyfit(current, voltage, 1)
slope, intercept = coefficients

# Calculate the fitted line
fitted_voltage = slope * current + intercept

# Print results
print(f"Slope (estimated internal resistance): {slope:.6f} Ω")
print(f"Intercept (estimated electromotive force): {intercept:.6f} V")

# Plotting
plt.scatter(current, voltage, label='Data points')
plt.plot(current, fitted_voltage, color='red', label='Least squares fit')
plt.xlabel('Current I [mA]')
plt.ylabel('Terminal Voltage V [V]')
plt.title('Voltage vs. Current with Linear Fit')
plt.legend()
plt.grid(True)
plt.show()

問題34

import numpy as np

# Constants
mu_0 = 4 * np.pi * 1e-7  # Vacuum permeability [H/m]
I = 1.0  # Current [A]
l = 1.0  # Length of wire segment [m]
r = 0.05  # Horizontal separation between wires [m]
d = 0.03  # Vertical separation between wires [m]

# Magnetic field component from one wire at the position of the other
H1 = I / (2 * np.pi * np.sqrt(r**2 + d**2))

# y-direction component of the total magnetic field from two wires
Hc = 2 * H1 * (r / np.sqrt(r**2 + d**2))

# Magnetic force on the wire
Fc = mu_0 * I * Hc * l

# Alternatively, direct formula from the text:
Fc_direct = (mu_0 * I**2 * r * l) / (np.pi * (r**2 + d**2))

# Output results
print(f"Magnetic field H1: {H1:.4e} A/m")
print(f"Composite magnetic field Hc: {Hc:.4e} A/m")
print(f"Electromagnetic force Fc (from Hc): {Fc:.4e} N")
print(f"Electromagnetic force Fc (direct formula): {Fc_direct:.4e} N")

問題37：電磁誘導

import math

def compute_voltage(M, g, R, B, l, mu0, theta_deg):
    theta = math.radians(theta_deg)
    numerator = M * g * R * (math.sin(theta) + mu0 * math.cos(theta))
    denominator = B * l * (math.cos(theta) - mu0 * math.sin(theta))
    V0 = numerator / denominator
    return V0

# 使用例
M = 0.5         # [kg]
g = 9.8         # [m/s^2]
R = 2.0         # [Ω]
B = 0.4         # [T]
l = 0.3         # [m]
mu0 = 0.3       # 摩擦係数
theta_deg = 30  # 角度[度]

V0 = compute_voltage(M, g, R, B, l, mu0, theta_deg)
print(f"V₀ = {V0:.3f} V")

問題52

import numpy as np

# 定数設定
q = 1.6e-19     # 電荷 [C]
m = 9.1e-31     # 質量 [kg]
E = 1e4         # 電場 [V/m]
B = 1.0         # 磁場 [T]
L = 0.05        # 距離 [m]
v = 1e6         # 初速度 [m/s]

# --------------------
# 式 (1) y = qEL^2 / (2mv^2)
# --------------------
y = q * E * L**2 / (2 * m * v**2)

# --------------------
# 式 (2) x ≒ qB L^2 / (2m v)
# --------------------
x = q * B * L**2 / (2 * m * v)

# 結果出力
print(f"Displacement in z-direction (y): {y:.4e} m")
print(f"Displacement in x-direction (x): {x:.4e} m")

問題53

import numpy as np

# 定数
q = 1.6e-19       # 電荷 [C]
B = 1.0           # 磁場 [T]
M = 1.67e-27      # 質量 [kg]
V0 = 5.0          # 振幅 [V]

# 周波数と周期の計算
f = q * B / (2 * np.pi * M)  # [Hz]
T = 1 / f                    # [s]

# 時間配列（2周期分）
t = np.linspace(0, 2*T, 5)  # 5点だけ表示（多すぎるとprintが重くなるため）

# 振動の式と各点での値を表示
print("振動の式: V(t) = V0 * sin(2πft)")
print(f"振幅 V0 = {V0} V")
print(f"電荷 q = {q} C")
print(f"磁場 B = {B} T")
print(f"質量 M = {M} kg")
print(f"周波数 f = {f:.3e} Hz")
print(f"周期 T = {T:.3e} s")

print("\n時刻と対応する振幅 V(t):")
for ti in t:
    vt = V0 * np.sin(2 * np.pi * f * ti)
    print(f"t = {ti:.3e} s -> V(t) = {vt:.3f} V")

問題54

# 定数
W = 2.3 * 1.6e-19          # 光電子の仕事関数 [J]
nu_0 = 5.6e14              # 限界振動数 [Hz]
c = 3.0e8                  # 光速 [m/s]
lambda_light = 6.3e-7      # 波長 [m]
W_prime = 1.9 * 1.6e-19    # Csにおける仕事関数 [J]

# プランク定数の計算
h = W / nu_0

# 光の振動数の計算
nu = c / lambda_light

# 限界振動数（Cs）の計算
nu0_cs = W_prime / h

# 結果出力
print(f"Plank定数 h: {h:.2e} J·s")
print(f"光の振動数 ν: {nu:.2e} Hz")
print(f"Csの限界振動数 ν0: {nu0_cs:.2e} Hz")

# 光電効果が起こるかどうか
if nu > nu0_cs:
    print("光電効果は起こる。")
else:
    print("光電効果は起こらない。")

問題55

import numpy as np
import matplotlib.pyplot as plt

# θ in radians (from 0 to π)
theta = np.linspace(0, np.pi, 100)
one_minus_cos = 1 - np.cos(theta)

# Constant: electron rest energy [MeV]
m_e_c2 = 0.511

# Assume incident photon energy is 0.662 MeV (e.g., from Cs-137)
E0_inv = 1 / 0.662  # inverse of initial energy [1/MeV]
inv_energy_cs = E0_inv + one_minus_cos / m_e_c2  # Theoretical line for Cs
inv_energy_mn = E0_inv + 0.8 * one_minus_cos / m_e_c2  # Slightly different slope for Mn

# Plot
plt.figure(figsize=(6, 6))
plt.plot(one_minus_cos, inv_energy_cs, label='Cs (experimental)', color='black')
plt.plot(one_minus_cos, inv_energy_mn, label='Mn (example)', color='red', linestyle='--')
plt.xlabel('1 - cos(θ)')
plt.ylabel('1 / hν [1/MeV]')
plt.title('Compton Scattering Plot')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題56

import numpy as np

# Constants
c = 3e8                 # Speed of light [m/s]
h = 6.626e-34           # Planck constant [J·s]
nu = 5e14               # Photon frequency [Hz]
L = 1.0                 # Length of cavity side [m]
t = 1.0                 # Time [s]
N = 1e20                # Number of photons

# (1) Number of collisions in time t
n = (c * t) / (2 * L)
print("Number of collisions per photon (n):", n)

# (3) Change in momentum per collision
delta_p = 2 * h * nu / c
print("Momentum change per collision Δp:", delta_p, "kg·m/s")

# Impulse by 1 photon in time t
I1 = delta_p * n
print("Impulse from one photon I₁:", I1, "N·s")

# Total impulse from all photons
I_total = N * I1
print("Total impulse from N photons:", I_total, "N·s")

# (5) Force = total impulse / time
F = I_total / t
print("Force F:", F, "N")

# (6) Pressure = Force / Area (Area = L²)
P = F / L**2
print("Radiation pressure P:", P, "Pa")

# (7) Total energy U of N photons
U = N * h * nu
print("Total energy U:", U, "J")

# Volume of the cube
V = L**3

# Energy density = U / V
energy_density = U / V
print("Energy density (U/V):", energy_density, "J/m³")

# Radiation pressure from energy: P = (1/3) * energy_density
P_from_energy = (1/3) * energy_density
print("Radiation pressure from energy (1/3 U/V):", P_from_energy, "Pa")

問題57

import numpy as np

# Constants
h = 6.6e-34        # Planck's constant [J·s]
m_e = 9.0e-31      # Electron mass [kg]
e = 1.6e-19        # Elementary charge [C]
d = 2.8e-10        # Atomic spacing [m]
theta_deg = 30     # Bragg angle in degrees
theta_rad = np.radians(theta_deg)  # Convert to radians

# (1) Bragg condition: 2d sinθ = nλ, solve for d when n = 4, λ = 7.0e-11
lam = 7.0e-11
n_given = 4
d_calc = (n_given * lam) / (2 * np.sin(theta_rad))
print(f"(1) Calculated d for n=4: {d_calc:.2e} m")

# (2) Max order n from inequality: 0 < nλ/(2d) ≤ 1
n_max = 2 * d / lam
print(f"(2) Max integer n satisfying Bragg condition: {int(n_max)}")

# (3) Electron velocity accelerated by voltage V
V = 1000  # volts
v = np.sqrt(2 * e * V / m_e)
print(f"(3) Electron speed from V = {V} V: {v:.2e} m/s")

# (4) De Broglie wavelength and Bragg condition with electrons
n_range = np.arange(1, 13)
V_list = (h**2 * n_range**2) / (2 * m_e * e * d**2)
print("(4) Voltage range satisfying Bragg condition:")
for n_val, v_val in zip(n_range, V_list):
    print(f"  n = {n_val}, V = {v_val:.2f} V")

# (5) Solve inequality: 1000 ≤ V ≤ 2000 =⇒ find range for n
lower = 1000
upper = 2000
V_expr = lambda n: (h**2 * n**2) / (2 * m_e * e * d**2)
n_lower = np.sqrt(lower * 2 * m_e * e * d**2 / h**2)
n_upper = np.sqrt(upper * 2 * m_e * e * d**2 / h**2)
print(f"(5) Integer n satisfying 1000 ≤ V ≤ 2000: {n_lower:.1f} ≤ n ≤ {n_upper:.1f}")

問題61

import numpy as np
import matplotlib.pyplot as plt

# 定数
Mp_c2 = 1000  # 陽子の静止エネルギー [MeV]

# 陽子が得た運動エネルギー K を 0〜50 MeV の範囲でプロット
K_values = np.linspace(0.1, 50, 500)  # K ≠ 0 にするため 0.1 から
E_values = 0.5 * (np.sqrt(2 * Mp_c2 * K_values) + K_values)

# グラフ描画
plt.figure(figsize=(8, 5))
plt.plot(K_values, E_values, label=r'$E = \frac{1}{2}(\sqrt{2M_pc^2K} + K)$', color='red')
plt.xlabel('Proton Kinetic Energy K [MeV]')
plt.ylabel('Gamma Ray Energy E [MeV]')
plt.title('Gamma Energy Required vs. Proton Kinetic Energy')
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

問題67：原子核

import math

# --- Constants from the problem ---
v = 1.9e6         # Speed [m/s]
d = 0.16e-3       # Distance [m]
survival_fraction_1 = 0.25     # From 25/100 (1/4)
survival_fraction_2 = 0.05     # From 5/100 (1/20)
log10_2 = math.log10(2)

# --- Step (4): Calculate half-life T ---
T = d / (2 * v)                # T = d / (2v)
print(f"T = {T:.2e} [s]")      # Should be ~4.2e-11 s

# --- Step (5): Calculate time t based on decay equation ---
t = ((1 + math.log10(2)) / log10_2) * T
print(f"t = {t:.2e} [s]")      # Time required for 5% survival

# --- Calculate distance ---
d2 = v * t
print(f"d = {d2:.2e} [m] = {d2 * 1000:.2f} [mm]")  # Convert to mm

import numpy as np

# 定数
e = 1.6e-19           # 電子の電荷 [C]
m = 8.8e-26           # 粒子の質量 [kg]
c = 3.0e8             # 光速 [m/s]
E_MeV = 1.0           # 運動エネルギー [MeV]
E_J = E_MeV * 1e6 * e # MeV -> J 変換

# (1) 運動エネルギーから速度を求める
v = np.sqrt(2 * E_J / m)
print(f"Speed v = {v:.2e} m/s")  # => 約 1.9e6 m/s

# (2) ドップラー効果による周波数変化（近似式）
v0 = 1.795e21  # 元の周波数 [Hz] = 1795 × 10^18
nu = (1 + v / c) * v0
print(f"Doppler-shifted frequency ν = {nu:.2e} Hz")  # => 約 1.806e21 Hz

追加10

import numpy as np
import matplotlib.pyplot as plt

# 定数
k = 9.0e9         # クーロン定数 [N·m²/C²]
Q = 1.0e-6        # 電荷 [C]
l = 0.1           # 距離 [m]
E = 1.0e5         # 一様電場 [V/m]

# x軸（lと−lを除外）
x = np.linspace(-0.099, 0.099, 1000)
x = x[x != l]     # 分母が0になる点を除外
x = x[x != -l]

# 電位 V(x)
V = k*Q/(l + x) - k*Q/(l - x) + E*x

# 電場 E(x) = -dV/dx
# 微分の解析式：
# dV/dx = -kQ/(l+x)^2 - kQ/(l−x)^2 + E
dVdx = -k*Q / (l + x)**2 - k*Q / (l - x)**2 + E
E_field = -dVdx  # 電場は V(x) のマイナス微分

# 描画
plt.figure(figsize=(10, 6))

# 電位グラフ
plt.subplot(2, 1, 1)
plt.plot(x, V, label='Potential V(x)', color='blue')
plt.ylabel('V(x) [V]')
plt.title('Electric Potential and Electric Field')
plt.grid(True)
plt.legend()

# 電場グラフ
plt.subplot(2, 1, 2)
plt.plot(x, E_field, label='Electric Field E(x)', color='red')
plt.xlabel('x [m]')
plt.ylabel('E(x) [V/m]')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

追加44

import numpy as np
import matplotlib.pyplot as plt

# 時間範囲
t = np.linspace(0, 2, 500)  # 2周期分

# パラメータ設定
A = 1.0       # 振幅
T = 1.0       # 周期
theta = np.pi / 4  # 例：45度
B = 1.0
d = 1.0
M = 1.0
L = 1.0
g = 9.8

# x(t) = A(1 - cos(2πt/T))
x = A * (1 - np.cos(2 * np.pi * t / T))

# I(t) = (Mg / Bd) * sinθ * (1 - cos(Bd * t / sqrt(ML)))
omega = B * d / np.sqrt(M * L)
I = (M * g / (B * d)) * np.sin(theta) * (1 - np.cos(omega * t))

# グラフ描画
plt.figure(figsize=(10, 6))

# x(t)
plt.subplot(2, 1, 1)
plt.plot(t, x, label='x(t)', color='blue')
plt.ylabel('x(t)')
plt.title('Displacement and Current')
plt.grid(True)
plt.legend()

# I(t)
plt.subplot(2, 1, 2)
plt.plot(t, I, label='I(t)', color='red')
plt.xlabel('t')
plt.ylabel('I(t)')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

追加52


import numpy as np
import matplotlib.pyplot as plt

# Constants for the equation y = A * x^2
# Based on: y = (2mE / (q B^2 L^2)) * x^2
A = 1.0  # You can adjust this based on the physical constants

# x range from the diagram (e.g., -90 to 10)
x_vals = np.linspace(-90, 10, 500)
y_vals = A * x_vals**2

# Plot
plt.figure(figsize=(8, 6))
plt.plot(x_vals, y_vals, label=r"$y = A x^2$", color="red")

# Plot key points like α (x=0, y=2), β (x=-60√2, y=-4), γ (x=-30√2, y=1)
x_alpha, y_alpha = 0, 2
x_beta, y_beta = -60 * np.sqrt(2), -4
x_gamma, y_gamma = -30 * np.sqrt(2), 1

plt.scatter([x_alpha, x_beta, x_gamma], [y_alpha, y_beta, y_gamma], color="black", zorder=5)
plt.text(x_alpha, y_alpha + 0.3, 'α', ha='center')
plt.text(x_beta, y_beta - 0.5, 'β', ha='center')
plt.text(x_gamma, y_gamma + 0.3, 'γ', ha='center')

# Axes settings
plt.axhline(0, color='gray', linewidth=0.5)
plt.axvline(0, color='gray', linewidth=0.5)
plt.xlabel("x")
plt.ylabel("y")
plt.title("Parabolic Trajectory of Charged Particle")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Pythonとコピペと生成AIで学ぶ 名問の森

はじめに

◆ 学習の進め方

問題1：放物運動

問題2：放物運動

問題3：放物運動

問題4：剛体の釣り合い

問題5：剛体の釣り合い

問題7：運動方程式

問題13：運動量保存則

【説明】

問題14：運動量保存則

問題15：保存則

問題17：保存則

問題18：保存則

問題27：等速円運動

問題30：円運動

問題37：単振動

問題38

問題39：単振動

(1) 最大摩擦力に達したときの静止摩擦係数の導出

(2) 力の式と単振動の中心

(3) 両端 $x_0, x_1$ の平均が中心

(4) 摩擦熱を考慮したエネルギー保存

(5) 左右に滑るときの振動中心

(6) 減衰の終わりの位置まで

(7) グラフ

各点の時間と位置のまとめ

問題41：単振動

問題の概要

微分方程式

解の形

問題43

問題44

問題45

問題46：分子運動

(1) 分子が壁と衝突する際の運動量の変化

(2) 衝突回数（時間あたり）

(3) 分子1個あたりの力積と全体の力

(4) 圧力 $P$ の導出

(5) 理想気体の状態方程式と温度の関係

(6) 内部エネルギー $U$

問題47：分子運動

(1) 分子がピストンに衝突する際の速度変化

(2) 運動エネルギーの減少

(3) 単位時間あたりの衝突回数

(4) 単位時間あたりのエネルギー減少量

(5) 容器の体積変化 $\Delta V$

(6) 多数の分子の平均化（全体のエネルギー変化）

(7) 三平方より：

(8) 内部エネルギーの変化

(9) エネルギーと温度の関係

微分方程式と断熱変化の式導出（続き）

(10) 微分形式での温度変化：

(11) 両辺を積分：

(12) よって、断熱変化の式：

問題49：熱力学

問題50

問題52：熱力学

問題53

問題54

問題55：熱力学

(1) 真空への膨張（温度一定）

(2) 温度変化の確認（断熱膨張）

(3) 圧力一定の定圧変化による温度変化

問題57

問題62

問題66

問題70

問題72

問題73

追加47

追加73

問題1：波の干渉

問題2：波動

問題3：光の干渉

問題6：波動

問題7：光の干渉

問題9

問題13：静電気単振動

Pythonとコピペと生成AIで学ぶ名問の森