阻尼牛顿法.docx

阻尼牛顿法求解二维函数极小值的程序说明
一 题目
利用阻尼牛顿法求函数的极小值点（迭代两次，一维搜索任选一种方法）。
二 阻尼牛顿法基本思想:
1) 给定初始点，收敛精度，置。
2) 计算、 、和
3) 求，其中为沿进行一维搜索的最佳步长。
4) 检查收敛精度。若，则，停机；否则置，返回步骤2，继续进行进行搜索。
改进后的阻尼牛顿法程序框图如下：

三 用阻尼牛顿法求函数
程序如下：
// 阻尼牛顿法 .cpp : Defines the entry point for the console application.
//

#include
#include
#include
#include
double fun1(double q1,double q2)
{
return(pow((q1-2),4)+pow((q1-2*q2),2)); //修改函数f(x1,x2)=(x1-2)*(x1-2)*(x1-2)*(x1-2)+(x1-2*x2)*(x1-2*x2)
}
double fun2(double g,double x,double y,double r1,double r2)
{
return (pow((x+g*y-2),4)+pow((x+g*y-2*(r1+g*r2)),2));//关于阻尼因子的函数
}
void main()
{
double A[2][1],B[2][2],C[2][1],D[2][1],X[2][1];
double E[2][1]={4,3};//迭代的初始点x0
int t=0,i=0,j=0;
double E0,x1,x2,x3,h(0.1);
double y1,y2,y3,m;
double a,b,k=0.618,a1,a2,f1,f2;
printf("输入收敛精度：");//输入标准收敛精度
std::cin>>E0;
do
{
D[0][0]=E[0][0];
D[1][0]=E[1][0];
A[0][0]=4*(D[0][0]-2)*(D[0][0]-2)*(D[0][0]-2)+2*D[0][0]-4*D[1][0];
A[1][0]=-4*(D[0][0]-2*D[1][0]);//A[0][0],A[1][0]为原函数梯度的各项
B[0][0]=1.0/(12.0*(D[0][0]-2)*(D[0][0]-2));
B[0][1]=1.0/(24.0*(D[0][0]-2)*(D[0][0]-2));
B[1][0]=1.0/(24.0*(D[0][0]-2)*(D[0][0]-2));
B[1][1]=(6.0*(D[0][0]-2)*(D[0][0]-2)+1)/(48.0*(D[0][0]-2)*(D[0][0]-2));//B[0][0],B[0][1],B[1][0],B[1][1]分别代表原函数的海赛矩阵的逆阵的各项
C[0][0]=-(B[0][0]*A[0][0]+B[0][1]*A[1][0]);
C[1][0]=-(B[1][0]*A[0][0]+B[1][1]*A[1][0]);//C[0][0],C[1][0]为搜索方向dk的各项
//下面利用外推法寻找函数2的区间，找单谷区间
x1=0;
x2=x1+h;
y1=fun2(x1,D[0][0],C[0][0],D[1][0],C[1][0]);
y2=fun2(x2,D[0][0],C[0][0],D[1][0],C[1][0]);
if(y2>y1)
{
h=-h;
x3=x1,y3=y1;
x1=x2,y1=y2;
x2=x3,y2=y3;
}
x3=x2+h;
y3=fun2(x3,D[0][0],C[0][0],D[1][0],C[1][0]);
while(y3 x2=x3,y2=y3;
x3=x2+h;
y3=fun2(x3,D[0][0],C[0][0],D[1][0],C[1][0]);
i++;
}
//下面利用黄金分割法寻找函数2极值
a=x1;
b=x3;
a1=b-k*(b-a);
a2=a+k*(b-a);
f1=fun2(a1,D[0][0],C[0][0],D[1][0],C[1][0]);
f2=fun2(a2,D[0][0],C[0][0],D[1][0],C[1][0]);
do
{
if(f1>=f2)
{
a=a1;
a1=a2;
f1=f2;
a2=a+k*(b-a);
f2=fun2(a2,D[0][0],C[0][0],D[1][0],C[1][0]);
}
else
{
b=a2;
a2=a1;
f2=f1;
a1=b-k*(b-a);
f1=fun2(a1,D[0][0],C[0][0],D[1][0],C[1][0]);
}
j++;
}while(fabs((b-a)/b)>=E0&&fabs((f2-f1)/f2)>=E0);
m=0.5*(a+b);//m为阻尼因子
E[0][0]=D[0][0]+m*C[0][0];
E[1][0]=D[1][0]+m*C[1][0];
printf("%d%15f10%15f10\n",t,E[0][0],E[1][0],fun1(E[0][0],E[1][0]));
t++;
}while(fabs(E[0][0]-D[0][0])>=E0&&fabs(E[1][0]-D[1][0])>=E0);
X[0][0]=E[0][0];
X[1][0]=E[1][0];
printf("迭代了%d次\n",t);
printf("极小点(x1,x2)=(%f10,%f10)\n",X[0][0],X[1][0]);
printf("极小值f(x1,x2)=%f10\n",fun1(X[0][0],X[1][0]));
}
程序运行结果：

四 结论
由该程序的运行结果可知，要求迭代两次后函数的极小值点在（2.52939210,1.46320410）处，经验证，运算结果完全正确。验证了该程序的可行性。