馬 旭 易 俗

(遼寧大學(xué)創(chuàng)新創(chuàng)業(yè)學(xué)院 遼寧 沈陽 110036)

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

馬 旭 易 俗

(遼寧大學(xué)創(chuàng)新創(chuàng)業(yè)學(xué)院 遼寧 沈陽 110036)

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

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

0 引 言

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

本文將就懸鏈線應(yīng)用中的多種實用形式，結(jié)合C語言描述，分別給出相關(guān)函數(shù)及應(yīng)用程序，求解懸鏈系數(shù)、弛垂度、懸索鏈長及懸點水平坐標(biāo)，并給出工程應(yīng)用實例。

1 兩懸點高度相同的情況

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

圖1 兩懸點等高的懸鏈線坐標(biāo)圖

懸鏈線的標(biāo)準(zhǔn)方程為：

(1)

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

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

1.1.1 常規(guī)數(shù)學(xué)方式求解弛垂度

弛垂度的標(biāo)準(zhǔn)方程為：

(2)

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

(3)

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

設(shè)：

(4)

(5)

則：

(6)

(7)

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

1.1.2 二分搜索法求解弛垂度

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

設(shè)：

(8)

(9)

化簡后可得：

(10)

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

double funz( double x )

{ double y;

y=sinh(x)/x;

return y;

}

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

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;

}

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

(11)

1.1.3 實 例

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

#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 數(shù)學(xué)求解與程序求解的比較

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

表1 二分搜索法與常規(guī)數(shù)學(xué)方法求解馳垂度對比

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

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

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

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

1.2.1 二分搜索法求鏈長

設(shè)：

(12)

(13)

化簡后可得：

(14)

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

double funv(double x)

{ double y;

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

return y;

}

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

(15)

1.2.2 實 例

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

#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]);

}

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

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

參考表3可見，搜索程序精度設(shè)為10-3時，誤差已經(jīng)可以在許多實踐環(huán)境中被接受，當(dāng)搜索程序精度設(shè)為10-5時，誤差就更小了。

2 兩懸點高度不同的情況

2.1 二分搜索法求懸點橫坐標(biāo)及弛垂度

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

圖2 兩懸點不等高的懸鏈線坐標(biāo)圖

由懸鏈線標(biāo)準(zhǔn)方程，可知：

(16)

(17)

(18)

(19)

(20)

(21)

即:

(22)

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

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

設(shè)：

(23)

化簡上式可得：

(24)

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

(25)

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

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;

}

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

2.2 實 例

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

#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節(jié)中描述的數(shù)學(xué)方法求解懸鏈系數(shù)a。利用二分搜索程序在δ<10-2和δ<10-3精度要求下求解的懸鏈系數(shù)a2和a3。與MATLAB求解的懸鏈系數(shù)a1對比如表4所示。

表4 二分搜索法與常規(guī)數(shù)學(xué)法求解懸鏈系數(shù)對比

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

3 懸鏈線應(yīng)用實例

3.1 懸點受力分析

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

圖3 懸點B受力分析

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

(26)

(27)

(28)

(29)

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

(30)

則：

(31)

(32)

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

3.2 懸鏈線應(yīng)用實例—吊桿架設(shè)

在特定的地區(qū)、地勢、地理等多種環(huán)境下，無法使用現(xiàn)代化的塔吊設(shè)施，完成貨物的定點吊運只能采用傳統(tǒng)的架設(shè)吊桿方法。

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

圖4 吊桿架設(shè)平面圖及受力分析

實踐中常規(guī)的吊桿架設(shè)方式可有如下3種：

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

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

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

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

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

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

Ll=L+HL×cos(α)

Lr=L-HL×cos(α)

H=HL×sinα

兩側(cè)懸鏈線系數(shù)分別為al、ar，懸點C在兩側(cè)懸鏈線上的水平坐標(biāo)為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(α)

下面程序中設(shè)置：

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

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

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

圖5 吊桿角度分析圖

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

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

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

圖6 例4主函數(shù)流程圖

例4程序描述如下：

/* a、s、l、x各數(shù)組0、1元素(例如：a[0]、a[1]、s[0]、s[1]…)分別表示左[右]側(cè)懸鏈線對應(yīng)值，h為兩條懸鏈線共同的懸點高度差，函數(shù)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(chǔ)、B兩點距離SS=100 m，吊桿長度HL=40 m、吊桿重量HG=200 kg、吊桿線密度P=10 kg/m，實驗數(shù)據(jù)如表5所示。

表5 吊桿架設(shè)實驗數(shù)據(jù)

實踐中亦可以構(gòu)造更詳細(xì)的表5，根據(jù)貨物重量的不同，查表調(diào)整兩根鏈索的長度，實現(xiàn)吊桿的架設(shè)。

4 結(jié) 語

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

[1] 白興蘭,段夢蘭,李強(qiáng).基于整體分析的鋼懸鏈線立管觸地點動力響應(yīng)分析[J].工程力學(xué),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] 李傳習(xí),陳洪林,柯紅軍.空間纜索懸索橋水平母線鞍座設(shè)計位置的確定[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] 馬旭,楊文美.醫(yī)療特征圖像邊緣曲率分析實例[J].醫(yī)學(xué)影像學(xué)雜志,2010,20(11):1713-1715.

[7] 顏玉嬌,尚偉偉.6自由度繩索牽引并聯(lián)機(jī)器人的懸鏈線建模與動力學(xué)分析[J].中國科學(xué)技術(shù)大學(xué)學(xué)報,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] 王丹,劉家新.一般狀態(tài)下懸鏈線方程的應(yīng)用[J].船海工程,2007,36(3):26-28.

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

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

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

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

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

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

[16] 林貴瑜,馮優(yōu)達(dá),苑登波.鋼絲繩拋物線理論與懸鏈線理論的計算誤差分析[J].建筑機(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。國家自然科學(xué)基金項目(11304135)；遼寧省自然科學(xué)基金項目(201602345)；沈陽市應(yīng)用基礎(chǔ)研究專項項目(F15-199-1-04)。馬旭，高級實驗師，主研領(lǐng)域：計算機(jī)應(yīng)用，人工智能。易俗，實驗師。