Pythonコードで学ぶ物理のエッセンス(物理の周辺)



## はじめに  

物理の学びは「式を覚える」ことよりも、「意味をつかむ」ことに価値があります。  
しかし、数式の意味や導出がピンとこないまま、暗記に走ってしまう経験、ありませんか？

このシリーズでは、物理や数学のエッセンスを**Pythonコードで“見て、動かして、確かめる”**ことを通して、理解を深めていきます。




## Physics Math Curriculum

このカリキュラムは、高校物理・数学の基本概念を、Pythonを使って視覚的・体験的に学ぶことを目的としています。誤差や単位変換から始まり、ベクトル・微積分・回路解析、そして三角関数・近似計算までを統一的に扱い、物理の背後にある数理的な構造を深く理解します。



---

### 1. 物理量のスケール感と単位系：SI接頭語とその変換
```python
def convert_physical_quantity(value, from_unit, to_unit):
    """
    SI接頭語付きの単位を変換する関数
    Convert from one SI-prefixed unit to another.

    例: convert_physical_quantity(1000, 'mm', 'm') → 1.0
    """
    prefixes = {
        'G': 1e9, 'M': 1e6, 'k': 1e3,
        '': 1,
        'm': 1e-3, 'μ': 1e-6, 'n': 1e-9
    }

    def extract_prefix(unit):
        for p in prefixes:
            if unit.startswith(p) and unit[len(p):] in ['m', 'g', 's', 'A', 'F', 'V']:
                return p, unit[len(p):]
        raise ValueError(f"Unknown unit format: {unit}")

    prefix_from, base_from = extract_prefix(from_unit)
    prefix_to, base_to = extract_prefix(to_unit)

    if base_from != base_to:
        raise ValueError(f"Incompatible units: {from_unit} vs {to_unit}")

    base_value = value * prefixes[prefix_from]
    return base_value / prefixes[prefix_to]

# 使用例 / Example
print(convert_physical_quantity(1000, 'mm', 'm'))      # 1.0
print(convert_physical_quantity(3.3, 'kV', 'V'))       # 3300.0
print(convert_physical_quantity(500, 'μs', 'ms'))      # 0.5

```

### 2. 測定と物理量の整合性：次元解析と有効数字の基礎
```python
# SymPyによる次元解析
from sympy import symbols
from sympy.physics.units import speed, acceleration
from sympy.physics.units.systems.si import dimsys_SI

v = speed
a = acceleration
print("Speed dimension:", dimsys_SI.get_dimensional_dependencies(v))
print("Acceleration dimension:", dimsys_SI.get_dimensional_dependencies(a))
```
```python
# 有効数字の丸めと評価
import math

def round_sig(x, sig=3):
    if x == 0:
        return 0
    return round(x, sig - int(math.floor(math.log10(abs(x)))) - 1)

x = 12345.6789
print("Original:", x)
print("3 significant figures:", round_sig(x, 3))
```
```python
# sin(x) と近似の比較（Taylor展開）
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-1, 1, 100)
true = np.sin(x)
approx = x - (x**3)/6

plt.plot(x, true, label='sin(x)')
plt.plot(x, approx, '--', label='Approx: x - x³/6')
plt.title('Taylor Approximation of sin(x)')
plt.xlabel('x')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.show()
```

---

