Pythonで学ぶ新応用数学改訂版(高専用大日本図書)1章 ベクトル解析

1・1 空間のベクトル

空間の基本ベクトルを i, j, k で表す。

ベクトル a を原点Oを始点として表したときの終点の座標を
$(a_x, a_y, a_z)$ とすると、aは次のように基本ベクトルの線形結合で表すことができる。

$$
a = a_x i + a_y j + a_z k
$$

ここで $a_x, a_y, a_z$ をそれぞれ x成分, y成分, z成分 といい、これをベクトルの成分表示という。

$$
a = (a_x, a_y, a_z)
$$

---

工学応用例

• 力ベクトル
  $\mathbf{F} = (F_x, F_y, F_z)$ [N]
  各成分は x, y, z方向の力。
• 電場ベクトル
  $\mathbf{E} = (E_x, E_y, E_z)$ [V/m]
  空間の電界の分布を表す。

---

大きさ（ノルム）

$$
|a| = \sqrt{a_x^2 + a_y^2 + a_z^2}
$$

---

工学応用例

• 速度ベクトルの大きさ = 速さ
  $|\mathbf{v}| = \sqrt{v_x^2 + v_y^2 + v_z^2}$ [m/s]
• 電流密度ベクトルの大きさ
  $|\mathbf{J}| = \sqrt{J_x^2 + J_y^2 + J_z^2}$ [A/m²]

---

ベクトルの和・差

$$
a \pm b = (a_x \pm b_x, ; a_y \pm b_y, ; a_z \pm b_z)
$$

---

工学応用例

• 並列に作用する力の合成
  $\mathbf{F}_{合} = \mathbf{F}_1 + \mathbf{F}_2$ [N]
• 電圧ベクトルの重ね合わせ
  $\mathbf{E}_{合} = \mathbf{E}_1 + \mathbf{E}_2$ [V/m]

---

スカラー倍

$$
ma = (ma_x, ma_y, ma_z) \quad (mは実数)
$$

---

工学応用例

• 電圧や力のスケーリング
  例えば、ギア比 m によるトルク変換：
  $\mathbf{T}{出力} = m \mathbf{T}{入力}$ [N·m]

---

位置ベクトル

点 P(x, y, z) の位置ベクトルは

$$
r = \overrightarrow{OP} = (x, y, z) \quad [m]
$$

---

工学応用例

• ロボット工学
  アーム先端の位置を $(x,y,z)$ [m] で表す。

• 電磁気学
  位置ベクトル r [m] でクーロンの法則

  $$
  \mathbf{E} = \frac{1}{4\pi\varepsilon_0}\frac{q}{|\mathbf{r}|^2}\hat{\mathbf{r}} \quad [V/m]
  $$

---

単位ベクトル

大きさ1のベクトルを単位ベクトルという。

$$
e = \pm \frac{a}{|a|}
$$

---

工学応用例

• 方向の指示
  力の方向だけを表す場合、単位ベクトル $\hat{F}$ を使う。
• 電場の方向
  $\mathbf{E} = |\mathbf{E}| \hat{e}$ と分解し、大きさ[V/m]と方向を分離。

---

内積（スカラー積）

2つのベクトルa, bのなす角をθとするとき、

$$
a \cdot b = |a||b|\cos \theta
$$

成分表示：

$$
a \cdot b = a_x b_x + a_y b_y + a_z b_z
$$

---

工学応用例

• 力と変位の仕事

  $$
  W = \mathbf{F} \cdot \Delta \mathbf{r} \quad [J = N \cdot m]
  $$

• 相関（信号処理）
  2つの信号ベクトルの内積 = 相関度（単位なし、正規化後）。

---

外積（ベクトル積）

$$
|a \times b| = |a||b|\sin \theta
$$

方向は右ねじの法則。

成分表示：

$$
a \times b = (a_y b_z - a_z b_y,; a_z b_x - a_x b_z,; a_x b_y - a_y b_x)
$$

---

工学応用例

• トルク

  $$
  \mathbf{T} = \mathbf{r} \times \mathbf{F} \quad [N \cdot m]
  $$

• ローレンツ力の磁場成分

  $$
  \mathbf{F}_B = q(\mathbf{v} \times \mathbf{B}) \quad [N]
  $$

  （q: 電荷[C], v: 速度[m/s], B: 磁束密度[T = N·s/(C·m)])

---

Python実装例

import numpy as np

# ベクトルの定義
a = np.array([1, 2, 3])  # ベクトル a = (1, 2, 3)
b = np.array([4, 5, 6])  # ベクトル b = (4, 5, 6)

# 1. ベクトルの大きさ（ノルム）
norm_a = np.linalg.norm(a)
print(f"ベクトル a の大きさ: {norm_a}")  # √(1² + 2² + 3²) = √14

# 2. ベクトルの和・差
sum_ab = a + b
diff_ab = a - b
print(f"ベクトルの和 a + b: {sum_ab}")  # (1+4, 2+5, 3+6) = (5, 7, 9)
print(f"ベクトルの差 a - b: {diff_ab}")  # (1-4, 2-5, 3-6) = (-3, -3, -3)

# 3. スカラー倍
m = 2
scalar_mult = m * a
print(f"スカラー倍 2a: {scalar_mult}")  # (2*1, 2*2, 2*3) = (2, 4, 6)

# 4. 単位ベクトル
unit_a = a / norm_a
print(f"ベクトル a の単位ベクトル: {unit_a}")  # a / |a|

# 5. 内積
dot_product = np.dot(a, b)
print(f"内積 a・b: {dot_product}")  # 1*4 + 2*5 + 3*6 = 32

# 6. 外積
cross_product = np.cross(a, b)
print(f"外積 a×b: {cross_product}")  # (a_y*b_z - a_z*b_y, a_z*b_x - a_x*b_z, a_x*b_y - a_y*b_x)

# 7. 工学応用例: 力と変位の仕事
F = np.array([10, 0, 0])  # 力ベクトル [N]（x方向に10N）
r = np.array([2, 3, 0])   # 変位ベクトル [m]
work = np.dot(F, r)
print(f"仕事 W = F・r: {work} [J]")  # 10*2 + 0*3 + 0*0 = 20 J

# 8. 工学応用例: トルク
r_arm = np.array([0, 0, 1])  # 力の作用点 [m]（z軸上）
F = np.array([10, 0, 0])     # 力ベクトル [N]（x方向）
torque = np.cross(r_arm, F)
print(f"トルク T = r×F: {torque} [N・m]")  # (0, -10, 0)

外積の導出

• $a, b$ を原点 $O$ を始点として描き，$a=\overrightarrow{OA},, b=\overrightarrow{OB}$ とおく。
  点 $A, B$ の $xy$ 平面への正射影をそれぞれ $A'(a_x,a_y,0),,B'(b_x,b_y,0)$ とし，$\overrightarrow{OA},\overrightarrow{OB}$ が $x$ 軸の正方向となす角をそれぞれ $\alpha,\beta$ とする。
  また $a$ と $b$ を軸の正の方向へ回したときの角をそれぞれ $\alpha,\beta$ とする（図は $0<\gamma<\tfrac{\pi}{2},\ 0<\beta-\alpha<\pi$ の場合）。

• このとき $a\times b$ の $z$ 成分 は

  $$
  |a\times b|\cos\gamma
  =2\triangle OAB\cdot \cos\gamma
  =2\triangle OA'B'
  =OA'\cdot OB'\sin(\beta-\alpha)
  $$

  $$
  =OA'\cdot OB'(\sin\beta\cos\alpha-\cos\beta\sin\alpha)
  $$

  $$
  =OA'_x\cdot OB'_y - OA'_y\cdot OB'_x
  =a_x b_y - a_y b_x .
  $$

• 他の軸でも同様に求め，成分をまとめると（外積の成分表示）：

  $$
  \boxed{,a\times b
  =(a_y b_z-a_z b_y),{\bf i}
  +(a_z b_x-a_x b_z),{\bf j}
  +(a_x b_y-a_y b_x),{\bf k},} \tag{4}
  $$

• （注）外積は行列式の展開公式で

  $$
  a\times b=
  \begin{vmatrix}
  {\bf i}&{\bf j}&{\bf k}\
  a_x&a_y&a_z\
  b_x&b_y&b_z
  \end{vmatrix}
  $$

  と書ける（行列式は形式的表現）。

• 例4　$\overrightarrow{OA}=(2,2,1),\ \overrightarrow{OB}=(1,2,3)$ のとき

  $$
  \overrightarrow{OA}\times\overrightarrow{OB}

  \begin{vmatrix}
  {\bf i}&{\bf j}&{\bf k}\
  2&2&1\
  1&2&3
  \end{vmatrix}
  =4{\bf i}-5{\bf j}+2{\bf k}=(4,-5,2).
  $$

外積の性質

(I) $ \ b\times a=-(a\times b)$（反交換則）
(II) $ a\times(b+c)=a\times b+a\times c$（分配則）
(III) $(ma)\times b=a\times(mb)=m(a\times b)$（$m$ は実数）
(IV) $a\neq0,\ b\neq0$ のとき $a\parallel b \Longleftrightarrow a\times b=0$.

（下に (II) の証明：行列式の線形性を使って
$\ a\times(b+c)=\begin{vmatrix} {\bf i}&{\bf j}&{\bf k}\ a_x&a_y&a_z\ b_x+c_x&b_y+c_y&b_z+c_z \end{vmatrix} =\begin{vmatrix}\cdots\ b\end{vmatrix}
+\begin{vmatrix}\cdots\ c\end{vmatrix} =a\times b + a\times c.$）

演習（抜粋）

• 問4 $a=(-1,3,2),\ b=(0,-1,0)$ のとき $a\times b$ を求めよ。
• 問5 空間内の三点 $A(2,1,3),,B(1,-1,2),,C(2,2,1)$。
  $AB\times AC$ を求め，三角形 $ABC$ の面積を求めよ。
• 問6 基本ベクトル ${\bf i},{\bf j}$ について $({\bf i}\times{\bf i})\times{\bf j}$ および ${\bf i}\times({\bf i}\times{\bf j})$ を求めよ。

---

工学的解説（単位つき）

平行四辺形面積・法線（面積ベクトル）

$$
{\bf S}=a\times b,\quad |{\bf S}|=|a||b|\sin\theta\ [\mathrm{m}^2],\quad
\hat{\bf n}=\dfrac{{\bf S}}{|{\bf S}|}\ (\text{無次元})
$$

• 例4の面積：$|(4,-5,2)|=3\sqrt5\ \mathrm{m}^2$。
  単位法線：$\pm\dfrac{1}{3\sqrt5}(4,-5,2)$（無次元）。
• 用途：CAD/メッシュ法線，電磁気の面積分（$\mathrm{d}{\bf S}=\hat{\bf n},\mathrm{d}S,[\mathrm{m}^2]$）。

トルク（モーメント）

$$
{\bf T}={\bf r}\times{\bf F}\quad[\mathrm{N\cdot m}]
$$

