1、高等流體力學(xué)班 級： 機(jī)設(shè)15學(xué)碩班 學(xué) 號： 2015200813 姓 名： 張湘楠 授課老師： 畢新勝 日 期： 2016年7月 1日 一、研究報告內(nèi)容：1、0.618的證明、一維搜索程序作業(yè)；2、單位矩陣程序作業(yè)；3、連桿機(jī)構(gòu)問題 自行選擇小型機(jī)械設(shè)計問題或其他工程優(yōu)化問題；（1）分析優(yōu)化對象，根據(jù)設(shè)計問題的要求，選擇設(shè)計變量，確立約束條件，建立目標(biāo)函數(shù)，建立優(yōu)化設(shè)計的數(shù)學(xué)模型并編制問題程序；（2）選擇適當(dāng)?shù)膬?yōu)化方法，簡述方法原理，進(jìn)行優(yōu)化計算；（3）進(jìn)行結(jié)果分析，并加以說明。4、寫出課程實踐心得體會，附列程序文本。5、為響應(yīng)學(xué)校2014年度教學(xué)工作會議的改革要求，探索新的課程考核評價方

2、法，特探索性設(shè)立一開放式考核項目，占總成績的5%。 試用您自己認(rèn)為合適的方式（書面）表達(dá)您在本門課程學(xué)習(xí)方面的努力、進(jìn)步與收獲。（考評將重點關(guān)注您的獨創(chuàng)性、簡潔性與可驗證性）。二、研究報告要求1、報告命名規(guī)則：學(xué)號姓名機(jī)械優(yōu)化設(shè)計課程實踐報告.doc2、報告提交郵址：（收到回復(fù)，可視為提交成功）。追 求：問題的工程性，格式的完美性，報告的完整性。不 追 求：問題的復(fù)雜性，方法的惟一性。評判準(zhǔn)則：獨一是好，先交為好；切勿拷貝。目錄：0.618的證明、一維搜索程序作業(yè) 關(guān)于的證明4 一維搜索的作業(yè) 采用matlab進(jìn)行編程 5 采用C語言進(jìn)行編程 7單位矩陣程

3、序作業(yè) 采用matlab的編程 9 采用c語言進(jìn)行編程 9機(jī)械優(yōu)化工程實例 連桿機(jī)構(gòu)11 自選機(jī)構(gòu)16課程實踐心得 20附列程序文本 21進(jìn)步，努力，建議25一、0.618的證明、一維搜索程序作業(yè)關(guān)于的證明黃金分割法要求插入點，的位置相對于區(qū)間兩端具有對稱性，即其中為待定常數(shù)。此外，黃金分割法還要求在保留下來的區(qū)間內(nèi)再插入一點所形成的區(qū)間新三段，與原來的區(qū)間三段具有相同的比例分布。黃金分割法還要求在保留下來的區(qū)間內(nèi)再插一點所形成的區(qū)間新三段，與原來的區(qū)間三段有相同的比例分布。b1211-圖一設(shè)原區(qū)間的長度為1，如圖一所示，保留下來的區(qū)間長度為，區(qū)間縮短率為。為了保持相同的分別比例。插入新點應(yīng)在

4、位置上，在原區(qū)間的1位置應(yīng)相當(dāng)于在保留區(qū)間的位置。故有： a2231(1-)2圖二解得編寫0.618的程序，并計算下列問題程序框圖（1）采用MATLAB進(jìn)行編程%fun.m:黃金分割法求極值點%輸入數(shù)據(jù)% a 搜索區(qū)間下限% b 搜索區(qū)間上限% e 精度%輸出數(shù)據(jù)：% x 極小值點%其他常量：% c1,c2,c3 區(qū)間% r 黃金分割比例0.618%定義函數(shù)function x=fun(a,b,e) r=0.618; c1=b-r*(b-a);c2=a+r*(b-a); y1=f(c1);y2=f(c2); while (abs(b-a)/b)e)&(abs(y2-y1)/y2)e) if y

5、1=y2 a=c1;c1=c2;y1=y2; c2=a+r*(b-a);y2=f(c2); else b=c2;c2=c1;y2=y1; c1=b-r*(b-a);y1=f(c1); end end x=0.5*(a+b);end 對f函數(shù)的確立10 function y=f(x) y=(x-2)2+3;end11 function y=f(x) y=cosxend如果要計算y=(x-2)2+3;的黃金分割法，則需要將圖所示的f腳本中的函數(shù)寫成如所示如果要計算y=cosx；黃金分割法，則需要將圖所示的f腳本中的函數(shù)寫成如所示fun函數(shù)表示對matlab的主程序語言。教學(xué)f函數(shù)NO.2運算結(jié)果：

