科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

エアコン室外機の騒音対策は、ナベヤ(株)のビルトイン防振マウントが良いみたい

ナベヤ(株)のビルトイン防振マウント(軽量)を使ったエアコン室外機の防振対策を紹介します。

f:id:sken20k:20190810210529p:plain

目標

目標は、低い周波数の振動 (~30Hz, ピアノの一番低い音くらい)がベランダの床から部屋の壁に伝わるのを防ぐことです。 周波数 f に対する伝達率は固有振動の周波数を fn を使って、f=3*fn の時 0.13 になるので、 fn を 10 Hz 程度まで小さくできれば、伝わる振動(30Hz )を一桁小さくできるはずです。

選定

以下の条件を確認しつながら防振マウントを選定します。 固有振動数( or バネ定数) / 許容加重 / 取り付け(ボルト直径, 安全性) / 耐候性 / 価格

取り付け

BBC35 シリーズより小さい防振マウントだと、取り付けネジが M6 になってエアコンには小さすぎるので BBC35 から選びます。BBC35 の取り付けは M8 なので、ワッシャを使えば不安なく取り付けられます。日曜大工では、こういう扱いやすさがなにより重要だったりします。

許容加重と固有振動数(バネ定数)

室外機は左右で重さが違うので、とりあえず、4 つの足に 7 + 7 + 3.5 + 3.5 = 21 kg の荷重を仮定してマウントを選定します(7kg はコンプレッサー側)。 BBC35D007 は荷重 5 kgf に対して固有振動数が 11 Hz、BBC35D013 は荷重10kgf に対して固有振動数が 11 Hz です。これらを組み合わせれば、室外機 (約 20kg) を保持できそうです。荷重が均等であれば BBC35D007 4つで対応できますが、室外機は左右の荷重に差があるので、許容加重の大きい BBC35D013 と組み合わせます。 コンプレッサー側は D013 x 2, ファン側に D007 x 2 を使うことにします。 カタログの固有振動数から求めたバネ定数は、D007 が 25 N/mm, D013 が 48 N/mm になりました。

耐候性

メーカーのカタログに以下のような記述があり、耐候性にもすぐれているようです。

エーテル系ポリウレタンで、水や紫外線に強く、ゴムのような劣化がほとんどありません。

価格

一個あたり 1500円くらいなので、4つで 6000円になりますが、防振性能、取り付けの簡便さ、安全性などを考えて、この製品におちつきました。

防振マウントの外観

外観については、このページの先頭の写真を参照してください。 BBC35D007 は手で簡単につぶれるくらいやわらかい素材です。 もちろん、手をはなせばきちんと復元します。 雨にもつよいし、耐候性もあるようです。金属部分がステンレスであることも安心材料です。

取り付け

f:id:sken20k:20190810210903j:plain 手元のエアコンの台(プラロックというらしい)は下向きのボルトに対応していなかったので、M8 の穴をあけました。プライヤーで下側のステンレスのリムを回して固定します。M8 だと溝の中でナットが回ってしまうので、 細めのドライバーを挿してナットの回転を防ぎながら締めました。スプリングワッシャーを入れる厚みはなかったので、下の溝の中は平ワッシャーのみです。

f:id:sken20k:20190810210944j:plain 4つ取り付けたところ。室外機を取り付けるときは、ホースに無理な力が加わらないように気をつける必要があるので、できれば二人以上で行った方が良いでしょう。

結果 (沈み込み量の計測)

f:id:sken20k:20190810092305p:plain 曲線は固有振動数(11Hz)が最大荷重に対応するとして荷重と固有振動数の関係をグラフにしました。 曲線の色付きの部分は製品の許容加重です。荷重が小さい側にのばしました(灰色の線)。 曲線上の数値は沈み込み量 [mm] です。 A~D は 4つの足の沈み込み量(ノギスで計測)によるプロットです。

記号 説明 沈み込み [mm] 加重 [kgf]
A コンプレッサー側 前面 1.4 6.8
B コンプレッサー側 後面 2.0 9.7
C ファン側 前面 0.9 2.3
D ファン側 後面 1.5 3.9

