高校物理の単振動の式 PythonとMATLAB

参考にした記事

import numpy as np
import matplotlib.pyplot as plt

# Define constants
m = 1.0        # Mass (kg)
k = 1.0        # Spring constant (N/m)
omega = np.sqrt(k / m)  # Angular frequency (rad/s)
period = 2 * np.pi / omega  # Period (s)

# Initial conditions
x0 = 1.0       # Initial displacement (m)
v0 = 0.0       # Initial velocity (m/s)

# Time array
t = np.linspace(0, 2 * period, 500)  # Simulate for two periods

# Equations of motion
position = x0 * np.cos(omega * t) + (v0 / omega) * np.sin(omega * t)  # Position
velocity = -x0 * omega * np.sin(omega * t) + v0 * np.cos(omega * t)   # Velocity
acceleration = -x0 * omega**2 * np.cos(omega * t) - v0 * omega * np.sin(omega * t)  # Acceleration

# Plot results
plt.figure(figsize=(10, 8))

# Position plot
plt.subplot(3, 1, 1)
plt.plot(t, position, color='blue', label='Position (X)')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.title(f'Harmonic Oscillator: Position, Velocity, and Acceleration\nPeriod (T) = {period:.4f} s')
plt.ylabel('Position (X) [m]')
plt.legend()
plt.grid()

# Velocity plot
plt.subplot(3, 1, 2)
plt.plot(t, velocity, color='green', label='Velocity (dX/dt)')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.ylabel('Velocity (dX/dt) [m/s]')
plt.legend()
plt.grid()

# Acceleration plot
plt.subplot(3, 1, 3)
plt.plot(t, acceleration, color='red', label='Acceleration (d²X/dt²)')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.xlabel('Time (t) [s]')
plt.ylabel('Acceleration (d²X/dt²) [m/s²]')
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

function [position, velocity, acceleration, period] = harmonic_oscillator(t, x0, v0, m, k)
% Harmonic Oscillator Simulation
% Inputs:
% t   - Time vector (s)
% x0  - Initial displacement (m)
% v0  - Initial velocity (m/s)
% m   - Mass (kg)
% k   - Spring constant (N/m)
%
% Outputs:
% position     - Displacement over time (m)
% velocity     - Velocity over time (m/s)
% acceleration - Acceleration over time (m/s^2)
% period       - Period of oscillation (s)

    % Calculate angular frequency and period
    omega = sqrt(k / m);           % Angular frequency (rad/s)
    period = 2 * pi / omega;       % Oscillation period (s)
    
    % Calculate position, velocity, and acceleration
    position = x0 * cos(omega * t) + (v0 / omega) * sin(omega * t);
    velocity = -x0 * omega * sin(omega * t) + v0 * cos(omega * t);
    acceleration = -x0 * omega^2 * cos(omega * t) - v0 * omega * sin(omega * t);
end

参考文献
https://www.akita-pu.ac.jp/reccs/umiyamoto/high.pdf

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
l = 1.0  # 振り子の長さ (m)
g = 9.8  # 重力加速度 (m/s^2)
v0 = 1.0  # 初期速度 (m/s)

omega = np.sqrt(g / l)  # 角振動数
period = 2 * np.pi / omega  # 振動の周期

# 時間の設定
t = np.linspace(0, 2 * period, 500)  # 2周期分の時間

# 位置の計算
position = (v0 / omega) * np.sin(omega * t)

# 一回微分 (速度)
velocity = v0 * np.cos(omega * t)

# 二階微分 (加速度)
acceleration = -v0 * omega * np.sin(omega * t)

# 周期を表示
print(f"Angular frequency (omega): {omega:.4f} rad/s")
print(f"Period of oscillation (T): {period:.4f} s")

# プロット
plt.figure(figsize=(10, 8))

# 位置のプロット
plt.subplot(3, 1, 1)
plt.plot(t, position, label="Position $x(t)$", color="blue")
plt.axhline(0, color="black", linestyle="--", linewidth=0.5)
plt.ylabel("Position [m]")
plt.legend()
plt.grid()

# 速度のプロット
plt.subplot(3, 1, 2)
plt.plot(t, velocity, label="Velocity $v(t)$", color="green")
plt.axhline(0, color="black", linestyle="--", linewidth=0.5)
plt.ylabel("Velocity [m/s]")
plt.legend()
plt.grid()

plt.subplot(3, 1, 3)
plt.plot(t, acceleration, label="Acceleration $a(t)$", color="red")
plt.axhline(0, color="black", linestyle="--", linewidth=0.5)
plt.xlabel("Time [s]")
plt.ylabel("Acceleration [m/s²]")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

おまけ、導体のレール型微分方程式
参考文献
https://okazaki-h.aichi-c.ed.jp/cms/wp-content/uploads/2023/02/iAbutsuri2_kagakuhousokutobibun_text.pdf

import numpy as np
import matplotlib.pyplot as plt

# 定数の設定
v0 = 5.0   # 初速度 (m/s)
B = 0.1    # 磁束密度 (T)
l = 1.0    # 導体棒の長さ (m)
m = 1.0    # 導体棒の質量 (kg)
R = 2.0    # 抵抗 (Ω)

# 時間配列の作成
t = np.linspace(0, 2000, 50000)  
# 電流の計算
i_t = (v0 * B * l / R) * np.exp(-B**2 * l**2 * t / (m * R))

# 電力の計算
p_t = (v0**2 * B**2 * l**2 / R) * np.exp(-2 * B**2 * l**2 * t / (m * R))

# プロット
plt.figure(figsize=(12, 8))

# 電流プロット
plt.subplot(2, 1, 1)
plt.plot(t, i_t, label=r'Current $i(t)$', color='blue')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.title("Electric Current and Power over Time")
plt.xlabel("Time (s)")
plt.ylabel("Current (A)")
plt.legend()
plt.grid()

# 電力プロット
plt.subplot(2, 1, 2)
plt.plot(t, p_t, label=r'Power $p(t)$', color='red')
plt.axhline(0, color='black', linestyle='--', linewidth=0.5)
plt.xlabel("Time (s)")
plt.ylabel("Power (W)")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()