• ${\bf r},[\mathrm{m}],\ {\bf F},[\mathrm{N}]$。レンチ・ロータの設計に必須。

角運動量

$$
{\bf L}={\bf r}\times{\bf p}\quad[\mathrm{kg,m^2/s}],\quad {\bf p}=m{\bf v}
$$

• ジャイロ，衛星姿勢制御。

ローレンツ力（磁場成分）

$$
{\bf F}_B=q({\bf v}\times{\bf B})\quad[\mathrm{N}]
$$

• $q,[\mathrm{C}],\ {\bf v},[\mathrm{m/s}],\ {\bf B},[\mathrm{T}]$。モータ・発電機。

---

演習の解答

• 問4
  $a\times b=(a_y b_z-a_z b_y,\ a_z b_x-a_x b_z,\ a_x b_y-a_y b_x)$
  $=(3\cdot0-2\cdot(-1),\ 2\cdot0-(-1)\cdot0,\ -1\cdot(-1)-3\cdot0)=(2,0,1)$（無次元の幾何ベクトル）。
  もし $a,b$ を長さ [m] とみなせば，単位は $[\mathrm{m}^2]$。

• 問5
  $AB=B-A=(-1,-2,-1),\ AC=C-A=(0,1,-2)$。
  $AB\times AC=(5,-2,-1)$。
  三角形面積 $=\tfrac12\lVert(5,-2,-1)\rVert=\tfrac12\sqrt{30}\ [\mathrm{m}^2]$。

• 問6
  ${\bf i}\times{\bf i}=0 \Rightarrow ({\bf i}\times{\bf i})\times{\bf j}=0$。
  ${\bf i}\times{\bf j}={\bf k}\Rightarrow {\bf i}\times({\bf i}\times{\bf j})={\bf i}\times{\bf k}=-{\bf j}$（右ねじ則）。

---

Python実装

import numpy as np

# 外積を計算する関数
def cross_product(a, b):
    return np.cross(a, b)

# 問4: a = (-1, 3, 2), b = (0, -1, 0) の外積
a = np.array([-1, 3, 2])
b = np.array([0, -1, 0])
result_q4 = cross_product(a, b)
print("問4: a × b =", result_q4)  # 出力: [2, 0, 1]

# 問5: AB × AC と三角形 ABC の面積
A = np.array([2, 1, 3])
B = np.array([1, -1, 2])
C = np.array([2, 2, 1])
AB = B - A
AC = C - A
result_q5 = cross_product(AB, AC)
area = np.linalg.norm(result_q5) / 2
print("問5: AB × AC =", result_q5)  # 出力: [ 5, -2, -1]
print("問5: 三角形 ABC の面積 =", area, "m^2")  # 出力: sqrt(30)/2

# 問6: 基本ベクトルの外積
i = np.array([1, 0, 0])
j = np.array([0, 1, 0])
k = np.array([0, 0, 1])
# (i) (i × i) × j
result_q6_1 = cross_product(cross_product(i, i), j)
print("問6 (i): (i × i) × j =", result_q6_1)  # 出力: [0, 0, 0]
# (ii) i × (i × j)
result_q6_2 = cross_product(i, cross_product(i, j))
print("問6 (ii): i × (i × j) =", result_q6_2)  # 出力: [0, -1, 0]

1) ベクトル関数・極限・連続

定義

$$
\mathbf{a}(t)=(a_x(t),a_y(t),a_z(t))
$$

各成分は実関数。位置ベクトルは

$$
\mathbf{r}(t)=(x(t),y(t),z(t))\ [\mathrm{m}]
$$

極限

$$
\lim_{t\to t_0}\mathbf{a}(t)=\mathbf{b}
\iff
\begin{cases}
\lim a_x=b_x,\
\lim a_y=b_y,\
\lim a_z=b_z
\end{cases}
$$

連続
$\lim_{t\to t_0}\mathbf{a}(t)=\mathbf{a}(t_0)$（成分ごとに連続）。

---

2) 微分（成分形・単位）

定義

$$
\mathbf{a}'(t)=\frac{d\mathbf{a}}{dt}
=\lim_{\Delta t\to0}\frac{\mathbf{a}(t+\Delta t)-\mathbf{a}(t)}{\Delta t}
=\Big(\frac{da_x}{dt},\frac{da_y}{dt},\frac{da_z}{dt}\Big)
$$

2次微分：$\displaystyle
\mathbf{a}''(t)=\Big(\frac{d^2a_x}{dt^2},\frac{d^2a_y}{dt^2},\frac{d^2a_z}{dt^2}\Big)
$

例（本の例と同じ）
$\mathbf{a}(t)=(\cos t,\sin t,1)\Rightarrow \mathbf{a}'(t)=(-\sin t,\cos t,0)$。
$t=0$: $\mathbf{a}'(0)=(0,1,0)$。

単位の読み替え（$\mathbf{r}(t)$ を位置[m]と解釈）

$$
\mathbf{v}(t)=\frac{d\mathbf{r}}{dt}\ [\mathrm{m/s}],\quad
\mathbf{a}(t)=\frac{d^2\mathbf{r}}{dt^2}\ [\mathrm{m/s^2}],\quad
\mathbf{j}(t)=\frac{d^3\mathbf{r}}{dt^3}\ [\mathrm{m/s^3}]
$$

---

3) 微分の公式（工学必携）

• スカラー倍：$\displaystyle \frac{d}{dt}[m(t)\mathbf{a}]=m'\mathbf{a}+m\mathbf{a}'$（$m$ は無次元または任意の物理量）
• 内積：$\displaystyle \frac{d}{dt}(\mathbf{a}!\cdot!\mathbf{b})=\mathbf{a}'!\cdot!\mathbf{b}+\mathbf{a}!\cdot!\mathbf{b}'$
• 外積：$\displaystyle \frac{d}{dt}(\mathbf{a}!\times!\mathbf{b})=\mathbf{a}'!\times!\mathbf{b}+\mathbf{a}!\times!\mathbf{b}'$
• 大きさ：$\displaystyle \frac{d}{dt}|\mathbf{a}|=\frac{\mathbf{a}\cdot\mathbf{a}'}{|\mathbf{a}|}$（$|\mathbf{a}|\neq0$）
• 合成関数：$\displaystyle \frac{d}{dt}\mathbf{a}(g(t))=\frac{d\mathbf{a}}{dg},g'(t)$

単位チェック例
$\frac{d}{dt}(\mathbf{F}!\cdot!\mathbf{v})$ は $[\mathrm{N}]\cdot[\mathrm{m/s}]=[\mathrm{W}]$ の時間微分 → $[\mathrm{W/s}]$。

---

4) 幾何量（運動学）

• 速さ：$\displaystyle v=|\mathbf{r}'(t)|\ [\mathrm{m/s}]$
• 弧長：$\displaystyle s(t)=\int_{t_0}^{t}|\mathbf{r}'(\tau)|,d\tau\ [\mathrm{m}]$
• 接ベクトル：$\displaystyle \mathbf{T}=\frac{\mathbf{r}'}{|\mathbf{r}'|}$（無次元）
• 曲率：$\displaystyle \kappa=\Big|\frac{d\mathbf{T}}{ds}\Big|=\frac{|\mathbf{r}'\times\mathbf{r}''|}{|\mathbf{r}'|^{3}}\ [\mathrm{m^{-1}}]$
• 法線：$\displaystyle \mathbf{N}=\frac{1}{\kappa}\frac{d\mathbf{T}}{ds}$（無次元）
• 従法線：$\displaystyle \mathbf{B}=\mathbf{T}\times\mathbf{N}$（無次元）
• 加速度の分解

$$
\mathbf{a}=\underbrace{\frac{dv}{dt}}{[\mathrm{m/s^2}]}\mathbf{T}
+\underbrace{v^2\kappa}{[\mathrm{m/s^2}]}\mathbf{N}
$$

（接線加速度＋法線（向心）加速度）

---

5) 回転運動（剛体）

• 角速度：$\boldsymbol{\omega}\ [\mathrm{rad/s}]$（radは無次元として扱う）
• 角加速度：$\boldsymbol{\alpha}\ [\mathrm{rad/s^2}]$
• 純回転の速度：$\displaystyle \mathbf{v}=\boldsymbol{\omega}\times\mathbf{r}\ [\mathrm{m/s}]$
• 加速度：$\displaystyle \mathbf{a}=\boldsymbol{\alpha}\times\mathbf{r}+\boldsymbol{\omega}\times(\boldsymbol{\omega}\times\mathbf{r})\ [\mathrm{m/s^2}]$

仕事率（出力）

$$
P=\mathbf{F}\cdot\mathbf{v}\ [\mathrm{W}],\qquad
P=\boldsymbol{\tau}\cdot\boldsymbol{\omega}\ [\mathrm{W}]
$$

（$\boldsymbol{\tau}$ はトルク [N·m]）

---

6) 電磁気・流体（時間変化を伴う応用）

• ローレンツ力：$\mathbf{F}=q(\mathbf{E}+\mathbf{v}\times\mathbf{B})\ [\mathrm{N}]$
• 磁束：$\Phi(t)=\int_S \mathbf{B}(t)\cdot d\mathbf{S}\ [\mathrm{Wb=T\cdot m^2}]$
• ファラデーの法則：$\mathcal{E}=-\frac{d\Phi}{dt}\ [\mathrm{V}]$。
  面や法線が時間で動くなら
  $\displaystyle \frac{d}{dt}(\mathbf{B}\cdot\mathbf{S}) =\dot{\mathbf{B}}\cdot\mathbf{S}+\mathbf{B}\cdot\dot{\mathbf{S}}$
  を用いる（内積微分則）。
• 体積流量：$\displaystyle Q=\int_A \mathbf{v}\cdot d\mathbf{S}\ [\mathrm{m^3/s}]$

---

7) 具体例（数値・単位つき）

(a) らせん運動（半径 $R,[\mathrm{m}]$、角速度 $\omega,[\mathrm{s^{-1}}]$、軸方向速度 $v_z,[\mathrm{m/s}]$）

$$
\mathbf{r}(t)=(R\cos\omega t,\ R\sin\omega t,\ v_z t)\ [\mathrm{m}]
$$

$$
\mathbf{v}(t)=(-R\omega\sin\omega t,\ R\omega\cos\omega t,\ v_z)\ [\mathrm{m/s}]
$$

$$
\mathbf{a}(t)=(-R\omega^2\cos\omega t,\ -R\omega^2\sin\omega t,\ 0)\ [\mathrm{m/s^2}]
$$

速さ $v=\sqrt{(R\omega)^2+v_z^2}\ [\mathrm{m/s}]$、
曲率 $\displaystyle \kappa=\frac{R\omega^2}{\big((R\omega)^2+v_z^2\big)^{3/2}}\ [\mathrm{m^{-1}}]$、
法線加速度 $v^2\kappa=\frac{(R\omega)^2}{\sqrt{(R\omega)^2+v_z^2}}\ [\mathrm{m/s^2}]$。

(b) トルクと出力