ほぼ期待通りの結果が得られ、低い周波数の振動も気にならなくなりました。 防振マウントの性能(低周波の低減)を最大に活かすには、追加ウェイトによる調整が必用です。 A, B は 荷重が 10kg, C,D は荷重が 5 kg になるようにウェイトを追加します。 今回は、体感で振動を十分に低減できたので、追加ウェイトによる調整は見送りました。

今回の室外機は小型のエアコンのものですが、大きいエアコンに比べて騒音/振動が大きいようです。 個体差なのか製品の特性なのかわかりませんが、小型のエアコンほど安く仕上げるために材料を削る影響が大きくなって、振動が出やすくなっているのではないかと疑っています。

注意/制限事項

今回の防振対策では、かなり単純な理論式にそって選定を行いました。固有振動数の式はこんな感じ。

 \displaystyle fn = \frac{1}{2\pi}\sqrt{\frac{k g}{L}} = \frac{1}{2\pi}\sqrt{\frac{k g}{W g}}

fn [Hz] 固有振動数, k [N/m] ばね定数, g =9.8 [m/s2] 重力加速度, L [N] 加重, W [kg] 重量

実際は、バネ定数の動的/静的の違い、バネ定数の温度特性、減衰特性...などなど、いろいろ考慮する必用があるので、どこまで性能が出ているかよくわかりません。指でマウントに触ると、マウントの上と下とで振動が全然違う(下側の振動は分からない)ので、効果は体感できます。対策前後の振動を定量的に計測できれば良いのですが、低い周波数の計測は難しいみたいです。

もう少しお安く...

固有振動数の目標を 11Hz から 14Hz に緩和すれば、おなじナベヤ(株)の防振フット(軽量型)も良さそうです。例えば BFL20D030 (荷重が 3~7 kg) を 4 つ使って、4つの足すべての荷重が7kg になるようにウェイトを追加すれば、最大の効果が得られるはずです。こちらは 1 個 500円くらいなので 4 個で 2000円となります。ボルトが M5 なので、外形の大きい平ワッシャ(特寸?)を用意する必要がありそうです(エアコン側は M10で共通らしい)。我が家のケースでは、固有振動数 14Hz でも十分だったかな、と思っています。

他メーカーの防振ゴム

低い周波数まで対応するには、バネ定数を小さく抑える必要があります。その観点で、 最後まで迷ったのが、倉敷化工の軽量用防振ゴム KXA-25-25H でした。 防振性能(バネ定数=14.4 N/mm)は魅力的ですが、1つ当たりの支持荷重が 3 kgf なので小型の室外機でも 6~7 個使って支えることになります。取り付け難しそうなのであきらめました。こちらを使えば防振部品の値段は半分くらいになりそうです。ただし、取り付け板や加工などの手間を考えると、トータルの費用はあまり変わらないかもしれません。厚めの合板を塗装したプレート(室外機を置ける1枚板)を用意するとよい気もします。マウントを上下で挟むようにして、荷重の偏りに応じて防振ゴムを適宜配置すれば、今回のものよりも性能の高い防振装置を安価に実現できる可能性があります。プレートに室外機を置いただけだと、横滑りが怖いので、試すなら安全対策もお忘れなく(室外機をプレートに固定する)。

金属バネ

バネ定数と耐荷重でスプリングを検索して、ばね通販 スプリングネットの 11-2820 という製品にたどりつきました。 バネ定数が 4.9 N/mm と小さくて、全タワミ(28.9mm)まで使えるなら 1 つで 14kg まで対応できる。 実際は 10mm 程度のタワミで 1つあたり 5kg の荷重で使うのが良さそう。 スプリングの上に室外機をそのまま載せると、たぶん危なっかしい感じになりそうです。 やわらかなスプリングなのでシリンダー内に装着しないと危ないかもしれません。 その辺の安全確保にかんする知識がないので、スプリングの使用はあきらめました。

参考 URL

  1. エアコン室外機の振動対策 (猫下僕の初心者DIY 日記) url0
  2. エアコン室外機の音を静かにする方法 (創造の館) url1
  3. エコキュートの騒音対策 (noahnoah 研究所) url2
  4. NOK(株) 防振ゴムカタログ (選定手順やゴムの特性が詳しい) url3
  5. ナベヤ(株) ビルトイン防振マウント カタログ url4
  6. 倉敷化工(株) 軽量用防振ゴム KXA/KXB シリーズ url5
  7. サミニ(株) ばね通販 スプリングネット (型番 11-2820) url6

