カテナリー曲線を描く

Balloon Platform Defenseで風船どうしをつなぐ糸を作るとき、2点の間に垂れ下がる糸、より一般にはカテナリー曲線と呼ばれる形の描き方を考える必要がありました。ネットで見つけた例はどれも、2点が同じ高さにあるといった制約を問題に課しているか、糸の始まりの角度のように、この場合には分からない情報を前提にしていました。分かっているのは糸の始点と終点、そして糸の長さだけです。カテナリー曲線を表す式についてはWikipediaに多くの情報があり、それを出発点にしました。

風船をドラッグすると、糸の端を動かせます。

2つの風船の間に垂れ下がる糸。以下で説明する方法で求めています。どちらかの風船をドラッグするか、糸の長さを変えてみてください。薄い破線は、同じカテナリー曲線の残りの部分です。

カテナリー曲線の式は次のとおりです。

\displaystyle y=a\cosh \left (\frac{x}{a} \right )=\frac{a\left ( e^{\frac{x}{a}}+e^{-\frac{x}{a}} \right )}{2}

この式は、糸の最も低い点がy軸と交わる点にあると仮定しています。実際には、xとyにそれぞれ定数が加わり、そのどちらもまだ決まっていません。最も低い点がどこに来るか分からないからで、これが教科書のどの例よりも複雑になる理由の一つです。（現実の数学の問題は、いつも教科書の例よりずっと複雑だと思いませんか。解析的に解ける積分が出てくる現実の問題には、まだお目にかかったことがありません。）しかしそれは後回しです。まず、曲線の幅の狭さ・広さを実質的に決める式中の定数aの値を求めなければなりません。2点のx座標が同じ場合、aは無限大になるので、ここに来る前にその場合を処理しておく必要があります。2点が非常に近く、aが計算できないほど大きくなる場合も同様です。（ここで「計算できないほど大きい」とは、指数関数の引数に使うと、使っている浮動小数点数の精度ではNaNや無限大が返ってくるような数のことです。）また当然ながら、2点の間隔が糸の長さより大きい場合もここまで来ないようにしてください。ひどい結果になります。

スケール係数を求める

aを求めるために、Wikipediaは次の式を挙げています（coshが曲線の位置を与えるのに対し、sinhが曲線の長さを与えるというカテナリー曲線の性質に基づいています）。

\displaystyle \sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )

ここでsは糸の長さ、hとvは始点と終点の間の水平・垂直方向の距離（絶対値）です。どれも既知の値なので未知数は1つだけですが、数値的に解く必要があります。最初は初期値をあまり気にせずに試したところ、収束するかどうかは運次第でした。幸い、よい初期値は簡単に見つかります。双曲線正弦のテイラー級数は次のとおりです。

\displaystyle \sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots

u = \frac{1}{4a^{2}}、c = \sqrt{s^{2}-v^{2}}と置き換えると、解くべき式は次のようになります。

\displaystyle c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )

sinhのテイラー級数の最初の3項を取って整理すると、次のようになります。

\displaystyle c=\frac{1}{\sqrt{u}}\left (h\sqrt{u}+\frac{\left ( h\sqrt{u} \right )^{3}}{3!} +\frac{\left ( h\sqrt{u} \right )^{5}}{5!} \right )
\displaystyle c=\frac{1}{\sqrt{u}}\left (h\sqrt{u}+\frac{h^{3}u\sqrt{u}}{3!} +\frac{h^{5}u^{2}\sqrt{u}}{5!} \right )
\displaystyle c=h+\frac{h^{3}u}{3!} +\frac{h^{5}u^{2}}{5!}

これを変形すると：

\displaystyle \frac{h^{5}}{120}u^{2}+\frac{h^{3}}{6}u+\left ( h-c \right )=0

これはuについての単純な二次方程式なので、その係数を学校で習った二次方程式の解の公式

\displaystyle x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}

に代入すれば、uのよい初期値、ひいてはaのよい初期値が得られます。

\displaystyle u = \frac{-\frac{1}{6}h^{3}+ \sqrt{\frac{1}{36}h^{6}-\frac{1}{30}h^{5}\left ( h-c \right )}}{\frac{1}{60}h^{5}}\text{、}

ただしa=\dfrac{1}{2\sqrt{u}}

この初期値なら、ニュートン法で解を求めるのに十分近い値です。f\left ( a \right )=0の形にすると、f\left ( a \right )=2a \sinh \left ( \frac{h}{2a} \right )-cおよび{f}'\left ( a \right )=2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )となります。