レンチ長 $L=0.25,[\mathrm{m}]$、力 $F=80,[\mathrm{N}]$（直交）
→ $|\boldsymbol{\tau}|=LF=20,[\mathrm{N\cdot m}]$。
角速度 $\omega=10,[\mathrm{rad/s}]$
→ 出力 $P=\tau\omega=200,[\mathrm{W}]$。

(c) 誘導起電力

コイル面積 $S=0.01,[\mathrm{m^2}]$、磁束密度が $B(t)=0.5\sin(100t),[\mathrm{T}]$（法線に平行）
→ $\Phi(t)=BS=0.005\sin(100t),[\mathrm{Wb}]$
→ $\mathcal{E}=-\dot{\Phi}= -0.5\cos(100t),[\mathrm{V}]$。

---

Python実装

import numpy as np
import sympy as sp
from sympy.vector import CoordSys3D

# 座標系と変数の定義
C = CoordSys3D('C')  # 3D座標系
t = sp.symbols('t')  # 時間変数
R, omega, v_z = sp.symbols('R omega v_z')  # 定数

# --- 具体例 (a): らせん運動 ---
print("--- 具体例 (a): らせん運動 ---")
# 位置ベクトル r(t) = (R*cos(ωt), R*sin(ωt), v_z*t)
r = R * sp.cos(omega * t) * C.i + R * sp.sin(omega * t) * C.j + v_z * t * C.k

# 速度 v(t) = dr/dt
v = r.diff(t)
print("速度 v(t) =", v)

# 加速度 a(t) = dv/dt
a = v.diff(t)
print("加速度 a(t) =", a)

# 速さ v = ||v(t)||
v_norm = sp.sqrt(v.dot(v))
print("速さ v =", v_norm)

# 曲率 κ = ||r' × r''|| / ||r'||^3
r_prime = v  # r'=v
r_double_prime = a  # r''=a
cross_product = r_prime.cross(r_double_prime)
cross_norm = sp.sqrt(cross_product.dot(cross_product))
r_prime_norm = sp.sqrt(r_prime.dot(r_prime))
kappa = cross_norm / r_prime_norm**3
print("曲率 κ =", kappa)

# 法線加速度 v^2 * κ
normal_acc = v_norm**2 * kappa
print("法線加速度 v^2*κ =", normal_acc)

# 数値例（R=1 m, ω=2 s^-1, v_z=0.5 m/s, t=0.1）
subs_dict = {R: 1, omega: 2, v_z: 0.5, t: 0.1}
v_num = v.subs(subs_dict)
v_norm_num = v_norm.subs(subs_dict)
a_num = a.subs(subs_dict)
kappa_num = kappa.subs(subs_dict)
normal_acc_num = normal_acc.subs(subs_dict)
print("\n数値例 (t=0.1):")
print("速度 v(0.1) =", [float(v_num.dot(C.i)), float(v_num.dot(C.j)), float(v_num.dot(C.k))], "m/s")
print("速さ v =", float(v_norm_num), "m/s")
print("加速度 a(0.1) =", [float(a_num.dot(C.i)), float(a_num.dot(C.j)), float(a_num.dot(C.k))], "m/s^2")
print("曲率 κ =", float(kappa_num), "m^-1")
print("法線加速度 =", float(normal_acc_num), "m/s^2")

# --- 具体例 (b): トルクと出力 ---
print("\n--- 具体例 (b): トルクと出力 ---")
L = 0.25  # レンチ長 [m]
F = 80    # 力 [N]
omega_val = 10  # 角速度 [rad/s]
tau = L * F  # トルク |τ| = L*F
P = tau * omega_val  # 出力 P = τ*ω
print("トルク |τ| =", tau, "N·m")
print("出力 P =", P, "W")

# --- 具体例 (c): 誘導起電力 ---
print("\n--- 具体例 (c): 誘導起電力 ---")
S = 0.01  # 面積 [m^2]
B_t = 0.5 * sp.sin(100 * t)  # 磁束密度 [T]
Phi = B_t * S  # 磁束 Φ = B*S [Wb]
E = -Phi.diff(t)  # 誘導起電力 E = -dΦ/dt [V]
print("磁束 Φ(t) =", Phi, "Wb")
print("誘導起電力 E(t) =", E, "V")

# 数値例（t=0.1）
E_num = E.subs(t, 0.1)
print("誘導起電力 E(0.1) =", float(E_num), "V")

（曲線・接線ベクトル・単位接線・弧長）

1) 曲線と位置ベクトル

• 空間内を動く点 $P$ の位置ベクトル

  $$
  \mathbf{r}(t)=(x(t),y(t),z(t))\quad[\mathrm{m}]
  $$

  時刻 $t$ により軌跡 $C$（曲線）を描く。

2) 接線ベクトルの導入

• 微小時間 $\Delta t$ の位置変化

  $$
  \Delta \mathbf{r}=\mathbf{r}(t+\Delta t)-\mathbf{r}(t)\quad[\mathrm{m}]
  $$

• 割線ベクトル

  $$
  \frac{\Delta \mathbf{r}}{\Delta t}\quad[\mathrm{m/s}]
  $$

  2点 $P(t),P(t+\Delta t)$ を結ぶ直線に平行。

• $\Delta t\to 0$ で極限が存在すれば

  $$
  \boxed{\ \mathbf{r}'(t)=\frac{d\mathbf{r}}{dt}=(\dot x,\dot y,\dot z)\quad[\mathrm{m/s}] \ }
  $$

  これは曲線 $C$ の接線ベクトル（接線方向と同向）。

3) 単位接線ベクトル（速度の方向）

$$
\boxed{\ \mathbf{t}(t)=\frac{\mathbf{r}'(t)}{\lVert\mathbf{r}'(t)\rVert}\ }\quad(\lVert\mathbf{t}\rVert=1,\ \text{無次元})
$$

• 速さ（スカラー速度）

  $$
  v(t)=\lVert\mathbf{r}'(t)\rVert=\sqrt{\dot x^{,2}+\dot y^{,2}+\dot z^{,2}}\quad[\mathrm{m/s}]
  $$

例6（テキストと同じ：円柱らせん）

定数 $a>0,\ b>0$ とし

$$
\mathbf{r}(t)=(a\cos t,\ a\sin t,\ bt)\quad[\mathrm{m}]
$$

$$
\mathbf{r}'(t)=(-a\sin t,\ a\cos t,\ b)\quad[\mathrm{m/s}]
$$

$$
\lVert\mathbf{r}'(t)\rVert=\sqrt{(-a\sin t)^2+(a\cos t)^2+b^2}=\sqrt{a^2+b^2}\quad[\mathrm{m/s}]
$$

$$
\boxed{\ \mathbf{t}(t)=\frac{1}{\sqrt{a^2+b^2}};(-a\sin t,\ a\cos t,\ b)\ } \ (\text{無次元})
$$

工学メモ（らせんの幾何）：半径 $a,[\mathrm{m}]$、1回転あたりのピッチ $p=2\pi b,[\mathrm{m}]$。傾斜角 $\phi$ は $\tan\phi=b/a$。
曲率 $\kappa=\dfrac{a}{a^2+b^2},[\mathrm{m^{-1}}]$、ねじれ $\tau=\dfrac{b}{a^2+b^2},[\mathrm{m^{-1}}]$。

4) 弧長（曲線の長さ）

区間 $t\in[a,b]$ で曲線の長さ $s$ は

$$
\boxed{\ s=\int_a^b \left\lVert\frac{d\mathbf{r}}{dt}\right\rVert dt
=\int_a^b \sqrt{\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2+\left(\frac{dz}{dt}\right)^2};dt\quad[\mathrm{m}] \ }
$$

• 例6 のらせん：

  $$
  s=\int_{t=a_0}^{t=b_0}\sqrt{a^2+b^2};dt=\sqrt{a^2+b^2},(b_0-a_0)\quad[\mathrm{m}]
  $$

---

工学的応用（単位つき）

1. ロボット・CNC/3Dプリンタのツールパス

• 経路 $\mathbf{r}(t),[\mathrm{m}]$、速さ $v=|\mathbf{r}'|,[\mathrm{m/s}]$。
• 走行距離 $s,[\mathrm{m}]$、等速なら加工時間 $T=s/v,[\mathrm{s}]$。
• らせん送り：半径 $a$ の穴加工で 1回転の距離 $s_{\mathrm{1rev}}=\sqrt{(2\pi a)^2+p^2},[\mathrm{m}]$。

1. ケーブル・配管長の見積

• コイル/スプリング形状をらせんで近似 → 必要長 $s,[\mathrm{m}]$ を弧長で算出。

1. 動力学（速度・出力）

• 速度ベクトル $\mathbf{v}=\mathbf{r}',[\mathrm{m/s}]$。
• 仕事率 $P=\mathbf{F}\cdot\mathbf{v},[\mathrm{W}]$、総仕事 $W=\int \mathbf{F}\cdot d\mathbf{r},[\mathrm{J}]$。

1. 回転機械の幾何

• ねじ/ボールねじ：ピッチ $p=2\pi b,[\mathrm{m}]$。
• 糸巻き・巻線：層ごとの半径 $a$ とピッチ $p$ から総長 $s$ を積分。

1. 電磁気（可動コイル）

• コイル面法線 $\hat{\mathbf{n}}$ と運動 $\mathbf{r}(t)$ から磁束 $\Phi(t)=\int\mathbf{B}\cdot d\mathbf{S},[\mathrm{Wb}]$ を評価。
• $\mathcal{E}=-d\Phi/dt,[\mathrm{V}]$（運動誘導は面や法線の時間変化を考慮）。

---

# Program Name: curve_tangent_arclength.py
# Creation Date: 20250821
# Overview: Compute and visualize a 3D curve r(t), its tangent r'(t), unit tangent t_hat(t), speed, and arc length. Includes an engineering readout (toolpath length and constant-speed machining time).
# Usage: Run as-is. Tweak parameters in the "Parameters" section (e.g., A_RADIUS_M, B_RISE_PER_RAD_M, T_START, T_END, NUM_SAMPLES, FEED_SPEED_MPS). The script prints key results and makes plots with matplotlib.

!pip install numpy matplotlib -q

import numpy as np
import matplotlib.pyplot as plt

# =========================
# Parameters / パラメータ一元管理（すべてここで設定）
# =========================
A_RADIUS_M           = 0.020      # Helix radius [m] / らせん半径 [m]
B_RISE_PER_RAD_M     = 0.002      # z rise per rad (pitch/(2π)) [m/rad] / 1ラジアンあたりの上昇量 [m/rad]
T_START              = 0.0        # Start parameter t [rad] / 開始パラメータt [rad]
T_END                = 6*np.pi    # End parameter t [rad] / 終了パラメータt [rad]
NUM_SAMPLES          = 1000       # Number of samples / サンプル数
FEED_SPEED_MPS       = 0.05       # Constant feed speed for time estimate [m/s] / 等速送りの仮定 [m/s]
PRINT_DIGITS         = 6          # Print precision / 表示桁数
PLOT_3D              = True       # Make 3D curve plot / 3D曲線プロットを行う
PLOT_SPEED           = True       # Plot speed vs t / 速度-時間のプロット
PLOT_UNIT_TANGENT    = True       # Plot unit tangent components / 単位接線ベクトル成分のプロット
PLOT_ARC_LENGTH      = True       # Plot arc length vs t / 弧長-時間のプロット