楕円パラメータの相互変換 (幾何学的 ⇔ 代数的)

楕円の表現

この記事では、楕円を表す2種類の形式(方程式)と、 それらの間の相互変換について述べる。 2種類の形式は、一般に canonical form, general form と呼ばれることが多いが、 ここでは便宜上、「幾何学的な表示形式」、「代数的な表示形式」と呼ぶことにする。 wikipedia(英語版) の ellipse の項をみると相互変換に関する解析解が掲載されている。 それを使っても良いが、式の形や導出がけっこうややこしいので、 別解として、比較的簡単に数値的に求める方法を紹介する。

f:id:sken20k:20181105181127p:plain

楕円の幾何学的な表示形式

幾何学的に楕円を表現する時は、パラメータとして  x_0, y_0, \theta, a, b を選ぶことが多い。 ここで、 x_0, y_0 は楕円の中心、  \theta は楕円の傾き、  a, b は楕円の半径(横方向と縦方向)である。  \cos \theta,  \sin \theta c, s で表すと、楕円の方程式は、以下の通り。

 \displaystyle {
         \frac{ ( c(x - x _ 0) + s(y - y _ 0) )^ 2  }{a^ 2}
      + \frac{ ( -s(x - x _ 0) + c(y - y _ 0) )^ 2  }{b^ 2} = 1
}

楕円の代数的な表示形式

代数的に楕円を表現する時は、係数を x, y の冪乗ごとにまとめる。  x^{2} の係数が 1 になるように規格化し、残りの項の係数を  A,B,C,D,E とすれば、 楕円の方程式(一般的な二次曲線)は、以下の通り。

 \displaystyle { 
x^{2} + A x y + B y^{2} + C x + D y + E = 0
}

幾何学的な表示(x0,y0,θ,a,b) ⇒ 代数的な表示(A,B,C,D,E)

幾何学的な表現を展開したあと、 x, y の冪でまとめなおして係数を比較する。

 
\begin{aligned}
\alpha A &= \gamma \\
\alpha B &= \beta \\
\alpha C &= -2\alpha x_{0} - \gamma y_{0} \\
\alpha D &= -2\beta y_{0} - \gamma x_{0} \\
\alpha E  &= \alpha x_{0}^{2} + \beta y_{0}^{2} + \gamma x_{0} y_{0} - 1 \\
\end{aligned}

ここで、

 \displaystyle{
\alpha = \frac{c^{2}}{a^{2}} + \frac{s^{2}}{b^{2}}, \qquad
\beta = \frac{s^{2}}{a^{2}} + \frac{c^{2}}{b^{2}}, \qquad
\gamma = 2 c s \left( \frac{1}{a^{2}} - \frac{1}{b^{2}} \right)
}

である。 代数的なパラメータ A,B,C,D,E は、上の関係から直ちに求まる。

代数的な表示(A,B,C,D,E) ⇒ 幾何学的な表示(x0,y0,θ,a,b)

 \gamma, \beta を消去して、全体を  \alpha で割ると以下を得る。


\begin{aligned}
C &= -2 x_{0} - A y_{0} \\
D &= -2 B y_{0} - A x_{0} \\
E &= x_{0}^{2} + B y_{0}^{2} + A x_{0} y_{0} - 1/\alpha
\end{aligned}

 x_{0}, y_{0} は、最初の2式を連立して求める。行列形式で表せば、

 \displaystyle{
\begin{pmatrix}
C \\
D
\end{pmatrix}
=
\begin{pmatrix}
-2 & -A \\
-A & -2B
\end{pmatrix}
\begin{pmatrix}
x_{0} \\
y_{0}
\end{pmatrix}
}

 \theta は、以下の関係を使って、A,B から直接求める。

 \displaystyle{
\tan(2\theta)=\frac{2cs}{c^{2}-s^{2}}=\frac{\gamma}{\alpha - \beta}=\frac{A}{1 - B}
}

 \alpha は、 x_{0}, y_{0} を最後の式に代入して求める。

 \displaystyle{
1/\alpha = x_{0}^{2} + B y_{0}^{2} + A x_{0} y_{0} - E
}

 a, b は、 \theta \alpha, \beta の定義式に代入して、それらを連立して求める。

 \displaystyle{
\begin{pmatrix}
\alpha \\
\beta
\end{pmatrix}
=
\begin{pmatrix}
\alpha \\
\alpha B
\end{pmatrix}
=
\begin{pmatrix}
c^{2} & s^{2} \\
s^{2} & c^{2}
\end{pmatrix}
\begin{pmatrix}
1/a^{2} \\
1/b^{2}
\end{pmatrix}
}