### 3. 電流の流れを解き明かす：キルヒホッフ則と行列解法
```python
import numpy as np
import matplotlib.pyplot as plt

# --- 抵抗と電圧の設定 / Set up resistors and voltage ---
R1 = 100  # Resistance in Ω
R2 = 200
R3 = 300
V1 = 10   # Voltage in V

# --- 直列抵抗 / Series resistance ---
# 合成抵抗 = R1 + R2 + R3
R_series = R1 + R2 + R3
print("=== Series Resistance ===")
print(f"R_series = {R_series} Ω\n")

# --- 並列抵抗 / Parallel resistance ---
# 合成抵抗 = 1 / (1/R1 + 1/R2 + 1/R3)
R_parallel = 1 / (1/R1 + 1/R2 + 1/R3)
print("=== Parallel Resistance ===")
print(f"R_parallel = {R_parallel:.2f} Ω\n")

# --- 分圧 / Voltage Divider ---
# Vout = Vin * (Rx / (R1 + Rx))
Rx = R2  # 出力を取り出す抵抗（例としてR2を使用） / Select R2 for output
Vout = V1 * (Rx / (R1 + Rx))
print("=== Voltage Divider ===")
print(f"Vout across R2 = {Vout:.2f} V\n")

# --- 連立方程式の構築: A x = b / Define matrix A and vector b ---
A = np.array([[R1 + R3, -R3], [-R3, R2 + R3]])
b = np.array([V1, 0])

# --- 1. クラメルの公式 / Cramer's Rule ---
detA = np.linalg.det(A)
A1 = A.copy(); A1[:, 0] = b
A2 = A.copy(); A2[:, 1] = b
I1_cramer = np.linalg.det(A1) / detA
I2_cramer = np.linalg.det(A2) / detA

# --- 2. NumPyのsolveで解く / Using np.linalg.solve ---
currents_solve = np.linalg.solve(A, b)

# --- 3. 最小二乗法 / Least Squares Method ---
currents_lstsq, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)

# --- 4. 勾配法 / Gradient Descent Method ---
x = np.zeros(2)          # 初期値 / Initial guess
alpha = 0.00001          # 学習率 / Learning rate
for _ in range(10000):   # 反復ステップ / Iteration
    grad = 2 * A.T @ (A @ x - b)  # 勾配 / Gradient
    x -= alpha * grad            # 更新 / Update
currents_gd = x

# --- 結果出力 / Display results ---
print("=== 1. Cramer's Rule ===")
print(f"I1 = {I1_cramer:.4f} A")
print(f"I2 = {I2_cramer:.4f} A\n")

print("=== 2. np.linalg.solve ===")
print(f"I1 = {currents_solve[0]:.4f} A")
print(f"I2 = {currents_solve[1]:.4f} A\n")

print("=== 3. Least Squares ===")
print(f"I1 = {currents_lstsq[0]:.4f} A")
print(f"I2 = {currents_lstsq[1]:.4f} A\n")

print("=== 4. Gradient Descent ===")
print(f"I1 = {currents_gd[0]:.4f} A")
print(f"I2 = {currents_gd[1]:.4f} A\n")

# --- 視覚化（電流計算）/ Visualization of current methods ---
methods = ['Cramer', 'Solve', 'Lstsq', 'GradDesc']
I1_vals = [I1_cramer, currents_solve[0], currents_lstsq[0], currents_gd[0]]
I2_vals = [I2_cramer, currents_solve[1], currents_lstsq[1], currents_gd[1]]

x_pos = np.arange(len(methods))
width = 0.35

plt.figure(figsize=(8, 5))
plt.bar(x_pos - width/2, I1_vals, width, label='I1')
plt.bar(x_pos + width/2, I2_vals, width, label='I2')
plt.xticks(x_pos, methods)
plt.ylabel('Current [A]')
plt.title('Comparison of Current Solutions by Method')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()


```