# =========================
# Curve definition / 曲線定義（r(t) = (x(t), y(t), z(t))）
# Default: circular helix r(t) = (a cos t, a sin t, b t)
# 既定：円柱らせん r(t) = (a cos t, a sin t, b t)
# =========================
def r_of_t(t):
    # English: Position vector r(t) [m]
    # 日本語：位置ベクトル r(t) [m]
    a = A_RADIUS_M
    b = B_RISE_PER_RAD_M
    return np.stack((a*np.cos(t), a*np.sin(t), b*t), axis=-1)

# =========================
# Numerical differentiation helpers / 数値微分ヘルパ
# =========================
def derivative(f_vals, t):
    # English: First derivative using central differences (np.gradient)
    # 日本語：中央差分（np.gradient）による1階微分
    return np.gradient(f_vals, t, axis=0)

def speed_from_rdot(rdot):
    # English: Speed ||r'(t)|| [m/s]
    # 日本語：速さ ||r'(t)|| [m/s]
    return np.linalg.norm(rdot, axis=1)

def unit_tangent(rdot):
    # English: Unit tangent t_hat = r'(t)/||r'(t)||
    # 日本語：単位接線ベクトル t_hat = r'(t)/||r'(t)||
    v = speed_from_rdot(rdot)
    # Avoid division by zero / 0除算回避
    v_safe = np.where(v == 0, 1.0, v)
    return (rdot.T / v_safe).T

def arc_length(speed, t):
    # English: Arc length s = ∫ ||r'(t)|| dt via cumulative trapezoid
    # 日本語：弧長 s = ∫ ||r'(t)|| dt（台形則の累積）
    s = np.zeros_like(speed)
    s[1:] = np.cumsum(0.5*(speed[1:] + speed[:-1])*(t[1:] - t[:-1]))
    return s

# =========================
# Main computation / 主計算
# =========================
t = np.linspace(T_START, T_END, NUM_SAMPLES)                 # Parameter grid [rad] / パラメータ格子 [rad]
r = r_of_t(t)                                                # r(t) [m]
rdot = derivative(r, t)                                      # r'(t) [m/s]
v = speed_from_rdot(rdot)                                    # Speed [m/s]
that = unit_tangent(rdot)                                    # Unit tangent (dimensionless) / 単位接線（無次元）
s = arc_length(v, t)                                         # Arc length [m]

# Engineering readouts / エンジニアリング出力
pitch_per_turn_m = 2*np.pi*B_RISE_PER_RAD_M                  # Pitch per revolution [m] / 1回転ピッチ [m]
length_per_turn_m = np.sqrt((2*np.pi*A_RADIUS_M)**2 + pitch_per_turn_m**2)  # 1回転の弧長 [m]
total_length_m = float(s[-1])                                # Total path length [m] / 全経路長 [m]
time_const_speed_s = total_length_m / FEED_SPEED_MPS         # Time at constant speed [s] / 等速での所要時間 [s]
curvature_helix_per_m = A_RADIUS_M / (A_RADIUS_M**2 + B_RISE_PER_RAD_M**2)  # Helix curvature [1/m] / 曲率 [1/m]
torsion_helix_per_m = B_RISE_PER_RAD_M / (A_RADIUS_M**2 + B_RISE_PER_RAD_M**2)  # Helix torsion [1/m] / ねじれ [1/m]

# =========================
# Printing / 数値結果の表示
# =========================
np.set_printoptions(precision=PRINT_DIGITS, suppress=True)
print("# Results / 結果")
print(f"Radius a [m]                       : {A_RADIUS_M:.{PRINT_DIGITS}f}")
print(f"Rise per rad b [m/rad]             : {B_RISE_PER_RAD_M:.{PRINT_DIGITS}f}")
print(f"Pitch per turn 2πb [m]             : {pitch_per_turn_m:.{PRINT_DIGITS}f}")
print(f"Length per turn [m]                : {length_per_turn_m:.{PRINT_DIGITS}f}")
print(f"Total arc length s [m]             : {total_length_m:.{PRINT_DIGITS}f}")
print(f"Const feed speed [m/s]             : {FEED_SPEED_MPS:.{PRINT_DIGITS}f}")
print(f"Const-speed time s/v [s]           : {time_const_speed_s:.{PRINT_DIGITS}f}")
print(f"Helix curvature κ [1/m]            : {curvature_helix_per_m:.{PRINT_DIGITS}f}")
print(f"Helix torsion τ [1/m]              : {torsion_helix_per_m:.{PRINT_DIGITS}f}")

# =========================
# Plots / プロット（英語表記）
# 各図は単独（サブプロットなし）、色は指定しない
# =========================
if PLOT_3D:
    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
    fig = plt.figure(figsize=(7,5))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot(r[:,0], r[:,1], r[:,2])
    # Draw a few unit tangents / 単位接線ベクトルを数点描画
    pick = np.linspace(0, len(t)-1, 12, dtype=int)
    scale = max(A_RADIUS_M, B_RISE_PER_RAD_M*2*np.pi)*0.3  # Short arrow scale / 矢印の長さ
    for idx in pick:
        ax.quiver(r[idx,0], r[idx,1], r[idx,2],
                  that[idx,0], that[idx,1], that[idx,2],
                  length=scale, normalize=True)
    ax.set_title("3D Curve r(t) and Sample Unit Tangents")
    ax.set_xlabel("x [m]")
    ax.set_ylabel("y [m]")
    ax.set_zlabel("z [m]")
    plt.show()

if PLOT_SPEED:
    plt.figure(figsize=(7,4))
    plt.plot(t, v)
    plt.title("Speed |r'(t)| vs t")
    plt.xlabel("t [rad]")
    plt.ylabel("Speed [m/s]")
    plt.grid(True)
    plt.show()

if PLOT_UNIT_TANGENT:
    plt.figure(figsize=(7,4))
    plt.plot(t, that[:,0], label="t_x")
    plt.plot(t, that[:,1], label="t_y")
    plt.plot(t, that[:,2], label="t_z")
    plt.title("Unit Tangent Components")
    plt.xlabel("t [rad]")
    plt.ylabel("Component [—]")
    plt.legend()
    plt.grid(True)
    plt.show()

if PLOT_ARC_LENGTH:
    plt.figure(figsize=(7,4))
    plt.plot(t, s)
    plt.title("Arc Length s(t)")
    plt.xlabel("t [rad]")
    plt.ylabel("Arc Length [m]")
    plt.grid(True)
    plt.show()

1) 例題1：弧長の計算

曲線 $\displaystyle \mathbf{r}(t)=\bigl(t,\ t^{2},\ \tfrac{2}{3}t^{3}\bigr)$（座標を長さ[m]と解釈）
区間 $t:0\to1$ の弧長 $s,[\mathrm{m}]$。

$$
\frac{d\mathbf{r}}{dt}=(1,\ 2t,\ 2t^{2})\quad[\mathrm{m/s}]
$$

$$
\left|\frac{d\mathbf{r}}{dt}\right|
=\sqrt{1+(2t)^2+(2t^2)^2}
=\sqrt{1+4t^2+4t^4}
=\sqrt{(1+2t^2)^2}=1+2t^2\quad[\mathrm{m/s}]
$$

$$
\boxed{,s=\int_{0}^{1}!\left|\frac{d\mathbf{r}}{dt}\right|dt
=\int_{0}^{1}(1+2t^2),dt
=\Bigl[t+\tfrac{2}{3}t^{3}\Bigr]_{0}^{1}
=\tfrac{5}{3}\ \mathrm{m},}
$$

2) 問11：弧長（解）

$\displaystyle \mathbf{r}(t)=(\cos t,\ \sin t,\ t)$（座標[m]、$t:0\to2\pi$）

$$
\mathbf{r}'(t)=(-\sin t,\ \cos t,\ 1)\quad[\mathrm{m/s}],\quad
|\mathbf{r}'(t)|=\sqrt{\sin^2 t+\cos^2 t+1}=\sqrt{2}\ [\mathrm{m/s}]
$$

$$
\boxed{,s=\int_{0}^{2\pi}\sqrt{2},dt=2\pi\sqrt{2}\ \mathrm{m},}
$$

---

3) §1.5 曲面：2変数ベクトル関数と偏導関数

2変数のベクトル関数 $\displaystyle \mathbf{a}(u,v)=(a_x(u,v),a_y(u,v),a_z(u,v))$。
偏導関数（ベクトル）：

$$
\frac{\partial \mathbf{a}}{\partial u}
=\Bigl(\frac{\partial a_x}{\partial u},\frac{\partial a_y}{\partial u},\frac{\partial a_z}{\partial u}\Bigr),\quad
\frac{\partial \mathbf{a}}{\partial v}
=\Bigl(\frac{\partial a_x}{\partial v},\frac{\partial a_y}{\partial v},\frac{\partial a_z}{\partial v}\Bigr)
$$

（$\mathbf{a}$ が位置[$\mathrm{m}$]なら $\partial\mathbf{a}/\partial u$ の単位は $[\mathrm{m}/u]$）

点の位置を与える曲面は $\mathbf{r}(u,v)\ [\mathrm{m}]$。

$$
\mathbf{r}_u\equiv\frac{\partial\mathbf{r}}{\partial u},\quad
\mathbf{r}_v\equiv\frac{\partial\mathbf{r}}{\partial v}
$$

• $v$ を固定し $u$ を変える曲線の接ベクトルは $\mathbf{r}_u$（同様に $\mathbf{r}_v$）。
• 点 $P=\mathbf{r}(u_0,v_0)$ における接平面は $\mathbf{r}_u,\mathbf{r}_v$ を含む平面。
• 法線ベクトルの一つ：$\ \mathbf{r}_u\times\mathbf{r}_v$（両方に直交）。
• 単位法線（向きは $\pm$ の2通り）：

$$
\boxed{,\mathbf{n}=\pm\frac{\mathbf{r}_u\times\mathbf{r}_v}{|\mathbf{r}_u\times\mathbf{r}_v|},}\quad(\text{無次元})
$$

• 面素（面積要素）：$\displaystyle dS=|\mathbf{r}_u\times\mathbf{r}_v|,du,dv\ [\mathrm{m^2}]$

例：円柱側面 $ \mathbf{r}(u,v)=(a\cos u,\ a\sin u,\ v)$

$$
\mathbf{r}_u=(-a\sin u,\ a\cos u,0),\quad
\mathbf{r}_v=(0,0,1)
$$

$$
\mathbf{r}_u\times\mathbf{r}_v=(a\cos u,\ a\sin u,0),\
|\mathbf{r}_u\times\mathbf{r}_v|=a
$$

$$
\mathbf{n}=(\cos u,\ \sin u,\ 0)\quad(\text{半径方向}),\quad
dS=a,du,dv\ [\mathrm{m^2}]
$$

