二分搜索法求解懸鏈線問題

馬 旭 易 俗

(遼寧大學創新創業學院 遼寧 沈陽 110036)

二分搜索法求解懸鏈線問題

馬 旭 易 俗

(遼寧大學創新創業學院 遼寧 沈陽 110036)

懸鏈線問題是現代工程實踐中經常遇到的問題之一。分析懸點等高及不等高兩種情況下懸鏈系數、弛垂度、懸點水平坐標的常規數學求解方式，基于其運算量大，誤差難以把握的不足，提出應用二份搜索程序設計思想輔助求解的具體思路。分別給出3種不同應用的完整程序及數據，比較分析數學方式與程序設計求解方法。以吊桿架設工程應用為實例，結合完整C語言程序及其主函數流程圖 ，較詳細地描述了相關各參數的程序求解過程。上述方法已多次應用于實踐，是簡單可行的。

二分搜索法 懸鏈線 懸鏈系數 弛垂度 誤差分析 C程序設計

0 引 言

懸索吊橋、雙曲拱橋、架空電纜、索道滑車等諸多實踐應用中，都涉及到懸鏈線[1-6]。懸鏈線函數是超越函數，應用中可以通過查雙曲函數表、拋物線逼近、級數展開求解高次方程等數學方法近似求解[7-9]。在懸鏈線接近水平時，常規數學運算可以滿足一般工程的需要，但針對傾斜度較大的懸鏈線，數學運算計算繁雜，精度難以提高，更難以把握[10-12]。許多工程實踐中，在指定精度要求下，利用二分搜索程序輔助求解，更為簡單實用[13]。

本文將就懸鏈線應用中的多種實用形式，結合C語言描述，分別給出相關函數及應用程序，求解懸鏈系數、弛垂度、懸索鏈長及懸點水平坐標，并給出工程應用實例。

1 兩懸點高度相同的情況

如圖1所示，A、B兩懸點等高，水平距離為L，懸鏈線鏈長為S，懸點與最低點垂直距離(弛垂度)為D，懸鏈線最低點與橫軸垂直距離為a(懸鏈系數)，x1、x2為兩懸點的水平坐標。

圖1 兩懸點等高的懸鏈線坐標圖

懸鏈線的標準方程為：

(1)

實際應用中，兩懸點間水平距離一般可以預知，無論是已知鏈長求弛垂度，還是已知弛垂度求鏈長，還是在工程實踐中分析懸點的受力，首要解決的問題就是求懸鏈系數a。

1.1 已知鏈長及兩懸點間水平距離求弛垂度

1.1.1 常規數學方式求解弛垂度

弛垂度的標準方程為：

(2)

參考圖1求解D，利用曲線積分可以得出下列結果：

(3)

求解上述超越方程，可以經過冪級數展開、麥克勞林展開、麥克勞林系數求解等一系列代數運算，在忽略高于二次的級數項后，可以近似得出下列結果[14]。

設：

(4)

(5)

則：

(6)

(7)

此公式在常規實踐中是可以應用的，但在誤差精度要求更高，尤其在誤差精度要求隨應用而調整的情況下，就很難適用了。

1.1.2 二分搜索法求解弛垂度

設有函數y=f(x)，與數學求解函數的方式不同，利用二分搜索程序設計思想，預先編制搜索原函數double fun(double x)和二分搜索函數double dichotomy (double (*fp)(), double x1, double x2, double y)，在一定范圍(x1,x2)、一定精度要求下，即可通過二分搜索法搜索到指定函數值y對應的x(x=f-1(y))，可以有效地解決上述實際問題[15-16]。

設：

(8)

(9)

化簡后可得：

(10)

基于函數y=sinh(x)/x建立搜索原函數如下：

double funz( double x )

{ double y;

y=sinh(x)/x;

return y;

}

基于指定搜索原函數fp()、指定函數值y，在(x1,x2)區間，按指定精度要求AC，進行二分搜索，返回搜索點近似值x，二分搜索函數如下：

double dichotomy (double (*fp)(),double x1, double x2, double y)

