实验报告阻尼牛顿法Word下载.docx
《实验报告阻尼牛顿法Word下载.docx》由会员分享，可在线阅读，更多相关《实验报告阻尼牛顿法Word下载.docx（15页珍藏版）》请在冰豆网上搜索。
骤
、
实验原理：
实验步骤：
1，画流程图，编写程序；
2，将目标函数代入；
3，编译运行，将结果保存
结
果
分
析
*********阻尼牛顿法计算结果**********
++++++一维搜索方法：
黄金分割法++++++
初始坐标：
x（0）=[0.0000000,0.0000000],f（0）=16.0000000
迭代轮数k=1
x
（1）=[-1.1249988,0.7499992],f
（1）=9.8125000
迭代精度：
0.000009986
迭代轮数k=2
x
（2）=[-1.1250000,0.7500000],f
（2）=9.8125000
0.000000114
*********************
阻尼牛顿法法优化最优点及目标函数值为：
x（*）=[-1.1250000,0.7500000],f（*）=9.8125000
算法程序实现
/*csssqj.cpp*/
#include<
string.h>
stdio.h>
math.h>
stdlib.h>
time.h>
#defineN2/*优化设计维数*/
#defineEPSIN0.000001/*迭代精度*/
#defineH_QJ1.0/*初始区间搜索步长*/
#defineY_F1/*一维搜索方法选择：
1——黄金分割法*/
/*2——二次插值法*/
FILE*fp;
charoutname[50]="
阻尼牛顿法计算结果.txt"
;
/*计算结果输出文件*/
/*给出初始点坐标*/
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[1]-1.0）*（x[1]-1.0）
+x[0]+x[1]+10.0;
/*阻尼牛顿法*/
returnf;
/*函数的梯度*/
voidgread（doublex[],doublegf[]）
gf[0]=8.0*（x[0]+1.0）+1.0;
gf[1]=4.0*（x[1]-1.0）+1.0;
/*Heisen矩阵*/
voidheisen（doublex[],doublehei[][N]）
hei[0][0]=8.0;
hei[0][1]=0.0;
hei[1][0]=0.0;
hei[1][1]=4.0;
/*计算f（xk+as）*/
doublexkadd（doublex[],doubled[],doublea）
doublex1[N];
i++）
x1[i]=x[i]+a*d[i];
returnhanshu（x1）;
/*输出选定的一维迭代方法*/
voidywddf（intyw）
switch（yw）
{
case1:
fprintf（fp,"
黄金分割法++++++\n\n"
）;
break;
case2:
二次插值法++++++\n\n"
}
/*输出当前迭代点坐标及目标函数值*/
doublexfout（doublex[],intm）
intj;
f=hanshu（x）;
if（-1==m）
fprintf（fp,"
x（*）=["
else
x（%3d）=["
m）;
for（j=0;
j<
N-1;
j++）
%15.7lf,"
x[j]）;
%15.7lf],f（*）=%15.7lf\n"
x[N-1],f）;
%15.7lf],f（%3d）=%15.7lf\n"
x[N-1],m,f）;
/*初始搜索区间的确定*/
voidcsssqj（doublex[],doubled[],doubleh,doubleab[]）
doublea1,a2,a3,f1,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;
f1=f2;
f2=f3;
f3=f1;
h=-h;
do
a1=a2;
a3=a2+h;
f3=xkadd（x,d,a3）;
h=2*h;
}while（f3<
f2）;
if（h>
0.0）
ab[0]=a1;
ab[1]=a3;
ab[0]=a3;
ab[1]=a1;
/*黄金分割法*/
voidgoldcut（doublex[],doubled[],doubleh,doubleebsin）
doublea1,a2,f1,f2,a,b,ab[2];
csssqj（x,d,h,ab）;
a=ab[0];
b=ab[1];
a1=b-0.618*（b-a）;
f1=xkadd（x,d,a1）;
a2=a+0.618*（b-a）;
if（f1>
{
a=a1;
a1=a2;
f1=f2;
a2=a+0.618*（b-a）;
f2=xkadd（x,d,a2）;
}
else
b=a2;
a2=a1;
f2=f1;
a1=b-0.618*（b-a）;
f1=xkadd（x,d,a1）;
}while（b-a>
ebsin）;
x[i]+=（b+a）*d[i]/2;
/*二次插值法*/
intrccz（doublex[],doubled[],doubleh,doubleebsin）
doublea1,a2,a3,a4,f1,f2,f3,f4,c1,c2,ab[2];
inti,p=0,k=0;
a1=ab[0];
a3=ab[1];
a2=（a1+a3）/2;
while
（1）
c1=（f3-f1）/（a3-a1）;
c2=（（f2-f1）/（a2-a1）-c1）/（a2-a3）;
if（0==c2）
p=1;
a4=0.5*（a1+a3-c1/c2）;
if（（a4-a1）*（a3-a4）<
=0.0）
p=2;
f4=xkadd（x,d,a4）;
if（1==k&
&
fabs（a4-a2）<
=ebsin）
if（a4<
a2）
if（f2>
f4）
{
f1=f2;
a1=a2;
k=1;
a2=a4;
f2=f4;
}
else
a3=a4;
f3=f4;
f3=f2;
a3=a2;
a1=a4;
f1=f4;
if（0==p）
f1=f4-f2?
a4:
a2;
f1=a2;
x[i]+=f1*d[i];
returnp;
/*矩阵求逆*/
voidatoa_1（doublea[][N],doubleb[][N]）
inti,j,k,m;
doublefact,ebs=0.0;
for（j=0;
b[i][j]=0.0;
ebs+=a[i][j]*a[i][j];
ebs=sqrt（ebs）/（N*N*100）;
b[i][i]=1.0;
if（fabs（a[i][i]）<
ebs）
for（k=i;
k<
k++）
if（fabs（a[i][i]）>
{
for（m=0;
m<
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;
a[i][k]/=fact;
b[i][k]/=fact;
i;
fact=a[k][i];
for（m=0;
a[i][m]-=a[i][m]*fact;
b[i][m]-=b[i][m]*fact;
for（k=i+1;
/*阻尼牛顿法*/
voidnewtow_method（doublex[],doubleh,doubleebsin,intyw）
doubled[N],t,f,s[N],hei[N][N],hei_1[N][N];
inti,j,k=1;
fprintf（fp,"
**********阻尼牛顿法计算结果**********\n\n"
ywddf（yw）;
/*输出一维迭代方法*/
\n"
f=xfout（x,0）;
gread（x,d）;
heisen（x,hei）;
atoa_1（hei,hei_1）;
for（i=0;
s[i]=0.0;
for（j=0;
s[i]+=hei_1[i][j]*d[j];
switch（yw）
case1:
goldcut（x,s,h,ebsin）;
break;
case2:
j=rccz（x,s,h,ebsin）;
gread（x,d）;
迭代轮数k=%3d\n"
k）;
f=xfout（x,k）;
t=0;
t+=d[i]*d[i];
t=sqrt（t）;
%15.9lf\n"
t）;
k++;
}while（t>
\n*********************\n"
f=xfout（x,-1）;
"
main（）
doublex0[N],h,ebsin;
intyw;
csd_x（x0）;
/*初始坐标*/
h=H_QJ;
ebsin=EPSIN;
yw=Y_F;
fp=fopen（outname,"
w"
newtow_method（x0,h,ebsin,yw）;
/*调用阻尼牛顿法*/
fclose（fp）;
return0;
备注：
不交此报告者，本次实验为“不合格”。
升级会员 