半開区間 $u\in[0,2\pi),\ v\in[z_1,z_2]$ の側面積は $S=\int a,du,dv=a(2\pi)(z_2-z_1)=2\pi a L\ [\mathrm{m^2}]$。

---

4) 工学的応用（単位つき）

1. CAD/CG・メッシュ

• 法線 $\mathbf{n}$（無次元）は陰影・レンダリングに必須。
• 面素 $dS[\mathrm{m^2}]$ で表面積、膜の質量 $m=\rho \int dS\ [\mathrm{kg}]$。

1. 電磁気のフラックス

• 磁束 $\displaystyle \Phi=\int_S \mathbf{B}\cdot\mathbf{n},dS\ [\mathrm{Wb}]$。
• 例：円柱側面に軸方向磁場 $\mathbf{B}=B_z\hat{\mathbf{z}}$ なら
  $\mathbf{n}\perp\hat{\mathbf{z}}\Rightarrow \mathbf{B}\cdot\mathbf{n}=0$ → 側面の磁束は 0。

1. 流体の体積流量

• $\displaystyle Q=\int_A \mathbf{v}\cdot\mathbf{n},dS\ [\mathrm{m^3/s}]$。
• 配管壁（円柱側面）は $\mathbf{n}\perp$ 流れ ⇒ 漏れない理想なら $Q=0$。

1. 応力・圧力の合力

• 面の法線に沿う力：$\displaystyle \mathbf{F}=\int_A p,\mathbf{n},dS\ [\mathrm{N}]$（圧力 $p[\mathrm{Pa=N/m^2}]$）。

1. 表面積と材料見積

• コーティング量、放熱面積、ラッピング長：$\int dS\ [\mathrm{m^2}]$。

---

速見チェック

• 弧長：$s=\int|\mathbf{r}'|dt\ [\mathrm{m}]$
• 接平面：$\mathrm{span}{\mathbf{r}_u,\mathbf{r}_v}$
• 単位法線：$\mathbf{n}=\pm(\mathbf{r}_u\times\mathbf{r}_v)/|\cdot|$
• 面素：$dS=|\mathbf{r}_u\times\mathbf{r}_v|,du,dv\ [\mathrm{m^2}]$

曲面の単位法線・面素・面積：詳解（すべてSI単位付き）

1) パラメータ曲面と微分

• 曲面：$\mathbf{r}(u,v)=[x(u,v),y(u,v),z(u,v)]$（位置の物理量 → [m]）
• 偏導関数：$\mathbf{r}_u=\partial\mathbf{r}/\partial u$（[m/単位u]），$\mathbf{r}_v=\partial\mathbf{r}/\partial v$（[m/単位v]）
• 直交性：$\mathbf{r}_u,\mathbf{r}_v$ はその点の 接平面 の基底

2) 法線と面素

• 法線（向き付き面素）

  $$
  \mathbf{r}_u\times\mathbf{r}_v\quad[\mathrm{m^2/(u,v)}]
  $$

  これは $\mathbf{r}_u,\mathbf{r}_v$ の両方に直交。右ねじで向き決定。

• 単位法線（無次元）

  $$
  \boxed{\ \mathbf{n}=\pm\frac{\mathbf{r}_u\times\mathbf{r}_v}{|\mathbf{r}_u\times\mathbf{r}_v|}\ }
  $$

• 面素（微小面積）

  $$
  \boxed{\ dS=|\mathbf{r}_u\times\mathbf{r}_v|;du,dv\quad[\mathrm{m^2}] }
  $$

• 面積

  $$
  \boxed{\ S=\iint_D|\mathbf{r}_u\times\mathbf{r}_v|;du,dv\quad[\mathrm{m^2}] }
  $$

3) グラフ面 $z=f(x,y)$ の特化式

$\mathbf{r}(x,y)=(x,y,f(x,y))$ と置く（$u!=!x,\ v!=!y$ は [m]）

$$
\mathbf{r}_x=(1,0,f_x),\ \mathbf{r}_y=(0,1,f_y),\
\mathbf{r}_x\times\mathbf{r}_y=(-f_x,-f_y,1)
$$

$$
\boxed{\ \mathbf{n}=\pm\frac{(-f_x,-f_y,1)}{\sqrt{1+f_x^2+f_y^2}}\ }(\text{無次元}),\quad
\boxed{\ dS=\sqrt{1+f_x^2+f_y^2};dx,dy\ [\mathrm{m^2}] }
$$

• 傾きが小さいときの近似：$\sqrt{1+f_x^2+f_y^2}\approx 1+\tfrac12(f_x^2+f_y^2)$

4) 代表例と単位

(a) 円柱側面 $\mathbf{r}(u,v)=(a\cos u,,a\sin u,,v)$

• $a$[m], $u$[rad], $v$[m]
  $\mathbf{r}_u=(-a\sin u,a\cos u,0),\ \mathbf{r}_v=(0,0,1)$
  $|\mathbf{r}u\times\mathbf{r}v|=a$ → 面素 $dS=a,du,dv,[\mathrm{m^2}]$
  単位法線 $\mathbf{n}=(\cos u,\sin u,0)$（無次元）
  側面積 $S=\int{0}^{2\pi}!\int{z_1}^{z_2}a,dv,du=2\pi a L,[\mathrm{m^2}]$

(b) 球面 $\mathbf{r}(\phi,\theta)=(R\cos\phi\cos\theta,,R\cos\phi\sin\theta,,R\sin\phi)$

• $R$[m], $\phi\in[-\tfrac\pi2,\tfrac\pi2]$, $\theta\in[0,2\pi)$
  $|\mathbf{r}\phi\times\mathbf{r}\theta|=R^2\cos\phi$
  面素 $dS=R^2\cos\phi,d\phi,d\theta,[\mathrm{m^2}]$
  面積 $S=4\pi R^2,[\mathrm{m^2}]$
  単位法線 $\mathbf{n}=\mathbf{r}/R$（外向き，無次元）

(c) 回転面の一般式

$\mathbf{r}(u,v)=(x(u)\cos v,,x(u)\sin v,,z(u))$
$|\mathbf{r}_u\times\mathbf{r}_v|=x(u)\sqrt{x'(u)^2+z'(u)^2}$

$$
\boxed{\ S=2\pi\int x(u)\sqrt{x'(u)^2+z'(u)^2},du\ [\mathrm{m^2}] }
$$

5) 例題の計算（本のp.14–16に一致）

例題2 $\mathbf{r}=(u,,v,,f(u,v))$

$$
\mathbf{n}=\pm\frac{(-f_u,-f_v,1)}{\sqrt{1+f_u^2+f_v^2}}\ (\text{無次元})
$$

例題3 $\mathbf{r}(u,v)=(u\cos v,,u\sin v,,u^2),\ 0!\le!u!\le!1,\ 0!\le!v!\le!2\pi$

$$
|\mathbf{r}_u\times\mathbf{r}_v|=u\sqrt{4u^2+1}
\Rightarrow
S=\int_0^{2\pi}!!\int_0^1 u\sqrt{4u^2+1},du,dv
=\frac{\pi}{6}(5\sqrt5-1)\ \mathrm{[m^2]}
$$

問13 $\mathbf{r}=(\cos u,,\sin u,,v),\ 0!\le!u!\le!\pi,\ 0!\le!v!\le!2$

$|\mathbf{r}_u\times\mathbf{r}_v|=1\Rightarrow S=\int_0^\pi!!\int_0^2 1,dv,du=2\pi\ \mathrm{[m^2]}$

6) 工学応用（積分にそのまま代入）

• 磁束：$\displaystyle \Phi=\iint_S \mathbf{B}\cdot\mathbf{n},dS\ [\mathrm{Wb=T,m^2}]$
  例：円柱側面＋軸方向 $\mathbf{B}=B_z\hat{\mathbf{z}}$ → $\mathbf{n}\perp\hat{\mathbf{z}}$ で $\Phi=0$。

• 体積流量：$\displaystyle Q=\iint_A \mathbf{v}\cdot\mathbf{n},dS\ [\mathrm{m^3/s}]$
  ノズル・ダクトの入口面で $A$ と $\mathbf{n}$ を使用。

• 圧力合力：$\displaystyle \mathbf{F}=\iint_A p,\mathbf{n},dS\ [\mathrm{N}]$
  均一圧力なら投影面積の法則：$\mathbf{F}$ の大きさ $=p\times$ 投影面積 $[\mathrm{N}]$。

• 表面積見積：塗装量 $m=\rho_{\text{coat}},S\ [\mathrm{kg}]$、放熱 $Q=h,S,(T_s-T_\infty)\ [\mathrm{W}]$。

7) 単位の注意

• $\mathbf{r}$[m] → $\mathbf{r}_u,\mathbf{r}_v$[m/param] → $\mathbf{r}_u\times\mathbf{r}_v$[m²/(param²)]。
  $u,v$ が角（rad, 無次元）や長さ（m）でも 最終的に $dS$ は[m²] に必ずなる。
• $\mathbf{n}$ は比（ベクトル/その大きさ）なので 無次元。

# Program Name: surface_normals_area.py
# Creation Date: 20250821
# Overview: Compute unit normals, surface element dS, and total area for parametric surfaces r(u,v). Includes book Example 3 and common surfaces; plots the surface and sampled normals.
# Usage: Run as-is. Choose SURFACE and parameters in the "Parameters" section. The script prints area [m^2] and draws matplotlib figures.

!pip install numpy matplotlib -q

import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

# ============================================================
# Parameters / パラメータ（ここだけ変更すれば全体に反映） -------------- 日本語/English
# ============================================================
SURFACE = "book_example3"   # "book_example3" | "cylinder_side" | "sphere_patch" | "revolution" | "graph" | "custom"
NU = 200                    # u-direction samples / u方向サンプル数
NV = 300                    # v-direction samples / v方向サンプル数
PLOT_3D = True              # Draw 3D surface / 3D表示を行う
QUIVER_STRIDE = (20, 30)    # Normal arrows stride (nu_step, nv_step) / 法線矢印の間引き
# ---- Physical/constants for each surface / 各曲面の定数・物理量 ----
# Example 3: r(u,v)=(u cos v, u sin v, u^2), u∈[0,1], v∈[0,2π]
U_MIN, U_MAX = 0.0, 1.0
V_MIN, V_MAX = 0.0, 2.0*np.pi
# Cylinder side: r(u,v)=(a cos u, a sin u, v)
CYL_RADIUS_M = 1.0          # [m]
CYL_Z1_M, CYL_Z2_M = 0.0, 2.0  # [m]
# Sphere patch: r(φ,θ)=(R cosφ cosθ, R cosφ sinθ, R sinφ)
SPHERE_R_M = 1.0            # [m]
SPHERE_PHI_MIN, SPHERE_PHI_MAX = -np.pi/4, np.pi/4   # [rad]
SPHERE_THETA_MIN, SPHERE_THETA_MAX = 0.0, 2*np.pi     # [rad]
# Surface of revolution: r(u,v)=(x(u)cos v, x(u)sin v, z(u))
def X_of_u(u): return 1.0 + 0.2*np.cos(2*u)          # [m]
def Z_of_u(u): return 0.5*u                           # [m]
REV_U_MIN, REV_U_MAX = 0.0, 2*np.pi                   # [rad]
REV_V_MIN, REV_V_MAX = 0.0, 2*np.pi                   # [rad]
# Graph surface: r(u,v)=(u, v, f(u,v))
def f_graph(u, v): return 0.2*(u**2 + v**2)           # [m]
GRAPH_U_MIN, GRAPH_U_MAX = -1.0, 1.0                  # [m]
GRAPH_V_MIN, GRAPH_V_MAX = -1.0, 1.0                  # [m]
# Custom parametric surface r(u,v) user-defined / 任意のカスタム
def r_custom(u, v):
    # English: Example custom saddle surface
    # 日本語：鞍状の例
    a = 1.0
    return np.stack([u, v, (u**2 - v**2)/(2*a)], axis=-1)  # [m]

