Pythonで理解する物理のエッセンス 力学・波動 (河合塾シリーズ)

Anaconda のインストール
https://www.python.jp/install/anaconda/windows/install.html

解説動画が豊富

加法定理の逆

import numpy as np
import matplotlib.pyplot as plt

# aとbの値を指定
a = 3
b = 4

# √(a^2+b^2)を計算
magnitude = np.sqrt(a**2 + b**2)
print("Magnitude:", magnitude)

# arctan(b/a)を計算（弧度法と度数法）
theta_rad = np.arctan2(b, a)
theta_deg = np.degrees(theta_rad)
print("Angle (radians):", theta_rad)
print("Angle (degrees):", theta_deg)

# 関数のプロット
x = np.linspace(0, 2*np.pi, 100)
y = a * np.sin(x) + b * np.cos(x)
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('a*sin(x) + b*cos(x)')
plt.title('Plot of a*sin(x) + b*cos(x)')
plt.grid(True)
plt.show()

p8〜10
運動方程式
x(t)をtで微分するとv(t)となる
v(t)をtで微分するとa(t)となる

import numpy as np
import matplotlib.pyplot as plt

# 関数定義
def x(t, initial_velocity, acceleration):
    return initial_velocity * t + 0.5 * acceleration * (t**2)

# 時間の範囲を設定
t_values = np.linspace(0, 10, 100)  # 0から10の範囲を100分割

# 初速度と加速度を指定
initial_velocity = 5
acceleration = 9.8

# x(t)の計算
position = x(t_values, initial_velocity, acceleration)

# x(t)のプロット
plt.figure(figsize=(10, 5))
plt.plot(t_values, position, label='x(t)')
plt.xlabel('Time (t)')
plt.ylabel('Position (x)')
plt.title('Position vs Time')
plt.legend()
plt.grid(True)
plt.show()

# x(t)をtで微分した結果をプロット
velocity = initial_velocity + acceleration * t_values
plt.figure(figsize=(10, 5))
plt.plot(t_values, velocity, label='Velocity')
plt.xlabel('Time (t)')
plt.ylabel('Velocity')
plt.title('Velocity vs Time')
plt.legend()
plt.grid(True)
plt.show()

# x(t)をtで2回微分した結果をプロット
acceleration_values = np.full_like(t_values, acceleration)
plt.figure(figsize=(10, 5))
plt.plot(t_values, acceleration_values, label='Acceleration')
plt.xlabel('Time (t)')
plt.ylabel('Acceleration')
plt.title('Acceleration vs Time')
plt.legend()
plt.grid(True)
plt.show()

結果

p5問題5

import numpy as np
import matplotlib.pyplot as plt

# Parameters
g = 9.8  # 加速度
v0 = 10  # 初速度
H = 5    # 高さ

# 解の公式を使ってtを計算する
discriminant = np.sqrt(v0**2 + 8*g*H)
t1 = (v0 + discriminant) / (2*g)
t2 = (v0 - discriminant) / (2*g)

# tの範囲を決定
t_min = 0
t_max = max(t1, t2) * 2  # 描写のため、少し余裕をもたせる

# tの値の範囲を生成
t_values = np.linspace(t_min, t_max, 100)

# 放物線の式
y_values = g * t_values**2 - v0 * t_values - 2*H

# グラフ描写
plt.plot(t_values, y_values, label='gt^2 - v0t - 2H')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Height (m)')
plt.title('Parabolic Trajectory')
plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
plt.legend()
plt.show()

結果

p13
y(t)＝-（g/((2(v0)^2)(cosθ)^2)）x^2+(tanθ)x

import numpy as np
import matplotlib.pyplot as plt

# Parameters
g = 9.8  # 加速度
v0 = 10  # 初速度
theta = np.radians(45)  # 角度（ラジアンに変換）

# xの範囲を決定
x_min = 0
x_max = 2 * v0**2 * np.cos(theta)**2 / g  # 最大のx値を計算

# xの値の範囲を生成
x_values = np.linspace(x_min, x_max, 100)

# yの値を計算
y_values = -(g / (2 * (v0**2) * (np.cos(theta)**2))) * x_values**2 + np.tan(theta) * x_values

# グラフ描写
plt.plot(x_values, y_values, label=r'$y(t) = -\frac{g}{2(v_0)^2(\cos\theta)^2} x^2 + \tan\theta \cdot x$')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Plot of y(x)')
plt.grid(color = 'gray', linestyle = '--', linewidth = 0.5)
plt.legend()
plt.show()