備考

「超精度くりこみ法」というのを知ったのをきっかけに、 楕円のパラメーター推定をあれこれ試しています。その過程で、 楕円パラメータの相互変換が必要になったので、忘備録として記事を書きました。 部分的なアークから楕円を推定する方法は、 天文の分野でも応用が広そうですね。参考文献に挙げた 「3次元~」は、 各種アルゴリズムソースコードが出版社のページで公開されていてます。 コードを眺めるだけで動かしていませんが、こういうサービスはうれしいですね。

実装例

def geo_fm_ana5(A, B, C, D, E):
    """
    Parameters
    ----------
    A~E : 楕円方程式の係数 
                 (x^2 + A*x*y + B*y^2 + C*x + D*y + E = 0)
    Returns
    -------
    xc, yc : [-]   中心座標
    th     : [rad] 回転角
    ra, rb : [-]   横半径, 縦半径
    """
    xc, yc = np.linalg.solve([[-2, -A], [-A, -2*B]], [C, D])
    th = np.arctan2(A, 1 - B)/2
    c = np.cos(th)
    s = np.sin(th)
    k = 1/(xc*xc + B*yc*yc + A*xc*yc - E)
    a_inv2, b_inv2 = np.linalg.solve([[c*c, s*s], [s*s, c*c]], [k, k*B])
    ra = np.sqrt(1/a_inv2)
    rb = np.sqrt(1/b_inv2)
    return xc, yc, th, ra, rb

参考文献

  1. 最小二乗法による楕円近似 / url1
  2. 楕円当てはめの精度比較:最小二乗法から超精度くりこみ法まで url2
  3. 金谷健一(他), "3次元コンピュータビジョン計算ハンドブック", 2016, 森北出版
  4. はてなブログで数式を書くurl3

図形を使った解法は、見ためも安心なのです(改)

コンピューターグラフィックスのアイデアが、ピタッとと当てはまる問題というのがあります。 教科書に出てきそうな Bang-Bang 制御の切り替え線の判定も、そんな問題の一つです。 図形の論理演算(and/or)を使って、制御の ON/OFF を判定する例を紹介します。

Bang-Bang 制御

Bang-Bang 制御は、位置  x と速度  \dot{x} に応じて、 一定量の制御を加えるかどうかを判定します。 いろんなパターンがありますが、たとえばこんな感じ(位相面)。

f:id:sken20k:20180525173650j:plain:w200

水色の領域では下向きに、緑色の領域では上向きに制御を加えます(速度を変化させる)。 折れ線になっているので、判定はちょっと面倒ですね。

よくある解法

よくある解法は、直線の交点を求めて x との大小関係で比較直線を選ぶというものでしょう。 そのままロジックを組むと、たくさんの交点を求めるのも面倒だし、 判定も if else がいくつも連なって、ややこしくなりそうです。

図形を使った解法

図形を使った解法はこうです。 比較的単純な領域を定義して、それらの論理和としてもとの領域を表現します。

f:id:sken20k:20180525173956j:plain:w400

制御を加える領域は、上の3つの領域の和になっています。 仕様(図計)を、プログラムに直せば、こんな感じになります。

def judge(x, y):
    y1, y2, y3, y4, y5 = f1(x), f2(x), f3(x),  f4(x), f5(x)
    if (y > y1 and y > y2) or (y>y3 and y>y4) or (y > y5):
        output = -1
   else:
        output = 0
    return output

仕様(図形を使った表現)と実装の対応が良いのも魅力です。 レビューをしっかりすれば、ほぼ間違えることはないのだけれど、 この手のロジックを検証する場合は、作った関数を使って、元の切り替え線(判定領域)を 作図してみるのが確実です(位相面を2次元でサンプリングする)。

おまけ

