#include <iostream.h>
#include <math.h>
const int dim=2;
//求解函数值：
double f(double t,double x[],double a[],int y){
double z;
z=4*((x[0]+t*a[0])-5)*((x[0]+t*a[0])-5)+((x[1]+t*a[1])-6)*((x[1]+t*a[1])-6);
return z;
}
//求解梯度：
double f_1(double array[],int k){
double m;
m=8*(array[0]-5);
return m;
}
double f_2(double array[],int k){
double m;
m=2*(array[1]-6);
return m;
}
//求步长和极值点：
double r(double x[],double a[],int y){
double a1,a2,a3,h,h0,y1,y2,y3,m_num3,m_num4,m_middle,m_h;
double x_x[dim],x_a[dim];
for (int i=0;i<dim;i++){
x_x[i]=x[i];
x_a[i]=a[i];
}
// float f(float t,float x[],float a[],int y);
a1=0;
y1=f(a1,x_x,x_a,2);
h0=1;
h=h0;
a2=h; y2=f(a2,x_x,x_a,2);
if(y2>y1){
h=(-h);
a3=a1;y3=y1;
a1=a2;y1=y2;
a2=a3;y2=y3;
a3=a2+h;
y3=f(a3,x_x,x_a,2);
while(y3<y2){
h=2*h;
a1=a2;y1=y2;
a2=a3;y2=y3;
a3=a2+h;
y3=f(a3,x_x,x_a,2);
}
if(a1<a3){
m_num3=a1;m_num4=a3;
}
else{
m_num3=a3;m_num4=a1;
}
}
if(y1>=y2){
a3=a2+h; y3=f(a3,x_x,x_a,2);
while(y3<y2){
h=2*h;
a1=a2;y1=y2;
a2=a3;y2=y3;
a3=a2+h;y3=f(a3,x_x,x_a,2);
}
m_num3=a1,m_num4=a3;
m_middle=a2;m_h=h;
}
//求极值点：
double a_0,a_1,a_2,a_m,b,y_1,y_2,y_0,c,d;
if(m_num3<m_num4){
a_0=m_num3;b=m_num4;
}
else{
a_0=m_num4;b=m_num3;
}
c=0.618;
d=0.000001;
a_1=b-c*(b-a_0); y_1=f(a_1,x_x,x_a,2);
a_2=a_0+c*(b-a_0); y_2=f(a_2,x_x,x_a,2);
if(y_1>=y_2){
a_0=a_1,a_1=a_2,y_1=y_2;
a_2=a_0+c*(b-a_0),y_2=f(a_2,x_x,x_a,2);
}
else{
b=a_2,a_2=a_1,y_2=y_1;
a_1=b-c*(b-a_0),y_1=f(a_1,x_x,x_a,2);
}
while(fabs((b-a_0)/b)>=d){
if(y_1>=y_2){
a_0=a_1;a_1=a_2;y_1=y_2;
a_2=a_0+c*(b-a_0);y_2=f(a_2,x_x,x_a,2);
}
else{
b=a_2,a_2=a_1,y_2=y_1;
a_1=b-c*(b-a_0),y1=f(a_1,x_x,x_a,2);
}
}
a_m=(a_0+b)/2; double m_num6=a_m;
y_0=f(a_m,x_x,x_a,2); double m_num7=y_0;
return a_m;
}
//主函数：
void main(){
double x[dim][dim],h[dim][dim],d[dim],f_dh[dim],b[dim],f_de[dim],e;
double g[2],x_x[2],h_h[dim][dim],h_1[dim][dim],h_2[dim][dim],hg[dim],
hg_2[dim][dim],gh[2],ghg;
int k=0;
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
if(i==j)
h[i][j]=1;
else
h[i][j]=0;
}
}
cout<<"请输入初始点"<<endl;
cin>>x[0][0]>>x[0][1];
cout<<"请输入精度："<<endl;
cin>>e;
double f_1(double array[],int k);
double f_2(double array[],int k);
f_dh[0]=f_1(x[0],dim);
f_dh[1]=f_2(x[0],dim);
if(sqrt(f_dh[0]*f_dh[0]+f_dh[1]*f_dh[1])<e){
cout<<"极值点为："<<x[1][0]<<" "<<x[1][1]<<endl;
cout<<"极值为："<<f(0,x[0],h[0],2)<<endl;
}
else{
d[0]=h[0][0]*f_dh[0]+h[0][1]*f_dh[1];
d[1]=h[1][0]*f_dh[0]+h[1][1]*f_dh[1];
double r(double x[],double a[],int y);
b[0]=r(x[0],d,dim);
x[1][0]=x[0][0]+b[0]*d[0];
x[1][1]=x[0][1]+b[0]*d[1];
f_de[0]=f_1(x[1],dim);
f_de[1]=f_2(x[1],dim);
while(sqrt(f_de[0]*f_de[0]+f_de[1]*f_de[1])>e){
if(k<(dim-1)){
g[0]=f_de[0]-f_dh[0];
g[1]=f_de[1]-f_dh[1];
x_x[0]=x[1][0]-x[0][0];
x_x[1]=x[1][1]-x[0][1];
h_1[0][0]=x_x[0]*x_x[0]/(x_x[0]*g[0]+x_x[1]*g[1]);
h_1[0][1]=x_x[0]*x_x[1]/(x_x[0]*g[0]+x_x[1]*g[1]);
h_1[1][0]=x_x[1]*x_x[0]/(x_x[0]*g[0]+x_x[1]*g[1]);
h_1[1][1]=x_x[1]*x_x[1]/(x_x[0]*g[0]+x_x[1]*g[1]);
hg[0]=h[0][0]*g[0]+h[0][1]*g[1];
hg[1]=h[1][0]*g[0]+h[1][1]*g[1];
hg_2[0][0]=hg[0]*g[0];
hg_2[0][1]=hg[0]*g[1];
hg_2[1][0]=hg[1]*g[0];
hg_2[1][1]=hg[1]*g[1];
h_2[0][0]=hg_2[0][0]*h[0][0]+hg_2[0][1]*h[0][1];
h_2[0][1]=hg_2[0][0]*h[1][0]+hg_2[0][1]*h[1][1];
h_2[1][0]=hg_2[1][0]*h[0][0]+hg_2[1][1]*h[0][1];
h_2[1][1]=hg_2[1][1]*h[0][1]+hg_2[1][1]*h[1][1];
gh[0]=g[0]*h[0][0]+g[1]*h[1][0];
gh[1]=g[0]*h[0][1]+g[1]*h[1][1];
ghg=gh[0]*g[0]+gh[1]*g[1];
h_2[0][0]=h_2[0][0]/ghg;
h_2[0][1]=h_2[0][1]/ghg;
h_2[1][0]=h_2[1][0]/ghg;
h_2[1][1]=h_2[1][1]/ghg;
for(i=0;i<dim;i++){
for(int j=0;j<dim;j++){
h_h[i][j]=h[i][j]+h_1[i][j]-h_2[i][j];
}
k=k+1;
}
for(i=0;i<dim;i++){
for(int j=0;j<dim;j++){
h[i][j]=h_h[i][j];
}
}
x[0][0]=x[1][0];
x[0][1]=x[1][1];
f_dh[0]=f_1(x[0],dim);
f_dh[1]=f_2(x[0],dim);
d[0]=-(h[0][0]*f_dh[0]+h[0][1]*f_dh[1]);
d[1]=-(h[1][0]*f_dh[0]+h[1][1]*f_dh[1]);
b[0]=r(x[0],d,dim);
x[1][0]=x[0][0]+b[0]*d[0];
x[1][1]=x[0][1]+b[0]*d[1];
f_de[0]=f_1(x[1],dim);
f_de[1]=f_2(x[1],dim);
}
else{
x[0][0]=x[1][0];
x[0][1]=x[1][1];
for(i=0;i<dim;i++){
for(int j=0;j<dim;j++){
if(i==j)
h[i][j]=1;
else
h[i][j]=0;
}
}
k=0;
f_dh[0]=f_1(x[0],dim);
f_dh[1]=f_2(x[0],dim);
d[0]=h[0][0]*f_dh[0]+h[0][1]*f_dh[1];
d[1]=h[1][0]*f_dh[0]+h[1][1]*f_dh[1];
b[0]=r(x[0],d,dim);
x[1][0]=x[0][0]+b[0]*d[0];
x[1][1]=x[0][1]+b[0]*d[1];
f_de[0]=f_1(x[1],dim);
f_de[1]=f_2(x[1],dim);
}
}
cout<<"极值点为："<<x[1][0]<<" "<<x[1][1]<<endl;
cout<<"极值为："<<f(b[0],x[0],d,2)<<endl;
}
}