結果
θ=45度のとき

連立方程式と行列とクラメルの公式

p19

import matplotlib.pyplot as plt
import numpy as np

# 重力加速度
g = 9.81  # m/s^2

# 角度 (度数法)
theta_deg = 30

# 角度をラジアンに変換
theta_rad = np.deg2rad(theta_deg)

# FとNの計算
F = g * np.sin(theta_rad)
N = g * np.cos(theta_rad)

# 原点から(N,F)のベクトル
vector_NF = np.array([N, F])

# グラフの描画
plt.figure(figsize=(6, 6))
plt.quiver(0, 0, vector_NF[0], vector_NF[1], angles='xy', scale_units='xy', scale=1, color=['r', 'b'])
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.xlabel('N')
plt.ylabel('F')
plt.title('Vector from Origin to (N,F)')
plt.grid()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

結果

オイラーの公式を用いた表現

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
m = 5  # 任意の実数
g = 9.8  # 重力加速度
theta = np.linspace(0, 2*np.pi, 100)  # 0から2πまでの角度

# オイラーの公式を使って複素数を生成
z = m * g * np.exp(1j * theta)

# 実部と虚部を取得
x = np.real(z)
y = np.imag(z)

# プロット
plt.figure(figsize=(8, 6))
plt.plot(x, y)
plt.title("Euler's Formula: m * g * e^(j*theta)")
plt.xlabel("Real")
plt.ylabel("Imaginary")
plt.grid()
plt.axis('equal')  # アスペクト比を保つ
plt.show()

結果

P20

import numpy as np

# 係数行列
coefficients = np.array([[np.sqrt(3)/2, -1/np.sqrt(2)],
                         [1/2, 1/np.sqrt(2)]])
# 定数項ベクトル
mg = 9.8  # 重力加速度
constants = np.array([0, mg])

# 連立方程式を解く
solution = np.linalg.solve(coefficients, constants)
x = solution[0]
y = solution[1]

print("x =", x)
print("y =", y)

結果

p25
並列のバネ

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# k1とk2の範囲を設定します
k1_range = np.linspace(0, 10, 100)  # 0から10までの範囲で100点を生成
k2_range = np.linspace(0, 10, 100)  # 0から10までの範囲で100点を生成

# k1とk2のそれぞれの値に対して計算を行います
k1, k2 = np.meshgrid(k1_range, k2_range)
result = np.zeros_like(k1)
mask = (k1 + k2) != 0
result[mask] = (k1[mask] * k2[mask]) / (k1[mask] + k2[mask])

# プロットします
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# (k1, k2, (k1*k2) / (k1+k2)) の3次元プロットを作成します
ax.plot_surface(k1, k2, result, cmap='viridis')

# 軸ラベルを設定します
ax.set_xlabel('k1')
ax.set_ylabel('k2')
ax.set_zlabel('(k1*k2) / (k1+k2)')

plt.show()

結果

p30
mとMは0以上
x(m,M)＝((m+2M)l)/(2(2m+M))

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def x(m, M, l):
    return ((m + 2*M) * l) / (2 * (2*m + M))

# m, M の範囲を設定
m_values = np.linspace(0, 10, 100)
M_values = np.linspace(0, 10, 100)

# m_values, M_values のすべての組み合わせに対して x を計算
X, Y = np.meshgrid(m_values, M_values)
Z = x(X, Y, 10)  # l = 10 とします

# 3Dプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')

# 軸ラベルの設定
ax.set_xlabel('m')
ax.set_ylabel('M')
ax.set_zlabel('x')

plt.title('3D Plot of x(m, M)')
plt.show()

結果

p78　マスバネモデルの微分方程式

import numpy as np
import matplotlib.pyplot as plt

def x(t, c1, c2, omega):
    return c2 * np.sin(np.sqrt(omega) * t) + c1 * np.cos(np.sqrt(omega) * t)

# パラメータの設定
omega = 1  # 任意の角周波数
t = np.linspace(0, 10, 1000)

# c1が0の場合のプロット
c1_0 = 0
c2_0 = 1  # ここではc2を1としていますが、任意の値を指定できます
y_c1_0 = x(t, c1_0, c2_0, omega)
plt.plot(t, y_c1_0, label='c1=0, c2=1')

