微分方程式で学ぶ：自由落下・ばね・一次遅れ系の数理モデル高校物理入門


## ■ はじめに

物体の運動や制御系のふるまいを記述する上で、\*\*微分方程式（Differential Equation）\*\*は極めて重要です。本記事では、以下の3つの系に対して、**数理モデルの導出と解法**を順を追って解説します。

1. 自由落下運動
2. 単振動（ばねの運動）
3. 一次遅れ系
---

## 【1】自由落下運動：運動方程式の導出と一般解

### ◆ 運動モデル

鉛直下向きを正とし、空気抵抗を無視した質点の運動を考えると、加速度は一定で**重力加速度 g**になります：

$$
\frac{d^2 x(t)}{dt^2} = g
$$

### ◆ 初期条件付きの一般解

両辺を2回積分すれば、以下のように\*\*位置関数 x(t)\*\*を得ます：

$$
x(t) = x_0 + v_0 t + \frac{1}{2} g t^2
$$

* $x_0$：初期位置
* $v_0$：初速度
* $g$：重力加速度（例：9.8 m/s²）


---

## 【自由落下：位置と速度の時間関数】

### ◆ 位置の式（2回積分）

$$
x(t) = x_0 + v_0 t + \frac{1}{2} g t^2
$$

* $x_0$：初期位置
* $v_0$：初速度
* $g$：重力加速度（例：9.8 m/s²）

---

### ◆ 速度の式（1回積分）

$$
v(t) = \frac{dx(t)}{dt} = v_0 + g t
$$

* 初速度 $v_0$ に、加速度 $g$ が時間に比例して加わる直線的増加（または減少）

---

## 【補足】加速度の式

加速度は定数：

$$
a(t) = \frac{dv(t)}{dt} = g
$$

---

##  まとめ：自由落下の3段階

| 物理量        | 式                                 | 備考           |
| ---------- | --------------------------------- | ------------ |
| 位置 $x(t)$  | $x_0 + v_0 t + \frac{1}{2} g t^2$ | 2階微分の解（積分2回） |
| 速度 $v(t)$  | $v_0 + g t$                       | 1階微分の解（積分1回） |
| 加速度 $a(t)$ | $g$                               | 定数（時間に依存しない） |



---

## 【2】ばねの微分方程式と線形近似（テイラー展開）

### ◆ 単振動の方程式

バネ定数 $k$、質量 $m$ の物体がばねに接続されているとき、ニュートンの運動方程式より：

$$
m \frac{d^2 x(t)}{dt^2} = -k x(t)
\quad \Rightarrow \quad
\frac{d^2 x(t)}{dt^2} + \omega^2 x(t) = 0
$$

ここで、$\omega = \sqrt{\frac{k}{m}}$ は**角振動数**。

### ◆ 一般解と振動

この2階線形常微分方程式の一般解は：

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

初期条件 $x(0), x'(0)$ によって係数 $A, B$ が定まります。



---

## 【ステップ 1】位置の初期条件から $A$ を求める

$$
x(0) = A \cos(0) + B \sin(0) = A \cdot 1 + B \cdot 0 = A
$$

$$
\Rightarrow A = x(0)
$$

---

## 【ステップ 2】速度の初期条件から $B$ を求める

まず、速度関数を求めます：

$$
x'(t) = \frac{d}{dt} \left[ A \cos(\omega t) + B \sin(\omega t) \right] = -A \omega \sin(\omega t) + B \omega \cos(\omega t)
$$

初期時刻 $t = 0$ に代入：

$$
x'(0) = -A \omega \sin(0) + B \omega \cos(0) = 0 + B \omega \cdot 1 = B \omega
$$