6、函數(shù)NO.1運算結(jié)果：（2）采用C語言進(jìn)行編程#include #define K 0.618double f(double); /*函數(shù)值計算函數(shù)聲明*/void main(void)double a,b,size;double a1,a2;int I;printf(“請輸入?yún)^(qū)間兩端點（端點值應(yīng)大于0）:”);scanf(“%lf,%lf”,&a,&b); /*輸入端點值*/printf(“請輸入精度:”);scanf(“%lf”,&size); /*輸入精度*/printf(“區(qū)間為(%lf,%lf),精度為%lfn”,a,b,size);printf(“序號t a1tt a2tt f(a

7、1)tt f(a2)n”);for(i=0;isize) /*用精度控制循環(huán)次數(shù)*/i+; a1=b-K*(b-a); /*按0.618法插入兩點*/ a2=a+K*(b-a); printf(“%2d:t%fta%ft%ft%fn”,I,a1,a2,f(a1),f(a2); /*輸出每次計算后a1,a2,f(a1),f(a2)的值*/if(f(a1)=f(a2)a=a1;elseb=a2;printf(“所求極小值點為:x=%lft極小值f(x)=%fn”,a,f(b);double f(double x) /*函數(shù)值計算函數(shù)*/double f;f=(x-2)*(x-2)+3;return

8、(f);對于y=cosx，須在程序中加一個#include（math.h）頭程序，以示我要調(diào)用函數(shù)語句。再將倒數(shù)第三行函數(shù)改為y=cos（x），再次編譯運算即可。運行如圖所示我們可以看出，大約在x=2處取到極值。與理論相符。運行如圖所示我們可以看出，大約在處取到極值。與理論相符。二、用簡單的語句寫一個單位矩陣采用matlab構(gòu)造：Matlab構(gòu)造的矩陣最是簡便，因為本身就有一個單位矩陣的函數(shù)。由于matlab是基于C語言而設(shè)立的一個數(shù)學(xué)運用軟件，所以他的集成度非常高，而且具有很好的開放性，于是我們得知eye函數(shù)，將其構(gòu)造。Eg：若要構(gòu)造一個3維的單位矩陣，則輸入eye（3）若要構(gòu)造一個n為的單

9、位矩陣，這輸入eye（n）用C語言構(gòu)造程序如圖所示：#include int main() int I,j,n; while(scanf(“%d”,&n) = 1) for(I = 0; I n; +i) for(j = 0; j n; +j) printf(“%d “,I = j); printf(“n”); return 0;三、機(jī)械優(yōu)化設(shè)計工程實例 連桿機(jī)構(gòu)問題（1）連桿機(jī)構(gòu)問題描述圖 1 機(jī)構(gòu)簡圖設(shè)計一曲柄連桿搖桿機(jī)構(gòu)，要求曲柄從時，搖桿的轉(zhuǎn)角最佳再現(xiàn)已知的運動規(guī)律：且=1，=5，為極位角，其傳動角允許在范圍內(nèi)變化。（2）數(shù)學(xué)模型的建立設(shè)計變量：這里有兩個獨立參數(shù)和。因此設(shè)計變量為目標(biāo)

10、函數(shù)：將輸入角分成30等分，并用近似公式計算，可得目標(biāo)函數(shù)的表達(dá)式約束條件：GX(1)=-X(1)0GX(2)=-X(2) 0GX(3)=-(X(1)+X(2)+6.00GX(4)=-(X(2)+4.0)+X(1) 0GX(5)=-(4.0+X(1)+X(2) 0GX(6)=-(1.4142*X(1)*X(2)-X(1)*2-X(2)*2)-16.00GX(7)=-(X(1)*2+X(2)*2+1.4142*X(1)*X(2)+36.00（3）程序編制C=SUBROUTINE FFX(N,X,FX)C=DIMENSION X(N)COMMON /ONE/I1,I2,I3,I4,NFX,I6NF

11、X=NFX+1P0=ACOS(1.0+X(1)*2-X(2)*2+25.0)/(10.0*(1.0+X(1)Q0=ACOS(1.0+X(1)*2-X(2)*2-25.0)/(10.0*X(2)T=90.0*3.1415926/(180.0*30.0)FX=0.0DO 10 K=0,30 PI=P0+K*TQE=Q0+2.0*(PI-P0)*2/(3.0*3.1415926)D=SQRT(26.0-10.0*COS(PI)AL=ACOS(D*D+X(2)*X(2)-X(1)*X(1)/(2.0*D*X(2)BT=ACOS(D*D+24.0)/(10.0*D) IF(PI.GE.0.0 .AND.

12、 PI.LT.3.1415926) THEN QI=3.1415926-AL-BT ELSE QI=3.1415926-AL+BT END IF IF(K.NE.0 .OR. K.NE.30)THEN FX=FX+(QI-QE)*2*T ELSE FX=FX+(QI-QE)*2*T/2.0 END IF12 CONTINUERETURNENDC=SUBROUTINE GGX(N,KG,X,GX)C=DIMENSION X(N),GX(KG)GX(1)=-X(1)GX(2)=-X(2)GX(3)=-(X(1)+X(2)+6.0GX(4)=-(X(2)+4.0)+X(1)GX(5)=-(4.0+X

13、(1)+X(2)GX(6)=-(1.4142*X(1)*X(2)-X(1)*2-X(2)*2)-16.0GX(7)=-(X(1)*2+X(2)*2+1.4142*X(1)*X(2)+36.0RETURNENDC=SUBROUTINE HHX(N,KH,X,HX)C=DIMENSION X(N),HX(KH)X(1)=X(1)RETURNEND輸入數(shù)據(jù)如下：2,7,04.3,3.22.,0.2,0.01,0.00001,0.000010,1,03.5,2.8,5.5,4.5（4） 可執(zhí)行程序的生成13 將程序文本保存成文件t001.txt，存放在OPT_EXAMFORTRAN0-fortran目

14、錄下。14 使用DOS操作界面，OPT_EXAMFORTRAN0-fortran目錄下，輸入命令 asumt.for+t001.txt t001.for 回車。生成t001.for文件。15 for1 t001; 回車；for2 回車；link t001; 回車，即生成可執(zhí)行程序t001.exe。從t001.dat文件中加載數(shù)據(jù)，運行t001.exe。t001t001.rtf 回車得到結(jié)果。操作過程如圖：操作步驟（5）結(jié)果與分析 結(jié)果顯示如下：= PRIMARY DATA = N= 2 KG= 7 KH= 0 X : .4300000E+01 .3200000E+01 FX: .7847605

15、E-02 GX: -.4300000E+01 -.3200000E+01 -.1500000E+01 -.2900000E+01 -.5100000E+01 -.6729390E+01 -.1218939E+02 X : .4300000E+01 .3200000E+01 FX: .7847605E-02 GX: -.4300000E+01 -.3200000E+01 -.1500000E+01 -.2900000E+01 -.5100000E+01 -.6729390E+01 -.1218939E+02 HX: PEN = .3974390E+01 R = .2000000E+01 C =

16、.2000000E+00 T0= .1000000E-01 EPS1= .1000000E-04 EPS2= .1000000E-04= ITERATION COMPUTE =* IRC = 1 R = .2000000E+01 PEN= .3974390E+01* IRC = 2 R = .4000000E+00 PEN= .3119361E+01* IRC = 3 R = .8000001E-01 PEN= .6407908E+00* IRC = 4 R = .1600000E-01 PEN= .1439420E+00* IRC = 5 R = .3200000E-02 PEN= .373

17、3451E-01* IRC = 6 R = .6400000E-03 PEN= .8986510E-02* IRC = 7 R = .1280000E-03 PEN= .2346609E-02* IRC = 8 R = .2560000E-04 PEN= .9042178E-03* IRC = 9 R = .5120001E-05 PEN= .5396702E-03* IRC =10 R = .1024000E-05 PEN= .4410887E-03* IRC =11 R = .2048000E-06 PEN= .4128061E-03* IRC =12 R = .4096001E-07 P

18、EN= .4071001E-03* IRC =13 R = .8192003E-08 PEN= .4014154E-03* IRC =14 R = .1638401E-08 PEN= .3992769E-03* IRC =15 R = .3276801E-09 PEN= .3987240E-03* IRC =16 R = .6553602E-10 PEN= .3984505E-03* IRC =17 R = .1310720E-10 PEN= .3983411E-03* IRC =18 R = .2621441E-11 PEN= .3982889E-03= OPTIMUM SOLUTION =

19、 IRC= 18 ITE= 35 ILI= 94 NPE= 447 NFX= 482 NGR= 0 R= .2621441E-11 PEN= .3982696E-03 X : .4135127E+01 .2315381E+01 FX: .3982623E-03 GX: -.4135127E+01 -.2315381E+01 -.4505072E+00 -.2180254E+01 -.5819746E+01 -.7079845E+01 -.3633157E-03 HX: Stop - Program terminated.從結(jié)果中得知當(dāng)X1=4.135127，X2=2.315381時，目標(biāo)函數(shù)取

20、得最小值, 搖桿的轉(zhuǎn)角最佳再現(xiàn)已知的運動規(guī)律：。 自選結(jié)構(gòu)模型設(shè)計某帶式輸送機(jī)減速器的高速級齒輪傳動。已知高速級輸入功率P1 = 10kW,小齒輪轉(zhuǎn)速n1 =960 r /min,傳動比i = 3. 2。齒輪材料和熱處理:大齒輪45號鋼(調(diào)質(zhì))硬度為217255HBS,小齒輪40Cr(調(diào)質(zhì))硬度為241286HBS,工作壽命15 年,假設(shè)每年工作300天,兩班制,帶式輸送機(jī)工作平穩(wěn),轉(zhuǎn)向不變。常規(guī)設(shè)計方案采用直齒圓柱齒輪: m=2.5, z1=30, d=1。解：（1）設(shè)計變量，單級直齒圓柱齒輪傳動的中心距 : f1X=121+mz1齒寬：f2x=mz1d將m, z1, d作為設(shè)計變量，即：

21、X=m，z1，dT=x1，x2，x3T2 ）目標(biāo)函數(shù)根據(jù)多目標(biāo)優(yōu)化的線性加權(quán)法建立體積最小的目標(biāo)函數(shù): f ( x) =1f1 ( x) +2f2 ( x)=11+x1x22 +2x1x2x3其中:1 ,2 是加權(quán)系數(shù),且1 +2 = 1,分別根據(jù)設(shè)計時徑向和軸向安裝位置的要求設(shè)定;取1 = 1表示要求中心距最小,取2 = 1則表示要求齒寬最小。3 ）約束條件模數(shù)的限制:對于傳遞動力的齒輪,通常要求模數(shù)不少于1. 5-2,得約束條件: G1x=2-x10 G2x=x1-1.50（2）小齒輪齒數(shù)的限制:小齒輪齒數(shù)z1應(yīng)不大于產(chǎn)生根切的最小齒數(shù)17 ,得約束條件: G3x=17-x20（3）齒寬系

22、數(shù)的限制:由于dmin d dmax ,約束條件為: G4x=x3-1.150 G5x=0.7-x30（4）齒面接觸強(qiáng)度的限制，根據(jù)公式并查表得約束條件： G6x=1.04107-2.916105x13x23x330 G7x=1.04107-8.95106x13x23x330（5）齒根彎曲強(qiáng)度的限制，根據(jù)公式查表得約束條件： G8x=1.51106-303.57x13x23x330 G9x=1.42106-2445.92x13x23x3304 建立數(shù)學(xué)優(yōu)化模型高速級齒輪傳動多目標(biāo)優(yōu)化設(shè)計的數(shù)學(xué)模型為:（1 取0.6，2取0.4）Fun(x)=min1 1+3.2x1x22 +2 x1x2x3X

23、=m，z1，dT=x1，x2，x3TG1x=2-x10G2x=x1-1.50G3x=17-x20G4x=x3-1.150G5x=0.7-x30G6x=1.04107-2.916105x13x23x330G7x=1.04107-8.95106x13x23x330G8x=1.51106-303.57x13x23x330G9x=1.42106-2445.92x13x23x3305編寫程序并運行結(jié)果目標(biāo)函數(shù)M文件：function f=zhwm(x)f=0.6*2.1*x(1)*x(2)+0.4*x(1)*x(2)*x(3);約束函數(shù)M文件：function c ceq=zhwy(x)c(1)=1.0

24、4*107-2.916*105*(x(1)*x(2)3*x(3);c(2)=1.04*107-8.95*106*(x(1)*x(2)3*x(3);c(3)=1.51*106-303.57*x(1)3*x(2)2*x(3);c(4)=1.42*106-2445.92*x(1)3*x(2)2*x(3);ceq=;優(yōu)化函數(shù)M文件：x0=2 32 1;lb=1.5 17 0.7;ub=2 inf 1.15;u=;x,fval=fmincon(zhwm,x0,lb,ub,zhwy)約束函數(shù)目標(biāo)函數(shù)經(jīng)過Matlab優(yōu)化并圓整后的齒輪參數(shù)如下：設(shè)計參數(shù)模數(shù)m/mm齒數(shù)z1齒寬系數(shù)d中心距/mm普通設(shè)計2.5

25、301157.5MATLAB優(yōu)化1.791127.43771.149984.5254圓整后參數(shù)1.75301117.6經(jīng)過計算，最小體積為87.15。四、課程學(xué)習(xí)心得從懶懶散散的寒假到匆匆忙忙的考試周，仿佛只是在彈指一揮間。然而我們就在這段時間內(nèi)，又學(xué)了幾門課程，而這些課程中，最有趣，最好玩的當(dāng)然就是機(jī)械優(yōu)化設(shè)計了機(jī)械優(yōu)化設(shè)計是一門理論性非常強(qiáng)的一門課，剛開始的時候不得要領(lǐng)，但是在王衛(wèi)榮老師的帶領(lǐng)下，逐漸登堂入室，漸漸有所了解。首先，這門課所依托的最強(qiáng)武器便是數(shù)學(xué)。數(shù)學(xué)即是方法，思維，和邏輯。其次所依托的便是算法程序語言。歸根結(jié)底的還是數(shù)學(xué)，因為數(shù)學(xué)賦予了算法以邏輯的力量，分析的過程，以及所能

26、解決的一切問題。還好我的數(shù)學(xué)基礎(chǔ)還算可以，邏輯思維還算通順，再加上老師所給予的方法，學(xué)習(xí)這門課程，一路上雖跌跌撞撞，但也有所收獲。所謂收獲，從方法論的角度上講，即時通過一個通用模型，解決一系列的問題，這是數(shù)學(xué)建模給我們的要求嗎，同時也是優(yōu)化這門課所賦予的核心思想。大二的時候，特別想搞建模，可是由于種種原因，其中最重要的原因就是懶惰，與之失之交臂。優(yōu)化這門課則彌補(bǔ)了我的部分遺憾。因為沒有建模的要求高，于是上手比較容易，同時又與我們自己的機(jī)械方面的課程緊密地結(jié)合在一起，過程輕車熟路。于是優(yōu)化真是一門寓學(xué)于樂的科目，只恨課程太多時間太少，不然真要好好地研究一番。說道大作業(yè)，那可是真的難，第一問和第二

27、問還好一點，第三問就不行了。說來慚愧，連桿機(jī)構(gòu)問題是借鑒了上屆學(xué)長的做法才學(xué)會的。因為太難了，所以就沒有用fortran語言做這些問題，全部換成了matlab。真要感謝優(yōu)化設(shè)計這門課程，如果不是這門課程的話，即使MATLAB已經(jīng)學(xué)過了，我估計也不會回過頭來復(fù)習(xí)，是優(yōu)化設(shè)計讓我重新拾起matlab，并且我可以說，凡是在優(yōu)化大作業(yè)中使用過的算法，程序，我已經(jīng)忘不掉了。因為fortran語言是一種相對陌生的語言，所以使用fortran變成的話有些力不從心，而我又是一個刨根究底的人，對于0.618的算法和單位矩陣的生成，我都是用的兩種算法，這樣保證思維的發(fā)散性。Matlab真的是一款很強(qiáng)大的軟體，他講

28、機(jī)會用到的所有常見數(shù)學(xué)模型全部都涵括了，所以，我們，親切地稱它為，萬能數(shù)學(xué)工具。當(dāng)學(xué)習(xí)變成一種樂趣的時候，每當(dāng)我們想做這種游戲的時候，時間，精力，甚至思路，仿佛全部有了，這就是積極所帶來的意義。學(xué)習(xí)優(yōu)化設(shè)計，傳達(dá)出這樣的信號：重要的不是考試的結(jié)果，而是解決問題的過程。所以，我覺得優(yōu)化設(shè)計這門課的評價方式，在我進(jìn)入大學(xué)以來，是第一次遇到，大作業(yè)，小作業(yè)的比值，竟然超過了考試成績的比值。換句話說，老師認(rèn)為，即使考試能考的很好，但是大作業(yè)，小作業(yè)不行，說明這只是一個會考試的工具，對于解決實際問題的能力還是有所欠缺。對！這樣的評價很接地氣。優(yōu)化設(shè)計這門課已經(jīng)結(jié)課了，考試已經(jīng)結(jié)束了。把大作業(yè)交上去，仿佛

29、這門課與我們就沒有什么瓜葛了一樣。不，絕對不是。解決問題的能力才剛剛開始。五、附列程序文本黃金分割法求極值點的matlab程序%fun.m:黃金分割法求極值點%輸入數(shù)據(jù)% a 搜索區(qū)間下限% b 搜索區(qū)間上限% e 精度%輸出數(shù)據(jù)：% x 極小值點%其他常量：% c1,c2,c3 區(qū)間% r 黃金分割比例0.618%定義函數(shù)function x=fun(a,b,e) r=0.618; c1=b-r*(b-a);c2=a+r*(b-a); y1=f(c1);y2=f(c2); while (abs(b-a)/b)e)&(abs(y2-y1)/y2)e) if y1=y2 a=c1;c1=c2;y

30、1=y2; c2=a+r*(b-a);y2=f(c2); else b=c2;c2=c1;y2=y1; c1=b-r*(b-a);y1=f(c1); end end x=0.5*(a+b);end 對f函數(shù)的確立function y=f(x) y=(x-2)2+3;endfunction y=f(x) y=cosxend黃金分割法求極值點的C語言程序#include #define K 0.618double f(double); /*函數(shù)值計算函數(shù)聲明*/void main(void)double a,b,size;double a1,a2;int i;printf(請輸入?yún)^(qū)間兩端點（端點值

31、應(yīng)大于0）:);scanf(%lf,%lf,&a,&b); /*輸入端點值*/printf(請輸入精度:);scanf(%lf,&size); /*輸入精度*/printf(區(qū)間為(%lf,%lf),精度為%lfn,a,b,size);printf(序號t a1tt a2tt f(a1)tt f(a2)n);for(i=0;isize) /*用精度控制循環(huán)次數(shù)*/i+; a1=b-K*(b-a); /*按0.618法插入兩點*/ a2=a+K*(b-a); printf(%2d:t%fta%ft%ft%fn,i,a1,a2,f(a1),f(a2); /*輸出每次計算后a1,a2,f(a1),f

32、(a2)的值*/if(f(a1)=f(a2)a=a1;elseb=a2;printf(所求極小值點為:x=%lft極小值f(x)=%fn,a,f(b);double f(double x) /*函數(shù)值計算函數(shù)*/double f;f=(x-2)*(x-2)+3;return(f);對于y=cosx，須在程序中加一個#include（math.h）頭程序，以示我要調(diào)用函數(shù)語句。再將倒數(shù)第三行函數(shù)改為y=cos（x），再次編譯運算即可。單位矩陣的matlab程序eye（x）單位矩陣的C語言程序#include int main() int I,j,n; while(scanf(“%d”,&n) =

33、 1) for(I = 0; I n; +i) for(j = 0; j n; +j) printf(“%d “,I = j); printf(“n”); return 0;連桿機(jī)構(gòu)的fortran語言C=SUBROUTINE FFX(N,X,FX)C=DIMENSION X(N)COMMON /ONE/I1,I2,I3,I4,NFX,I6NFX=NFX+1P0=ACOS(1.0+X(1)*2-X(2)*2+25.0)/(10.0*(1.0+X(1)Q0=ACOS(1.0+X(1)*2-X(2)*2-25.0)/(10.0*X(2)T=90.0*3.1415926/(180.0*30.0)FX

34、=0.0DO 10 K=0,30 PI=P0+K*TQE=Q0+2.0*(PI-P0)*2/(3.0*3.1415926)D=SQRT(26.0-10.0*COS(PI)AL=ACOS(D*D+X(2)*X(2)-X(1)*X(1)/(2.0*D*X(2)BT=ACOS(D*D+24.0)/(10.0*D) IF(PI.GE.0.0 .AND. PI.LT.3.1415926) THEN QI=3.1415926-AL-BT ELSE QI=3.1415926-AL+BT END IF IF(K.NE.0 .OR. K.NE.30)THEN FX=FX+(QI-QE)*2*T ELSE FX=

35、FX+(QI-QE)*2*T/2.0 END IF10 CONTINUERETURNENDC=SUBROUTINE GGX(N,KG,X,GX)C=DIMENSION X(N),GX(KG)GX(1)=-X(1)GX(2)=-X(2)GX(3)=-(X(1)+X(2)+6.0GX(4)=-(X(2)+4.0)+X(1)GX(5)=-(4.0+X(1)+X(2)GX(6)=-(1.4142*X(1)*X(2)-X(1)*2-X(2)*2)-16.0GX(7)=-(X(1)*2+X(2)*2+1.4142*X(1)*X(2)+36.0RETURNENDC=SUBROUTINE HHX(N,KH,X,HX)C=DIMENSION X