# c2が0の場合のプロット
c1_1 = 1  # ここではc1を1としていますが、任意の値を指定できます
c2_1 = 0
y_c2_0 = x(t, c1_1, c2_1, omega)
plt.plot(t, y_c2_0, label='c1=1, c2=0')

# グラフの設定
plt.xlabel('Time (t)')
plt.ylabel('x(t)')
plt.title('Plot of x(t)')
plt.grid(True)
plt.legend()
plt.show()

結果

import numpy as np
import matplotlib.pyplot as plt

# 関数 x(t) = A * sin(omega * t + B)
def x(t, A, omega, B):
    return A * np.sin(omega * t + B)

# 導関数 x'(t)
def dx(t, A, omega, B):
    return A * omega * np.cos(omega * t + B)

# 二階導関数 x''(t)
def d2x(t, A, omega, B):
    return -A * omega**2 * np.sin(omega * t + B)

# パラメータ設定
A = 1
omega = 1
B = 0

# 時間の範囲を設定
t = np.linspace(0, 10, 1000)

# x(t)、x'(t)、x''(t)の計算
xt = x(t, A, omega, B)
xt_prime = dx(t, A, omega, B)
xt_double_prime = d2x(t, A, omega, B)

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

plt.subplot(3, 1, 1)
plt.plot(t, xt, label='$x(t)$')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, xt_prime, label="$x'(t)$")
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(t, xt_double_prime, label="$x''(t)$")
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

結果

p85

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
ω = 1.0
λ = 0.1

# 時間の範囲
t = np.linspace(0, 10, 1000)

# 関数
x = np.cos(ω * t) * np.exp(-λ * t)

# 極大値と極小値のインデックスを見つける
max_indices = np.where((np.diff(np.sign(np.diff(x))) < 0) & (np.diff(np.sign(np.diff(x))) != 0))[0] + 1
min_indices = np.where((np.diff(np.sign(np.diff(x))) > 0) & (np.diff(np.sign(np.diff(x))) != 0))[0] + 1

# プロット
plt.plot(t, x, label='x = (cos(ωt))e^(-λt)')
plt.scatter(t[max_indices], x[max_indices], color='red', label='極大値')
plt.scatter(t[min_indices], x[min_indices], color='green', label='極小値')
plt.xlabel('t')
plt.ylabel('x')
plt.legend()
plt.grid(True)
plt.show()

結果　2π／ωずつ　極小値、極大値を繰り返す

p92
楕円運動の力学的エネルギー保存則の式

def kinetic_energy(m, v):
    return 0.5 * m * v**2

def potential_energy(G, M, m, r):
    return - G * M * m / r

def total_energy(kinetic, potential):
    return kinetic + potential

# パラメータ
m = 1.0  # 質点の質量
v = 10.0  # 速度
G = 6.67430e-11  # 万有引力定数
M = 5.972e24  # 中心天体の質量
r = 6371000.0  # 質点の中心天体からの距離（地球の半径）

# 運動エネルギーとポテンシャルエネルギーを計算
T = kinetic_energy(m, v)
U = potential_energy(G, M, m, r)

# 合計エネルギーを計算
E_total = total_energy(T, U)

print("運動エネルギー:", T)
print("ポテンシャルエネルギー:", U)
print("合計エネルギー:", E_total)

結果

p104　波の式

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータの設定
A = 1  # 振幅
T = 1  # 時間周期
λ = 1  # 波長

# x, t の範囲を設定
x = np.linspace(-5, 5, 100)
t = np.linspace(0, 10, 100)

# x, t のメッシュグリッドを作成
X, T = np.meshgrid(x, t)

# y(x,t)の計算
Y = A * np.sin((T/T) - (X/λ))

# 3次元プロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, T, Y, cmap='viridis')

ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('y(x,t)')

plt.show()

結果

p106
定常波のプロット 2つの波

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# パラメータ設定
x = np.linspace(0, 2*np.pi, 100)
y = np.linspace(0, 2*np.pi, 100)
X, Y = np.meshgrid(x, y)

# 2つの波の定義
wave1 = np.sin(X) * np.cos(Y)
wave2 = np.sin(X) * np.sin(Y)

# 波のプロット
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 波1のプロット
ax.plot_surface(X, Y, wave1, cmap='viridis')

# 波2のプロット
ax.plot_surface(X, Y, wave2, cmap='plasma')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Amplitude')

plt.show()

結果

p135
近似式

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 関数の定義
def func1(x, a):
    return (1 + x) ** a