# ============================================================
# Surface selector / 曲面の選択と微分（解析式） ------------------- 日本語/English
# ============================================================
def r_uv(u, v):
    # Return r(u,v) in [m] / 位置ベクトル [m]
    if SURFACE == "book_example3":
        return np.stack([u*np.cos(v), u*np.sin(v), u**2], axis=-1)
    if SURFACE == "cylinder_side":
        a = CYL_RADIUS_M
        return np.stack([a*np.cos(u), a*np.sin(u), v], axis=-1)
    if SURFACE == "sphere_patch":
        R = SPHERE_R_M
        x = R*np.cos(u)*np.cos(v)
        y = R*np.cos(u)*np.sin(v)
        z = R*np.sin(u)
        return np.stack([x, y, z], axis=-1)
    if SURFACE == "revolution":
        x = X_of_u(u)
        z = Z_of_u(u)
        return np.stack([x*np.cos(v), x*np.sin(v), z], axis=-1)
    if SURFACE == "graph":
        return np.stack([u, v, f_graph(u, v)], axis=-1)
    if SURFACE == "custom":
        return r_custom(u, v)
    raise ValueError("Unknown SURFACE")

def ru(u, v):
    # ∂r/∂u [m/ unit(u)] / 解析微分
    if SURFACE == "book_example3":
        return np.stack([np.cos(v), np.sin(v), 2*u], axis=-1)
    if SURFACE == "cylinder_side":
        a = CYL_RADIUS_M
        return np.stack([-a*np.sin(u), a*np.cos(u), np.zeros_like(u)], axis=-1)
    if SURFACE == "sphere_patch":
        R = SPHERE_R_M
        return np.stack([-R*np.sin(u)*np.cos(v),
                         -R*np.sin(u)*np.sin(v),
                          R*np.cos(u)], axis=-1)
    if SURFACE == "revolution":
        x = X_of_u(u); z = Z_of_u(u)
        dx = (X_of_u(u + 1e-6) - X_of_u(u - 1e-6)) / (2e-6)  # safe numeric derivative
        dz = (Z_of_u(u + 1e-6) - Z_of_u(u - 1e-6)) / (2e-6)
        return np.stack([dx*np.cos(v), dx*np.sin(v), dz], axis=-1)
    if SURFACE == "graph":
        # fx = ∂f/∂u, fy = ∂f/∂v at fixed v / 数値微分
        du = 1e-6
        fx = (f_graph(u+du, v) - f_graph(u-du, v))/(2*du)
        return np.stack([np.ones_like(u), np.zeros_like(u), fx], axis=-1)
    if SURFACE == "custom":
        # Central difference generic / 一般の中央差分
        du = 1e-6
        return (r_custom(u+du, v) - r_custom(u-du, v))/(2*du)
    raise ValueError("Unknown SURFACE")

def rv(u, v):
    # ∂r/∂v [m/ unit(v)] / 解析微分
    if SURFACE == "book_example3":
        return np.stack([-u*np.sin(v), u*np.cos(v), np.zeros_like(u)], axis=-1)
    if SURFACE == "cylinder_side":
        return np.stack([np.zeros_like(u), np.zeros_like(u), np.ones_like(u)], axis=-1)
    if SURFACE == "sphere_patch":
        R = SPHERE_R_M
        return np.stack([-R*np.cos(u)*np.sin(v),
                          R*np.cos(u)*np.cos(v),
                          np.zeros_like(u)], axis=-1)
    if SURFACE == "revolution":
        x = X_of_u(u)
        return np.stack([-x*np.sin(v), x*np.cos(v), np.zeros_like(u)], axis=-1)
    if SURFACE == "graph":
        dv = 1e-6
        fy = (f_graph(u, v+dv) - f_graph(u, v-dv))/(2*dv)
        return np.stack([np.zeros_like(u), np.ones_like(u), fy], axis=-1)
    if SURFACE == "custom":
        dv = 1e-6
        return (r_custom(u, v+dv) - r_custom(u, v-dv))/(2*dv)
    raise ValueError("Unknown SURFACE")

def domain():
    # Parameter domain / パラメータ領域
    if SURFACE == "book_example3":
        return (U_MIN, U_MAX, V_MIN, V_MAX)
    if SURFACE == "cylinder_side":
        return (0.0, 2*np.pi, CYL_Z1_M, CYL_Z2_M)
    if SURFACE == "sphere_patch":
        return (SPHERE_PHI_MIN, SPHERE_PHI_MAX, SPHERE_THETA_MIN, SPHERE_THETA_MAX)
    if SURFACE == "revolution":
        return (REV_U_MIN, REV_U_MAX, REV_V_MIN, REV_V_MAX)
    if SURFACE == "graph":
        return (GRAPH_U_MIN, GRAPH_U_MAX, GRAPH_V_MIN, GRAPH_V_MAX)
    if SURFACE == "custom":
        return (-1.0, 1.0, -1.0, 1.0)
    raise ValueError("Unknown SURFACE")

# ============================================================
# Discretization / 離散化と面素・法線の計算 --------------------- 日本語/English
# ============================================================
u0, u1, v0, v1 = domain()
uu = np.linspace(u0, u1, NU)
vv = np.linspace(v0, v1, NV)
U, V = np.meshgrid(uu, vv, indexing="ij")  # U: (NU,NV), V: (NU,NV)

R = r_uv(U, V)                      # r(u,v) [m]
Ru = ru(U, V)                       # ∂r/∂u
Rv = rv(U, V)                       # ∂r/∂v
N = np.cross(Ru, Rv)                # Oriented area vector / 面素の向き [m^2/(u v)]
dS_density = np.linalg.norm(N, axis=2)   # ||ru×rv|| / 面素密度
# Unit normal (avoid 0) / 単位法線（0除算回避）
n_eps = np.where(dS_density==0, 1.0, dS_density)
Nhat = (N.transpose(2,0,1) / n_eps).transpose(1,2,0)  # 無次元

# Total area via Riemann sum / 面積（長方形近似）
du = (u1 - u0)/(NU - 1)
dv = (v1 - v0)/(NV - 1)
AREA_m2 = float(dS_density.sum() * du * dv)

# Closed-form reference for book example 3 / 例題3の厳密値
REF_m2 = None
if SURFACE == "book_example3":
    REF_m2 = (math.pi/6.0)*(5.0*math.sqrt(5.0) - 1.0)  # [m^2]

# ============================================================
# Print results / 結果の表示（物理単位を明記） -------------------- 日本語/English
# ============================================================
print("# Surface area result / 曲面積の結果")
print(f"S_numeric [m^2] : {AREA_m2:.10f}")
if REF_m2 is not None:
    err = AREA_m2 - REF_m2
    rel = err/REF_m2
    print(f"S_exact   [m^2] : {REF_m2:.10f}")
    print(f"Abs error [m^2] : {err:.3e}")
    print(f"Rel error   [-] : {rel:.3e}")

# ============================================================
# Plotting / 可視化（英語表記、単図、色指定なし） ---------------- 日本語/English
# ============================================================
if PLOT_3D:
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(R[:,:,0], R[:,:,1], R[:,:,2], alpha=0.7, linewidth=0, antialiased=True)
    # Quiver samples / 法線の矢印を間引いて描画
    step_u, step_v = QUIVER_STRIDE
    Uq = R[::step_u, ::step_v, 0]
    Vq = R[::step_u, ::step_v, 1]
    Zq = R[::step_u, ::step_v, 2]
    Nx = Nhat[::step_u, ::step_v, 0]
    Ny = Nhat[::step_u, ::step_v, 1]
    Nz = Nhat[::step_u, ::step_v, 2]
    # Scale of arrows / 矢印スケール
    span = max(R[:,:,0].ptp(), R[:,:,1].ptp(), R[:,:,2].ptp())
    ax.quiver(Uq, Vq, Zq, Nx, Ny, Nz, length=0.1*span, normalize=True)
    ax.set_title("Parametric Surface r(u,v) with Unit Normals")
    ax.set_xlabel("x [m]")
    ax.set_ylabel("y [m]")
    ax.set_zlabel("z [m]")
    plt.show()

# ============================================================
# Notes / 補足（日本語と英語）
# - dS_density = ||∂r/∂u × ∂r/∂v|| gives the surface element per (du,dv).
# - 面素密度 dS_density は「(du,dv) あたりの面積」。面積 S は Σ dS_density * du * dv。
# - Unit normal Nhat は法線の向きのみ（無次元）。物理積分では n*dS を使う。
# - 例題3の厳密面積 S = (π/6)(5√5 − 1) [m^2] と数値結果の相対誤差を表示。
# ============================================================

練習問題1（p.17）詳解

前提：座標を長さと解釈すれば $ \mathbf{r}$[m]、角度は無次元（rad）です。面積は $[\mathrm{m^2}]$。

---

1.

$\mathbf{a}+\mathbf{b}+\mathbf{c}=0 \Rightarrow \mathbf{a}=-(\mathbf{b}+\mathbf{c})$

$$
\mathbf{a}\times\mathbf{b}=(-\mathbf{b}-\mathbf{c})\times\mathbf{b}
=-(\mathbf{b}\times\mathbf{b})-(\mathbf{c}\times\mathbf{b})
=\mathbf{b}\times\mathbf{c}
$$

同様に $\mathbf{c}\times\mathbf{a}=\mathbf{b}\times\mathbf{c}$。
$\boxed{\mathbf{a}\times\mathbf{b}=\mathbf{b}\times\mathbf{c}=\mathbf{c}\times\mathbf{a}}$（無次元）

---

2.

$\mathbf{a}=(a_1,a_2,a_3),\ \mathbf{b}=(b_1,b_2,b_3),\ \mathbf{c}=(c_1,c_2,c_3)$

(1)

$$
\mathbf{a}\cdot(\mathbf{b}\times\mathbf{c})
=\begin{vmatrix}
a_1&a_2&a_3\
b_1&b_2&b_3\
c_1&c_2&c_3
\end{vmatrix}
$$

