1.h

非线性方程组解法牛顿下山法,求解非线性方程组

#include <math.h>
class newton_diedai
{private:
double xi,yi,zi;//迭代初始向量
double x[8],y[8],z[8];//8个节点坐标
int count;//  数据样本数
int index[5];//记录乔里斯基分解中的行交换号
double A[4][4];
double B[4]; //A*X=B,A、B矩阵
double b[4];//b为迭代向量，后返回结果
double xx,yy,zz;
double Ni[8];//节点分配系数
public:
void input(double *X,double *Y,double *Z,double *x0,double *y0,double *z0)
{for(int i=0;i<8;i++)
 {x[i]=X[i];
  y[i]=Y[i];
  z[i]=Z[i];
 }
xi=1;yi=0;zi=0;
xx=*x0;yy=*y0;zz=*z0;
//cout<<xx<<yy<<zz<<endl;
}
void calculate_matrix_f()//计算系数矩阵
{A[1][1]=((1+yi)*(1+zi)*x[0]+(1-yi)*(1+zi)*x[1]+(1-yi)*(1-zi)*x[2]+(1+yi)*(1-zi)*x[3]
        -(1+yi)*(1+zi)*x[4]-(1-yi)*(1+zi)*x[5]-(1-yi)*(1-zi)*x[6]-(1+yi)*(1-zi)*x[7])/8;
 A[1][2]=((1+xi)*(1+zi)*x[0]-(1+xi)*(1+zi)*x[1]-(1+xi)*(1-zi)*x[2]+(1+xi)*(1-zi)*x[3]
        +(1-xi)*(1+zi)*x[4]-(1-xi)*(1+zi)*x[5]-(1-xi)*(1-zi)*x[6]+(1-xi)*(1-zi)*x[7])/8; 
 A[1][3]=((1+xi)*(1+yi)*x[0]+(1+xi)*(1-yi)*x[1]-(1+xi)*(1-yi)*x[2]-(1+xi)*(1+yi)*x[3]
        +(1-xi)*(1+yi)*x[4]+(1-xi)*(1-yi)*x[5]-(1-xi)*(1-yi)*x[6]-(1-xi)*(1+yi)*x[7])/8;
 A[2][1]=((1+yi)*(1+zi)*y[0]+(1-yi)*(1+zi)*y[1]+(1-yi)*(1-zi)*y[2]+(1+yi)*(1-zi)*y[3]
        -(1+yi)*(1+zi)*y[4]-(1-yi)*(1+zi)*y[5]-(1-yi)*(1-zi)*y[6]-(1+yi)*(1-zi)*y[7])/8;
 A[2][2]=((1+xi)*(1+zi)*y[0]-(1+xi)*(1+zi)*y[1]-(1+xi)*(1-zi)*y[2]+(1+xi)*(1-zi)*y[3]
        +(1-xi)*(1+zi)*y[4]-(1-xi)*(1+zi)*y[5]-(1-xi)*(1-zi)*y[6]+(1-xi)*(1-zi)*y[7])/8;
 A[2][3]=((1+xi)*(1+yi)*y[0]+(1+xi)*(1-yi)*y[1]-(1+xi)*(1-yi)*y[2]-(1+xi)*(1+yi)*y[3]
        +(1-xi)*(1+yi)*y[4]+(1-xi)*(1-yi)*y[5]-(1-xi)*(1-yi)*y[6]-(1-xi)*(1+yi)*y[7])/8;
 A[3][1]=((1+yi)*(1+zi)*z[0]+(1-yi)*(1+zi)*z[1]+(1-yi)*(1-zi)*z[2]+(1+yi)*(1-zi)*z[3]
        -(1+yi)*(1+zi)*z[4]-(1-yi)*(1+zi)*z[5]-(1-yi)*(1-zi)*z[6]-(1+yi)*(1-zi)*z[7])/8;
 A[3][2]=((1+xi)*(1+zi)*z[0]-(1+xi)*(1+zi)*z[1]-(1+xi)*(1-zi)*z[2]+(1+xi)*(1-zi)*z[3]
        +(1-xi)*(1+zi)*z[4]-(1-xi)*(1+zi)*z[5]-(1-xi)*(1-zi)*z[6]+(1-xi)*(1-zi)*z[7])/8;
 A[3][3]=((1+xi)*(1+yi)*z[0]+(1+xi)*(1-yi)*z[1]-(1+xi)*(1-yi)*z[2]-(1+xi)*(1+yi)*z[3]
        +(1-xi)*(1+yi)*z[4]+(1-xi)*(1-yi)*z[5]-(1-xi)*(1-yi)*z[6]-(1-xi)*(1+yi)*z[7])/8;
/*for(int i=1;i<=3;i++)
 {for(int j=1;j<=3;j++)
{cout<<A[i][j]<<endl;
}
}*/
}

// 杜利特尔分解
void cludcmp()
{
 int n=0;
 int i,imax,j,k;
 double big,dum,sum,temp,vv[4],d=1.0;
 for(i=1;i<=3;i++)
    {big=0.0;
     for(j=1;j<=3;j++)
        {if((temp=fabs(A[i][j]))>big) big=temp;
        }
     vv[i]=1.0/big;
     }
 for(j=1;j<=3;j++)
    {for(i=1;i<j;i++)
        {sum=A[i][j];
         for(k=1;k<i;k++) sum-=A[i][k]*A[k][j];
         A[i][j]=sum;
        }
     big=0.0;
     for(i=j;i<=3;i++)
        {sum=A[i][j];
         for(k=1;k<j;k++) sum-=A[i][k]*A[k][j];
         A[i][j]=sum;
         if((dum=vv[i]*fabs(sum))>=big) {big=dum;imax=i;}
        }
     if(j!=imax)
       {for(k=1;k<=3;k++)
           {dum=A[imax][k];A[imax][k]=A[j][k];A[j][k]=dum;
           }
        d=-d;
        vv[imax]=vv[j];
        }
      index[j]=imax;
      if(j!=n)
        {dum=1.0/A[j][j];
         for(i=j+1;i<=3;i++) A[i][j]*=dum;
        }
     }
 for( i=1;i<=3;i++)
{for( j=1;j<=3;j++)
{cout<<A[i][j]<<endl;
}
}
}
void clufbsb()
 {
//  int n;
  int i,ii=0,ip,j;
  double sum;
  for(i=1;i<=3;i++)
     {ip=index[i];sum=B[ip];B[ip]=B[i];
      if(ii) {for(j=ii;j<=i-1;j++) sum-=A[i][j]*B[j];}
      else if(sum) ii=i;
      B[i]=sum;
     }
  for(i=3;i>=1;i--)
     {sum=B[i];
      for(j=i+1;j<=3;j++) sum-=A[i][j]*B[j];
      B[i]=sum/A[i][i];
     }
//  cout<<B[1]<<"  "<<B[2]<<" "<<B[3]<<endl;  
 }
void calculate()
{calculate_matrix_f();
  cludcmp();
double max;
 b[1]=xi;b[2]=yi;b[3]=zi;
// cout<<B[1]<<"  "<<B[2]<<" "<<B[3]<<endl; 
for(;;)
 {B[1]=((1+b[1])*(1+b[2])*(1+b[3])*x[0]+(1+b[1])*(1-b[2])*(1+b[3])*x[1]+(1+b[1])*(1-b[2])*(1-b[3])*x[2]+(1+b[1])*(1+b[2])*(1-b[3])*x[3]
      +(1-b[1])*(1+b[2])*(1+b[3])*x[4]+(1-b[1])*(1-b[2])*(1+b[3])*x[5]+(1-b[1])*(1-b[2])*(1-b[3])*x[6]+(1-b[1])*(1+b[2])*(1-b[3])*x[7])/8-xx;
  B[2]=((1+b[1])*(1+b[2])*(1+b[3])*y[0]+(1+b[1])*(1-b[2])*(1+b[3])*y[1]+(1+b[1])*(1-b[2])*(1-b[3])*y[2]+(1+b[1])*(1+b[2])*(1-b[3])*y[3]
      +(1-b[1])*(1+b[2])*(1+b[3])*y[4]+(1-b[1])*(1-b[2])*(1+b[3])*y[5]+(1-b[1])*(1-b[2])*(1-b[3])*y[6]+(1-b[1])*(1+b[2])*(1-b[3])*y[7])/8-yy;
  B[3]=((1+b[1])*(1+b[2])*(1+b[3])*z[0]+(1+b[1])*(1-b[2])*(1+b[3])*z[1]+(1+b[1])*(1-b[2])*(1-b[3])*z[2]+(1+b[1])*(1+b[2])*(1-b[3])*z[3]
      +(1-b[1])*(1+b[2])*(1+b[3])*z[4]+(1-b[1])*(1-b[2])*(1+b[3])*z[5]+(1-b[1])*(1-b[2])*(1-b[3])*z[6]+(1-b[1])*(1+b[2])*(1-b[3])*z[7])/8-zz;
  B[1]=-B[1];B[2]=-B[2];B[3]=-B[3];
  // cout<<B[1]<<"  "<<B[2]<<" "<<B[3]<<endl;
//calculate_matrix_f();
//cludcmp();
  clufbsb();
max=fabs(B[1])>fabs(B[2])?fabs(B[1]):fabs(B[2]);max=max>fabs(B[3])?max:fabs(B[3]);  
b[1]+=B[1];b[2]+=B[2];b[3]+=B[3];
//xi=b[1];yi=b[2];zi=b[3];
if(max<1e-5) break;
}
//******************************8
//求节点分配系数
Ni[0]=(1+b[1])*(1+b[2])*(1+b[3])/8;
Ni[1]=(1+b[1])*(1-b[2])*(1+b[3])/8;
Ni[2]=(1+b[1])*(1-b[2])*(1-b[3])/8;
Ni[3]=(1+b[1])*(1+b[2])*(1-b[3])/8;
Ni[4]=(1-b[1])*(1+b[2])*(1+b[3])/8;
Ni[5]=(1-b[1])*(1-b[2])*(1+b[3])/8;
Ni[6]=(1-b[1])*(1-b[2])*(1-b[3])/8;
Ni[7]=(1-b[1])*(1+b[2])*(1-b[3])/8;
//*************************
cout<<" "<<b[1]<<" "<<b[2]<<" "<<b[3]<<endl;
for(int i=0;i<8;i++)
{cout<<"   "<<Ni[i]<<endl;
//max+=Ni[i];Ni[i]分配系数的和是1
}
cout<<"sfsdf   "<<max<<endl;
}
};