吳 珂, 高守棟

(1.安徽建筑大學(xué)土木工程學(xué)院，安徽 合肥 230601；2.廣州大學(xué)土木工程學(xué)院，廣東 廣州 510006)

0 引 言

對受壓構(gòu)件穩(wěn)定性分析方法有：中性平衡法、微分方程法、能量法、變分法、勢能駐值法、差分法、有限單元法等[1]。

差分法又稱差商法，是一種用差商近似導(dǎo)數(shù)的方法，分為向前差商型微分公式、向后差商型微分公式和中心差商型微分公式[2]。將一端固定一端鉸接的軸心受壓桿件劃分為n段，運用差分公式使得微分方程轉(zhuǎn)化為代數(shù)方程的形式，從而將微分方程的問題簡化為代數(shù)方程組的問題。但是，該方程組為含有一個未知量的代數(shù)方程組，若使用傳統(tǒng)方法，求解代數(shù)方程組工作量大，求解一元高次多項式的根也很困難[1]。

受一篇文章的啟發(fā)[3]，在差分法的基礎(chǔ)上，先使用MATLAB求出代數(shù)方程組的解，該解是含有一個未知量的高次多項式。再使用人工智能算法中的遺傳算法求解一元高次多項式的最小正根，從而得到一端固定一端鉸接的軸心受壓桿件的臨界荷載Pcr。差分法結(jié)合遺傳算法求Pcr，求解思路清晰，過程簡單。

1 代數(shù)方程組的求解[1]

1.1 n=2時，代數(shù)方程組的求解

函數(shù)y=f(x)的圖形見圖1，縱坐標(biāo)y(i-2h)、y(i-h)、y(i)、y(i+h)、y(i+2h)分別表示 函數(shù)f(x)在i-2h、i-h、i、i+h、i+2h處的函數(shù)值，h為各相鄰點間的距離，于是有二階導(dǎo)數(shù)的差分公式

(1)

四階導(dǎo)數(shù)的差分公式為

(2)

下端固定、上端鉸接的軸心受壓桿件圖見圖2，其四階微分方程為

(3)

將公式(2)和(3)代入(1)得到任一點i的差分方程為

ui+2h-4ui+h+6ui-4ui-h+ui-2h+

(4)

將桿件分為兩段，其圖形見圖3，此時n=2，i=1，h=L/n=L/2，由公式(4)和邊界條件可得方程組

u0=0

u2=0

u-1-u1=0

u1+u3=0

(5)

(6)

令A(yù)x=b，其中

A是一個5階方陣，且含有未知量T。為了使得Ax=b有非零解，則必有A的行列式為0，det(A)=0 ?6-2T=0。

圖1

1.2 n>2時，代數(shù)方程組的求解

當(dāng)n>2，表示桿件劃分的段數(shù)≥3，此時矩陣A的階數(shù)≥6，如果采用手算則計算量很大，且容易出錯。采用MATLAB中的函數(shù)det(A)[6]，會返回A的行列式的值，能夠快速準(zhǔn)確的計算A的行列式。調(diào)用的代碼為f=det(A)，即將函數(shù)det(A)返回的值賦予f。n=3時，方程組Ax=b的形式為：

(7)

圖2

圖3

表1 矩陣A行列式的值

圖4

2 遺傳算法求解一元高次多項式的根

2.1 遺傳算法簡介

1975年，Michigan大學(xué)的J. H. Holland對遺傳算法進(jìn)行了系統(tǒng)的研究，并發(fā)表論文[7]闡述了遺傳算法的基本理論和基本方法。De Jong是Holland的學(xué)生，他在其老師研究的基礎(chǔ)上，使用計算機做了大量數(shù)值實驗并得出具有指導(dǎo)意義的結(jié)論[8]。在計算機技術(shù)快速進(jìn)步的80年代，遺傳算法在工程優(yōu)化領(lǐng)域、經(jīng)濟預(yù)測領(lǐng)域、機器學(xué)習(xí)等領(lǐng)域得以成功應(yīng)用。90年代至今，遺傳算法已經(jīng)形成一套完整的體系，并隨著人工智能技術(shù)的發(fā)展不斷完善。

2.2 運用遺傳算法求T[9]

通過GUI使用MATLAB遺傳算法工具箱[10]，版本為2014A，對表1中關(guān)于T的多項式進(jìn)行求解。

當(dāng)n=2時，f(T)=6-2T。由于ga函數(shù)只解決最小值的問題，所以直接使用f(T)=6-2T作為目標(biāo)函數(shù)求T是不可取的，必須要對f(T)作適當(dāng)變換。這里采用加絕對值的方法，將求解f(T)=6-2T=0的根的問題轉(zhuǎn)化為求g(T)=|f(T|)=|6-2T|的最小值的問題，求得T的值為3。

當(dāng)n=3時，f(T)=-3T2+16T-19,g(T)=|-3T2+16T-19|。函數(shù)f(T)有兩個根，使用ga函數(shù)能夠把所有的根求出。隨著n的增大，f(T)的根會越來越多，要求的是f(T)的最小正根，因此需要對ga函數(shù)施加約束。約束條件為T2，求解的結(jié)果T=1.78。需要注意的是，n不同則約束條件也不同，這里用g(T)的函數(shù)圖像來修正約束條件。n=2，3，…，10時T的值見表2。

表2 T值

3 臨界荷載Pcr及誤差計算[2]

表3 Pcr值及誤差值

4 里查森外推法[1]

分析表3可知，當(dāng)n大于5時，|er|的收斂速度顯著降低，出現(xiàn)工作量增加但效果不佳的情況。為了減少計算工作量的同時快速減小誤差，使用里查森外推法來計算Pcr。里查森外推公式為

(8)

式中n1為第一次的分段數(shù)，所得到的近似解為β1，n2為第二次的分段數(shù)，近似解為β2。使用里查森外推公式需要用到兩次計算得到的近似解，這部分工作已經(jīng)在前面做好了。

綜合數(shù)據(jù)和圖4可知：(1)誤差由1.93%降低到0.540%，說明里查森外推法在減少計算工作量的同時能夠快速減小誤差。(2)兩條曲線的斜率在n=10處均發(fā)生了突變，說明里查森外推法加快了求解的收斂速度。

5 結(jié) 論

(1)將差分法同遺傳算法相結(jié)合應(yīng)用于求解軸心受壓桿件的臨界荷載是一次全新的嘗試，計算了n=2，…，10時的臨界荷載值，使用里查森進(jìn)行了一次外推，并將計算結(jié)果與精確解進(jìn)行比較分析，分析的結(jié)果表明這種新的嘗試得到的結(jié)果是可靠的。

(2)使用遺傳算法在求解方程根的過程中需要對目標(biāo)函數(shù)施加約束條件，且約束條件是隨著方程的改變而變化的，這給求解過程帶來了較大障礙。如何更好的解決這個問題，將是下一步研究的工作。