はじめに

スプライン曲線の区分多項式の次数が３次である理由について、物理量の最適化の観点から導出されることを記載してみる。

今回の結論

　　データ点を通る曲線を物理的な連続体であるとした場合、

　　　　　曲げによる弾性エネルギーを最小にする曲線

　　　は ３次多項式による曲線 である。

設定課題式の導出

連続体に蓄えられる曲げによる弾性エネルギー（概略）

連続体を「梁」とする。

梁全体の弾性エネルギー

梁に荷重がかかった時、荷重による微小区間の中立軸から $y$ 離れた場所のひずみ $\varepsilon$ は

\begin{eqnarray}
dx &=& \rho d \theta \\
\varDelta &=& (\rho + y)d \theta - \rho d\theta = y d\theta \\
\\
\varepsilon &=& \frac{\varDelta}{dx} = \frac{yd\theta}{\rho d\theta} = \frac{y}{\rho}

\end{eqnarray}

応力 $\sigma$ は微小区間に与えられたモーメント $M$ と比例するので、断面二次モーメントを $I$ として位置 y の応力は

\sigma = \frac{M}{I}y

応力とひずみに対するフックの法則は

\sigma = E \varepsilon

したがって、

\varepsilon = \frac{\sigma}{E} = \frac{M}{EI}y

２つの $\varepsilon$ の式より

\frac{1}{\rho} = \frac{M}{EI}

これらより、

d\theta = \frac{dx}{\rho} = \frac{M}{EI}dx

荷重による外力が加える微小区間に対する仕事は、すなわちこの微小区間に蓄えられる弾性エネルギーなので

\begin{eqnarray}
dU &=& \frac{1}{2}M d\theta \\
   &=& \frac{M^2}{2EI}
\end{eqnarray}

これを梁全体に対して積分すると、梁に蓄えられる弾性エネルギーが求まる。

U = \int_L dU = \int_L \frac{M^2}{2EI} dx

弾性曲線方程式

たわみ曲線の微小区間 $ds$ と曲率半径 $\rho$、微小角 $d\theta$ の関係は

ds = \rho (-d\theta)

マイナス（－）は反時計回りに $d\theta$ が増加することによる。
たわみ曲線の傾きは

\frac{dy}{dx} = tan \theta

たわみによる水平線からのたわみ角は微小であるとして、以下の近似を使用する。

\begin{eqnarray}
ds &\fallingdotseq& dx \\
tan \theta &\fallingdotseq& \theta
\end{eqnarray}

これらより、

\begin{eqnarray}
\kappa = \frac{1}{\rho} 
 &=& - \frac{d \theta}{ds} \\
 &\fallingdotseq& - \frac{d}{dx} tan(\theta) \\
 &=& - \frac{d}{dx} \left (\frac{d y}{dx} \right ) \\
 &=& - \frac{d^2 y}{dx^2}

\end{eqnarray}

先の結果とあわせて、

\frac{d^2 y}{dx^2} = - \frac{M}{EI}

（別導出）

\begin{eqnarray}
\frac{1}{\rho} &=& - \frac{d\theta}{ds} = - \frac{d\theta}{dx} \frac{dx}{ds} \\
\end{eqnarray}

\begin{eqnarray}
tan(\theta) &=& \frac{dy}{dx} \\
\\
\therefore \quad \theta &=& tan^{-1} \left ( \frac{dy}{dx} \right ) \\
\\
\frac{d\theta}{dx} &=& \frac{\frac{d^2 y}{dx^2}}{1+ {\left( \frac{dy}{dx} \right )}^2}
\end{eqnarray}

\begin{eqnarray}

\frac{dx}{ds} = cos(\theta) = \frac{1}{\sqrt{1+tan^{2}(\theta)}} = \frac{1}{\sqrt{1+ {\left ( \frac{dy}{dx} \right )}^2}}

\end{eqnarray}

\begin{eqnarray}

\frac{1}{\rho} &=& - \frac{\frac{d^2 y}{dx^2}}{1+ {\left( \frac{dy}{dx} \right )}^2} \cdot \frac{1}{\sqrt{1+ {\left ( \frac{dy}{dx} \right )}^2}} \\

&=& - \frac{\frac{d^2 y}{dx^2}}{\left \{ 1+ {\left( \frac{dy}{dx} \right )}^2 \right \}^{\frac{3}{2}}} 

\end{eqnarray}

ここで、たわみ角 $\theta$ が小さいとすると、

\frac{dy}{dx} = tan(\theta) \quad \thickapprox \quad 0

したがって、

\begin{eqnarray}
\frac{1}{\rho} &\fallingdotseq& - \frac{d^2y}{dx^2} \\
\\
\therefore \quad \frac{d^2 y}{dx^2} &=& - \frac{M}{EI}
\end{eqnarray}

【今回の課題】最小化する対象式

梁の弾性エネルギー式に弾性曲線方程式を代入すると、

\begin{eqnarray}
U &=& \int_L \frac{M^2}{2EI} dx \\
  &=& \int_L \left \{ \frac{1}{2EI} \cdot E^2 I^2 \left ( \frac{d^2y}{dx^2} \right )^2 \right \} dx \\
  &=& \frac{EI}{2} \int_L \left ( \frac{d^2y}{dx^2} \right )^2 dx