（混合積＝3×3行列式）

(2)

$$
\boxed{\ \mathbf{a}\cdot(\mathbf{b}\times\mathbf{c})
=\mathbf{b}\cdot(\mathbf{c}\times\mathbf{a})
=\mathbf{c}\cdot(\mathbf{a}\times\mathbf{b})\ }
$$

（循環置換で不変）

---

3. A(0,1,2), B(-1,2,3), C(2,-2,2)

$$
\overrightarrow{AB}=B-A=(-1,1,1),\quad
\overrightarrow{AC}=C-A=(2,-3,0)
$$

$$
\overrightarrow{AB}\times\overrightarrow{AC}
=\begin{vmatrix}
\mathbf{i}&\mathbf{j}&\mathbf{k}\
-1&1&1\
2&-3&0
\end{vmatrix}
=(3,2,1)
$$

(1) 平行四辺形の面積

$$
\boxed{S=|\overrightarrow{AB}\times\overrightarrow{AC}|
=\sqrt{3^2+2^2+1^2}=\sqrt{14}}\ [\mathrm{m^2}]
$$

(2) 両方に垂直な単位ベクトル

$$
\boxed{\ \pm\frac{(3,2,1)}{\sqrt{14}}\ }\ (\text{無次元})
$$

---

4. $\mathbf{r}(t)=(e^t,e^{-t},\sqrt2,t)$

$$
\mathbf{r}'(t)=(e^t,-e^{-t},\sqrt2),\quad
|\mathbf{r}'(t)|=\sqrt{e^{2t}+e^{-2t}+2}=\boxed{e^t+e^{-t}}\ [\mathrm{m/s}]
$$

(1) 単位接線

$$
\boxed{\ \mathbf{t}(t)=\dfrac{(e^t,-e^{-t},\sqrt2)}{e^t+e^{-t}}\ }\ (\text{無次元})
$$

(2) 弧長 $t:0\to1$

$$
\boxed{\ s=\int_0^1 (e^t+e^{-t}),dt
=\bigl[e^t-e^{-t}\bigr]_0^1
=e-\tfrac1e\ }\ [\mathrm{m}]
$$

---

5. $\mathbf{r}(t)=(2t,\ \sqrt3,t^{2},\ t^{3})$（$t\ge0$）

$$
\mathbf{r}'(t)=(2,\ 2\sqrt3,t,\ 3t^{2}),\quad
|\mathbf{r}'(t)|=\sqrt{4+12t^{2}+9t^{4}}=\boxed{,2+3t^{2},}\ [\mathrm{m/s}]
$$

弧長 $t:0\to3$

$$
\boxed{\ s=\int_0^3 (2+3t^{2}),dt
=\bigl[2t+t^{3}\bigr]_0^3
=33 }\ [\mathrm{m}]
$$

---

6. $a>0$, $\mathbf{r}(u,v)=(a\cos u\sin v,\ a\sin u\sin v,\ a\cos v)$,

$D:\ 0\le u\le \pi/2,\ \pi/3\le v\le \pi/2$（半径 $a$ の球の一部）

(1) 単位法線
球面では外向き法線 $=\mathbf{r}/a$

$$
\boxed{\ \mathbf{n}(u,v)=\bigl(\cos u\sin v,\ \sin u\sin v,\ \cos v\bigr)\ }\ (\text{無次元})
$$

(2) 面積要素 $|\mathbf{r}_u\times\mathbf{r}_v|=a^{2}\sin v$ より

$$
S=\int_{0}^{\pi/2}!!\int_{\pi/3}^{\pi/2} a^{2}\sin v,dv,du
=a^{2}\Bigl(\tfrac{\pi}{2}\Bigr)\Bigl[-\cos v\Bigr]_{\pi/3}^{\pi/2}
=\boxed{\ \tfrac{\pi}{4},a^{2}\ }\ [\mathrm{m^2}]
$$

---

§2 スカラー場とベクトル場：勾配を“徹底的に” （単位つき）

1) ナブラと勾配の本質

• スカラー場 $\phi(x,y,z)$（例：温度 $T[\mathrm{K}]$、圧力 $p[\mathrm{Pa}]$、電位 $V[\mathrm{V}]$）

• ベクトル場 $\mathbf{a}(x,y,z)$（例：流速 $\mathbf{u}[\mathrm{m/s}]$、電場 $\mathbf{E}[\mathrm{V/m}]$）

• ナブラ

  $$
  \nabla=\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z}\right)
  $$

• 勾配（グラディエント）

  $$
  \boxed{\nabla\phi=\left(\frac{\partial\phi}{\partial x},\frac{\partial\phi}{\partial y},\frac{\partial\phi}{\partial z}\right)}
  \quad[\ \text{単位}=\phi/\mathrm{m}\ ]
  $$

“最急上昇”の意味

任意の単位ベクトル $\mathbf{e}$ 方向の変化率（方向微分）

$$
\boxed{D_{\mathbf{e}}\phi=\lim_{h\to0}\frac{\phi(\mathbf{r}+h\mathbf{e})-\phi(\mathbf{r})}{h}
=\nabla\phi\cdot\mathbf{e}}
$$

$\mathbf{e}$ と $\nabla\phi$ のなす角を $\theta$ とすると

$$
D_{\mathbf{e}}\phi=|\nabla\phi|\cos\theta
$$

ゆえに最大増加率は $|\nabla\phi|$（方向は $\nabla\phi$ ）、最小は $-|\nabla\phi|$（方向 $-\nabla\phi$）。

2) 等位面（等高面）と法線

等位面 $\phi(x,y,z)=c$ 上の曲線 $\mathbf{r}(t)$ を考える：

$$
\frac{d}{dt}\phi(\mathbf{r}(t))=\nabla\phi\cdot\frac{d\mathbf{r}}{dt}=0
\Rightarrow \boxed{\ \nabla\phi\ \perp\ \text{接ベクトル}\ }
$$

したがって $\nabla\phi$ は等位面の法線ベクトル（単位は $\phi/\mathrm{m}$）。
点 $\mathbf{r}_0$ での接平面は

$$
\nabla\phi(\mathbf{r}_0)\cdot(\mathbf{r}-\mathbf{r}_0)=0
$$

全微分（線形近似）

$$
d\phi=\nabla\phi\cdot d\mathbf{r}
\quad[\ \phi\ ]
$$

3) 代表座標系での勾配

• 円柱座標 $(r,\theta,z)$：

  $$
  \boxed{\ \nabla\phi=\frac{\partial\phi}{\partial r},\mathbf{e}r
  +\frac{1}{r}\frac{\partial\phi}{\partial\theta},\mathbf{e}\theta
  +\frac{\partial\phi}{\partial z},\mathbf{e}_z\ }
  $$

• 球座標 $(r,\theta,\varphi)$（$\theta$：極角，$\varphi$：方位角）：

  $$
  \boxed{\ \nabla\phi=\frac{\partial\phi}{\partial r},\mathbf{e}r
  +\frac{1}{r}\frac{\partial\phi}{\partial \theta},\mathbf{e}\theta
  +\frac{1}{r\sin\theta}\frac{\partial\phi}{\partial \varphi},\mathbf{e}_\varphi\ }
  $$

いずれも単位は $\phi/\mathrm{m}$。

4) 二階微分：ヘッセ行列と曲率

$$
H(\phi)=\begin{bmatrix}
\phi_{xx}&\phi_{xy}&\phi_{xz}\
\phi_{yx}&\phi_{yy}&\phi_{yz}\
\phi_{zx}&\phi_{zy}&\phi_{zz}
\end{bmatrix}\quad[\ \phi/\mathrm{m^2}\ ]
$$

方向 $\mathbf{e}$ での“二階の変化率”（方向2階微分）

$$
D_{\mathbf{e}}^{(2)}\phi=\mathbf{e}^{\mathsf T}H(\phi),\mathbf{e}\quad[\ \phi/\mathrm{m^2}\ ]
$$

極大・極小・鞍点判定や等位面の局所曲率評価に使う。

5) 物理・工学へのダイレクト接続（式と単位）

• 熱伝導（フーリエ）：$\mathbf{q}=-k,\nabla T$
  $\mathbf{q}[\mathrm{W/m^2}],\ k[\mathrm{W/(m,K)}],\ \nabla T[\mathrm{K/m}]$
• 静電場：$\mathbf{E}=-\nabla V$（$[\mathrm{V/m}]$）
• 流体の圧力勾配力：体積力密度 $-\nabla p,[\mathrm{N/m^3}]$
• 重力ポテンシャル：$\mathbf{g}=-\nabla\Phi$（$[\mathrm{m/s^2}]$）
• 最適化・機械学習：勾配降下
  $\mathbf{x}_{k+1}=\mathbf{x}_k-\eta,\nabla f(\mathbf{x}_k)$（$\eta$：学習率）
• 拡散方程式：$\displaystyle \frac{\partial \phi}{\partial t}=D,\Delta\phi$, $\Delta=\nabla\cdot\nabla$（ラプラシアン $[\phi/\mathrm{m^2}]$）

6) “勾配の線積分” と保存力

$$
\int_{\mathcal{C}}\nabla\phi\cdot d\mathbf{r}=\phi(\mathbf{r}_B)-\phi(\mathbf{r}_A)
$$

→ 経路に依らず始終点だけで決まる（勾配場は保存場）。
ベクトル場 $\mathbf{F}$ が $\mathbf{F}=\nabla\phi$ を持つ ⇔ $\nabla\times\mathbf{F}=\mathbf{0}$（単連結領域）。

7) 手を動かすミニ例（単位そろえ）

例1：温度場と熱流

$T(x,y,z)=300+5x-2z,[\mathrm{K}]$（$x,z[\mathrm{m}]$）

$$
\nabla T=(5,0,-2)\ [\mathrm{K/m}]
$$

銅の $k=200,[\mathrm{W/(m,K)}]$ ⇒
$\mathbf{q}=-k\nabla T=(-1000,0,400),[\mathrm{W/m^2}]$

例2：等位面の接平面

$\phi=x^2+y^2+z$（無次元）・点 $(1,2,3),[\mathrm{m}]$

$$
\nabla\phi=(2x,2y,1)\Big|_{(1,2,3)}=(2,4,1)
$$

接平面：$2(x-1)+4(y-2)+1(z-3)=0$

例3：方向微分と最急上昇

上の点で $\mathbf{e}=\frac{1}{\sqrt2}(1,1,0)$ 方向

$$
D_{\mathbf{e}}\phi=\nabla\phi\cdot\mathbf{e}=\frac{2+4}{\sqrt2}=\frac{6}{\sqrt2}\ [1/\mathrm{m}]
$$

最大増加率は $|\nabla\phi|=\sqrt{2^2+4^2+1^2}=\sqrt{21}\ [1/\mathrm{m}]$。

8) ありがちな落とし穴

• スケールの混在：座標の単位を揃えないと $\nabla\phi$ の単位が崩れる。
• 非直交座標の基底：円柱・球座標では $\mathbf{e}\theta,\mathbf{e}\varphi$ が位置依存。式をそのまま使う。
• ノイズ：数値微分はノイズを増幅。スムージング or 物理モデルの正則化（例：$\Delta\phi$ 罰則）を併用。