{ double x, p;

do { x=(x1+x2)/2;

p=fp(x);

if (p>y) x2=x;

if (p

} while ( fabs(p-y) > AC );

return x;

}

由懸鏈線實踐確定的數據范圍可預先設定3>t>0，設Z=S/L，調用二分搜索函數dichotomy( funz, 0, 3, Z )可以搜索求得原函數Z=sinh(t)/t中t值，再利用下面公式即可求得懸鏈系數a及垂弛度D：

(11)

1.1.3 實 例

例1：給定一組鏈長S、懸點間距離L，利用二分搜索法，在誤差δ<10-3條件下，求懸鏈系數a及弛垂度D，主函數程序描述如下：

#include ″stdio.h″

#include ″math.h″

#define AC 1e-3

main()

{ double s[5]={25,35,45,55,65},

l[5]={20,30,40,50,60},

a[5],d[5];

int i;

for ( i=0; i<5; i++)

{ double t,p;

t=dichotomy( funz, 0, 3, s[i]/l[i] );

a[i]=l[i]/(2*t);

p=l[i]/(2*a[i]);

d[i]=a[i]*cosh(p)-a[i];

}

for ( i=0; i<5; i++)

printf(″ L=%6.3f, S=%6.3f, a=%6.3f, D=%6.3f ″, l[i],s[i],a[i],d[i]);

}

1.1.4 數學求解與程序求解的比較

指定精度要求為δ<10-5、δ<10-4，分別利用1.1.3節中程序求解得懸鏈線的弛垂度為D2、D3，通過MATLAB按1.1.1節中數學方法求得懸鏈線的弛垂度為D1，它們對比如表1所示。

表1 二分搜索法與常規數學方法求解馳垂度對比

參考表1可見，利用1.1.1節中數學解法求得的弛垂度D介于二分搜索程序誤差設置在10-5與10-4之間。

實踐中，可以在指定精度要求下，預先確定鏈長S與懸點距離L的比值表。具體工程中即可通過查表方式確定t值，再通過上面所述的公式求系數a、弛垂度D。查表求解更為簡單實用。表2為δ<10-6時鏈長S與懸點距離L的比值表。

表2 鏈長S與懸點距離L的比值表

1.2 已知弛垂度及兩懸點間水平距離求鏈長

1.2.1 二分搜索法求鏈長

設：

(12)

(13)

化簡后可得：

(14)

基于函數y=(cosh(x)-1)/x建立搜索原函數如下：

double funv(double x)

{ double y;

y=(cosh(x)-1)/x;

return y;

}

參照例1，在3>t>0范圍內，設V=D/(L/2)，調用二分搜索函數dichotomy(funv,0,3,V)搜索求解原函數V=(cosh(t)-1)/t中t值，再利用下面公式可以求得懸鏈系數a和鏈長S：

(15)

1.2.2 實 例

例2：給定一組弛垂度D，懸點距離L。利用二分搜索法，在誤差δ<10-3條件下，求懸鏈系數a及鏈長S，主函數程序描述如下：

#include ″stdio.h″

#include ″math.h″

#define CC 1e-3

main()

{ double l[5]={20.000,30.000,40.000, 50.000,60.000},

d[5]={6.636,7.924,9.030,10.015,10.911},a[5], s[5];

int i;

for ( i=0; i<5; i++)

{ double t,v,p,y;

v=2*d[i]/l[i];

t=dichotomy( funv,0,3,v );

a[i]=l[i]/(2*t);

p=l[i]/(2*a[i]);

s[i]=2*a[i]*sinh(p);

}

for ( i=0; i<5; i++)

printf(″L=%6.3f, D=%6.3f, a=%6.3f, S=%6.3f ″, l[i],d[i],a[i],s[i]);

}

上述程序設定的初始弛垂度D為表1中的D2。指定不同精度要求(δ<10-3，δ<10-5)利用二分搜索程序求解的懸鏈線鏈長(S1，S2)與表1中初始S對比如表3所示。

表3 二分搜索法所求鏈長S與真實鏈長對比

參考表3可見，搜索程序精度設為10-3時，誤差已經可以在許多實踐環境中被接受，當搜索程序精度設為10-5時，誤差就更小了。

2 兩懸點高度不同的情況

2.1 二分搜索法求懸點橫坐標及弛垂度

如圖2所示，A、B兩懸點水平距離為L，懸鏈線鏈長為S。較低的懸點A與最低點的垂直距離(弛垂度)為D，兩懸點的垂直距離為H，懸鏈線最低點與橫軸垂直距離為a，x1、x2為兩懸點的水平坐標。

圖2 兩懸點不等高的懸鏈線坐標圖

由懸鏈線標準方程，可知：

(16)

(17)

(18)

(19)

(20)

(21)

即:

(22)

通過上式可知，在給定L、S、H時，可以通過將sinh(x)級數展開，再解高次方程求得懸鏈系數a，進而求得x1及弛垂度D(公式推導略)。

考慮到高次方程運算的復雜度限制，在級數展開時通常只能選取前2至3項。上述解法不但誤差不好調整，在一些實際操作中也不夠靈活、方便。

設：

(23)

化簡上式可得：

(24)

設Z=sqrt(S×S+H×H)/L，調用二分搜索函數dichotomy(funz,0,3, Z)可搜索求解t，從而求得懸鏈系數a(函數funz()定義參照1.1.2節)。已知懸點垂直距離H推導公式：

(25)

可以參見例1、例2，建立搜索函數如下：

double funh(double x)

{ double y1,y2,x1,x2,t;

x1=(x+l[i]) / a[i];

x2=x/a[i];

t=a[i]*( cosh(x1)-cosh(x2) );

return t;

}

調用進行二分搜索dichotomy(funh,-3×a,3×a,H)，即可求得x1，進而弛垂度D可求。

2.2 實 例

例3：給定一組鏈長S，懸點距離L，懸點高度差H，求解懸鏈系數a、弛垂度D與低點的橫坐標X1，主函數程序描述如下:

#include ″stdio.h″

#include ″math.h″

#define CC 1e-3

double s[5]={ 30.000, 30.000, 30.000, 90.000, 90.000 },l[5]={ 20.000, 20.000, 20.000, 60.000, 60.000 },h[5]={ 0.000, 5.000, 10.000, 0.000, 30.000 },a[5], x[5], d[5];

int i;

main()

{ double t,z;

for ( i=0; i<5; i++)

{ z=sqrt( s[i]*s[i]-h[i]*h[i] ) / l[i];

t=dichotomy( funz,0,3,z );

a[i]=l[i]/(2*t);

x[i]=dichotomy( funh, -3*a[i], 3*a[i], h[i] );

d[i]=a[i]*( cosh(x[i]/a[i])-1 );

}

for ( i=0; i<5; i++)

printf(″ L=%6.3f, S=%6.3f, H=%6.3f, a=%6.3f, x=%6.3f, D=%6.3f ″,l[i], s[i], h[i], a[i], x[i], d[i]);

}

可以借助MATLAB，按2.1節中描述的數學方法求解懸鏈系數a。利用二分搜索程序在δ<10-2和δ<10-3精度要求下求解的懸鏈系數a2和a3。與MATLAB求解的懸鏈系數a1對比如表4所示。

表4 二分搜索法與常規數學法求解懸鏈系數對比

參考表4可見，利用數學解法求得的懸鏈系數a介于二分搜索程序誤差設置在10-2與10-3之間。表中省略了弛垂度D與低點的橫坐標X1的輸出。

3 懸鏈線應用實例

3.1 懸點受力分析

懸鏈線應用實踐中，懸點的受力情況是首要考慮的問題，如圖3所示。

圖3 懸點B受力分析

設鏈索線密度為ρ，則懸點B水平、垂直及切線方向受力Tx、Ty、T可表示如下：

(26)

(27)

(28)

(29)

在式(26)、式(28)中代入鏈長：

(30)

則：

(31)

(32)

式(26)、式(28)、式(29)、式(31)和式(32)經常在懸鏈線應用中被選擇使用。

3.2 懸鏈線應用實例—吊桿架設

在特定的地區、地勢、地理等多種環境下，無法使用現代化的塔吊設施，完成貨物的定點吊運只能采用傳統的架設吊桿方法。

已知預吊起的貨物重量WG，已具備的吊桿長度HL，吊桿重量為HG，兩根鏈索的線密度ρ，怎樣架設吊桿才能準確地將貨物吊運到指定位置P處。工程實踐中，可以通過多種方式將吊桿固定在兩根鏈索的同一平面上，本例只探討吊桿及兩根鏈索所在平面的架設問題，圖4為吊桿架設的平面圖。

圖4 吊桿架設平面圖及受力分析

實踐中常規的吊桿架設方式可有如下3種：

(1) 已知左右兩根鏈索的長度Sl與Sr，AB兩點已選定，水平距離2L，吊桿架設的支點O為AB中點，求解吊桿的水平夾角α，從而確定吊桿基點O與指定吊點P的水平距離OP，可以驗證鏈索長度是否合適。

(2) 已知左右兩根鏈索的長度Sl與Sr，吊桿的基點位置O已選定，求解AB兩點的水平距離2L，從而確定AB兩點的固定位置，在選定處固定鏈索即可。

(3) 吊桿的基點位置O已固定，AB兩點也已選定，已知一條鏈索的長度，求解另一條鏈索的長度，按要求配置鏈索即可。

下面僅就上述第一種架設方法給以實例說明，另外兩種架設方法與之相似，本文略。

例4：已知左右兩側鏈索線密度ρ，A、B兩點的水平距離2L，架設的吊桿長度HL，重HG。輸入兩根鏈索長度Sl、Sr，預吊起的貨物重WG，求吊桿水平夾角α及貨物吊點與吊桿支點水平距離PX。

分析圖4可知，左右兩側懸鏈線懸點水平距離Ll、Lr及懸點垂直距離H分別為：

Ll=L+HL×cos(α)

Lr=L-HL×cos(α)

H=HL×sinα

兩側懸鏈線系數分別為al、ar，懸點C在兩側懸鏈線上的水平坐標為Xl、Xr，則懸點C水平及垂直作用力Tl、Tr、Vl、Vr分別為：

Tl=al×ρTr=ar×ρ

懸點C處順時針力矩MS與逆時針力矩MN分別為：

MS=Tr×HL×sinα+(Vl+Vr+0.5×HG+WG)×

HL×cos(α)

MN=Tl×HL×sin(α)

吊桿平衡時MS=MN，可以整理等式為：

(Vl+Vr+0.5×HG+WG)×HL×cos(α)=

(Tl-Tr)×HL×sin(α)

下面程序中設置：

M1=(Tl-Tr)×HL×sin(α)

M2=(Vl+Vr+0.5×HG+WG)×HL×cos(α)

圖5為兩根鏈索分別拉直時的吊桿平面圖。

圖5 吊桿角度分析圖

通過余弦定理可求吊桿極限角度α1與α2：

吊桿水平角度α實際取值范圍應在α1與α2之間，即α1<α<α2。

例4采用二分搜索法在α1與α2之間搜索α(本例中為變量c)，函數funz()、funh()、dichotomy()可以參見前文，圖6為主函數程序流程圖。

圖6 例4主函數流程圖

例4程序描述如下：

/* a、s、l、x各數組0、1元素(例如：a[0]、a[1]、s[0]、s[1]…)分別表示左[右]側懸鏈線對應值，h為兩條懸鏈線共同的懸點高度差，函數funz()與funh()定義參照例1與例3。 */

#include ″stdio.h″

#include ″math.h″

#define CC 1e-3

double s[2], h, l[2], a[2], x[2];

int i;

double funz(double x) /*略*/

double funh(double x) /*略*/

double dichotomy(double (*fp)(),double x1,double x2,double p1) /*略*/

main()

{ double c1, c2, c, ss, p, wg, hl, hg, z, h, m1, m2, px;

double t[2],v[2];

ss=100;

hl=40, hg=200, p=10;

printf(″Enter SL & SR :″); scanf(″%lf %lf″,&s[0],&s[1]);

printf(″Enter WG= ″);

scanf(″%lf″,&wg);

c1=3.14-acos( (hl*hl+ss*ss/4-s[0]*s[0]) / (hl*ss) );

c2=acos( (hl*hl+ss*ss/4-s[1]*s[1]) / (hl*ss) );

do { double temp;

c=(c1+c2)/2;

l[0]=ss/2+hl*cos(c);

l[1]=ss/2-hl*cos(c);

h=hl*sin(c);

for(i=0;i<2;i++)

{ z=sqrt(s[i]*s[i]-h*h)/l[i];

temp=dichotomy(funz,0,3,z);

a[i]=l[i]/(2*temp);

x[i]=dichotomy(funh, -3*a[i], 3*a[i],h);

t[i]=a[i]*p;

v[i]=temp*sinh((x[i]+l[i]) / a[i]);

}

m1=(t[0]-t[1])*hl*sin(c);

m2=(v[0]+v[1]+hg/2+wg)*hl*cos(c);

if (m1>m2) c1=c;

else c2=c;

} while ( fabs(c1-c2)>1e-6);

px=hl*cos(c);

printf( ″C=%6.2f,PX=%6.2f″, c/3.14*180,px);

}

【程序輸入】

Enter SL & SR ：70 60

Enter WG ：500

【程序輸出】

C=78.55 PX=7.97

程序中已賦值A、B兩點距離SS=100 m，吊桿長度HL=40 m、吊桿重量HG=200 kg、吊桿線密度P=10 kg/m，實驗數據如表5所示。

表5 吊桿架設實驗數據

實踐中亦可以構造更詳細的表5，根據貨物重量的不同，查表調整兩根鏈索的長度，實現吊桿的架設。

4 結 語

懸鏈線問題核心是超越函數的求解。無論是怎樣的代數轉換，運算量和精度要求都難以把握。實踐證明，通過二分搜索程序設計輔助求解可以有效地解決上述問題，在實踐中是簡單可行的。

[1] 白興蘭,段夢蘭,李強.基于整體分析的鋼懸鏈線立管觸地點動力響應分析[J].工程力學,2014,31(12):249-256.

[2] Thethi R,Moros T.Soil interaction effects on simple catenary riser response[R].Deepwater Pipeline & Riser Technology Conference,Houston,Texas,2011:1-24.

[3] 李傳習,陳洪林,柯紅軍.空間纜索懸索橋水平母線鞍座設計位置的確定[J].中外公路,2015,25(6):90-94.

[4] Nader M,Maroney B.The New San Francisco-Oakland Bay Bridge[C]//Structures Congress,2013:588-598.

[5] 邱為鋼.懸鏈線的幾何特征[J].物理與工程,2015,25(3):28-34.

[6] 馬旭,楊文美.醫療特征圖像邊緣曲率分析實例[J].醫學影像學雜志,2010,20(11):1713-1715.

[7] 顏玉嬌,尚偉偉.6自由度繩索牽引并聯機器人的懸鏈線建模與動力學分析[J].中國科學技術大學學報,2015,45(7):546-554.

[8] Zi B,Ding H F,Cao J B,et al.Integrated mechanism design and control for completely restrained hybrid-driven based cable paralled manipulators[J].Journal of Intellignent and Robotic Systems,2014,74(34):643-661.

[9] 王丹,劉家新.一般狀態下懸鏈線方程的應用[J].船海工程,2007,36(3):26-28.

[10] 王建國,逄煥平,李雪峰.懸索橋線形分析的懸鏈線單元法[J].應用力學學報,2008,25(4):626-631.

[11] 王長昌.基于懸鏈線方程的跨接電纜長度計算[J].鐵道車輛,2013,51(6) 4-6,33.

[12] 陳常松,陳政清,顏東煌.懸索橋主纜初始位形的懸鏈線方程精細迭代分析法[J].工程力學,2006,23(8):62-68.

[13] 馬旭,張壽鵬.一種基于弧長與弦高的半徑函數逼近算法[J].遼寧大學學報,2016,43(2):130-133.

[14] 邢富沖.懸鏈線弛垂度的計算方法[J].數學的實踐與認識,2004,34(11):98-101.

[15] 馬旭,王大勇.趣味智能模擬程序設計實例[J].計算機與數字工程,2016,44(5):979-982.

[16] 林貴瑜,馮優達,苑登波.鋼絲繩拋物線理論與懸鏈線理論的計算誤差分析[J].建筑機械化,2010(3):42-43,55.

BINARYSEARCHMETHODFORSOLVINGCATENARYPROBLEM

Ma Xu Yi Su

(SchoolofInnovationandEntrepreneurship,LiaoningUniversity,Shenyang110036,Liaoning,China)

Catenary problem is one of the problems often encountered in modern engineering practice. In the two cases that two hanging points at the same height and different heights, the paper analyzes the conventional mathematical methods for solving the catenary coefficient, sag and the abscissa of the hanging points. Due to the large amount of calculation and the error is difficult to grasp, the paper proposes the binary search programming specific ideas to assist in solving, and offers three complete programs and data for different applications, compares and analyzes the mathematical methods and programming methods. Taking the boom erection project as an example, the paper provides a complete C program and the main function flowchart, and describes in more detail the procedures to solve the relevant parameters. The above method has been applied to practice many times, and it is simple and feasible.

Binary search method Catenary Catenary coefficient Sag Error analysis C programming

TP39

A

10.3969/j.issn.1000-386x.2017.09.011

2016-11-28。國家自然科學基金項目(11304135)；遼寧省自然科學基金項目(201602345)；沈陽市應用基礎研究專項項目(F15-199-1-04)。馬旭，高級實驗師，主研領域：計算機應用，人工智能。易俗，實驗師。