最優(yōu)化方法課程實驗報告

項目維搜索算法(一)［實驗目的］編寫加步探索法、對分法、Newton法的程序。［實驗準備］1.掌握一維收搜索中搜索區(qū)間的加步探索法的思想及迭代步驟；2.掌握對分法的思想及迭代步驟；3.掌握Newton法的思想及迭代步驟。［實驗內容及步驟］編程解決以下問題：1.用加步探索法確定一維最優(yōu)化問題min9(t)=13—2t+1t>0的搜索區(qū)間，要求選取10=°，h0=1,a=2.加步探索法算法的計算步驟:⑴選取初始點t￡［0，+8)(或，￡U［°，t］)計算甲=9(t)給出初始步長h〉0I/J>1-1八、、00max,I牙00?l口a*h0A0,加步系數以>1,令k=0。比較目標函數值.令t=t+h,計算9=9(t)若9<9，轉(3)，否k+1kkk+1k+1k+1k貝膊專(4)。(o\加十捽步中甘Ah=ah同時At=t,t=t,k=k+1姓加大探索步長?令k+1k，同時，令k'kk+1,,轉(2)。反向探索.若k=0,轉換探索方向,令氣=-氣,t=\，轉(2)。否則，停止迭代，令a=代，令a=min{t，t}，b=max{t，t}加步探索法算法的計算框圖程序清單加步探索法算法程序見附錄1實驗結果運行結果為：I'D:\Pnogr3mFile&(croso-HVisualStudio\k/lyPr該問題的根的搜索空間是：［^31PressanyIce-tocontinue.2.用對分法求解min中(t)=t(t+2),已知初始單谷區(qū)間［a，b］=［—3,5］,要求按精度￡=。.3,￡=0.001分別計算.對分法迭代的計算步驟：確定初始搜索區(qū)間［a,b］,要求中'(a)<0,中啊>0。計算［a,b］的中點c=1(a+b)2'若甲'(c)V0，則a=c,轉(4);若甲'(c)=0，則t*=c,轉(5);若甲'(c)>0，則b=c，轉(4).

1/八(4)右|a-bl<e，則t*=-(a+b),轉(5);否則轉(2).^2⑸打印t*，結束程序清單對分法程序見附錄2實驗結果運行結果為：I'D:\PnogramFiles(x86AMicrc￡oftVisualStudic\MyProject世解的結果是：TPressanykmytocontinue3.用Newton法求解min中(t)=13—2t+1/已知初始單谷區(qū)間[a，b]=[0，1]，要求精度8=0.01.Newton法的計算步驟

⑴確定初始搜索區(qū)間［。,仞，要求甲’(a)<0，平’(b)>0(2)選定t0⑶計算t='?！准?(t0)若11-10l^e，則t0=t，轉(3);否則轉(5).打印t，中(t)，結束.Newton法的計算框圖程序清單Newton法程序見附錄3實驗結果運行結果為:項目二一維搜索算法(二)［實驗目的］編寫黃金分割法、拋物線插值法的程序。［實驗準備］掌握黃金分割法的思想及迭代步驟；掌握拋物線插值法的思想及迭代步驟。［實驗內容及步驟］編程解決以下問題：1.用黃金分割法求解min中(t)=t(t+2),已知初始單谷區(qū)間［a，b］=［—3,5］,要求精度8=0.001.黃金分割法迭代步驟:確定中(t)的初始搜索區(qū)間［a，b］.計算12=a+0.382(b-a)計算ti=a+0.618(b-a)若1ti-t2l<e,則打印t*=1^2，結束；否則轉(5).判別是否滿足甲V甲：若滿足，則置a=t，t=t，甲=平1222121然后轉(3);否則，置b=t，t=t，甲=甲，t=a+。(b—a)，甲=甲(t)11212222然后轉(4).黃金分割法的計算框圖：程序清單黃金分割法程序見附錄4實驗結果運行結果為：2.用拋物線插值法求解minf(x)=8x3—2x2—7x+3已知初始單谷區(qū)間［",』］-心二頃.拋物線插值法的計算步驟：甲(t)〈甲(t),所以相對t來說t是好點，故劃掉區(qū)間［t,t］,保留［t,t］為新區(qū)TOC\o"1-5"\h\z000210間,故置t=t，甲(t)=甲(t),t=t，甲(t)=甲(t),t保持不變;2020001甲(t)＞甲(t),所以相對t來說t是好點，故劃掉區(qū)間［t,t］,保留［t,t］為新區(qū)間，0012故置t1=t，^()=中(t1),t0與12保持不變；程序清單拋物線插值法程序見附錄5實驗結果運行結果為:項目三常用無約束最優(yōu)化方法(一)［實驗目的］編寫最速下降法、Newton法(修正Newton法)的程序。［實驗準備］.掌握最速下降法的思想及迭代步驟。.掌握Newton法的思想及迭代步驟;.掌握修正Newton法的思想及迭代步驟。［實驗內容及步驟］編程解決以下問題：用最速下降法求minf(X)=x2+25x;，X0=［2,2］r，￡=0.01.最速下降法計算步驟(1)取初始點X(0)，容許誤差(精度)￡>0，令k=0(2)計算p(k)=-Vf(X(k))(3)檢驗||p(k<e?若是迭代終止,取X*=X(k),否則轉4。(4)求最優(yōu)步長孔：minf(X(k)+Xp(k))=f(X(k)+Rp(k))(一維搜索)令X(k+i)=X(k)+Xp(k),令k:=k+1,轉20k最速下降法的計算框圖最速下降法程序見附錄6實驗結果運行結果為:2.用Newton法求minf(X)二60一10x-4x+X2+%2一%%初始點X0=[0，0]T，8=0.01Newton法的計算步驟(1)給定初始點x(0),及精度8>0，令k=0；(2)若HVf(X(k))1~8,停止，極小點為％(k)，否則轉步驟(3);(3)計算N2f(X(k))]T,令s(k)=-[H(X(k))]TVf(X(k))；令%(k+i)=%(k)+s(k),k=k+1,轉步驟(2)。程序清單Newton法程序見附錄7實驗結果運行結果為:Frcssanykeytocontinuc3.用修正Newton求minf(X)=4(%+1)2+2(%一1)2+%+%+10,

1212初始點X0=[0，0]T，8=0.01.修正Newton的計算步驟(1)給定初始點%(0),及精度8>0,令k=0;(2)若HVf(X(k叫~8,停止，極小點為％(k)，否則轉步驟(3);

(3)計算N^f(X(人))J',令S(人)=-[H(X(k))J1巧(X(k))；(4)用一維搜索法求a,使得f(X(k)+a(Sk滬minfX^a)Sk,)令a>0X(k+1)=X(k)+a(k)S(k),k=k+1,轉步驟(2)。程序清單修正Newton程序見附錄8實驗結果運行結果為：T'DeprogramFiles；x8c)\Micro=oftVis-uaJStudicVMyPrcyects^zuiyciijXC^bLg^zryuan.exe'向量相乘，兩矩陣相乘結果為、二狀函藪FM湖農頊系數'S-iS^—農頊系數'-12.375常數項系數?16椽盅newt<m法函數ESL或〉T就』*1十2<X2-1)A2#S1■*?2*10一的極小點為i自堡童取值為:向量相乘，兩矩陣相乘結果為、>rl=-lA25直代次數為：1Pi'bisean事keytocointinLLB項目四常用無約束最優(yōu)化方法（二）實驗目的編寫共軛梯度法、變尺度法（DFP法和BFGS法）程序。實驗準備1.掌握共軛方向法的思路及迭代過程；掌握共軛梯度法的思想及迭代步驟；3.掌握DFP法和BFGS法的思想及迭代步驟。實驗內容及步驟編程解決以下問題：1.用共軛梯度法求得min(時+4弋),取初始點％=[1，1]r,e=0.01.共輛梯度法算法步驟(1)給定初始點x(0),及精度e>0;(2)若|阿(x(0))1<￡,停止,極小值點為X(0),否則轉步驟(3);(3)取p(0)=—Nf(x(0)),且置k=0；(4)用一維搜索法求t，使得f(xk)+tpk)=minfx+tp)令，X(k+1)=x(k)+tp(k),轉步驟5；k(5)若||Vf(x(k+1))|<e,停止,極小值點為x(k+1),否則轉步驟(6);(6)若k+1=n，令x(0)=x(n),轉步驟(3),否則轉步驟(7);Vf(x(k+1))12(7)令p(k+1)=—Vf(x(k+1))+人p(k),人=,置k=k+1,轉步驟kk||Vf(x(k))||2。程序清單共軛梯度法程序見附錄9實驗結果運行結果為：

D:\ProgramFiles(x86)\IMicro5oftVisualStudio\MyPrqject5\zuiyou\Debug\ziyuan.exe"阿輸入函數的元數值忡請輸入初始值：11=-B.000000,pl=fl.013131304^p2=B_ra.b：l=I-4.500BB0.1-500000]=B.4?4238xl=^0.060000,x2-000000,pi=0?000090^p2=-0.000005[a^b]-[-4.566660,1.500000]pl-0-666606,p2—6.600005020913最代解為xi-g.GeGeGo^xE-e.eoeoea最教色匡蘇值為白.GGGOOOrressanjpkeytocontimiE：—2.用共軛梯度法求min/（X）=2^+X2-x1七，自定初始點，e=0.01.程序清單共軛梯度法程序見附錄9實驗結果運行結果為:rr"D:\ProgramiFiles(K86)\MicrD5aftVisualStudioVMyPT^jertsXzui^ouXDebugXziyuan.exe8清輸A函分的打數道n=W_.j—一—清輸A函分的打數道n=W睛輸入初始值:Y1=n一MMWMMMrY2=H.1買=-M.MWMMMMP=-rt.MMMM11Pt=rt.僵技由函數0-B300G0Fressany虹母tncontinue.3.用DFP法求minf(X)=4(%—5)2+(氣-6)2,初始點X0=[8,9]…0.01.DFP法的具體迭代步驟如下:(1)給定初始點對),迭代精度，，維數n(2)置0—k,單位矩陣I—,計算矽?。(3)計算搜索方向-?(4)進行一維搜索求危，使得迭代新點濟3氣舟"7(5)檢驗是否滿足迭代終止條件H'l\lII?若滿足,停止迭代,輸出最優(yōu)解野T),_廠，；否則進行下一步。(6)檢查迭代次數，若k=n，則置x^-xin),轉向步驟(2);若k<n，則進行下一步。(7)計算:-必'TA'然后，置k+1-k,轉向步驟（3）。DFP法的計算框圖程序清單DFP法程序見附錄10實驗結果運行結果為:DFF法函數”心“2〉=4牌1-5〉七+牌2-6〉七.的極小點為:f=0刑得板小點為日變量卅伯.為，xl=5^2=6點目nukeutncritinu皓一項目五常用約束最優(yōu)化方法［實驗目的］編寫外點罰函數法、外點罰函數法的程序。［實驗準備］掌握外點罰函數法的思想及迭代步驟；掌握內點罰函數法的思想及迭代步驟。［實驗內容及步驟］編程解決以下問題：1.用外點罰函數法編程計算minf(X)=-尤+尤，Jg(X)=ln尤>0,.［.(X)=氣+%-1=0，精度e=10-5.外點法的計算步驟：(1)給定初始點x(0),初始罰因子從1),放大系數c>1;允許誤差e>0,設置k=1;(2)以x(k-1)作為搜索初始點，求解無約束規(guī)劃問題minf(x)+人P(x)，令X(k)為所求極小點。(3)當人P(x(們)<°,則停止計算,得到點x(k)；否則，令Xk+1)=ck(k),返回(2)執(zhí)行。程序清單外點罰函數法程序見附錄11實驗結果實驗結果為運行結果為：：s(x86hMicroso-HVisualStudio\MyProjects^zuiyou\Debug\ziyLjan.exe'擔外點法求解：tiinf<X>=-X1+X2gl<X>=lnX2>=0s.t.g2<X>=Xl+X2-l=0結果是：x*=7.0567e-0171f<x*>=1Pressanykeytocontinue