\end{eqnarray}

この $U$ を最小化する関数 $y(x)$ を求めるのが今回の課題。

変分法による解法

最小化する最適関数 $y(x)$ は大きく３つの範囲に分けられる。
データ点を $p_0 \thicksim p_n$ とすると、
　
　・データ点 $p_0$ に到達するまで
　・データ点範囲内 $p_0 \thicksim p_n$
　・データ点 $p_n$ から先

の３つの範囲である。
x 軸の範囲で示せば、

( - \infty, x_0) \\
[x_0, x_n] \\
(x_n, \infty)

である。

１．データ点範囲内の関数

データ点 $p_0 \thicksim p_n$ の範囲（ x の範囲 $[x_0, x_n]$）において、エネルギー $U$ を最小化する関数 $y(x)$ を求める。
スプライン曲線は、データ点の間、区分ごとに異なる多項式を用いる。
よって、全区間を一括して表す単一関数 $y(x)$ は存在しない。
区分ごとに多項式を求めてつなぎ合わせる。
したがって、エネルギー $U$ を最小化する関数 $y(x)$ も区分ごとに違い、それぞれで求める事になる。

関数 $y_0(x)$ がエネルギー $U$ を最小化する関数であった場合、この関数 $y_0(x)$ が微小に変化しても（わずかに違う関数に変わっても）エネルギー $U$ はほとんど変化しないはずである（停留点）。

関数 $y_0(x)$ をわずかに変化させた関数を

y(x) = y_0(x) + \varepsilon \eta(x)

とする。
関数の両端での $y_0(x)$ の値、$y_0(x_j), y_0(x_{j+1})$ が変わってしまうと意味がなくなるので、$\eta(x)$ の境界条件を以下とする。

\eta(x_j) = \eta(x_{j+1}) = 0 \\
\eta'(x_j) = \eta'(x_{j+1}) = 0

この条件により、$y_0(x)$ に足しても $x_j, x_{j+1}$ では関数の値は変化しない。
また、この両端において接線の傾きも変化しない。
この境界条件を満たす限り、$\eta(x)$ は任意である。
$\varepsilon$ は任意の小さい定数である。

次に、

U = \frac{EI}{2} \int_{x_j}^{x_{j+1}} \left ( \frac{d^2y}{dx^2} \right )^2 dx

において、

f(y'') = \left ( \frac{d^2y}{dx^2} \right )^2 = (y'')^2

とおき、

I[y] = \int_{x_j}^{x_{j+1}} f(y'') dx

とすると、$U$ を最小化することと $I$ を最小化することは同じである。
$I[y]$ は汎関数である。

汎関数 $I[y]$ の微分を考える。

y'' = y_0'' + \varepsilon \eta''

なので、

\begin{eqnarray}
\frac{dI}{d\varepsilon} &=& \int_{x_j}^{x_{j+1}} \frac{\partial f}{\partial y''} \frac{d y''}{d \varepsilon} dx \\
&=& \int_{x_j}^{x_{j+1}}  \frac{\partial f}{\partial y''} \eta''(x) dx \\
&=& \left [ \frac{\partial f}{\partial y''} \eta'(x)  \right ]_{x_j}^{x_{j+1}} - \int_{x_j}^{x_{j+1}} \frac{d}{dx} \left ( \frac{\partial f}{\partial y''} \right ) \eta'(x) dx \\
&=& \left [ \frac{\partial f}{\partial y''} \eta'(x)  \right ]_{x_j}^{x_{j+1}}
 - \left [ \frac{d}{dx} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x)  \right ]_{x_j}^{x_{j+1}}
+ \int_{x_j}^{x_{j+1}} \frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x) \; dx \\

\end{eqnarray}

$\eta(x)$ の境界条件を適用すると前２項がゼロになるので

\frac{dI}{d\varepsilon} = \int_{x_j}^{x_{j+1}} \frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x) \; dx

これをゼロにするような $y(x)$ を求めればよい。

\begin{eqnarray}
\frac{dI}{d\varepsilon} &=& \int_{x_j}^{x_{j+1}} \frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x) \; dx \\
&=& 0
\end{eqnarray}

$\eta(x)$ は境界条件以外任意なので、これが恒等的にゼロになるためには、

\frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) = 0

である必要がある。
これを解けばエネルギー $U$ を最小にする関数 $y_0(x)$ が求まる。

これを積分すると、

\begin{eqnarray}
\frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) &=& 0 \\
\frac{d}{dx} \left ( \frac{\partial f}{\partial y''} \right ) &=& C_1 \\
\left ( \frac{\partial f}{\partial y''} \right ) &=& C_1x + C_2 \\
\end{eqnarray}

ここで $f(y'') = (y'')^2$ を元に戻すと、

\begin{eqnarray}
\frac{d}{d y''}(y'')^2 &=& C_1x + C_2 \\
2y'' &=& C_1x + C_2 \\
\frac{d^2 y}{dx^2} &=& C_3 x + C_4 \\
\frac{d y}{dx} &=& C_5 x^2 + C_6 x + C_7 \\
\\
y(x) &=& C_8 x^3 + C_9 x^2 + C_{10} x + C_{11} \\
\end{eqnarray}