def func2(x, a):
    return (1 + a * x)

def func3(x):
    return np.sin(x)

def func4(x):
    return np.cos(x)

def func5(x):
    return np.log(1 + x)

# データの生成
x = np.linspace(-5, 5, 100)
a = np.linspace(-2, 2, 100)
X, A = np.meshgrid(x, a)
Y1 = func1(X, A)
Y2 = func2(X, A)
Y3 = func3(x)
Y4 = func4(x)
Y5 = func5(x)

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

ax1 = fig.add_subplot(231, projection='3d')
ax1.plot_surface(X, A, Y1, cmap='viridis')
ax1.set_title('(1+x)^a vs (1+ax)')

ax2 = fig.add_subplot(232, projection='3d')
ax2.plot_surface(X, A, Y2, cmap='viridis')
ax2.set_title('(1-x)^a vs (1-ax)')

ax3 = fig.add_subplot(233)
ax3.plot(x, Y3, label='sin(x)')
ax3.plot(x, x, label='x')
ax3.legend()
ax3.set_title('sin(x) vs x')

ax4 = fig.add_subplot(234)
ax4.plot(x, Y4, label='cos(x)')
ax4.plot(x, np.ones_like(x), label='1')
ax4.legend()
ax4.set_title('cos(x) vs 1')

ax5 = fig.add_subplot(235)
ax5.plot(x, Y5, label='log(1+x)')
ax5.plot(x, x, label='x')
ax5.legend()
ax5.set_title('log(1+x) vs x')

plt.show()

結果

p134

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def intensity(x, y, d, wavelength, D, I0):
    r = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    phi = (np.pi * d / wavelength) * np.sin(theta)
    return I0 * (np.cos(phi))**2 * (np.sinc((d / wavelength) * np.sin(theta)))**2 * (D / (D**2 + r**2))

# パラメータの設定
d = 0.1  # スリット間の距離 (m)
wavelength = 0.0005  # 光の波長 (m)
D = 1  # スリットからスクリーンまでの距離 (m)
I0 = 1  # 中央の明るさ

# スクリーンの座標を生成
x = np.linspace(-0.02, 0.02, 100)
y = np.linspace(-0.02, 0.02, 100)
X, Y = np.meshgrid(x, y)

# 光の強度を計算
I = intensity(X, Y, d, wavelength, D, I0)

# 3Dプロット
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, I, cmap='hot', rstride=5, cstride=5)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_zlabel('Intensity')
ax.set_title('Intensity Distribution on the Screen')
plt.show()

結果

import matplotlib.pyplot as plt

# パラメーターの値を設定
l = 5  # m
d = 0.001  # m
x = 0.02  # m

# 各点の座標を計算
points = [
    (0, 0),
    (l, 0),
    (l, d/2),
    (l, d/2 + x),
    (0, d)
]

# プロットする点の座標を抽出
x_values = [point[0] for point in points]
y_values = [point[1] for point in points]

# プロット
plt.plot(x_values, y_values, 'ro')  # 'ro' は赤い点を意味する
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Points Plot')
plt.grid(True)
plt.show()

結果

p161 万有引力の位置エネルギー

import numpy as np
import matplotlib.pyplot as plt

def gravitational_force(m1, m2, r):
    """
    二つの物体間の万有引力を計算する関数。

    Parameters:
        m1 (float): 一つ目の物体の質量（キログラム単位）。
        m2 (float): 二つ目の物体の質量（キログラム単位）。
        r (float or array-like): 二つの物体間の距離（メートル単位）。

    Returns:
        float or array-like: 万有引力（ニュートン単位）。
    """
    G = 6.67430e-11  # 万有引力定数（m^3/kg/s^2）
    force = G * m1 * m2 / (r ** 2)
    return force

# 距離の範囲を定義
distances = np.linspace(1e9, 2e11, 100)  # 100万kmから2億kmまでの範囲

# 各距離での引力を計算
mass_earth = 5.972e24  # 地球の質量（キログラム単位）
mass_sun = 1.989e30  # 太陽の質量（キログラム単位）
forces = gravitational_force(mass_earth, mass_sun, distances)

# グラフ描画
plt.figure(figsize=(10, 6))
plt.plot(distances, forces, color='blue')
plt.title('太陽と地球の間の万有引力')
plt.xlabel('距離（メートル）')
plt.ylabel('万有引力（ニュートン）')
plt.grid(True)
plt.show()

結果