3次元CG のモデリング(画像生成)では、基本的な形状の組みあわせによって 複雑な形状を表現することがあります。 例えば、レンズの形は、球と球の and 領域として定義できます。 上の、Bang-Bang 制御の判定ロジックは、そこからアイデアをいただきました。 むかし CGや図形処理で遊んでいましたが、そういう経験が結構役立つ場面は多い気がします。 図形処理は、応用範囲が広いということですかね。

多変数最適化の Nelder-Mead は、なかなか使い勝手が良いのです

多変数最適化のアルゴリズムの一種である Nelder-Mead は、なかなか使い勝手の良いアルゴリズムですね。 汎用的に使える(対象が最小二乗に限定されない)ことに加えて、 微分情報を使わないところも魅力です。 評価関数の性格が多少悪くても、上手に最適をみつけてくれる、便利屋さんといったところでしょうか。

Nelder-Mead のお勧めバージョン

Nelder-Mead には、いろんな改良/拡張が存在しますが、 場合分けの見た目の安心感? から、以下の 5 種類の動作を含んだ設計がお勧めです。

  1. reflection
  2. expansion
  3. inside-contraction
  4. outside-contraction
  5. shrinking

細かく言うと、最悪, 2番悪, 最良 の 3点の情報を使って 7 通りに場合分けしたあと、 動作が同じになるところをまとめると、最終的に 5 ケースとなります。 例えば http://adl.stanford.edu/aa222/lecture_notes_files/chapter6_gradfree.pdf

アルゴリズムの妥当性と試験について

この手のアルゴリズムは、経験則の寄せ集め的なところがあって、 絶対的な正解という考え方には馴染まない気がします。 適用限界を明示するのが難しい。 あと設計やコーディングに間違いを含んでいても、 なんとなく動いてしまって、しかも、本来と異なる経路で最適値に収束したりして、 外部仕様を定めにくいですね。つまり、ブラックボックス的な試験設計は難しい。 開発手法(管理)の観点からは、興味深い対象かもしれません(扱いが難しいということ)。 というわけで、ホワイトボックス的な試験に重点を置くのが良さそうです。とにかく全てのパスを通す。 このとき、見やすい図を作って直感的に動作を確認するのは、 そう悪くない方法かと思ったりもします。例を下図に示します。

f:id:sken20k:20180513141435j:plain

適当な評価関数(2次元)を用意して、乱数で初期値を選ぶようにして、 見やすい図が得られるまで乱数のシードをいろいろ試しました。。 スライダーなどの GUI で初期値を選べるようにすると、おもしろいデモが 作れそうですね。

パラメータの一部を固定して評価したいとき

最適化の問題をあつかっていると、一部のパラメータを固定して推定したい。 という場面に出くわすことがあります。 パラメータの組み合わせ毎に最適化対象の関数を定義しはじめると、 評価関数だらけになってしまって、何をやってるのか分からなくなりがちです。 最適化のルーチンの I/F に手を加えて、要素毎に推定する/しないを指定するフラグを追加 しても良いのですが、 python のようなインタープリターを使うのであれば、 固定パラメータを補う関数を動的に生成するのも悪くないですね。 こうすれば、ベーシックな最適化ルーチンに手を加えずシンプルに保ったまま、 目的を達することができます。

落穂ひろい

Nelder-Mead 以外で、実用の見地から押さえておきたい最適化アルゴリズムは、 やっぱり Levenberg-Marquardt でしょう(LM法)。 適用できる問題は最小二乗形式に限定されますが、 変数の数が多めのデリケートな問題でも、良好な結果が得られることが多いですね。 推定誤差や correlation を評価がしやすいところもポイントが高いです。 注意が必要なのは、ダンピングファクターの初期値の調整くらいでしょうか。 実装にもよりますが、ダンピングファクターと収束判定の関係を理解していないと、 使いこなせない場面が、時々あります。 最適化ルーチンをブラックボックス的に利用するのは良いのですが、こういう 特性を理解した上でないと、使いこなせないこともある。 ブラックボックスの表面に、そういう細かな特性を認識できているかどうかで話が大分違ってきそうです。 個々のアルゴリズムや理論の前提条件や適用限界に注意を促したり、それらを身に付けるための方法論が 求められる。そう思う次第です。

剛体回転に関するオイラーの定理は、絵的に理解したいのです

「剛体回転におけるオイラーの定理」について、 天球図を使った直感的な説明を紹介します。 代数的な説明だけだと、なんとなく分かった気がしない...という人におすすめです。