### 4. 物理定数の計算精度：平方根をニュートン法で求める
```python
import numpy as np
import matplotlib.pyplot as plt

# ================================
# 1. 近似値法 Approximation by trial
# ================================
def approx_sqrt(value):
    guess = 0.1
    while guess * guess < value:
        guess += 0.001
    return round(guess, 3)

# ================================
# 2. 分解法 Decomposition method
# ================================
def decomposition_sqrt(value):
    factor = 0.01
    inside = value / (factor ** 2)
    return factor * np.sqrt(inside)

# ================================
# 3. 加減法的に桁ごとに近づける方法の簡易版
# ================================
def digit_by_digit_sqrt(value, digits=3):
    # 小数点以下 digits 桁まで計算する
    from decimal import Decimal, getcontext
    getcontext().prec = digits + 5
    x = Decimal(value)
    guess = Decimal(0)
    step = Decimal(1)
    for _ in range(digits + 2):
        while (guess + step) ** 2 <= x:
            guess += step
        step /= 10
    return float(round(guess, digits))

# ================================
# 4. ニュートン法 Newton's method
# ================================
def f(x, A): return x**2 - A
def df(x): return 2 * x

def newton_method(A, x0=1.0, tol=1e-10, max_iter=100):
    x = x0
    for _ in range(max_iter):
        fx = f(x, A)
        dfx = df(x)
        if dfx == 0:
            return None
        x_new = x - fx / dfx
        if abs(x_new - x) < tol:
            return x_new
        x = x_new
    return x

# ================================
# 比較と可視化 / Comparison & Visualization
# ================================
target = 0.0145

approx = approx_sqrt(target)
decomp = decomposition_sqrt(target)
digit = digit_by_digit_sqrt(target)
newton = newton_method(target, x0=0.1)
actual = np.sqrt(target)

print("=== Square Root of 0.0145 ===")
print(f"近似値法 Approximation      : {approx}")
print(f"分解法 Decomposition       : {decomp}")
print(f"加減法（桁ごとに近づける）: {digit}")
print(f"ニュートン法 Newton Method : {newton}")
print(f"math.sqrt() 正確な値       : {actual}")

# 可視化 / Plot
x_vals = np.linspace(0, 0.2, 400)
y_vals = f(x_vals, target)

plt.plot(x_vals, y_vals, label=f'f(x) = x² - {target}')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(actual, color='black', linestyle=':', label=f'√{target} (actual)')
plt.axvline(approx, color='blue', linestyle='--', label='Approximation')
plt.axvline(decomp, color='green', linestyle='--', label='Decomposition')
plt.axvline(digit, color='orange', linestyle='--', label='Digit-by-digit')
plt.axvline(newton, color='red', linestyle='--', label='Newton Method')
plt.title(f'Comparison of Square Root Methods (√{target})')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.grid(True)
plt.show()

```





---

### 5. 放物運動と二次方程式の解
```python
import numpy as np
import matplotlib.pyplot as plt
import cmath

# --- 初期条件 --- #
g = 9.8          # 重力加速度 [m/s²]
v0 = 10.0        # 初速度 [m/s]
a = -0.5 * g     # 二次項の係数
b = v0           # 一次項の係数
c = 0            # 定数項（地面の高さと一致）

# --- 解の公式による時間の計算 --- #
D = b**2 - 4*a*c  # 判別式
t1 = (-b + cmath.sqrt(D)) / (2*a)
t2 = (-b - cmath.sqrt(D)) / (2*a)

print(f"判別式 D = {D}")
print(f"落下する時間の解 t1 = {t1:.4f} 秒")
print(f"落下する時間の解 t2 = {t2:.4f} 秒")

# --- 時間範囲の設定（プロット用） --- #
t_vals = np.linspace(0, float(t2.real), 300)
y_vals = v0 * t_vals - 0.5 * g * t_vals**2

# --- グラフ描画 --- #
plt.figure(figsize=(8, 5))
plt.plot(t_vals, y_vals, label='y(t) = v₀t - ½gt²')
plt.axhline(0, color='gray', linestyle='--', linewidth=0.7)
plt.axvline(float(t2.real), color='red', linestyle=':', label=f'落下時間 ≈ {t2.real:.2f} s')
plt.title('Projectile Motion: y(t) = v₀t - ½gt²')
plt.xlabel('Time [s]')
plt.ylabel('Height [m]')
plt.grid(True)
plt.legend()
plt.show()

```

---