\displaystyle a_{n+1}=a_{n}-\frac{f\left ( a_{n} \right )}{{f}'\left ( a_{n} \right )}

なので、

\displaystyle \begin{aligned} a_{n+1} &= a_{n}-\frac{2a \sinh \left ( \frac{h}{2a} \right )-c}{2 \sinh \left ( \frac{h}{2a} \right )-\frac{h}{a}\cosh \left ( \frac{h}{2a} \right )} \\ &= a_{n}-\frac{a \sinh \left ( \frac{h}{2a} \right )-0.5c}{ \sinh \left ( \frac{h}{2a} \right )-\frac{h}{2a}\cosh \left ( \frac{h}{2a} \right )} \end{aligned}

現在は、数列が前の値から0.001以内に収束したかどうかを調べています。たいていは2〜4回の反復で収束しますが、x座標どうしが近いと10回以上かかることもあります。計算の大部分は双曲線関数の計算であり、関数をさらに微分してもその倍数が増えるだけのはずなので、高次のHouseholder法も計算量をあまり増やさずに使えるでしょう。これでスケールが分かったので、平行移動を求めます。

平行移動を求める

標準形の式y=a\cosh \left (\frac{x}{a} \right )は、カテナリー曲線がy軸に関して対称で、最下点が(0,a)にあると仮定していますが、私たちが描く曲線は画面上のどこにあってもかまいません。そこで、正しい位置に移す平行移動(p,q)を求める必要があり、式は次のようになります。

\displaystyle y-q=a\cosh \left (\frac{x-p}{a} \right )

（念のため：この時点でaは既知、xとyは変数なので、未知なのはpとqだけです。）最初に考えたのは、この式にxとyの既知の値、つまり糸の両端の座標を代入することでした。理論上はこれでうまくいきますが、できる限り式を簡単にしても、計算の途中で倍精度浮動小数点数では扱えない値がたびたび出てきました。結局別の方法を見つけることになりましたが、そちらのほうが効率的でもありました。

カテナリー曲線の式でcoshをsinhに置き換えると、曲線の長さが得られるのでした。この形に左右の端点のx座標と、描きたい糸の既知の長さ（s）を代入すると、次の式が得られます。

\displaystyle \begin{aligned} & a\sinh \left (\frac{x_{right}-p}{a} \right ) \\ & \quad -a\sinh \left (\frac{x_{left}-p}{a} \right )=s \end{aligned}

双曲線関数の恒等式の一つに、次のものがあります。

\displaystyle \begin{aligned} & \sinh x - \sinh y \\ & \quad =2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right ) \end{aligned}

したがって：

\displaystyle \begin{aligned} & \sinh \left (\frac{x_{right}-p}{a} \right )-\sinh \left (\frac{x_{left}-p}{a} \right ) \\ & \quad =2\cosh\left ( \frac{x_{right}+x_{left}-2p}{2a} \right ) \\ & \qquad \times \sinh \left ( \frac{x_{right}-x_{left}}{2a} \right ) \end{aligned}

coshの引数では左右のx座標を足して2で割っているので、中点の座標になります。そこで、それで置き換えてしまいましょう。sinhの引数の引き算は、先ほどhと呼んだ値と同じなので、このsinhも既知の量として計算でき、元の式は次のように変形できます。

\displaystyle \cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}

これで逆双曲線余弦関数（arcosh）を使ってpを求められますが、注意が必要です。coshは偶関数なので、返ってくる答えは常に正ですが、負の答えもあり、pの候補は2つあります。

\displaystyle p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )

どちらの値が必要かは、端点のy座標を見れば分かります。ひもを使って試してみると分かるように、糸の長さとx座標が同じなら、右のy座標が左のy座標より低いとき、曲線の最下点は（そもそも両端のx座標の間にあるなら）中央より右に来ます。つまり正のarcoshから得られるほうです。逆の場合はその反対です。

風船をドラッグすると、糸の端を動かせます。

標準形のカテナリー曲線y = a cosh(x/a)（薄い線、最下点は(0, a)）を(p, q)だけ平行移動して、風船の間に垂らしたもの。破線の曲線は、arcoshから得られるpのもう一方の値を使ったもので、左端からqを求めると、右端を通りません。風船を上下に動かして、正しい値が中央のどちら側にあるかを確かめてみてください。

もう一つ起こりうる問題として、y座標が等しいとき、曲線の最下点はちょうど中央に来るはずなので、coshの引数は0、arcoshを取る値は1になるはずです。ところが浮動小数点数の誤差のため、1よりわずかに小さくなることが多く、arcoshを取るときにエラーになります。この計算の前にy座標が等しい（またはほぼ等しい）かどうかを確認しておくのが賢明です。計算も省けてこのエラーも防げます。その場合は、pがちょうど両端点の中間にあるとみなせばよいのです。

aとpが分かったので、qは、端点の一方のxとyの値を次の式に代入し、

\displaystyle y-q=a\cosh \left (\frac{x-p}{a} \right )

次のように変形するだけで求められます。

\displaystyle q=y-a\cosh \left (\frac{x-p}{a} \right )

これでqを求めれば、曲線を描くための最終的な式が得られます。

\displaystyle y=a\cosh \left (\frac{x-p}{a} \right )+q

あとは曲線を描くだけです。パラメータを求める計算のほとんどの段階では、誤差がどう伝わるかを見積もるのが難しかったので、双曲線関数をできるだけ正確に計算しました。しかし計算の大部分を占めるこの最後の段階では、使う引数の大きさも許容できる誤差も分かっているので、ある次数までのテイラー級数で十分な精度が出るかどうかを予測できます。私の経験では、たいていの場合テイラー級数を使う価値があります。

← すべての記事