オイラーの定理

オイラーの名の付く定理や公式は、いくつもあって紛らわしいのですが、 その中の一つに「剛体回転におけるオイラーの定理」というのがあります。

【定理】 「点のまわりの回転変位は, その点を通る1つの軸のまわりの回転によって達せられる」

3次元座標の座標回転の言葉で言い直すと、

「原点まわりの任意の座標回転は、ある方向ベクトルを回転の軸とする1回の回転で達成できる」

1回の回転で変換を達成する回転の軸を「オイラー軸」あるいは「オイラー回転軸」と呼び、 その時の回転の角度を「オイラー回転角」と呼んだりします。

オイラー角」と名前が似てますが別物です。「オイラー角」が座標軸回りの3回の回転の組み合わせであるのに対して、「オイラー回転角」は任意方向回りの1回の回転です。

幾何学的な説明

f:id:sken20k:20180427211221j:plain

上の図は、座標系 X-Y-Z から X'-Y'-Z' への回転を表しています。 ひとまず Z を除いて、X と Y に注目した上で、回転軸の成分が保存される (回転軸とXの成す角は、回転軸とX'の成す角と等しい)ことに注目すれば、 回転軸の方向は、 X, X' を底辺とする二等辺(球面)三角形の頂点に存在するはずです。つまり、 X を X' に移す回転軸は X と X' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 同様に、Y を Y' に移す回転軸は Y と Y' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 これらを同時に満たす点は2つの大円の交点として求まるので、これを P とする。 (逆方向の解も存在するが、どちらを選んでも議論の一般性は失われないので見やすい方を選ぶ)。

これで、オイラー回転軸の候補 P が手に入りましたが、 X', Y' を同時に満たす回転(角度)が存在するかどうか分かりません。 以下、それを確認します。

XをX'に移す回転角をφx、YをY'に移す回転角をφyとする。 オイラーの定理が成り立つためには、φx とφy が一致すればよい。

Pの成分が座標回転に対して不変なので、 球面三角 P-X-Y と P-X'-Y' は合同であることが分かる。 Pを頂点とする角度から重なり部分を除いた角度がφxとφyなので、 φxとφyは確かに等しい。

Z, Z' についても同様の議論が成り立つ。

これで、オイラーの定例が成り立つことが確認できました。

代数的な説明

先に幾何学的な説明をしましたが、 オイラーの定理は代数的に説明されることの方が多いですね。 座標回転を表す行列の固有値固有ベクトルを使った議論は、 あいまいさなく簡潔に証明できるところが利点と思います。 こちらは、どの教科書にも載っているので説明は省略します。

おまけ

オイラーの定理の絵的な理解は、実践的(応用範囲が広い)知見と思います。 図解した教科書が少ない(気がするので)ここに紹介しました。

補間と言えば akima-spline は知らなきゃもぐりだろう、と言ってみたいのです

スプライン補間でもっともメジャーなのは cubic-spline でしょうか。 ある種の理解しやすさや微分も一緒に求まるという使い良さもあるのですが、 データにギャップがある場合に、補間結果が振動するという やっかいな特性を持っています。 下の図の赤線は、サンプル点(黒点)を cubic-spline で補間した結果です。 ギャップの前後で補間値があばれているのが分かります。

同じ図に、akima-spline で補間した結果を青線で示します。 ギャップの前後も補間値があばれることなく、スムースな補間結果が得られことが分かります。

f:id:sken20k:20180420234011j:plain

参考までに、別の例も載せておきます。 f:id:sken20k:20180420235709j:plain

補間方法には、それぞれ利点/欠点はあるので、 akima-spline が万能というわけではありませんが、 区間ごとに推定した軌道や姿勢をつなぐ場面では、つなぎ目に微小なギャップが 現れることが多いので、そういうとき akima-spline のなめらかさ(オーバーシュートの無さ)は とてもうれしい特性なのでした。

akima は、秋間浩さんという日本人の名前に由来します。

高級な補間公式が、万能とは限らないのです。

