0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

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

Last updated at Posted at 2025-03-28

はじめに

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

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

Physics Math Curriculum

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


1. 物理量のスケール感と単位系:SI接頭語とその変換

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. 測定と物理量の整合性:次元解析と有効数字の基礎

# 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))
# 有効数字の丸めと評価
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))
# 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. 電流の流れを解き明かす:キルヒホッフ則と行列解法

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. 物理定数の計算精度:平方根をニュートン法で求める

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. 放物運動と二次方程式の解

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. ベクトルの内積・外積と投影

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. 加速度・速度・位置の関係:微積分で運動をとらえる

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を体感


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.加法定理と余弦定理の活用

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. 和積公式

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. 振幅合成と位相:三角関数合成による波形の変換

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. 物理測定の限界:誤差・有効数字とその取り扱い

# 絶対誤差と相対誤差
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でシミュレーション
# 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
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?