实验报告阻尼牛顿法

太原理工高校机械学院机测系课程上机试验报告课程名称：机械优化设计班级日期成果评定姓名试验室老师签名试验名称用阻尼牛顿法程序解题所用软件C++DEV实验目的及内容试验目的：L把握并能够建立最优化基本类型问题的数学模型。.把握最优化方法的基本概念、基本理论和基本方法，奠定最优化的理论基础。.能够娴熟编制和调试最优化方法的程序，奠定解决实际中的优化问题的基础试验内容：理解阻尼牛顿法并编写相关程序求其最优解。牛顿法程序考核题于(X)=4(玉+1)2+2(x2-1)2+.+(2+10初始点：X°=[0,0]T,梯度精度£=0.01ebs=sqrt(ebs)/(N*N*100);for(i=0;i<N;i++)b[i][i]=1.0;for(i=0;i<N;i++)(if(fabs(a[i][i])<ebs)(for(k=i;k<N;k++)(if(fabs(a[i][i])>ebs)(for(m=0;m<N;m++)(fact=a[i][m];a[i][m]=a[k][m];a[k][m]=fact;fact=b[i][m];b[i][m]=b[k][m];b[k][m]=fact;)break;}})fact=a[i][i];for(k=0;k<N;k++)(a[i][k]/=fact;b[i][k]/=fact;)for(k=0;k<i;k++)(fact=a[k][i];for(m=0;m<N;m++)(a[i][m]-=a[i]fm]*fact;b[i][m]-=b[i][m]*fact;})for(k=i+l；k<N;k++)(fact=a[k][i];for(m=0;m<N;m++)(a[i][m]-=a[i][m]*fact;return;)/*阻尼牛顿法*/voidnewtow_method(doublex[],doubleh,doubleebsin,intyw)(doubled[N],t,f,s[N],hei[N][N],hei_l[N][N];inti,j,k=l;fprintf(fp,”**********阻尼牛顿法计算结果**********\n\n");ywddf(yw);/*输出一维迭代方法*/fprintf(通”初始坐标:\n");f=xfout(x,0);fprintf(fp,n\n");gread(x,d);do(heisen(x,hei);atoa_l(hei,hei_l);for(i=0;i<N;i++)(s[i]=0.0;for(j=0;j<N;j++)s[i]+=hei_l「皿*d「];)switch(yw)(goldcut(x,s,h,ebsin);break;j=rccz(x,s,h,ebsin);break;)gread(x,d)；fprintf(fp,”迭代轮数k=%3d\nM,k);f=xfout(x,k);t=0;for(i=0;i<N;i++)t+=d[i]*d[i];t=sqrt(t);fprintf(fp,"k++;}while(t>ebsin);fprintf(fp"\n*********************\n")•fprintf（fp,”阻尼牛顿法法优化最优点及目标函数值为：\n）f=xfout(x,-l);fprintfCfp/迭代精度：”)；fprintf(fp,"％15.9mn”,t);return;main()(doublexO[N],h,ebsin;intyw;csd_x(xO);/*初始坐标*/h=H_QJ;ebsin=EPSIN;yw=Y_F;fp=fopen(outname,"wn);newtow_method(xO,h,ebsin,yw);/*调用阻尼牛顿法*/fclose(fp);return0;备注：不交此报告者，本次试验为“不合格”。试验原理:实验原理步骤试验步骤：1,画流程图，编写程序；2,将目标函数代入；3,编译运行，将结果保存*********阻尼牛顿法计算结果**********++++++一维搜寻方法：黄金分割法+++++Tf(0)=X(0)=[16.00000000.0000000,迭代轮数k=1x(1)=[-1.1249988,f(1)二9.8125000迭代精度：0.000009986迭代轮数k=2x(2)=[-1.1250000,f(2)二9.8125000迭代精度：0.000000114实*1*7，*1**1**1**1**1*rTwrTweiw•卜•卜•卜.卜.卜"卜4、rTwrTw*Tw一卜验阻尼牛顿法法优化最优点及目标函数值为：x(*)=[-1.1250000,结f(*)=9.8125000迭代精度：0.000000114果初始坐标:0.0000000],0.7499992],0.7500000],0.7500000],算法程序实现/*csssqj.cpp*/#include<string.h>#include<stdio.h>#include<math.h>#include<stdlib.h>#include<time.h>#defineN2#defineEPSIN0.000001#defineH_QJ1.0#defineY_F1/*#defineN2#defineEPSIN0.000001#defineH_QJ1.0#defineY_F1/*初始区间搜寻步长*//*一维搜寻方法选择：1——黄金分割法*//*2——二次插值法*/FILE*fp;charoutname[50]="阻尼牛顿法计算结果/*计算结果输出文件*//东给出初始点坐标*/voidcsd_x(doublex0[])(inti;for(i=0;i<N;i++)/*初始点为坐标原点的状况*/x0[i]=0.0;return;)/*目标函数号doublehanshu(doublex[|)(doublef;/*/*阻尼/*阻尼f=4.0*(x[0]+1.0)*(x[0]+1.0)+2,0*(x[l]-1.0)*(x[l]-1.0)+x[0]+x[l]+10.0;/*阻尼returnf;)/*函数的梯度*/voidgread(doublex[],doublegff])(/*阻尼牛顿法*/gf[0]=8.0*(x[0]+1.0)+1.0;gf[l]=4.0*(x[l]-1.0)+1.0;return;)/*Heisen矩阵*/voidheisen(doublex[],doubleheiflfN])(hei[0][0]=8.0;hei[OJ[l]=O.O;hei[l|[0]=0.0;hei[l][l]=4.0;return;)/*计算f(xk+as)*/doublexkadd(doublex[],doubled[],doublea)(inti;doublexl|N];for(i=0;i〈N;i++)xl[i]=x[i]4-a*d[i];returnhanshu(xl);)/*输出选定的一维迭代方法*/voidywddf(intyw)(switch(yw)(黄金分割法++++++\n\n");二次插值法++++++\n\n”);fprintf(fp,"++++++黄金分割法++++++\n\n");二次插值法++++++\n\n”);fprintf(fp,”++++++一维搜寻方法:break;)return;)/*输出当前迭代点坐标及目标函数值*/doublexfout(doublex[],intm)(intj;doublef;f=hanshu(x);if(-l=m)fprintf(fp；'x(*)=「)；elsefprintf(fp,"x(%3d)=[n,m);for(j=0;j<N-l;j++)fprintf(fp,"%15.71f,n,x[jj);if(-l==m)fprintf(fp；'%15.71f],f(*)=%15.71f\nH,x[N-l],f);elsefprintf(fp,"%15.71f],f(%3d)=%15.71f^nH,x[N-l],m,0;returnf;/*初始搜寻区间的确定*/voidcsssqj(doublex[],doubled[],doubleh,doubleab[])(doubleal,a2,a3,fl,f2,f3;a2=0.0;a3=a2+h;f2=xkadd(x,d,a2);f3=xkadd(x,d,a3);if(f3>f2)(a2=a3;a3=0.0;fl=f2;f2=f3;f3=fl;h=・h;)do(al=a2;a2=a3;fl=f2;f2=f3;a3=a2+h;f3=xkadd(x,d,a3);h=2*h;}while(f3<f2);if(h>0.0)ab[O]=al;ab[l]=a3;)else(ab[0]=a3;ab|l]=al;)return;)/*黄金分割法*/voidgoldcut(doublex[],doubledf],doubleh,doubleebsin)(doublea1,a2,f1,f2,a,b,ab[2];inti;csssqj(x,d,h,ab);a=ab[0];b=ab[l];al=b-0.618*(b-a);fl=xkadd(x,d,al);a2=a+0.618*(b-a);f2=xkadd(x,d,a2);do(if(fl>f2)(a=al;al=a2;fl=f2;a2=a+0.618*(b-a);f2=xkadd(x,d,a2);)else(b=a2;a2=al;f2=fl;al=b-0.618*(b-a);fl=xkadd(x,d,al);}while(b-a>ebsin);for(i=0;i<N;i++)x[i]+=(b+a)*d[i]/2;return;)/*二次插值法*/intrccz(doublex[],doubled[],doubleh,doubleebsin)(doublea1,a2,a3,a4,fl,f2,f3,f4,cl,c2,ab[2];inti,p=0,k=0;csssqj(x,d,h,ab);al=ab[0];a3=ab[l];a2=(al+a3)/2;fl=xkadd(x,d,al);f2=xkadd(x,d,a2);f3=xkadd(x,d,a3);while(l)(cl=(f3-fl)/(a3-al);c2=((f2-fl)/(a2-al)-cl)/(a2-a3);if(0==c2)(P=l；break;}a4=0.5*(a1+a3-c1/c2);if((a4-al)*(a3-a4)<=0.0)(P=2；break;)f4=xkadd(x,d,a4);if(l==k&&fabs(a4-a2)<=ebsin)break;if(a4<a2)(if(f2>f4)(fl=f2;al=a2;k=l;a2=a4;f2=f4;)else(a3=a4;f3=f4;))else(if(f2>f4)(f3=f2;a3=a