補間公式には、ラグランジュ、エルミート、スプライン...と、各種ありますが、 私がよく使うのは、線形補間(直線補間)とAkima-Spline の 2 種類くらいです(JPL の惑星歴のチェビシェフは除く)。 線形補間(直線補間)とは、随分と原始的なものを使うのですね...と思う人がいるかもしれませんが、 高級な補間公式で安全性を保障するのは、結構むずかしかったりします。

高級な補間公式の弱点

高級な補間公式は、(使い方を間違わなければ) 少ないサンプル点で高い精度を実現できますが、 それらは以下のような欠点ももっています。

  • 単調性が保証されない
  • 高階差分を適用するとき桁落ちが発生しやすい

これらの欠点について注意が払われることは、少ないように思います。 数値計算の教科書も、精度の観点から議論が多くを占めて、 単調性が保証されない点に関しては、たとえ記述されていたとしても さらりと一言述べられるだけの事が多いようです。 記述量が少ないため、読み飛ばされることもしばしば。 実運用や解析ツールの設計においては、とても大切な特性ですが、 痛い目に会わないと、なかなか身に付かない知見の一つかと思い、 ここに紹介します。

単調性の保証が大事な時

例えば、軌道情報が時系列データ(点列)として与えられたときに、 距離 R がある閾値を下回るタイミングを求めたいとする。 このとき、単調性が保証されない補間公式を使うと、 偽のイベントがジッターのように出現する可能性がある。 次の図では、本当のタイミングは time が 1.1 を超えあたりにありそうだけど、 高次の補間を使って求めようとしたために 0.9 ~ 1.0 の間に偽の事象が出現しています。

f:id:sken20k:20180417185817j:plain

補間の時に最後の Time = 1.2 の点を使わないのが原因ではないか? と思う人がいるかもしれませんが、Time = 1.0 までの点が揃ったところで、 直近のN 点から事象を判定する場面は、そう無理な仮定ではなかろうと思います。 高次の補間公式で単調性を保証するのは、いろいろとむずかしいので、 とにかく安全を保証したい場面でであれば、単調性が保証できる線形補間を使うことを お勧めする次第です。

上の図は、単調性が保たれない箇所を拡大しています。元のスケールで補間の様子全体を 表示すると以下のようになります(6点を使ったラグランジュ補間)。 f:id:sken20k:20180420183939j:plain

高階差分の桁落ち

高次の補間は、単調性の話以外にも精度の面で問題が生じる場合があります。 具体的には、補間公式には高階の差分演算が現れるため、似たような数値の減算の繰り返し が多くなり、次数が高くなるほど、つまり階数が深くなるほど、桁落ちが生じやすくなるからです。 これに対抗するには、最初のデータの有効桁を増やすしかありません。 時刻を例にとれば、準ユリウス日ユリウス日を時刻に使っていたりすると、 高階差分のために思わぬ誤差が生じることがあるので注意が必要です。 最近では、2000年からの時刻が使われたりすることもありますが、高い精 度が求められる場合には、解析期間の開始時刻を時刻原点に採用するのが無難です。 普段、あまり次数の高い補間公式を使わないようにしていても、DOP853 という 高精度の数値積分のルーチンに組み込まれた densout (稠密出力)などを利用する 場面で、あまり意識することなく、高階の差分を使ってしまうこともあります。

落穂ひろい

高次の補間公式が振動する話は有名ですが、単調性に焦点を当てた議論というのは 少ないように思います。教科書を書いている数値計算の専門家にとっては あたりまえのことであるので、ページ数が限られている教科書においては、 話が精度に偏ってしまうのも仕方ないことかもしれません。あるいは、 数値計算の専門家(研究者)には、実運用のシステムを設計/運用する経験はほとんどないと 思われるので、適用場面で重用になる特性に興味が薄いことも一因かと疑ったりしています。 というわけで、たとえあたりまえの話であっても、ある程度のボリュームを割いて、議論する意味 というのはあるのかなと思う次第です。もちろん、知識の効率的な吸収という観点とは相反するわけだけれども、 ある種の無駄は必要だろうという立場をとります。

少し一般的な話をすると、およそ公式や手法には、それぞれ利点だけでなく欠点もあるものです。 使える場面や条件を明晰に意識できるかどうかが、応用が利くかどうかの分かれ目ですね。 利点だけでなく欠点や制約を2つ3つ列挙できるかどうかが、その手法の理解度を測るバロメータ だと言う人もいます。

Akima Spline については、別記事にて...