### 6. ベクトルの内積・外積と投影
```python
import numpy as np
import matplotlib.pyplot as plt

A = np.array([3, 2])
F = np.array([2, 3])
dot_product = np.dot(A, F)
cross_product = np.cross(A, F)
proj_F_on_A = (np.dot(F, A) / np.dot(A, A)) * A

print(f"内積 A・F = {dot_product}")
print(f"外積 A×F = {cross_product}")
print(f"F を A 方向に投影したベクトル = {proj_F_on_A}")

fig, ax = plt.subplots()
ax.quiver(0, 0, A[0], A[1], angles='xy', scale_units='xy', scale=1, color='blue', label='A')
ax.quiver(0, 0, F[0], F[1], angles='xy', scale_units='xy', scale=1, color='red', label='F')
ax.quiver(0, 0, proj_F_on_A[0], proj_F_on_A[1], angles='xy', scale_units='xy', scale=1, color='green', label='Projection')
ax.set_xlim(-1, 5)
ax.set_ylim(-1, 5)
ax.set_aspect('equal')
ax.set_title('Dot Product and Projection')
ax.grid(True)
ax.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()
```



### 7. 加速度・速度・位置の関係：微積分で運動をとらえる
```python
import numpy as np
import matplotlib.pyplot as plt

# --- 時間設定 --- #
t = np.linspace(0, 10, 1000)
dt = t[1] - t[0]

# --- 力学：重力加速度による運動 --- #
g = 9.8  # 重力加速度 [m/s^2]
a = np.full_like(t, g)  # a(t) = g （一定）

def numerical_integrate(y, dt, initial=0):
    integ = [initial]
    for i in range(1, len(y)):
        area = 0.5 * (y[i] + y[i - 1]) * dt
        integ.append(integ[-1] + area)
    return np.array(integ)

def numerical_derivative(y, dt):
    return np.gradient(y, dt)

v = numerical_integrate(a, dt, initial=0)  # 速度 v(t)
x = numerical_integrate(v, dt, initial=0)  # 位置 x(t)

# --- 電気：電圧波形と微分・積分 --- #
v_elec = np.sin(t)                          # 電圧 v(t) = sin波
i_elec = numerical_derivative(v_elec, dt)   # 電流 i(t) = dv/dt
q_elec = numerical_integrate(i_elec, dt)    # 電荷 q(t) = ∫i dt

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

# 力学
plt.subplot(3, 2, 1)
plt.plot(t, x, label='Position x(t)', color='blue')
plt.ylabel('x(t) [m]')
plt.title('Position (x) under Gravity')
plt.grid(True)
plt.legend()

plt.subplot(3, 2, 3)
plt.plot(t, v, label='Velocity v(t)', color='green')
plt.ylabel('v(t) [m/s]')
plt.title('Velocity (v) under Gravity')
plt.grid(True)
plt.legend()

plt.subplot(3, 2, 5)
plt.plot(t, a, label='Acceleration a(t)', color='red')
plt.ylabel('a(t) [m/s²]')
plt.xlabel('Time [s]')
plt.title('Constant Acceleration (a = g)')
plt.grid(True)
plt.legend()

# 電気
plt.subplot(3, 2, 2)
plt.plot(t, q_elec, label='Charge q(t)', color='purple')
plt.ylabel('q(t) [C]')
plt.title('Charge (q) from Sin Voltage')
plt.grid(True)
plt.legend()

plt.subplot(3, 2, 4)
plt.plot(t, i_elec, label='Current i(t)', color='orange')
plt.ylabel('i(t) [A]')
plt.title('Current (i) = dv/dt')
plt.grid(True)
plt.legend()

plt.subplot(3, 2, 6)
plt.plot(t, v_elec, label='Voltage v(t)', color='brown')
plt.ylabel('v(t) [V]')
plt.xlabel('Time [s]')
plt.title('Voltage (v) = sin(t)')
plt.grid(True)
plt.legend()

plt.suptitle('Kinematics under Constant Gravity vs Electrical Analogues', fontsize=14)
plt.tight_layout()
plt.show()

```

---

