最优化方法课程实验报告

项目一一维搜索算法（一）［实验目的］编写加步探索法、对分法、Newton法的程序。［实验准备］掌握一维收搜索中搜索区间的加步探索法的思想及迭代步骤;掌握对分法的思想及迭代步骤；掌握Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：用加步探索法确定一维最优化问题min中（t）=13—2t+1t>0的搜索区间，要求选取‘。=°，h0=1，口=2.加步探索法算法的计算步骤:t£［。，+8）（或t£U［0，t］）「+曾甲=甲（t）莎中初力4,中k*h>0（1）选取仞始点0 0 max，计异0 0・给口出仞如步长h。>0，加步系数以>L令上=0。⑵比较目标函数值.令tk+1=tk+«，计算％+1=甲（七+1），若甲*〈气，转（3）,否贝膊叙4）。 * * * * * */c、―■［•穴rh\y人h=ahe人t=t，t=t，k=k+1七卜/o\⑶加大探索步长.令k+1k，同时，令kk k+1 ，转（2）。⑷反向探索.若k=0，转换探索方向，令hk=~hk，t=tk+1，转（2）。否则，停止迭代，a=min{a=min{t,t},b=max{t,t}令程序清单加步探索法算法程序见附录1实验结果运行结果为：rI'D:\ProgramFiles［k86)\IV1icrosoftVisuaIStcidio\MyPr；该间题的根的腰素空间是■:[0^1Pi*essanyl<ei^tocontinue用对分法求解min中(t)=t(t+2),已知初始单谷区间［a，b］=［-3,5］，要求按精度8=0.3,8=0.001分别计算.对分法迭代的计算步骤：确定初始搜索区间［a,b］,要求中'(a)<0，中S>0。， 1，八计算［a,b］的中点c=-(a+b).若甲'(c)<0,则a=c，转(4)；若甲'(c)=0，则t*=c，转(5)；若甲'(c)>0，则b=c，转(4).1若Ia-b1<8，则t*=f(a+b)，转(5)；否则转(2).打印t*，结束对分法的计算框图对分法程序见附录2实验结果运行结果为：I'D:\PnogramFiles[x86AMicrcEoftVi&ualStudic\MyProject世解的结果是:T _Pressanyk=ytocontinue用Newton法求解min中(t)=13—2t+1,已知初始单谷区间[M]=[0,1]，要求精度e=0.01.Newton法的计算步骤⑴确定初始搜索区间[a,b],要求平'(a)<0，甲'(b)>0选定10计算t=t0一甲'(t0)/甲"(t0)若11—101>£，则10=t，转(3)；否则转(5).打印'，中(t)，结束.Newton法的计算框图

程序清单Newton法程序见附录3实验结果运行结果为：项目二一维搜索算法(二)［实验目的］编写黄金分割法、抛物线插值法的程序。［实验准备］掌握黄金分割法的思想及迭代步骤；掌握抛物线插值法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：用黄金分割法求解min中(t)=t(t+2)已知初始单谷区间［a，b］=［—3,5］，要求精度8=0.001黄金分割法迭代步骤:确定中(t)的初始搜索区间[a,b].计算12=a+0.382(b-a)计算'=a+0.618(b一a)若11-11<£，则打印t*=七;，结束；否则转(5).TOC\o"1-5"\h\z1 2 2判别是否满足中<^:若满足，则置a=t,t=t,中=中1 2 22 1 2 1然后转(3)；否则，置b=t,t=t，中=中,t=以+。(b—a),中=中(t)11 2 1 22 2 2然后转(4).黄金分割法的计算框图：程序清单黄金分割法程序见附录4实验结果运行结果为：''D:\ProgramFiles(xBGJA.MicrosoftVrsualStudi-o\MyPrqjecTts^zuiyou\_Debug\ziyuan.exe,功丁迭代欢数为蚌］迭代皓巢为root=0-B56函数值近似为:fCroot^=0.115136Pressa_nykeytocontinue.用抛物线插值法求解minf(x)=8x3—2x2—7x+3已知初始单谷区间M=［°，W=0-001.抛物线插值法的计算步骤：中(t)v中(t)，所以相对t来说t是好点，故划掉区间［t,t］，保留［t,t］为新区间，0 0 0 2 1 0故置t=t，中(t)=中(t),t=t，甲(t)=甲(t),t保持不变；2 0 2 0 0 0 1甲(t)＞平财,所以相对t来说10是好点，故划掉区间［ti，t］,保留［t,t2］为新区间，故置匕=t,中(t)=中(t「，t°与t2保持不变；程序清单抛物线插值法程序见附录5实验结果运行结果为:项目三常用无约束最优化方法(一)［实验目的］编写最速下降法、Newton法（修正Newton法）的程序。［实验准备］掌握最速下降法的思想及迭代步骤。掌握Newton法的思想及迭代步骤；掌握修正Newton法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：用最速下降法求minf（X）=x2+25x;,X0=［2,2］tt=0.01最速下降法计算步骤（1） 取初始点X（0）,容许误差（精度）8>0,令一=0（2） 计算p（k）=-Vf（X（k））（3） 检验|p（』|<£?若是迭代终止,取X*=X（k）,否则转4。（4） 求最优步长人：minf（X（k）+Xp（k））=f（X（k）+1p（k））（一维搜索）令X（k+i）=X（k）+1p（k）,令k:=k+1,转20k最速下降法的计算框图程序清单最速下降法程序见附录6实验结果运行结果为：用Newton法求TOC\o"1-5"\h\zminf(X)二60一10x—4x+x2+%2一xx

1 2 1 2 12'初始点X0=[0，0]t，£=0.01.Newton法的计算步骤给定初始点x(0)，及精度e>0，令k=0；Vf(X(k)) … …若 ，停止，极小点为x(k)，否则转步骤(3)；、■曾[V2f(X(k))]TaS(k)=一[H(X(k))]TVf(X(k))计算 ，令 ;令x(k+D=x(k)+s(k)，k=k+1，转步骤⑵。程序清单Newton法程序见附录7实验结果运行结果为:rrcss-ankeytoconiinLie用修正Newton求minf(X)=4(x+1)2+2(x-1)2+x+x+101 2 1 2 ,初始点X0=[0，0]r，8=0.01.修正Newton的计算步骤给定初始点x(0)，及精度8>0，令k=0；若'W(X以))"'8，停止，极小点为x(k)，否则转步骤(3)；、4曾[?2f(X(k))]-1人5(k)=-[H(X(k))]-1Vf(X(k))计算 ，令 ;(4)用一维搜索法求侦，使得f(X(k)+a(k)S(k))=minf(X(k)+aS(k))，令O>0X(k+1)=X(k)+a(k)s(k)，k=k+1，转步骤(2)。程序清单修正Newton程序见附录8实验结果运行结果为:T"ID^ProgramFiles-x3oro=oftVisuadStud<3\MyPrpjects\zuiyciU\DebLg^1ziyijan.eMe'h—m阵的逆矩阵为=0,125 H0__0.25程嘴矗锚鑫向rra乘，两矩腓乘结殿：1.125-0.75二:J数ninE咨心的二次顶系数;6.Itt7b—次顶系数;-12.375常被攻秀'挫；16芋引了，I犒度"的士=■旨三nHwmu注函数Etxl.展>70(」*1>*2+2<»21r2*Kl的极小点为,F-9,9125B箜星取值为;xl—L.125以一L鄙芯:坦「1P«eckeytocontintie项目四常用无约束最优化方法（二）实验目的编写共轭梯度法、变尺度法（DFP法和BFGS法）程序。实验准备掌握共轭方向法的思路及迭代过程；掌握共轭梯度法的思想及迭代步骤；掌握DFP法和BFGS法的思想及迭代步骤。实验内容及步骤编程解决以下问题：1.用共轭梯度法求得min(x1+4xP，取初始点％=[1，中，£=0.01.共辆梯度法算法步骤给定初始点X(o)，及精度£>0；若时(x(0))1<8，停止，极小值点为x(0)，否则转步骤(3)；取p(0)=—Nf(x(0))，且置k=0；用一维搜索法求七，使得f(x(k)+t^p(k))=minf(x(k)+tp(k))，令，x(k+1)=x(k)+tp(k)，转步骤5；

k若|(x(k+1))1<8，停止，极小值点为x(k+1)，否则转步骤(6)；若k+1=n，令x(0)=x(n)，转步骤(3)，否则转步骤(7)；令p(k+1)=-Vf(x(k+1))+人p(k)，人= ，置k=k+1，转步骤(4)。kk ||Vf(x(k))||2程序清单共轭梯度法程序见附录9实验结果运行结果为：V"D:\Pragr3mFiles(x86)\MicrosoftVisualStudio\MyProjects\zuiyou\Debug\ziyuan.exep阿输入函数的元数值n=2请输入初始值：11El=-0_66GG02#x2=-0-000390,.pl=3_SQQQS4ll-p2=0-000008a.h]=1-4.E00000.1.590339]1=G_G66664.p2=0.083003,fc=0_494238m±=-G.GGGGBfl^x2=0.809933_pl=0_ =-0.000005[a.hJ=[-4.S00000.1.500339]pl=G,6SS66Q^p2=-0.000005^t=0_32SQ1Q最加解去冲=-0-008000『x2=-0-0BBBBB晟弩的函数值为0-0O0000Pi'esEJanytocontinue.2.用共轭梯度法求minf(X)=2%2+号一，自定初始点，£=O.°l程序清单共轭梯度法程序见附录9实验结果运行结果为:si=0_aesseo,^2=a.eseeoB,pi=-e-eaiaBsa,p2=-&.Baaaii[a,508003,1-560000]pl--0.0B0000,p2--0.060011Pt=fl„020000耄佑解为型-0-丽羽如，盗-0・09BS00最终由函数值为0.羽配如Pessens蝠9tocentinme.用DFP法求minf(X)=4(\一5)2+(七-6)2，初始点X°=[8,9]e=0.01.DFP法的具体迭代步骤如下：给定初始点，：'”，迭代精度'，维数n置0—k，单位矩阵I一顼；，计算可"-..。计算搜索方向a进行一维搜索求"，使.赤•'••.•一®叫一*.也山F'ij得迭代新点&EWtF（5） 检验是否满足迭代终止条件「MFE"wE若满足，停止迭代，输出最优解为T星，；否则进行下一步。（6） 检查迭代次数，若k=n，则置”7:"，转向步骤（2）；若k<n，则进行下一步。（7） 计算：』*「，寸;]■卢上笊「二心、心、Rr：l-h7a','3"K:，qjt然后，置k+1-k，转向步骤（3）。DFP法的计算框图程序清单DFP法程序见附录10实验结果

运行结果为:(xl,x2>=4(Xl-5>"2t<X2-6>"2.的极』点为:f=0弟.■得板小度为日变量那伯.为，xl=5x2=6注en"keutnentinli巳一项目五常用约束最优化方法［实验目的］编写外点罚函数法、外点罚函数法的程序。［实验准备］掌握外点罚函数法的思想及迭代步骤；掌握内点罚函数法的思想及迭代步骤。［实验内容及步骤］编程解决以下问题：用外点罚函数法编程计算minf(X)=-尤+尤,Jg(X)=ln尤>0,

..］々(X)=气+%-1=0，精度e=10-5.外点法的计算步骤：给定初始点X(o),初始罚因子人⑴，放大系数c>1；允许误差e>0，设置k=1；以x(k-1)作为搜索初始点，求解无约束规划问题min f(x)+XP(x)，令x(k)为所求极小点。(3)当人P(x(k))<e,则停止计算，得到点x(k)；否则，令九k+D=ck(k),返回(2)执行。程序清单外点罚函数法程序见附录11实验结果实验结果为

运行结果为：家'、•：】=":上一：二Wd二一:.二一上：:一阳外点法求解:minf<X>=-X1+X2yl<K>-lnK2>-0s.t.g2<X>=X1+X2-1=0结果是：x*=7.0567e-0171f<x*>=lPressanykeytocontinue用内点罚函数法编程计算min1(x+1)3+x3 1 2x-1>0,s.t.<1x>0.初始点取为X0=[3,4]T，初始障碍因子取u1=10，缩小系数取C=0-1内点罚步骤:(1)给定初始内点x(0)eS,允许误差e>0,障碍参数Y⑴，缩小系数bG(0,1)，置k=1；(2)以(2)以x(k-1)为初始点，求解下列规划问题:minf(minf(x)+y(k)B(x)s.t.xeS，令x(k)为所求极小点(3)如果y(k)(3)如果y(k)B(x(k))<e则停止计算，得到结果x(k)，(4)否则令y(k+i)=by(k)置k=k+1,返回(2)。内点罚计算框图程序清单内点罚函数法程序见附录12实验结果运行结果为:■"D:\PrcigramFiles&B6}\MicrosoftVisuialSidi^oukDebug\2iyuan-e>e'用内具法求斛“inFCH'TKIH颛ql<H>-Xli>-9附录1#include<>#include<>#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=;doubleb=;StepAdding_Search(a,b);cout<<"该问题的根的搜索空间是：["<<a<<'',''<<b<<'']\n”;return0;}附录2#include<>#include<>#defineepsdoublefun(doublex){returnx*x+2*x;}doublediff_fun(doublex){return2*x+2;}doubledichotomy(doublea,doubleb){doublec=;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=;doubleb=;doubleresult=dichotomy(a,b);cout<<"求解的结果是:"<<result<<endl;return0;}附录3#include<>#include<>#defineepsdoublefun(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=;while(true){t=result-(diff_fun(result)/diff_2_fun(result));if(fabs((t-result))<eps){returnresult;}else{result=t;}}}intmain(){doublea=;doubleb=;doublex=newton(a,b);cout<<"求解的结果是x="<<x<<endl;cout<<"此时f(x)="<<fun(x)<<endl;return0;}附录4#include<iostream>#include<>usingnamespacestd;doublefun(doublet){return(t*t+2*t);}voidGoldensection(double(*pfun)(doublet)){intmaxflag=1000,k=1;doublea=-3,b=5,err=,t,t1,t2;do{t1=a+*(b-a);t2=*(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="<<kT<<endl;cout<<"迭代结果：近似根为root="<<t<<endl;cout<<"函数值近似为：f(root)="<<fun(t)<<endl;}}intmain(){Goldensection(fun);return0;}附录5#include<iostream>#include<>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=,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<>#include<>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<<"最速下降法求解最优解程序运行结果"<<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[1]=-p[1];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<<"满足精度的最优近似结果x[1],x[2]分别为："<<endl;cout<<"x[1]="<<x[0]<<endl;cout<<"x[2]="<<x[1]<<endl;cout<<endl;cout<<"满足进度要求所得的最优值为："<<endl;cout<<"minf(x)="<<y<<endl;}附录7#include<>#include""极小点为："<<"f="<<f0<<endl;cout<<"自变量取值为："<<endl;for(i=0;i<n;i++){cout<<"x"<<i+1<<"="<<x[i]<<endl;}cout<<"迭代次数为："<<tap<<endl;}附录9#include<>#include<>#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=;/*初始值*/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<t1t:t1;*b=t>t1t: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);/*调用进退法搜索区间*/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)/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最优解为x1=%lf,x2=%lf”,x[0],x[1]);printf("\n最终的函数值为％lf”,f(x,g,t));}main(){gtd();}附录10#include<>#include<>#include<>intconstn=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])<<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])voidHesse_fun(doublex[n],doubleHesse[n][n]);<<endl;;cout<<"g2(X)=X2>=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;{g0[0]=pow(x[0]+1,2)-Rk/pow(x[1]-1,2);g0[1]=1-Rk/pow(x[1],2);}voidHesse_fun(doublex[n],doubleHesse[n][n])//求出与Rk对应的Hesse矩阵Hesse[0][0]=2*(x[0]+1)+2*Rk/pow(x[1]-1,3);Hesse[0][1]=Hesse[1][0]=0;Hesse[1][1]=2*Rk/pow(x[1],3);}intH(doubleg0[n],doublec){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