grads.cpp

最优化理论与算法（第2版）这本书中的课后作业。用C++实现的一些具体算法。

#include<iostream.h>
#include<math.h>

void grads(double A[2][2],double b[2],double c,double x0[2]);

void main()
{
	double A1[2][2]={{1,0},{0,2}};
	double b1[2]={0,0};
	double c1=0;
	double x01[2]={4,4};

	cout<<"共轭梯度法求第(1)小题的最小值："<<endl;
	grads(A1,b1,c1,x01);

	double A2[2][2]={{2,-4},{0,4}};
	double b2[2]={0,2};
	double c2=0;
	double x02[2]={0,0};

	cout<<"共轭梯度法求第(2)小题的最小值："<<endl;
	grads(A2,b2,c2,x02);


	double A3[2][2]={{2,0},{0,4}};
	double b3[2]={-4,-4};
	double c3=6;
	double x03[2]={1,3};

	cout<<"共轭梯度法求第(3)小题的最小值："<<endl;
	grads(A3,b3,c3,x03);

	double A4[2][2]={{4,4},{0,2}};
	double b4[2]={3,-4};
	double c4=0;
	double x04[2]={3,4};

	cout<<"共轭梯度法求第(4)小题的最小值："<<endl;
	grads(A4,b4,c4,x04);

	double A5[2][2]={{4,4},{0,10}};
	double b5[2]={0,0};
	double c5=0;
	double x05[2]={2,-2};

	cout<<"共轭梯度法求第(5)小题的最小值："<<endl;
	grads(A5,b5,c5,x05);
}

void grads(double A[2][2],double b[2],double c,double x0[2])
{
    double fk[2];
	double gk[2],g;
	double dk[2];
	double langta;
	double x1[2];
	double betak;
	int n=50;
	double y1,y2,y3,y;

	fk[1]=A[1][1]*x0[1]+A[1][2]*x0[2]+b[1];
	fk[2]=A[2][1]*x0[1]+A[2][2]*x0[2]+b[2];

	//double gk[2],g;
	gk[1]=fk[1];
	gk[2]=fk[2];
	//g=sqrt(pow(gk[1],2)+pow(gk[2],2));

	//double dk[2];
	dk[1]=-gk[1];
	dk[2]=-gk[2];

	//double langta;
	langta=-((gk[1]*dk[1]+gk[2]*dk[2])/((dk[1]*A[1][1]+dk[2]*A[2][1])*dk[1]+(dk[1]*A[1][2]+dk[2]*A[2][2])*dk[2]));
	
	//double x1[2];
	x1[1]=x0[1]+langta*dk[1];
	x1[2]=x0[2]+langta*dk[2];

	//double betak;
	//int n=50;

	for(int k=1;k<n;k++)
	{
	gk[1]=A[1][1]*x1[1]+A[1][2]*x1[2]+b[1];
	gk[2]=A[2][1]*x1[1]+A[2][2]*x1[2]+b[2];
	g=sqrt(pow(gk[1],2)+pow(gk[2],2));

	if(g==0) break;

	betak=(pow(gk[1],2)+pow(gk[2],2))/(pow(dk[1],2)+pow(dk[2],2));
	dk[1]=betak*dk[1]-gk[1];
	dk[2]=betak*dk[2]-gk[2];

	langta=-((gk[1]*dk[1]+gk[2]*dk[2])/((dk[1]*A[1][1]+dk[2]*A[2][1])*dk[1]+(dk[1]*A[1][2]+dk[2]*A[2][2])*dk[2]));
	x1[1]=x1[1]+langta*dk[1];
	x1[2]=x1[2]+langta*dk[2];

	dk[1]=-gk[1];
	dk[2]=-gk[2];
	}

	//double y1,y2,y3,y;
	y1=(x1[1]*A[1][1]+x1[2]*A[2][1])*x1[1]+(x1[1]*A[1][2]+x1[2]*A[2][2])*x1[2];
	y2=b[1]*x1[1]+b[2]*x1[2];
	y3=c;
	y=y1+y2+y3;

	cout<<"\t\t\t\t"<<y<<endl;
	cout<<x1[1]<<endl;
	cout<<x1[2]<<endl;
	//return y;
}