### 8. 回転運動と三角関数：単位円でsin, cosを体感
```python

import numpy as np
import matplotlib.pyplot as plt

# θの範囲を定義（0 〜 2π）
theta = np.linspace(0, 2*np.pi, 360)
x = np.cos(theta)
y = np.sin(theta)

# 単位円の描画
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(x, y, label='Unit Circle')
ax.axhline(0, color='gray', linewidth=0.5)
ax.axvline(0, color='gray', linewidth=0.5)

# 三角比の例：60度（π/3ラジアン）
theta_deg = 60
theta_rad = np.deg2rad(theta_deg)
x_val = np.cos(theta_rad)
y_val = np.sin(theta_rad)

# cos成分（実部）、sin成分（虚部）、複素数ベクトルの描画
ax.plot([0, x_val], [0, y_val], color='red', label=r'$e^{i\theta}$ vector')
ax.plot([0, x_val], [0, 0], '--', color='green', label='cos θ (Re)')
ax.plot([x_val, x_val], [0, y_val], '--', color='blue', label='sin θ (Im)')

# 矢印とアノテーション
ax.annotate('Re', xy=(1.1, 0), ha='center')
ax.annotate('Im', xy=(0, 1.1), va='center')
ax.annotate(r'$e^{i\theta}$', xy=(x_val, y_val), xytext=(x_val + 0.1, y_val + 0.1),
            arrowprops=dict(arrowstyle='->', color='red'), color='red')

# フォーマット設定
ax.set_aspect('equal')
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2, 1.2)
ax.set_title("Trigonometric Ratios and Euler's Formula on the Unit Circle")
ax.legend()
plt.grid(True)
plt.show()

```

---

### 9.加法定理と余弦定理の活用
```python
import numpy as np
import matplotlib.pyplot as plt

# 加法定理：sin(x + α) = sin x cos α + cos x sin α
x = np.linspace(0, 2 * np.pi, 500)
alpha = np.pi / 4  # 45度
lhs = np.sin(x + alpha)
rhs = np.sin(x) * np.cos(alpha) + np.cos(x) * np.sin(alpha)

plt.plot(x, lhs, label=r'sin(x + π/4)', linewidth=2)
plt.plot(x, rhs, '--', label=r'sin(x)cos(π/4) + cos(x)sin(π/4)', linewidth=2)
plt.title("Sine Addition Formula")
plt.xlabel("x")
plt.ylabel("Value")
plt.legend()
plt.grid(True)
plt.show()

# 余弦定理の可視化：c² = a² + b² - 2ab cos(C)
a = 5
b = 7
C = np.deg2rad(60)  # 角度60度をラジアンに変換
c = np.sqrt(a**2 + b**2 - 2*a*b*np.cos(C))
print(f"辺cの長さ: {c:.3f}")

# 三角形の座標を描く
A = np.array([0, 0])
B = np.array([b, 0])
C_x = a * np.cos(C)
C_y = a * np.sin(C)
C_point = np.array([C_x, C_y])

fig, ax = plt.subplots()
triangle_x = [A[0], B[0], C_point[0], A[0]]
triangle_y = [A[1], B[1], C_point[1], A[1]]

ax.plot(triangle_x, triangle_y, marker='o')
ax.text(*A, "A", fontsize=12, ha='right')
ax.text(*B, "B", fontsize=12, ha='left')
ax.text(*C_point, "C", fontsize=12, va='bottom')
ax.set_title("Visualization of the Cosine Law")
ax.grid(True)
ax.set_aspect('equal')
plt.show()
```

---

### 10. 和積公式
```python
import numpy as np
import matplotlib.pyplot as plt

# 時間 t を 0 から 2秒まで 500点生成
# Generate time values from 0 to 2 seconds
t = np.linspace(0, 2, 500)

# 角周波数 ω [rad/s] を設定（例：2π rad/s → 周波数1Hz）
# Set angular frequency ω (e.g., 2π rad/s for 1 Hz)
omega = 2 * np.pi

# 波の重ね合わせ sin(ωt) * cos(ωt)
# Wave superposition using product of sine and cosine
product_wave = np.sin(omega * t) * np.cos(omega * t)

# 和積公式で変換：sin(ωt)cos(ωt) = 1/2 sin(2ωt)
# Use trigonometric identity: sin(ωt)cos(ωt) = 1/2 sin(2ωt)
sum_wave = 0.5 * np.sin(2 * omega * t)

# グラフ描画
# Plotting the waves
plt.plot(t, product_wave, label=r'$\sin(\omega t)\cos(\omega t)$')
plt.plot(t, sum_wave, '--', label=r'$\frac{1}{2}\sin(2\omega t)$')

# タイトルとラベルの設定
# Set title and axis labels
plt.title('Wave Superposition in Time Domain via Product-to-Sum Identity')
plt.xlabel('Time t [s]')
plt.ylabel('Amplitude')

# 凡例とグリッド表示
# Show legend and grid
plt.legend()
plt.grid(True)

# グラフ表示
# Display the plot
plt.show()

```