よって、$y(x)$ は３次多項式である。

以上により、各区分において

　エネルギー $U$ を最小にする区分多項式 $y(x)$ は３次式

であることが分かった。

２．データ点`p_0`までの範囲の関数

x の範囲 $(-\infty, x_0)$ において、エネルギー $U$ を最小にする関数 $y(x)$ を求める。
同じく、$y_0(x)$ がエネルギー $U$ を最小にする関数であるとして、関数を微小変化させる。
今回は範囲の左端が $-\infty$ であるので境界条件を課すのは $x_0$ の点においてだけであり、

\eta(x_0) = 0 \\
\eta'(x_0) = 0

とする。
先と同じく汎関数 $I[y]$ の微分を考える。

\begin{eqnarray}
\frac{dI}{d\varepsilon} &=& \int_{- \infty}^{x_0} \frac{\partial f}{\partial y''} \frac{d y''}{d \varepsilon} dx \\
&=& \left [ \frac{\partial f}{\partial y''} \eta'(x)  \right ]_{- \infty}^{x_0}
 - \left [ \frac{d}{dx} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x)  \right ]_{- \infty}^{x_0}
+ \int_{- \infty}^{x_0} \frac{d^2}{dx^2} \left ( \frac{\partial f}{\partial y''} \right ) \eta(x) \; dx \\

\end{eqnarray}

ここで今回の境界条件を用いてもどの項もゼロにはならない。
$- \infty$ で $\eta(x)$ がゼロでないからである。
よって、この微分が恒等的にゼロになるためには、

\begin{eqnarray}
\frac{\partial f}{\partial y''} = 0
\end{eqnarray}

である必要がある。

$f(y'') = (y'')^2$ を戻すと、

\begin{eqnarray}
\frac{d}{d y''} (y'')^2 &=& 0 \\
2(y'') &=& 0 \\
y' &=& C_1 \\
\\
y(x) &=& C_1 x + C_2
\end{eqnarray}

よって、$y(x)$ は１次多項式である。

以上により、最初のデータ点 $p_0$ まで（$(- \infty, x_0)$ の範囲）において、

　エネルギー $U$ を最小にする区分多項式 $y(x)$ は１次式

であることが分かった。

３．データ点`p_n`から先の範囲の関数

x の範囲 $(x_n, \infty)$ において、エネルギー $U$ を最小にする関数 $y(x)$ を求める。
同じく、$y_0(x)$ がエネルギー $U$ を最小にする関数であるとして、関数を微小変化させる。
今回は範囲の右端が $\infty$ であるので境界条件を課すのは $x_n$ の点においてだけであり、

\eta(x_n) = 0 \\
\eta'(x_n) = 0

とする。

$(-\infty, x_0)$ の時と同じく、汎関数 $I[y]$ の微分においてどの項もゼロにならないので、恒等的にゼロになるためには、

\begin{eqnarray}
\frac{\partial f}{\partial y''} = 0
\end{eqnarray}

である必要がある。
したがって結論は同じく、$y(x)$ は１次多項式である。

以上により、データ点 $p_n$ から先 $(x_n, \infty)$ の範囲において、

　エネルギー $U$ を最小にする区分多項式 $y(x)$ は１次式

であることが分かった。

スプライン曲線全体

３つの区間

　・データ点 $p_0$ に到達するまで
　・データ点範囲内 $p_0 \thicksim p_n$
　・データ点 $p_n$ から先

それぞれで弾性エネルギーを最小にする曲線を求めつなぎ合わせれば、全体においても弾性エネルギーを最小にする曲線となる。
※ 尚、接合点において曲率を合わせない場合はその限りでない。

結論

データ点を通る曲線を物理的な連続体であるとした場合、曲げによる弾性エネルギーを最小にする曲線は

　〇　範囲 $(-\infty, x_0)$ においては、　$y(x) = a + bx$
　〇　範囲 $[x_0, x_n]$ においては、　$y(x) = a + bx + cx^2 + dx^3$
　〇　範囲 $(x_n, \infty)$ においては、　$y(x) = a + bx$

である。
各範囲において、境界点で微分値が連続するように関数を求めつなぎ合わせることで、曲げの弾性エネルギーが最小となるスプライン曲線が求められる。

留意点

その１

汎関数 $I[y]$ を微分しゼロになる点を求めることで、弾性エネルギーを最小とする曲線（関数） $y(x)$ を求めた。
しかし、汎関数の微分値がゼロになる点は一般的に最小を与える点とは限らず、最大、極大、極小を与える場合もある。
今回は最小を与えるが、その理由は今回触れていない。

その２

弾性曲線方程式の導出において、「たわみ角が小さい」という条件下で導いた。
しかし一般にデータ点間をスプライン補間する場合において、たわみ曲線のたわみ角が小さいとは限らない。

前提条件下でない場合においての理論の正当性は別途検証が必要と思われる。

おわり

以上です。
間違いなどありましたらご指摘いただければ幸いです。