lwq.cpp

规划数学 数学专业的人应该能用的到的一个小 程序。

#include<iostream.h>
#include<math.h>
#include<stdlib.h>
const double dps=1e-5;//算法中用到的精度
const double dps1=1e-15;//求方程根（得到λ）用到的精度。这个精度比较高为了确保λ比较精确
const int max=100;//解方程求λ时为了判断无正根而设置的
void make2darray(double **&x,int rows,int cols)//动态申请2维数组
{ x=new double*[rows];
  for(int i=0;i<rows;i++) x[i]=new double[cols];
}
void delete2darray(double ** &x,int rows)//收回分配的2维数组
{ for(int i=0;i<rows;i++)
   delete[] x[i];
  delete []x;
  x=0;
}
void calp(double **f,double *p,double *x,int xnum,int step)//step为最高次幂,xnum变元数量
{//第一步骤(1)
	for(int i=0;i<xnum;i++)
	{
		double t=1.0;
		 p[i]=0;
		  for(int j=0;j<step;j++)
		  {
			  p[i]+=f[i][j+1]*(j+1)*t;
			   t*=x[i];
		  }
		  p[i]=-p[i];
	}
}
double value(double *p,int xnum)//求||P||，这里没有加平方根
{
	double s=0;
	 for(int i=0;i<xnum;i++)
		 s+=p[i]*p[i];
	   return s;
}
void changex(double *x,double *p,double num,int xnum)//第（5）步
{
	for(int i=0;i<xnum;i++)
		x[i]+=num*p[i];
}
double calpx(double *px,int step,double x)//计算函数的值，在求λ解方程时用到。
{ 
	double sum=0,t=1.0;
	 for(int i=0;i<=step;i++)
	 {
		 sum+=px[i]*t;
		  t*=x;
	 }
	  return sum;
}
double getfvalue(double **f,int step,double *x,int xnum)//计算f(X)
{
	double sum=0,t;
	 for(int i=0;i<xnum;i++)
	 {
		 t=1.0;
		 for(int j=0;j<=step;j++)
		 {
			 sum+=f[i][j]*t;
			  t*=x[i];
		 }
	 }
	 return sum;
}
double getsolution(double *px,int step)//解方程，得到一个正根,用二分法
{
	double a=0,b=1.21,mid;    
	 int flag=0;
	  while(calpx(px,step,a)*calpx(px,step,b)>0&&flag<max)
	  {  
		  flag++;
		   b*=2;
	  }
	  if(flag>=max)
	  {
		  cout<<"找不到大于0的λ!"<<endl;
		   system("pause");
		    exit(0);
	  }
	  while(fabs(b-a)>dps1)
	  {  
		  mid=(a+b)/2;
		   if(calpx(px,step,mid)*calpx(px,step,a)>0)
			   a=mid;
		    else
			   b=mid;
	  }        
	  return (a+b)/2;
}
int getc(int n,int i)//计算组合数
{
   if(i<=0) return 1;
	int up=1,down=1;
	 for(int j=0;j<n;j++)
		 up*=(j+1);
	  for(j=0;j<i;j++)
		  down*=(j+1);
	   for(j=0;j<n-i;j++)
		   down*=(j+1);
	    return up/down;
}
double cimi(double x,int n)//计算x的n次幂
{
	double s=1.0;
	 for(int i=0;i<n;i++)
		 s*=x;
	   return s;
}
double getnum(double **f,double *p,double *x,int xnum,int step)//求得满足要求的λ
{
	double **f1,*temp=new double[step+1],answer;
	 make2darray(f1,xnum,step+1);
	  for(int i=0;i<xnum;i++)
		  for(int j=0;j<=step;j++)
			  f1[i][j]=f[i][j];                 
		  for(i=0;i<xnum;i++)
			  for(int j=0;j<=step;j++)
			  {
				  double bak=f1[i][j];
				   f1[i][j]=0;
				    for(int k=0;k<=j;k++)
					   f1[i][k]+=bak*getc(j,k)*cimi(x[i],j-k)*cimi(p[i],k);
			  }          
				  for(i=0;i<step;i++)
				  {
		             temp[i]=0;
		              for(int j=0;j<xnum;j++)
			               temp[i]+=f1[j][i+1]*(i+1);
				  }
				     temp[step]=0;
	  answer=getsolution(temp,step);
	   delete2darray(f1,xnum);
	    delete []temp;
		 return answer;
}
void show(double *p,int  xnum,int type,int flag)//type==0显示x，==1显示p
{
	if(type==0)
		cout<<"X"<<flag<<"=(";
	 else
		 cout<<"P"<<flag<<"=(";
	  for(int i=0;i<xnum-1;i++)
		  cout<<p[i]<<",";
	        cout<<p[i]<<")"<<endl;
}
void main()
{
	double **f,*p,*x;
	 int xnum,step,flag=0;
	  cout<<"请输入变元的个数!"<<endl;
	   cin>>xnum;
	    cout<<"请输入变元的最高次幂!"<<endl;
		 cin>>step;
	  make2darray(f,xnum,step+1);
	   p=new double[xnum];
	    x=new double[xnum];
		 for(int i=0;i<xnum;i++)
		 {
			cout<<"请输入第"<<i+1<<"个变量的系数，从低到高!("<<step+1<<"个)"<<endl;
			 for(int j=0;j<=step;j++)
				 cin>>f[i][j];
		 }
		 cout<<"请输入初始点!从低到高!("<<xnum<<"个)"<<endl;
		  for(i=0;i<xnum;i++)
			  cin>>x[i];   
		 calp(f,p,x,xnum,step);    
		   while(value(p,xnum)>dps)
		   { 
			  cout<<"第"<<flag<<"次迭代--------------------"<<endl;
			   show(p,xnum,1,flag);
			    show(x,xnum,0,flag);
				 cout<<"f(";
				 for(i=0;i<xnum;i++)
					 cout<<"X"<<i<<",";
				      cout<<'\b'<<")="<<getfvalue(f,step,x,xnum)<<endl;
			     double num=getnum(f,p,x,xnum,step);
			      cout<<"λ="<<num<<endl;
			       changex(x,p,num,xnum);
				    calp(f,p,x,xnum,step);
					 flag++;
		   }
		   for(i=0;i<xnum;i++)
			   for(int j=0;j<=step;j++)
			   {
				   if(fabs(f[i][j])>dps1&&j!=0)
				     cout<<f[i][j]<<"X"<<i<<"(^"<<j<<")+";
				    else
						if(fabs(f[i][j])>dps1)
							cout<<f[i][j]<<"+";
			   }
			         cout<<'\b'<<"=F(";
					  for(i=0;i<xnum;i++)
						  cout<<"X"<<i<<",";
					        cout<<'\b'<<")的最小值解为："<<endl;
		      show(x,xnum,0,flag);
			    cout<<"f(";
				  for(i=0;i<xnum;i++)
					 cout<<"X"<<i<<",";
				      cout<<'\b'<<")min="<<getfvalue(f,step,x,xnum)<<endl;
			   system("pause");
		          delete2darray(f,xnum);
				   delete []p;
				    delete []x;
		  
}