---

### 11. 振幅合成と位相：三角関数合成による波形の変換
```python
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 2*np.pi, 500)
a = 2
b = 3
y = a*np.sin(x) + b*np.cos(x)

R = np.sqrt(a**2 + b**2)
phi = np.arctan2(b, a)
y_composite = R * np.sin(x + phi)

plt.plot(x, y, label='a sin(x) + b cos(x)')
plt.plot(x, y_composite, '--', label='R sin(x + φ)')
plt.title('Sine and Cosine Combination')
plt.xlabel('x')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.show()
```

---

### 12. 物理測定の限界：誤差・有効数字とその取り扱い
```python
# 絶対誤差と相対誤差
true_value = 43.20
measured_value = 43.15

absolute_error = abs(measured_value - true_value)
relative_error = absolute_error / true_value

print(f"絶対誤差: {absolute_error:.5f} g")
print(f"相対誤差: {relative_error*100:.3f} %")
```
```python
# 有効数字に関する物理工学的理解と演算をPythonでシミュレーション
# Understanding and simulation of significant figures using Python

import numpy as np
import math

# 有効数字の丸め関数
# Function to round to specified significant figures
def round_sig(x, sig):
    if x == 0:
        return 0
    return round(x, sig - int(math.floor(math.log10(abs(x)))) - 1)

# 加算・減算：有効数字ではなく「末位をそろえる」処理
# For addition/subtraction, match the decimal place of the least precise value
def add_with_sig_figs(a, b):
    # 小数点以下の桁数の取得
    a_decimals = len(str(a).split(".")[1]) if '.' in str(a) else 0
    b_decimals = len(str(b).split(".")[1]) if '.' in str(b) else 0
    min_decimals = min(a_decimals, b_decimals)
    result = a + b
    return round(result, min_decimals)

# 乗算・除算：有効数字の少ない方に合わせる
# For multiplication/division, match the number of significant figures

def mul_with_sig_figs(a, b):
    a_sig = len(str(a).replace(".", "").lstrip("0"))
    b_sig = len(str(b).replace(".", "").lstrip("0"))
    min_sig = min(a_sig, b_sig)
    return round_sig(a * b, min_sig)

def div_with_sig_figs(a, b):
    a_sig = len(str(a).replace(".", "").lstrip("0"))
    b_sig = len(str(b).replace(".", "").lstrip("0"))
    min_sig = min(a_sig, b_sig)
    return round_sig(a / b, min_sig)

# === テスト用の例 ===
# Example for addition
print("加算（2.34 + 5.6）：", add_with_sig_figs(2.34, 5.6))  # -> 7.9
print("加算（12.345 + 6.7）：", add_with_sig_figs(12.345, 6.7))  # -> 19.0

# Example for multiplication
print("乗算（1.23 × 4.5）：", mul_with_sig_figs(1.23, 4.5))  # -> 5.5
print("乗算（1.234 × 4.56）：", mul_with_sig_figs(1.234, 4.56))  # -> 5.63
print("乗算（1.2345 × 4.567）：", mul_with_sig_figs(1.2345, 4.567))  # -> 5.638

# Example for division
print("除算（4.321 / 1.23）：", div_with_sig_figs(4.321, 1.23))  # -> 約3.51
print("除算（1.23 / 4.321）：", div_with_sig_figs(1.23, 4.321))  # -> 約0.285
```