2.用內點罰函數法編程計算min[3(七+1)3+%」，0.1JX-1>0,

X>0.0.1初始點取為X0=[3,4H，初始障礙因子取"1=10，縮小系數取C內點罰步驟：給定初始內點X(0)GS，允許誤差e>0,障礙參數Y⑴，縮小系數bG(0,1)，置k=1；以X(k-1)為初始點，求解下列規(guī)劃問題：minf(x)+y(k)B(x)人，令X(k)為所求極小點st.XGS(3)如果y(k)B(X(k))<e，則停止計算，得到結果X(k),(4)否則令y(k+1)=by(k),置k=k+1,返回(2)。內點罰計算框圖

程序清單內點罰函數法程序見附錄12實驗結果運行結果為:附錄1#include<iostream.h>#include<LIMITS.H>#defineMAX20480doublefun(doublex){returnx*x*x-2*x+1;}doubleMax_Num(doublea,doubleb){if(a>b)returna;elsereturnb;}doubleMin_Num(doublea,doubleb){if(a<b)returna;elsereturnb;}voidStepAdding_Search(double&a,double&b){doublet[MAX]={0};doubleh[MAX]={0};doublef[MAX]={0};doubleresult=0;doublep=2;t[0]=0;h[0]=1;f[0]=fun(t[0]);for(intk=0;k<MAX-1;k++){t[k+1]=t[k]+h[k];f[k+1]=fun(t[k+1]);if(f[k+1]<f[k]){h[k+1]=p*h[k];result=t[k];t[k]=t[k+1];f[k]=f[k+1];}elseif(k==0){h[k]=-h[k];result=t[k+1];}else{a=Min_Num(result,t[k+1]);b=Max_Num(result,t[k+1]);}}}intmain(){doublea=0.0;doubleb=0.0;StepAdding_Search(a,b);cout<<"該問題的根的搜索空間是：["<<a<<”,"<<b<<”]\n”;return0;}附錄2#include<iostream.h>#include<math.h>#defineeps0.001doublefun(doublex){returnx*x+2*x;}doublediff_fun(doublex){return2*x+2;}doubledichotomy(doublea,doubleb){doublec=0.0;if((diff_fun(a)<0)&&(diff_fun(b)>0)){while(true){c=(a+b)/2;if(diff_fun(c)<0){a=c;if(fabs((a-b))<eps){return(a+b)/2;}}elseif(diff_fun(c)==0){returnc;}else{b=c;if(fabs((a-b))<eps){return(a+b)/2;}}}}}intmain(){doublea=-3.0;doubleb=5.0;doubleresult=dichotomy(a,b);cout<<"求解的結果是:"<<result<<endl;return0;}附錄3#include<iostream.h>#include<math.h>#defineeps0.01doublefun(doublex){returnx*x*x-2*x+1;}doublediff_fun(doublex){return3*x*x-2;}doublediff_2_fun(doublex){return6*x;}doublenewton(doublea,doubleb){doubleresult=(a+b)/2;doublet=0.0;while(true){t=result-(diff_fun(result)/diff_2_fun(result));if(fabs((t-result))<eps){returnresult;}else{result=t;}}}intmain(){doublea=0.0;doubleb=3.0;doublex=newton(a,b);cout<<"求解的結果是x="<<x<<endl;cout<<"此時f(x)="<<fun(x)<<endl;return0;}附錄4#include<iostream>#include<math.h>usingnamespacestd;doublefun(doublet){return(t*t+2*t);}voidGoldensection(double(*pfun)(doublet)){intmaxflag=1000,k=1;doublea=-3,b=5,err=0.001,t,t1,t2;do{t1=a+0.618*(b-a);t2=b-0.618*(b-a);if(t1=t2){a=t2;b=t1;}else{if(t1<t2){a=t2;}else{b=t1;}}k++;}while(fabs(a-b)>err&&k<maxflag);if(k>=maxflag)cout<<endl<<"黃金分割法迭代失敗!迭代次數為k="<<k<<endl;else{t=(a+b)/2;cout<<endl<<"黃金分割法迭代成功!迭代次數為k="<<k-1<<endl;cout<<"迭代結果：近似根為root="<<t<<endl;cout<<"函數值近似為：f(root)="<<fun(t)<<endl;}}intmain(){Goldensection(fun);return0;}附錄5#include<iostream>#include<math.h>usingnamespacestd;doublefun(doublex){return(8*x*x*x-2*x*x-7*x+3);}doublemax(doublea,doubleb){if(a>b)returna;elsereturnb;}doublemin(doublea,doubleb){if(a>b)returnb;elsereturna;}voidParainterpolation(double(*pfun)(doublex)){doublea=0,b=2,err=0.001,x=0,x0=1,f,f0;do{x=((x0*x0-b*b)*fun(a)+(b*b-a*a)*fun(x0)+(a*a-x0*x0)*fun(b))/(2*((x0-b)*fun(a)+(b-a)*fun(x0)+(a-x0)*fun(b)));f0=fun(x0);f=fun(x);if(f=f0){a=min(x,x0);b=max(x,x0);x0=(a+b)/2;}else{if((fun(x)-f0)*(x-x0)>0){b=max(x,x0);x0=min(x,x0);}else{a=min(x,x0);x0=max(x,x0);}}}while(fabs(x-x0)>err);x=(x+x0)/2;cout<<"迭代結果："<<endl;cout<<"近似根:"<<x<<endl;cout<<"函數值近似為:"<<fun(x)<<endl;}intmain(){Parainterpolation(fun);return0;}附錄6#include<iostream.h>#include<math.h>doublelamda(doublex[2],doublep[2],doublea[2]){doublelam1,lam2;lam1=(pow(a[0],3)*x[0]*x[0]+pow(a[1],3)*x[1]*x[1]);lam2=-(pow(a[0]*x[0],2)+pow(a[1]*x[1],2));doubles;s=-lam2/(2*lam1);returns;}voidmain(){cout<<"最速下降法求解最優(yōu)解程序運行結果"<<endl;cout<<endl;doublelamd,x[3],a[6];doublep[2],g[2],e,y,m,n;inti=0;cout<<"請輸入精度e"<<endl;cin>>e;cout<<"請輸入初始點x[0],x[1]的值:\n"<<endl;cin>>m;cin>>n;x[0]=m;x[1]=n;cout<<"函數通式為f(x)=a[0]x1*x1+a[1]x2*x2+a[2]x1*x2+a[3]x1+a[4]x2+a[5]"<<endl;cout<<"請依次輸入函數的系數：a[0]、a[1]、a[2]、a[3]、a[4]、a[5]:"<<endl;for(i=0;i<6;i++)cin>>a[i];p[0]=(2*a[0]*x[0]+a[2]*x[1]+a[3]);p[1]=(2*a[1]*x[1]+a[2]*x[0]+a[4]);g[0]=-p[0];g[1]=-p[1];i=0;cout<<endl;while(sqrt(g[0]*g[0]+g[1]*g[1])>e&&i<=200){lamd=lamda(x,g,a);x[0]=x[0]+lamd*g[0];x[1]=x[1]+lamd*g[1];p[0]=2*a[0]*x[0];p[1]=2*a[1]*x[1];

g[0]=-p[0];g[i]=-p[i];i++;cout<<******************************************<<endl;cout<<"第"<<i<<”次迭代結果："<<endl;cout<<"p的模為："<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<"x的值”<<x[0]<<""<<x[1]<<endl;cout<<******************************************<<endl;cout<<endl;}y=(a[0]*x[0]*x[0]+a[1]*x[0]*x[1]+a[2]*x[0]*x[1]+a[3]*x[0]+a[4]*x[1]+a[5]);cout<<"此時滿足精度要求的p的模為:"<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<endl;cout<<"滿足精度的最優(yōu)近似結果x[1],x[2]分別為cout<<******************************************<<endl;cout<<"第"<<i<<”次迭代結果："<<endl;cout<<"p的模為："<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<"x的值”<<x[0]<<""<<x[1]<<endl;cout<<******************************************<<endl;cout<<endl;}y=(a[0]*x[0]*x[0]+a[1]*x[0]*x[1]+a[2]*x[0]*x[1]+a[3]*x[0]+a[4]*x[1]+a[5]);cout<<"此時滿足精度要求的p的模為:"<<sqrt(g[0]*g[0]+g[1]*g[1])<<endl;cout<<endl;cout<<"滿足精度的最優(yōu)近似結果x[1],x[2]分別為："<<endl;cout<<"x[1]="<<x[0]<<endl;cout<<"x[2]="<<x[1]<<endl;cout<<endl;cout<<"滿足進度要求所得的最優(yōu)值為："<<endl;cout<<"minf(x)="<<y<<endl;}附錄7#include<iostream.h>#include"matrix.h"http://矩陣的有關運算intconstn1=2;doublef(double*x);//目標函數；voidg(double*x,double*y);//目標函數的梯度；voidmain(){inti;doublex0[n1],x1[n1],g0[n1],g1[n1],H[n1][n1],h[n1][n1],p[n1],f0,f1,temp,e;x0[0]=0;x0[1]=0;e=0.01;H[0][0]=2;H[0][1]=-1;H[1][0]=-1;H[1][1]=2;f0=f(x0);g(x0,g0);matrix_inv(H,h);matrix_mul(h,g0,p);//matrix_minus(x0,p,x1);for(i=0;i<n1;i++)x1[i]=x0[i]-p[i];g(x1,g1);f1=f(x1);temp=g1[0]*g1[0]+g1[1]*g1[1];while(temp>e){for(i=0;i<n1;i++)f0=f1;matrix_mul(h,g0,p);for(i=0;i<n1;i++)x1[i]=x0[i]-p[i];g(x1,g1);f1=f(x1);temp=g1[0]*g1[0]+g1[1]*g1[1];}cout<<"當X取X("<<x1[0]<<","<<x1[1]<<")時,Min(f(X)="<<f(x1)<<endl;}doublef(double*x){return60-10*x[0]-4*x[1]+x[0]*x[0]+x[1]*x[1]-x[0]*x[1];}voidg(double*x,double*y){y[0]=2*x[0]-x[1]-10;y[1]=2*x[1]-x[0]-4;}附錄8#include<iostream.h>#include<math.h>#include<iomanip.h>#definen2//正定二次函數的自變量個數doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2])//輸入變量為函數自變量初值{inti,j;doublef=0;for(i=0;i<n;i++)//計算X人2部分{f+=pow(x[i],2)*f_xs[i];}for(;i<2*n;i++)//計算X部分{f+=x[i%n]*f_xs[i];}f+=f_xs[i];//計算常數項部分for(i=0;i<n;i++)//計算XiXj部分{for(j=i+1;j<n;j++){f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];}}returnf;}voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1]){inti,j;for(i=0;i<n;i++){Q[i][i]=2*f_xs[i];}for(i=0;i<n;i++){Q[i][n]=f_xs[n+i];}for(i=0;i<n;i++){for(j=i+1;j<n;j++){Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];}}}voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n])//自變量初值，正定二次函數的各項系數，返回梯度的初值{inti,j;for(i=0;i<n;i++){g0[i]=0;//清零for(j=0;j<n;j++){if(i==j)g0[i]+=Q[i][j]*x[j];elseg0[i]+=Q[i][j]*x[j];}}for(i=0;i<n;i++){g0[i]+=Q[i][n];}cout<<endl<<"梯度g0的值="<<endl;for(i=0;i<n;i++)cout<<setw(10)<<g0[i]<<endl;cout<<endl;}voidHesse_fun(doubleQ[n][n+1],doubleHesse[n][n]){inti,j;for(i=0;i<n;i++){for(j=0;j<n;j++){Hesse[i][j]=Q[i][j];}}cout<<"Hesse矩陣:"<<endl;for(i=0;i<n;i++){for(j=0;j<n;j++){cout<<setw(10)<<Hesse[i][j]<<"";}cout<<endl;}}intH(doubleg0[n],doublec)//判別準則：返回1結束,返回0繼續(xù)迭代{doubles=0;for(inti=0;i<n;i++){s+=pow(g0[i],2);}if(sqrt(s)<c)return1;elsereturn0;}voidinv(doublea[n][n],doublec[n][n])//求Hesse矩陣的逆矩陣{doublem[n];//輔助乘數doubleT;//存儲換行臨時變量doubletemp=0;doublet;//最大列主元inttap1=0,tap2=0;//最大列主元下標for(inti=0;i<n;i++){for(intj=0;j<n;j++){if(i==j){c[i][j]=1;}else{c[i][j]=0;}}}for(ints=0;s<n;s++){t=a[s][s];//賦初值for(intp=s;p<n;p++)//選取最大列主元{if(fabs(a[p][s])>t){t=a[p][s];tap1=p;tap2=s;}}if(t==0){cout<<"列主元為零,失?。?<<endl;break;}if(tap1!=tap2)//換行{for(inti=0;i<n;i++)//s{T=a[tap1][i];a[tap1][i]=a[s][i];a[s][i]=T;T=c[tap1][i];//逆矩陣換行c[tap1][i]=c[s][i];c[s][i]=T;}}tap1=tap2=0;//置零for(intj=0;j<n;j++)//單位化,(注意：逆矩陣每列都單位化,而原矩陣則不受影響){a[s][j]=a[s][j]/t;c[s][j]=c[s][j]/t;}for(inti=0;i<n;i++)//消元{if(i!=s){m[i]=a[i][s]/a[s][s];for(intj=0;j<n;j++)//注意：逆矩陣每列都消元,而原矩陣則不受影響{a[i][j]=a[i][j]-m[i]*a[s][j];c[i][j]=c[i][j]-m[i]*c[s][j];}}}}cout<<endl<<endl;cout<<"Hesse矩陣的逆矩陣為:"<<endl;for(i=0;i<n;i++){for(intj=0;j<n;j++){cout<<setw(10)<<c[i][j]<<"";}cout<<endl;}}voidxiang_cheng(doublea[n][n],doubleb[n],doublec[n])//n*n矩陣乘以n*1矩陣{inti,j,k;if(n==n){cout<<"兩矩陣階數相符可以相乘！"<<endl;for(i=0;i<n;i++){c[i]=0;}for(k=0;k<n;k++)//c矩陣的第k行{for(j=0;j<n;j++){c[k]+=a[k][j]*b[j];}}cout<<"Hesse矩陣的逆矩陣和梯度向量相乘,兩矩陣相乘結果為:"<<endl;for(i=0;i<n;i++){cout<<setw(10)<<c[i]<<endl;}}elsecout<<"兩個矩陣階數不符,不能相乘"<<endl;}〃開始計算minf(Xk+tPk)時的步長t的值，由于這是n元二次函數所以minf(t)是關于t>0的二次函數,先求二次方程a,b,c系數，用一階導為零。voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3])//t[3]中返回的是a,b,c的系數值{inti,j;t[0]=t[1]=t[2]=0;t[0]=fun(p,f_xs)-f_xs[2*n];for(i=n;i<2*n;i++){t[0]-=f_xs[i]*p[i%n];}for(i=0;i<n;i++){t[1]+=2*f_xs[i]*x[i]*p[i];}for(;i<2*n;i++){t[1]+=f_xs[i]*p[i%n];}for(i=0;i<n;i++){for(j=i+1;j<n;j++){t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);}}t[2]=fun(x,f_xs);cout<<endl<<"二次函數minf(Xk)的二次項系數："<<t[0]<<"一次項系數："<<t[1]<<"常數項系數："<<t[2]<<endl<<endl;}〃二次函數一階導為零計算t的值，t>0doublet_c(doublet[3]){cout<<"用一階導求得步長為：t="<<-t[1]/t[0]/2<<endl;return-t[1]/t[0]/2;}voidmain(){doublef_xs[n+n+1+(n-1)*n/2]={4,2,9,-3,16,0};//正定二次函數的各項系數,第一個為：X1人2系數,第二個為：X2人2系數,第三個為：X1系數,第四個為:X2系數,第五個為：常數項,第六個為：X1X2系數;//n元的系數存放類推doublex[n]={0,0};//函數自變量初值doublef0;//函數值doubleg0[n];//梯度的值doubleQ[n][n+1];//正定二次函數對應的實對稱矩陣doublec=0.01;//H準則限值doubleHesse[n][n];//Hesse矩陣doubleinv_Hesse[n][n];//Hesse矩陣的逆doublet[3];//返回求minf()時t的二次函數的a,b,c的系數值doublet_bc;//步長doublep[n];//保存矩陣相乘結果一加負號后變成下降方向inti,tap=0;//迭代次數Q_fun(f_xs,Q);//計算正定二次函數對應的實對稱矩陣ok!f0=fun(x,f_xs);//求函數初值ok!D_fun(x,Q,g0);//返回梯度的初值while(H(g0,c)==0){Hesse_fun(Q,Hesse);//返回Hesse矩陣的值？？一個二次函數的Hesse矩陣是變的還是不變inv(Hesse,inv_Hesse);//返回Hesse矩陣的逆矩陣xiang_cheng(inv_Hesse,g0,p);for(i=0;i<n;i++)p[i]*=(-1);abc(x,p,f_xs,t);//開始計算minf(Xk+tPk)時的步長t的值,t_bc=t_c(t);//求一階導來計算tfor(i=0;i<n;i++){x[i]=x[i]+t_bc*p[i];}tap++;f0=fun(x,f_xs);D_fun(x,Q,g0);}cout<〈"修正newton法"<<endl;cout<<"函數f(x1,x2)=4(X1+1)人2+2(X2-1)人2+X1+X2+10.的極小點為："<<"f="<<f0<<endl;cout<<"自變量取值為："<<endl;for(i=0;i<n;i++){cout<<"x"<<i+1<<"="<<x[i]<<endl;}cout<<"迭代次數為："<<tap<<endl;}附錄9#include<stdio.h>#include<math.h>#defineN10#defineepspow(10,-6)doublef(doublex[],doublep[],doublet)(doubles;s=pow(x[0]+t*p[0],2)+4*pow(x[1]+t*p[1],2);returns;}voidsb(double*a,double*b,doublex[],doublep[]){doublet0,t1,t,h,alpha,f0,f1;intk=0;t0=2.5;/*初始值*/h=1;/*初始步長*/alpha=2;/*加步系數*/f0=f(x,p,t0);t1=t0+h;f1=f(x,p,t1);while(1){if(f1<f0){h=alpha*h;t=t0;t0=t1;f0=f1;k++;}else{if(k==0){h=-h;t=t1;}else{*a=t<t1?t:t1;*b=t>t1?t:t1;break;}}t1=t0+h;f1=f(x,p,t1);}}doublehjfg(doublex[],doublep[]){doublebeta,t1,t2,t;doublef1,f2;doublea=0,b=0;double*c,*d;c=&a,d=&b;sb(c,d,x,p);/*調用進退法搜索區(qū)間*/printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf”,x[0],x[1],p[0],p[1]);printf("\n[a,b]=[%lf,%lf]”,a,b);beta=(sqrt(5)-1.0)/2;t2=a+beta*(b-a);f2=f(x,p,t2);t1=a+b-t2;f1=f(x,p,t1);while(1){if(fabs(t1-t2)<eps)break;else{if(f1<f2){t=(t1+t2)/2;b=t2;t2=t1;f2=f1;t1=a+b-t2;f1=f(x,p,t1);}else{a=t1;t1=t2;f1=f2;t2=a+beta*(b-a);f2=f(x,p,t2);}}}t=(t1+t2)/2;returnt;}voidgtd(){doublex[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0;inti,k,n;printf("請輸入函數的元數值n=");scanf("%d",&n);printf("\n請輸入初始值卅n");for(i=0;i<n;i++)scanf("%lf",&x[i]);f0=f(x,g,t);g[0]=2*x[0];g[1]=50*x[1];mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的長度*/if(mod1>eps){p[0]=-g[0];p[1]=-g[1];k=0;while(1)(t=hjfg(x,p);/*調用黃金分割法求t的值*/printf("\np1=%lf,p2=%lf,t=%lf”,p[0],p[1],t);x[0]=x[0]+t*p[0];x[1]=x[1]+t*p[1];g[0]=2*x[0];g[1]=50*x[1];/*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf”,x[0],x[1],g[0],g[1]);*/mod2=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的長度*/if(mod2<=eps)break;else{if(k+1==n){g[0]=2*x[0];g[1]=50*x[1];p[0]=-g[0];p[1]=-g[1];k=0;}else{nanda=pow(mod2,2)/pow(mod1,2);printf("\nnanda=%lf,mod=%lf",nanda,mod2);p[0]=-g[0]+nanda*p[0];p[1]=-g[1]+nanda*p[1];mod1=mod2;k++;}}printf("\n");}}printf("\n最優(yōu)解為x1=%lf,x2=%lf",x[0],x[1]);printf("\n最終的函數值為%lf",f(x,g,t));}main(){gtd();}附錄10#include<iostream.h>#include<math.h>#include<iomanip.h>intconstn=2;//正定二次函數的自變量個數doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2]);//輸入變量為函數自變量初值voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1]);voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n]);intH(doubleg0[n],doublec);//判別準則：返回1結束,返回0繼續(xù)迭代voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3]);doublet_c(doublet[3]);//二次函數一階導為零計算t的值,t>0voidmain(){doublef_xs[n+n+1+(n-1)*n/2]={4,1,-40,T2,136,0};doublex[n]={8,9};//函數自變量初值doublef0;//函數值doubleg0[n];//梯度的值doubleQ[n][n+1];〃求梯度處值設置的中間變量，包含兩部分：正定二次函數對應的實對稱矩陣，還有一次項系數doublec=0.01;//H準則限值doublet[3];//返回求minf()時t的二次函數的a,b,c的系數值doublet_bc;//步長doublep[n];//保存下降方向doubleH0[n][n];//保存模擬Hesse矩陣的逆doubley[n];//y(k)=g0(k+1)-g0(k)doubles[n];//s(k)=X(k+1)-X(k)doubles_temp[n][n]={0};//計算保存矩陣doubles_temp2[n][n]={0};doubles_temp3[n][n]={0};doubles_tl[n]={0};doubletemp;//臨時值inti,j,k,flag=0,tap=0;//迭代次數Q_fun(f_xs,Q);//計算正定二次函數對應的實對稱矩陣f0=fun(x,f_xs);//求函數初值D_fun(x,Q,g0);//返回梯度的初值do{for(i=0;i<n;i++)//給H0[n][n]的處值賦單位矩陣{for(j=0;j<n;j++){if(i==j)H0[i][j]=1;elseH0[i][j]=0;}}for(i=0;i<n;i++)p[i]=(-1)*g0[i];k=0;//step2;do{abc(x,p,f_xs,t);//開始計算minf(Xk+tPk)時的步長t的值，t_bc=t_c(t);//求一階導來計算tfor(i=0;i<n;i++){x[i]=x[i]+t_bc*p[i];s[i]=t_bc*p[i];//保存計算之值X(k+1)-X(k)}for(i=0;i<n;i++)y[i]=g0[i];//保存之類的f0=fun(x,f_xs);D_fun(x,Q,g0);//step3;for(i=0;i<n;i++)y[i]=g0[i]-y[i];//保存計算g0(k+1)-g0(k)if(H(g0,c)==0)//即不滿足小于c{if(k!=n){//y//shavedone!temp=0;//初值for(i=0;i<n;i++)temp+=s[i]*y[i];for(i=0;i<n;i++){for(j=0;j<n;j++){s_temp[i][j]=s[i]*s[j]/temp;}}〃求出S(k)S(k)t/S(k)tykfor(i=0;i<n;i++)//初值s_tl[i]=0;for(i=0;i<n;i++){for(j=0;j<n;j++){s_tl[i]+=y[j]*H0[j][i];}}temp=0;//初值for(i=0;i<n;i++){temp+=s_tl[i]*y[i];}〃這時s_tl[n]和s_temp2[n][n]可以利用for(i=0;i<n;i++)//初值s_tl[i]=0;for(i=0;i<n;i++){for(j=0;j<n;j++){s_tl[i]+=H0[i][j]*y[j];}}for(i=0;i<n;i++)//初值for(j=0;j<n;j++)s_temp2[i][j]=0;for(i=0;i<n;i++){for(j=0;j<n;j++){s_temp2[i][j]=s_tl[i]*y[j];}}for(i=0;i<n;i++)//初值for(j=0;j<n;j++)s_temp3[i][j]=0;for(i=0;i<n;i++){for(j=0;j<n;j++){for(intii=0;ii<n;ii++){s_temp3[i][j]+=s_temp2[i][ii]*H0[ii][j];}}}for(i=0;i<n;i++){for(j=0;j<n;j++){H0[i][j]=H0[i][j]+s_temp[i][j]-s_temp3[i][j]/temp;}}for(i=0;i<n;i++){temp=0;for(j=0;j<n;j++){temp+=H0[i][j]*g0[j];}p[i]=(-1)*temp;}k++;}else{f0=fun(x,f_xs);break;}}else{flag=1;break;}tap++;}while(H(g0,c)==0);if(flag==1)//符合梯度準則跳出{flag=0;break;}}while(k==n);cout<<"DFP法"<<endl;cout<<"函數f(x1,x2)=4(X1-5)人2+(X2-6)人2.的極小點為："<<"f="<<f0<<endl;cout<<"取得極小點為自變量取值為："<<endl;for(i=0;i<n;i++){cout<<"x"<<i+1<<"="<<x[i]<<endl;}}doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2])//輸入變量為函數自變量初值{inti,j;doublef=0;for(i=0;i<n;i++)//計算X人2部分{f+=pow(x[i],2)*f_xs[i];}for(;i<2*n;i++)//計算X部分{f+=x[i%n]*f_xs[i];}f+=f_xs[i];〃計算常數項部分for(i=0;i<n;i++)//計算XiXj部分{for(j=i+1;j<n;j++){f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];}}returnf;}voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1])//Q[n][n]是二次函數的正定矩陣，但Q的第n+1列存儲一次項的系數{inti,j;for(i=0;i<n;i++){Q[i][i]=2*f_xs[i];}for(i=0;i<n;i++){Q[i][n]=f_xs[n+i];}for(i=0;i<n;i++){for(j=i+1;j<n;j++){Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];}}}voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n])//自變量初值，正定二次函數的各項系數，返回梯度的初值{inti,j;for(i=0;i<n;i++){g0[i]=0;//清零for(j=0;j<n;j++){if(i==j)g0[i]+=Q[i][j]*x[j]；elseg0[i]+=Q[i][j]*x[j]；}}for(i=0;i<n;i++){g0[i]+=Q[i][n];}cout<<endl<<"梯度g0的值="<<endl;for(i=0;i<n;i++)cout<<setw(10)<<g0[i]<<endl;cout<<endl;}intH(doubleg0[n],doublec)//判別準則：返回1結束,返回0繼續(xù)迭代{doubles=0;for(inti=0;i<n;i++){s+=pow(g0[i],2);}if(sqrt(s)<c)return1;elsereturn0;}voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3])//t[3]中返回的是a,b,c的系數值{inti,j;t[0]=t[1]=t[2]=0;t[0]=fun(p,f_xs)-f_xs[2*n];for(i=n;i<2*n;i++){t[0]-=f_xs[i]*p[i%n];}for(i=0;i<n;i++){t[1]+=2*f_xs[i]*x[i]*p[i];}for(;i<2*n;i++){t[1]+=f_xs[i]*p[i%n];}for(i=0;i<n;i++){for(j=i+1;j<n;j++){t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[i]);}}t[2]=fun(x,f_xs);cout<<endl<<"二次函數minf(Xk)的二次項系數："<<t[0]<<"一次項系數："<<t[1]<<"常數項系數："<<t[2]<<endl<<endl;}doublet_c(doublet[3]){cout<<"用一階導求得步長為：t="<<-t[1]/t[0]/2<<endl;return-t[1]/t[0]/2;}附錄11#include<iostream.h>#include<math.h>intconstn=2;staticdoubleMk=1;doublefun_original(doublex[n]);//既是主函數,也是F(X,Mk)之一doubleg_fun(doublex[n]);//約束函數doubleh_fun(doublex[n]);;//等式函數doublea_e_f(doublex[n]);doublefun(doublex[n]);//讓單純形法求解此函數最小值voidcompare(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3]);//在調用之前LGH先給x[n+1]賦值voidzx(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3],doublexn_1[n],doublexn_2[n]);doubleH(doublex[n+1][n],doubleLGH[3][n],doublec);//判別準則voiddanchun(doubleoriginal_point[n]);voidmain(){intk=1;inti;doubler=10;doublec=0.000001;doublex[n]={0,0.5};//給定的初始點在某個增廣目標函數的求解域之中do{danchun(x);//ln(x2)>=0時F(X,Mk)=-X1+X2;無極小點,不用管這個增廣目標函數if(Mk*a_e_f(x)<c)break;Mk=r*Mk;}while(1);cout<<"用外點法求解："<<endl;cout<<"minf(X)=-X1+X2"<<endl;cout<<"g1(X)=lnX2>=0"<<endl<<"s.t."<<endl;;cout<<"g2(X)=X1+X2-1=0"<<endl;cout<<"結果是："<<endl;cout<<"x*="<<endl;for(i=0;i<n;i++)cout<<x[i]<<"";cout<<endl;cout<<"f(x*)="<<fun_original(x)<<endl;}doublefun_original(doublex[n])〃既是主函數,也是F(X,Mk)之一{doubley;y=-x[0]+x[1];returny;}doubleg_fun(doublex[n])//約束函數{doubley;y=pow(log(x[1]),2);returny;}doubleh_fun(doublex[n])//等式函數{doubley;y=x[0]+x[1]-1;returny;}doublea_e_f(doublex[n]){doubley;y=Mk*(pow(h_fun(x),2)+pow(g_fun(x),2));returny;}doublefun(doublex[n])//讓單純形法求解此函數最小值{doubley;y=fun_original(x)+a_e_f(x);returny;}voidcompare(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3])//在調用之前LGH先給x[n+1]賦值{inti;doubley[n+1];for(i=0;i<n+1;i++)y[i]=fun(x[i]);//y[i]=fun(x[i][n]);//////////////////////////////////////////////////?doubletemp[2];temp[0]=y[0];for(i=0;i<n+1;i++)//temp[0]存儲最小的值XL,temp[1]存儲對應x的下標{if(y[i]<=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[0]=temp[1];//記錄標記for(i=0;i<n;i++)LGH[0][i]=x[int(temp[1])][i];for(i=0;i<n+1;i++)//temp[0]存儲最大的值XH,temp[1]存儲對應x的下標{if(y[i]>=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[2]=temp[1];//for(i=0;i<n;i++)LGH[2][i]=x[int(temp[1])][i];temp[0]=y[0];if(LGH_tap[2]==0){temp[0]=y[1];}for(i=0;i<n+1;i++)//尋找次大值{if(i==LGH_tap[2])//if(i==temp[1])continue;if(y[i]>=temp[0]){temp[0]=y[i];temp[1]=i;}}LGH_tap[1]=temp[1];for(i=0;i<n;i++)LGH[1][i]=x[int(temp[1])][i];}voidzx(doublex[n+1][n],doubleLGH[3][n],doubleLGH_tap[3],doublexn_1[n],doublexn_2[n]){inti,j;for(i=0;i<n;i++)xn_1[i]=xn_2[i]=0;for(i=0;i<n+1;i++)//返回重心{if(i==LGH_tap[2])//輪流加n次(去掉一次)continue;for(j=0;j<n;j++){xn_1[j]+=x[i][j];}}for(i=0;i<n;i++)xn_1[i]/=n;for(i=0;i<n;i++)//返回反射點xn_2[i]=2*xn_1[i]-LGH[2][i];}doubleH(doublex[n+1][n],doubleLGH[3][n],doublec)//判別準則{inti;doubles=0;for(i=0;i<n+1;i++)//小心i從0到n即(i=0;i<n+1;i++){s+=pow(fun(x[i])-fun(LGH[0]),2);//s+=pow(fun(x[i][n])-fun(LGH[0][n]),2);//////////////}if(s<=c)return1;elsereturn0;}voiddanchun(doubleoriginal_point[n]){doublex[n+1][n]={0};doublec=0.000001;doublea=2;doubleb=0.5;doubleLGH[3][n]={0};doubleLGH_tap[3]={0};doublexn_1[n]={0};doublexn_2[n]={0};doublexn_3[n]={0};doublexn_4[n]={0};inti,j,tap=0;for(i=0;i<n;i++)x[0][i]=original_point[i];for(i=1;i<n+1;i++)//假設只給出一個初始點X0(1,1),求另外兩個初始點{for(j=0;j<n;j++){if(j+1==i)x[i][j]=x[0][j]+1;elsex[i][j]=x[0][j];}}compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它們在x[n+1][n]中的下標zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射點X(n+2)while(H(x,LGH,c)==0){if(fun(xn_2)<fun(LGH[0]))//if(fun(xn_2)<fun(LGH[0][n]))//////////////////{for(i=0;i<n;i++){xn_3[i]=xn_1[i]+a*(xn_2[i]-xn_1[i]);//擴張}if(fun(xn_3)<fun(LGH[0])){for(i=0;i<n;i++){x[int(LGH_tap[2])][i]=LGH[2][i]=xn_3[i];}}else{for(i=0;i<n;i++){x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}}}else{if(fun(xn_2)<fun(LGH[1])){for(i=0;i<n;i++){x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}//gotostep2;}else{if(fun(xn_2)<fun(LGH[2]))//elseif(fun(xn_2)<fun(LGH[2][n])){for(i=0;i<n;i++){x[int(LGH_tap[2])][i]=LGH[2][i]=xn_2[i];}}else//means(fun(xn_2)>=fun(LGH[2][n])){for(i=0;i<n;i++){xn_4[i]=xn_1[i]+b*(xn_2[i]-xn_1[i]);//壓縮}if(fun(xn_4)<fun(LGH[2])){for(i=0;i<n;i++){x[int(LGH_tap[2])][i]=LGH[2][i]=xn_4[i];}//gotostep2;}else{for(i=0;i<n+1;i++){if(i==LGH_tap[0])continue;for(j=0;j<n;j++){x[i][j]=(LGH[0][j]+x[i][j])/2;}}}//gotostep2}}}//That'sit!!!compare(x,LGH,LGH_tap);//返回Xl,Xh,Xg及它們在x[n+1][n]中的下標zx(x,LGH,LGH_tap,xn_1,xn_2);//返回中心X(n+1)和反射點X(n+2)tap++;}for(i=0;i<n;i++)original_point[i]=LGH[0][i];//把在某個Mk值下增廣目標函數的最優(yōu)值返回}附錄12#include<iostream.h>#include<math.h>#include<iomanip.h>constintn=2;staticdoubleRk=10;//Rk是懲罰因子staticintp=1;doublefun_original(doublex[n]);//既是主函數,也是F(X,Rk)之一doublea_e_f(doublex[n]);doublefun(doublex[n]);//讓單純形法求解此函數最小值voidg_fun(doublex[n],doubleg0[n]);//牛頓法子函數開始,求出與Rk對應的梯度向量.voidHesse_fun(doublex[n],doubleHesse[n][n]);//求出與Rk對應的Hes