dfp.txt

本程序用于机械优化设计工程计算

#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;
	}
}