$$
\Rightarrow B = \frac{x'(0)}{\omega}
$$

---

##  結論：初期条件からの係数の決定式

$$
\boxed{A = x(0)}, \quad \boxed{B = \frac{x'(0)}{\omega}}
$$

---




### ◆ 線形近似とテイラー展開

例えば、非線形項 $\sin(x)$ を含む場合、$x \approx 0$ のとき：

$$
\sin(x) \approx x - \frac{x^3}{6} + \dots \approx x
$$

として、**線形化**すると振動方程式に帰着できます（振り子の近似など）。

---

## 【3】一次遅れ系の微分方程式とその応答

### ◆ 一次遅れ系とは

時定数 $\tau$ をもつ一次遅れシステムの微分方程式は：

$$
\tau \frac{dy(t)}{dt} + y(t) = Ku(t)
$$

* $y(t)$：出力
* $u(t)$：入力（例：ステップ関数）
* $K$：ゲイン

### ◆ 解法：定数入力 $u(t) = u_0$

定常入力 $u_0$ に対して、初期値 $y(0) = y_0$ の解は：

$$
y(t) = Ku_0 \left(1 - e^{-t/\tau}\right) + y_0 e^{-t/\tau}
$$

この式は、**指数関数的漸近応答**を示し、時間が経つにつれて $Ku_0$ に収束します。

---

## 【まとめ】

| 系       | 微分方程式                                | 一般解                                            |
| ------- | ------------------------------------ | ---------------------------------------------- |
| 自由落下    | $\frac{d^2x}{dt^2} = g$              | $x(t) = x_0 + v_0 t + \frac{1}{2}gt^2$         |
| 単振動（ばね） | $\frac{d^2x}{dt^2} + \omega^2 x = 0$ | $A \cos(\omega t) + B \sin(\omega t)$          |
| 一次遅れ系   | $\tau \frac{dy}{dt} + y = Ku$        | $y(t) = Ku(1 - e^{-t/\tau}) + y_0 e^{-t/\tau}$ |

---


以下に、自由落下運動の**微分方程式をラプラス変換し、伝達関数として表現**する手順を詳述します。

---

## 【1-B】自由落下運動のラプラス変換と伝達関数

### ◆ 微分方程式の再掲

自由落下運動は以下の2階微分方程式で記述されます：

$$
\frac{d^2 x(t)}{dt^2} = g
$$

### ◆ ラプラス変換（初期条件つき）

ラプラス変換を施すと：

$$
\mathcal{L} \left[ \frac{d^2 x(t)}{dt^2} \right] = s^2 X(s) - s x(0) - x'(0)
$$

よって、左辺：

$$
s^2 X(s) - s x_0 - v_0 = \mathcal{L}[g] = \frac{g}{s}
$$

---

### ◆ X(s) の解法

$$
s^2 X(s) = \frac{g}{s} + s x_0 + v_0
\quad \Rightarrow \quad
X(s) = \frac{g}{s^3} + \frac{x_0}{s} + \frac{v_0}{s^2}
$$

---

### ◆ 逆ラプラス変換で元の関数へ

$$
x(t) = x_0 + v_0 t + \frac{1}{2} g t^2
$$

これは、解析解と一致。

---


## 【2-B】ばねの微分方程式と伝達関数（単振動系）

### ◆ 微分方程式の再掲

$$
m \frac{d^2 x(t)}{dt^2} + k x(t) = 0
\quad \Rightarrow \quad
\frac{d^2 x(t)}{dt^2} + \omega^2 x(t) = 0,\quad \omega^2 = \frac{k}{m}
$$

ここでは、ばねに**外力 $f(t)$** を加えた場合を考えます：

$$
m \frac{d^2 x(t)}{dt^2} + k x(t) = f(t)
$$

---

### ◆ ラプラス変換による伝達関数の導出

両辺をラプラス変換（初期値ゼロと仮定）：

$$
m s^2 X(s) + k X(s) = F(s)
\quad \Rightarrow \quad
X(s) \left( m s^2 + k \right) = F(s)
$$

よって、伝達関数 $G(s) = \frac{X(s)}{F(s)}$ は：

$$
G(s) = \frac{1}{m s^2 + k} = \frac{1}{m(s^2 + \omega^2)}
$$

---

### ◆ 特徴

* **2次系（Second-order system）**
* 外力 $f(t)$ を加えたときの**位置応答**を記述
* 振動現象の基本モデル（調和振動子）

---

## 【3-B】一次遅れ系の伝達関数

### ◆ 微分方程式の再掲

$$
\tau \frac{dy(t)}{dt} + y(t) = Ku(t)
$$

---

### ◆ ラプラス変換

ラプラス変換（初期値ゼロ）を施す：

$$
\tau s Y(s) + Y(s) = K U(s)
\quad \Rightarrow \quad
Y(s)(\tau s + 1) = K U(s)
$$

---

### ◆ 伝達関数

$$
G(s) = \frac{Y(s)}{U(s)} = \frac{K}{\tau s + 1}
$$

---

### ◆ 特徴

* **1次系（First-order lag system）**
* ゲイン $K$、時定数 $\tau$
* ステップ応答は：$y(t) = Ku_0 \left(1 - e^{-t/\tau} \right)$
* 応答速度は $\tau$ によって決まる

---

## 📌 まとめ：伝達関数早見表

| システム名        | 微分方程式                   | 伝達関数 $G(s)$                   |
| ------------ | ----------------------- | ----------------------------- |
| 自由落下（加速度→位置） | $\frac{d^2x}{dt^2} = g$ | $\frac{1}{s^2}$               |
| 単振動（力→位置）    | $m\ddot{x} + kx = f(t)$ | $\frac{1}{m(s^2 + \omega^2)}$ |
| 一次遅れ系（入力→出力） | $\tau \dot{y} + y = Ku$ | $\frac{K}{\tau s + 1}$        |

---

```Python
import numpy as np
import matplotlib.pyplot as plt

# -------------------- Parameter Definitions --------------------
# 自由落下（Free Fall）
g = 9.8       # [m/s^2]
x0 = 10.0     # [m]
v0 = 5.0      # [m/s]

# 単振動（Spring Oscillation）
m = 1.0       # [kg]
k = 4.0       # [N/m]
omega = np.sqrt(k / m)        # [rad/s]
f = omega / (2 * np.pi)       # [Hz]
T = 1 / f                     # [s]
A = 2.0
B = 0.0

# 一次遅れ系（First-order lag system）
tau = 1.5     # [s]
K = 1.0       # unit gain
u0 = 1.0      # step input
y0 = 0.0

# Time vector
t = np.linspace(0, 10, 1000)

# -------------------- Equations --------------------
# 自由落下
x_free = x0 + v0 * t + 0.5 * g * t**2
v_free = v0 + g * t
a_free = np.full_like(t, g)

# 単振動
x_spring = A * np.cos(omega * t) + B * np.sin(omega * t)
v_spring = -A * omega * np.sin(omega * t) + B * omega * np.cos(omega * t)
a_spring = -omega**2 * x_spring

# 一次遅れ系
y_first_order = K * u0 * (1 - np.exp(-t / tau)) + y0 * np.exp(-t / tau)

# -------------------- Plotting --------------------
fig, axs = plt.subplots(3, 2, figsize=(12, 12))

# 自由落下
axs[0, 0].plot(t, x_free)
axs[0, 0].set_title("Free Fall: Position vs Time")
axs[0, 0].set_xlabel("Time [s]")
axs[0, 0].set_ylabel("Position [m]")
axs[0, 0].grid(True)

axs[1, 0].plot(t, v_free)
axs[1, 0].set_title("Free Fall: Velocity vs Time")
axs[1, 0].set_xlabel("Time [s]")
axs[1, 0].set_ylabel("Velocity [m/s]")
axs[1, 0].grid(True)

axs[2, 0].plot(t, a_free)
axs[2, 0].set_title("Free Fall: Acceleration vs Time")
axs[2, 0].set_xlabel("Time [s]")
axs[2, 0].set_ylabel("Acceleration [m/s^2]")
axs[2, 0].grid(True)

# 単振動と一次遅れ系
axs[0, 1].plot(t, x_spring)
axs[0, 1].set_title(f"Spring Oscillation: Position vs Time\nFrequency = {f:.2f} Hz, Period = {T:.2f} s")
axs[0, 1].set_xlabel("Time [s]")
axs[0, 1].set_ylabel("Position [m]")
axs[0, 1].grid(True)

axs[1, 1].plot(t, v_spring)
axs[1, 1].set_title("Spring Oscillation: Velocity vs Time")
axs[1, 1].set_xlabel("Time [s]")
axs[1, 1].set_ylabel("Velocity [m/s]")
axs[1, 1].grid(True)

axs[2, 1].plot(t, y_first_order)
axs[2, 1].set_title("First-order Lag System: Output vs Time")
axs[2, 1].set_xlabel("Time [s]")
axs[2, 1].set_ylabel("Output")
axs[2, 1].grid(True)

plt.tight_layout()
plt.show()

```