9) まとめ（現場で使う最小セット）

• 勾配：$\nabla\phi$（単位 $\phi/\mathrm{m}$）
• 法線：等位面に直交（接平面：$\nabla\phi(\mathbf{r}_0)\cdot(\mathbf{r}-\mathbf{r}_0)=0$）
• 方向微分：$D_{\mathbf{e}}\phi=\nabla\phi\cdot\mathbf{e}$
• 保存場：$\mathbf{F}=\nabla\phi \Rightarrow \int\mathbf{F}\cdot d\mathbf{r}=\phi(B)-\phi(A)$
• 物理：$\mathbf{q}=-k\nabla T,\ \mathbf{E}=-\nabla V,\ -\nabla p$ など

# Program Name: gradient_toolkit.py
# Creation Date: 20250821
# Overview: Compute ∇phi, directional derivative, tangent plane of level sets, and verify path-independence for gradient fields. Includes Cartesian/Cylindrical/Spherical formulae and simple plots.
# Usage: Run as-is. Adjust parameters in the "Parameters" section. The script prints gradients, directional derivatives, tangent-plane equation, and optionally draws quiver plots.

!pip install numpy sympy matplotlib -q

# ================================
# Parameters / パラメータ（一元管理）
# ================================
DO_PLOT = True                 # 2D可視化を行うか / whether to plot
GRID_N = 25                    # 可視化グリッド密度 / grid density for plotting
Z_SLICE = 0.0                  # z=定数の2Dスライス [m] / z-slice [m]
EXAMPLE = "poly3d"             # "poly3d" | "temperature" | "custom"
# 無次元多項式場 / dimensionless scalar field
#   phi(x,y,z) = x^2 + y^2 + z
# 温度場 / temperature field
#   T(x,y,z) = 300 + 5x - 2z [K]
# カスタムは下の custom_phi で定義 / custom field defined below
POINT = (1.0, 2.0, 3.0)        # 評価点 [m] / evaluation point [m]
E_DIR = (1/2**0.5, 1/2**0.5, 0.0)  # 方向微分の単位ベクトル / unit direction for D_e phi
DOMAIN_X = (-2.0, 2.0)         # x範囲 [m] / x-range [m]
DOMAIN_Y = (-2.0, 2.0)         # y範囲 [m] / y-range [m]

# =====================================
# Imports / インポート
# =====================================
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# =====================================
# Scalar fields / スカラー場の定義（関数）
# =====================================
# 日本語: カスタム場の定義（単位はユーザが意味づけ） / English: define your custom scalar
def custom_phi(x, y, z):
    # 例: ガウス分布（無次元） / Example: dimensionless Gaussian
    return sp.exp(-(x**2 + y**2 + z**2))

# 日本語: 場の選択に応じたシンボリック関数を返す / English: return sympy expression of phi
def select_phi(EXAMPLE):
    x, y, z = sp.symbols('x y z', real=True)
    if EXAMPLE == "poly3d":
        phi = x**2 + y**2 + z          # 無次元 / dimensionless
        unit = "1/m"                   # ∇phi の単位（phiが無次元→1/m） / grad unit
        return phi, unit
    if EXAMPLE == "temperature":
        # T = 300 + 5x - 2z [K]; x,z [m]
        phi = 300 + 5*x - 2*z
        unit = "K/m"                   # 勾配の単位 / grad unit
        return phi, unit
    if EXAMPLE == "custom":
        phi = custom_phi(x, y, z)
        unit = "1/m"
        return phi, unit
    raise ValueError("unknown EXAMPLE")

# =====================================
# Gradient, directional derivative, tangent plane
# 勾配・方向微分・接平面（シンボリック→数値）
# =====================================
x, y, z = sp.symbols('x y z', real=True)
phi_expr, grad_unit = select_phi(EXAMPLE)
grad_expr = (sp.diff(phi_expr, x), sp.diff(phi_expr, y), sp.diff(phi_expr, z))  # ∇phi

# 日本語: 数値評価関数 / English: numeric lambdas
phi_f = sp.lambdify((x, y, z), phi_expr, "numpy")
gx_f  = sp.lambdify((x, y, z), grad_expr[0], "numpy")
gy_f  = sp.lambdify((x, y, z), grad_expr[1], "numpy")
gz_f  = sp.lambdify((x, y, z), grad_expr[2], "numpy")

# =====================================
# Evaluate at a point / 点での評価
# =====================================
px, py, pz = POINT
gx = float(gx_f(px, py, pz))
gy = float(gy_f(px, py, pz))
gz = float(gz_f(px, py, pz))
grad_vec = np.array([gx, gy, gz], dtype=float)

# 日本語: 方向微分 D_e phi = ∇phi ⋅ e / English: directional derivative
e = np.array(E_DIR, dtype=float)
e = e / np.linalg.norm(e)  # 念のため正規化 / ensure unit
D_e = float(np.dot(grad_vec, e))

# 日本語: 接平面の式 ∇phi(r0)⋅(r - r0) = 0 / English: tangent plane to level surface
# ここで c = phi(r0) / level set value at the point
phi0 = float(phi_f(px, py, pz))
# 平面の標準形: gx*(x-px)+gy*(y-py)+gz*(z-pz)=0
plane_coeff = (gx, gy, gz, -(gx*px + gy*py + gz*pz))  # ax+by+cz+d=0 の (a,b,c,d)

# =====================================
# Print results / 結果表示（単位明記）
# =====================================
print("# Gradient / 勾配")
print(f"phi(x,y,z) = {sp.simplify(phi_expr)}")
print(f"∇phi = ({sp.simplify(grad_expr[0])}, {sp.simplify(grad_expr[1])}, {sp.simplify(grad_expr[2])})  [{grad_unit}]")
print(f"At P=({px:.3f},{py:.3f},{pz:.3f}) [m]: ∇phi(P) = ({gx:.6g}, {gy:.6g}, {gz:.6g}) [{grad_unit}]")

print("\n# Directional derivative / 方向微分")
print(f"e = ({e[0]:.6g}, {e[1]:.6g}, {e[2]:.6g}) [-]")
print(f"D_e phi(P) = {D_e:.6g} [{grad_unit}]")

print("\n# Tangent plane of level surface at P / 等位面の接平面（点P）")
print(f"phi(P) = {phi0:.6g}")
a,b,c,d = plane_coeff
print(f"Plane: {a:.6g}*x + {b:.6g}*y + {c:.6g}*z + {d:.6g} = 0  (normal = ∇phi(P))")

# =====================================
# Cylindrical and spherical gradient formulae (report only)
# 円柱・球座標の勾配（式を表示）
# =====================================
r, th, zz = sp.symbols('r theta z', real=True, positive=True)
phi_cyl = sp.Function('phi')(r, th, zz)
grad_cyl = (sp.diff(phi_cyl, r),
            (1/r)*sp.diff(phi_cyl, th),
            sp.diff(phi_cyl, zz))
print("\n# Cylindrical coordinates gradient / 円柱座標の勾配（基底 e_r,e_theta,e_z）")
print(f"∇phi = (∂phi/∂r) e_r + (1/r)(∂phi/∂θ) e_θ + (∂phi/∂z) e_z  [phi/m]")

rr, tht, ph = sp.symbols('r theta phi', real=True, positive=True)
phi_sph = sp.Function('phi')(rr, tht, ph)
grad_sph = (sp.diff(phi_sph, rr),
            (1/rr)*sp.diff(phi_sph, tht),
            (1/(rr*sp.sin(tht)))*sp.diff(phi_sph, ph))
print("# Spherical coordinates gradient / 球座標の勾配（基底 e_r,e_theta,e_phi）")
print(f"∇phi = (∂phi/∂r) e_r + (1/r)(∂phi/∂θ) e_θ + (1/(r sinθ))(∂phi/∂φ) e_φ  [phi/m]")

# =====================================
# Path-independence check for gradient field
# 勾配場の線積分（経路非依存の確認）
# =====================================
A = np.array([0.0, 0.0, 0.0])
B = np.array([1.0, 1.0, 0.0])

def line_integral_grad_along_path(param_to_xyz, t0, t1, N=1000):
    # 数値線積分 ∫ ∇phi · dr / numerical line integral
    ts = np.linspace(t0, t1, N+1)
    integral = 0.0
    for i in range(N):
        t_mid = 0.5*(ts[i]+ts[i+1])
        P = np.array(param_to_xyz(t_mid), dtype=float)
        gradP = np.array([gx_f(*P), gy_f(*P), gz_f(*P)], dtype=float)
        dP = np.array(param_to_xyz(ts[i+1]), dtype=float) - np.array(param_to_xyz(ts[i]), dtype=float)
        integral += float(np.dot(gradP, dP))
    return integral

# 経路1: 直線 / straight line
def path1(t): return A + t*(B - A)
# 経路2: 折れ線 (x方向→y方向) / polyline
def path2(t):
    if t <= 0.5:
        return np.array([2*t, 0.0, 0.0])
    else:
        s = 2*(t-0.5)
        return np.array([1.0, s, 0.0])

I1 = line_integral_grad_along_path(path1, 0.0, 1.0, N=2000)
I2 = line_integral_grad_along_path(path2, 0.0, 1.0, N=2000)
phiA = float(phi_f(*A)); phiB = float(phi_f(*B))

print("\n# Line integral of gradient / 勾配の線積分（経路非依存の確認）")
print(f"Path1 integral = {I1:.6g}  [phi]")
print(f"Path2 integral = {I2:.6g}  [phi]")
print(f"phi(B)-phi(A) = {phiB-phiA:.6g}  [phi]")

# =====================================
# 2D Quiver plot at z = Z_SLICE / zスライスでの勾配ベクトルの矢印図
# =====================================
if DO_PLOT:
    xs = np.linspace(DOMAIN_X[0], DOMAIN_X[1], GRID_N)
    ys = np.linspace(DOMAIN_Y[0], DOMAIN_Y[1], GRID_N)
    Xg, Yg = np.meshgrid(xs, ys, indexing="xy")
    Zg = np.full_like(Xg, Z_SLICE, dtype=float)
    Gx = gx_f(Xg, Yg, Zg)
    Gy = gy_f(Xg, Yg, Zg)

    plt.figure(figsize=(6, 5))
    plt.quiver(Xg, Yg, Gx, Gy, angles="xy", scale_units="xy", scale=1)
    plt.title("Gradient Field on z = const")
    plt.xlabel("x [m]")
    plt.ylabel("y [m]")
    plt.axis("equal")
    plt.tight_layout()
    plt.show()

# =====================================
# Notes / 補足（日本語・English）
# 日本語: ∇phi の単位は（phiの単位）/m。方向微分 D_e phi も同じ単位。接平面の法線は ∇phi(P)。
# English: Gradient unit is (unit of phi)/m. Directional derivative has same unit. Plane normal equals ∇phi at the point.
# =====================================