繪製懸鏈線

在為 Balloon Platform Defense 製作氣球之間的繩子時,我得弄清楚如何畫出懸掛在兩點之間的繩子,這種形狀通常稱為懸鏈線。我在網路上找到的所有例子,要嘛給問題加了限制,例如兩個點位於同一高度,要嘛仰賴一些已知條件,例如繩子起始時的角度,而這些在這裡都不知道——已知的只有繩子的起點、終點和繩長。維基百科提供了許多關於懸鏈線方程式的資訊,我就以此為起點。 懸鏈線的方程式為:
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 的值,也就是方程式中基本上決定曲線寬窄的常數。注意:如果兩點在 x 軸上的位置相同,a 就是無窮大,所以在走到這一步之前就應該處理好這種情況——甚至包括兩點非常接近、a 大到無法計算的情況。(這裡「大到無法計算」指的是:當作指數函數的參數時,在你所用的浮點數精度下會使它傳回 NaN 或無窮大的任何數。)當然,如果兩點之間的距離大於繩長,也不要讓程式碼走到這一步——結果會很糟糕。

計算縮放係數

為了計算 a,維基百科給出了下面的方程式(它基於懸鏈線的一個性質:cosh 給出曲線的位置,而 sinh 給出曲線的長度):
\sqrt{s^{2}-v^{2}}=2a\sinh \left ( \frac{h}{2a} \right )

其中 s 是繩長,h 和 v 是起點與終點之間水平和垂直方向的距離(取絕對值)。這些都是已知量,只剩一個未知數,但必須用數值方法求解。我起初沒怎麼講究初始值,結果收斂與否全憑運氣。幸好,找到一個好的初始值並不難。雙曲正弦的泰勒級數為:
\sinh x=x+\frac{x^{3}}{3!}+\frac{x^{5}}{5!}+\cdots

令 u = \frac{1}{4a^{2}},c = \sqrt{s^{2}-v^{2}},要解的方程式就變成:
c=\frac{1}{\sqrt{u}}\sinh \left ( h\sqrt{u} \right )

取 sinh 泰勒級數的前三項並化簡,得到:
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 )
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 )

c=h+\frac{h^{3}u}{3!} +\frac{h^{5}u^{2}}{5!}

整理後得到:

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

這是一個關於 u 的簡單的一元二次方程式,把它的係數代入中學學過的公式解
x = \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}
就能得到 u 的一個好的初始值,進而得到 a 的初始值。
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}},其中 a=\frac{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 )
由於
a_{n+1}=a_{n}-\frac{f\left ( a_{n} \right )}{{f}'\left ( a_{n} \right )}
,
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 )}
我目前檢查數列是否收斂到與前一個值相差不超過 0.001,這通常需要 2 到 4 次疊代,不過當 x 座標比較接近時,有時要十次以上。更高階的 Householder 法大概也可以用,而且不會增加太多計算量,因為大部分計算都花在求雙曲函數上,而函數的更高階導數應該只會多出它們的倍數。既然縮放已經知道了,就可以計算平移了。

計算平移

標準方程式
y=a\cosh \left (\frac{x}{a} \right )
假設懸鏈線對稱於 y 軸,最低點在 (0,a),而我們要畫的曲線可能在螢幕上的任何位置,所以必須求出一個把它移到正確位置的平移 (p,q),方程式變為:
y-q=a\cosh \left (\frac{x-p}{a} \right )
(提醒一下:此時 a 已知,x 和 y 是變數,所以只有 p 和 q 仍是未知數。)我的第一個想法是把 x 和 y 的已知值——也就是繩子的兩個端點——代入方程式。理論上這是可行的,但無論我怎樣盡力化簡,得到的方程式在計算過程中常常產生雙精度浮點數無法處理的中間值。最後我只好另找辦法,結果反而更有效率。 還記得嗎?把懸鏈線方程式中的 cosh 換成 sinh,就得到曲線的長度。把左、右端點的 x 座標以及要畫的繩子的已知長度(s)代入這個形式,得到方程式:
a\sinh \left (\frac{x_{right}-p}{a} \right )-a\sinh \left (\frac{x_{left}-p}{a} \right )=s
雙曲函數的一個恆等式是:
\sinh x - \sinh y=2\cosh\left ( \frac{x+y}{2} \right ) \sinh \left ( \frac{x-y}{2} \right )
因此:
\sinh \left (\frac{x_{right}-p}{a} \right )-\sinh \left (\frac{x_{left}-p}{a} \right )=2\cosh\left ( \frac{x_{right}+x_{left}-2p}{2a} \right ) \sinh \left ( \frac{x_{right}-x_{left}}{2a} \right )
在 cosh 的參數中,左右兩個 x 座標相加再除以 2,得到的是中點的座標,所以不妨直接用它代替。sinh 參數中的減法得到的正是前面稱為 h 的值,所以這個 sinh 也可以當作已知量算出,原方程式可以整理為:
\cosh\left ( \frac{x_{middle}-p}{a} \right ) = \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )}
現在可以用反雙曲餘弦(arcosh)解出 p,但需要小心。由於 cosh 是偶函數,我們得到的總是正的答案,但還有一個負的答案,所以 p 有兩個可能的值。
p =x_{middle}\pm a\cosh^{-1}\left ( \frac{s}{2a\sinh \left ( \frac{h}{2a}\right )} \right )
要判斷用哪一個,可以看端點的 y 座標。拿一段繩子試一試就會發現:在繩長和 x 座標都相同的情況下,如果右端的 y 座標比左端低,曲線的最低點(如果它確實落在兩個 x 座標之間)就在中點的右邊(因此對應正的 arcosh),反之亦然。 還可能出現另一個問題:當兩個 y 座標相等時,曲線的最低點應該正好在中間,所以 cosh 的參數應該為零,取 arcosh 的那個值相應地應該為一。由於浮點誤差,它常常算出來略小於一,在取 arcosh 時就會出錯。明智的做法是在開始這一步計算之前,先檢查兩個 y 座標是否相等(或幾乎相等),既能省去計算,又能避免這個錯誤;在這種情況下,直接認為 p 正好位於兩個端點的中間即可。 既然 a 和 p 都已知,只要把其中一個端點的 x 和 y 值代入方程式
y-q=a\cosh \left (\frac{x-p}{a} \right )
並整理為:
q=y-a\cosh \left (\frac{x-p}{a} \right )
就能算出 q。求出 q 之後,你就得到了用來畫曲線的最終方程式。
y=a\cosh \left (\frac{x-p}{a} \right )+q
現在畫曲線就很簡單了。在計算參數的大多數步驟中,我都盡可能精確地計算雙曲函數,因為很難判斷誤差會如何傳播。但在最後這一步——它佔了大部分計算量——你知道所用參數有多大、能容忍多大的誤差,可以預測展開到某一次冪的泰勒級數是否足夠精確。根據我的經驗,大多數時候都值得使用。