病态方程求根

    技术2022-05-19  20

    #include<stdio.h>#include<math.h>#include<stdlib.h>

    //连乘表达式double function(double x){ int i; double s=1; for(i=1;i<21;i++)  s*=(x-i); return(s);}

    double diff(double x){ int i; double s=0; for(i=1;i<21;i++)  s+=function(x)/(x-i); return(s);}

    int qr(long double *a,int n,long double *u,long double *v,long double eps,int itmax){  int i,j,k,ii,jj,kk; long double x,y,p,q,r;                          long double q00,q01,q02,q11,q12,q22;            /* 用于求Q*/ int is1,is2,n1;                             long double *a1;                                 long double b,c,s;                              if(itmax == 0)                            /* 已经不能再迭代*/ {   printf("fail/n");  return(0); } if(n==1)                                   /* 矩阵是1阶*/ {  u[0] = a[0];  v[0] = 0.0;  return(1); } if(n==2)                                  {   b = (a[0]+a[3]);  c = a[0]*a[3]-a[1]*a[2];  s = b*b-4.0*c;  y = sqrt(fabs(s));   if (s > 0.0)                             {    if(b > 0.0)    u[0] = (b+y)/2.0;   else    u[0]=(b-y)/2.0;     v[0]=0.0;    u[1]=c/u[0]; v[1]=0.0;  }  else                                     {    u[0] = b/2.0; v[0] = y/2.0;    u[1] = u[0]; v[1] = -v[0];     }  return(1); } is1 = 0; is2 = 0; while(is2 < n-1)                      {  is2++;  j = is2*n+is2;  if(fabs(a[j-1]) < eps*(fabs(a[j-n-1])+fabs(a[j])))   {   n1 = is2-is1;   a1 = (long double*)malloc(n1*n1*sizeof(long double));   for(i=0; i<n1; i++)    for(j=0; j<n1; j++)     a1[i*n1+j] = a[(i+is1)*n+j+is1];       qr(a1,n1,u+is1,v+is1,eps,itmax);       free(a1);    is1 = is2;  } } if(is1>0 && is1<n)                            {  n1=n-is1;  a1=(long double*)malloc(n1*n1*sizeof(long double));  for(i=0; i<n1; i++)   for(j=0; j<n1; j++)    a1[i*n1+j] = a[(i+is1)*n+j+is1];   qr(a1,n1,u+is1,v+is1,eps,itmax);        free(a1);   return(1); } else if(is1  == n)  return(1); for (k=0; k<n-1; k++) {  if(k==0)                                   {   x = a[(n-2)*n+n-2]+a[n*n-1];   y = a[(n-2)*n+n-2]*a[n*n-1]-a[(n-2)*n+n-1]*a[(n-1)*n+n-2];   p = a[0]*(a[0]-x)+a[1]*a[n]+y;   q = a[n]*(a[0]+a[n+1]-x);   r = a[n]*a[2*n+1];  }  else                       {    p = a[k*n+k-1];    q = a[(k+1)*n+k-1];   if(k != n-2)     r= a[(k+2)*n+k-1];   else    r= 0.0;  }  if ((fabs(p)+fabs(q)+fabs(r))!=0.0)      {   if(p<0.0)                              s = -sqrt(p*p+q*q+r*r);   else    s = sqrt(p*p+q*q+r*r);   if(k!=0)     a[k*n+k-1]=-s;   q00 = -p/s;                        q01 = -q/s;    q02 = -r/s;    q11 = -q00-q02*r/(p+s);   q12 = q01*r/(p+s);   q22 = -q00-q01*q/(p+s);   for(j=k; j<n; j++)             {    ii = k*n+j;     jj = (k+1)*n+j;    kk = (k+2)*n+j;    p = q00*a[ii]+q01*a[jj];    q = q01*a[ii]+q11*a[jj];    r = q02*a[ii]+q12*a[jj];    if (k!=n-2)                        {     p = p+q02*a[kk];     q = q+q12*a[kk];     r = r+q22*a[kk];      a[kk] = r;    }    a[ii] = p;    a[jj] = q;    }      j=k+3;                          if(j>=n-1)     j=n-1;   for(i=0; i<=j; i++)             {     ii = i*n+k;     jj = i*n+k+1;    p = q00*a[ii]+q01*a[jj];    q = q01*a[ii]+q11*a[jj];    r = q02*a[ii]+q12*a[jj];    if(k!=n-2)                   {      kk = i*n+k+2;     p = p+q02*a[kk];     q = q+q12*a[kk];     r = r+q22*a[kk];      a[kk]=r;    }    a[ii]=p;    a[jj]=q;    }  }  if(k > 0)                        {   a[(k+1)*n+k-1] = 0.0;   if(k != n-2)    a[(k+2)*n+k-1] = 0.0;  } }                              i = qr(a,n,u,v,eps,itmax-1);     return(i);}

    void main(){ int i,j,n=20; double x0,x[20]; long double eps=1.0e-12;      //eps 精度要求,用于判断元素是否为0 int itmax=50000;        itmax 最大迭代次数 long double u[20],v[20];      u   保存特征值的实部.             long double *H; long double b[21]={1,-2.09999999880791e2,20615,-1.256850e6,5.3327946e7,  -1.672280820e9,4.0171771630e10,-7.56111184500e11,1.1310276995381e13,  -1.35585182899530e14,1.307535010540400e15,-1.0142299865511500e16,6.3030812099294900e16,  -3.11333643161391000e17,1.206647803780370000e18,-3.599979517947610000e18,8.037811822645050000e18,  -1.2870931245151000000e19,1.3803759753640700000e19,-8.752948036761600000e18,2.432902008176640000e18};    for(i=0;i<20;i++) {  x[i]=i+1.1; } printf("原方程的根如下:/n"); for(j=0;j<20;j++) {  do  {   x0=x[j];   x[j]=x0-function(x0)/diff(x0);  }while(fabs(x[j]-x0)>1e-8);  printf("x[%d]=%lf/n",j,x[j]); } printf("方程的19次系数做微小扰动+pow(2,-23)后的解如下:/n");  H=(long double*)malloc(sizeof(long double)*n*n);    /* 生成的上H矩阵*/ for(i=0;i<n;i++)  H[i]=-1.0*b[i+1]/b[0];                /* 第一行*/ for(i=n;i<n*n;i++)  H[i]=0.0;                               /* 下面为0*/ for(i=1;i<n;i++)  H[i*n+i-1]=1.0;                         /* 次对角线为1*/ qr(H,n,u,v,eps,itmax);    for(j=19;j>=0;j--)       {  if(v[j]>0)   printf("x[%d]=%.12lf+%.12lfi/n",20-j,u[j],v[j]);    else if(v[j]<0)   printf("x[%d]=%.12lf%.12lfi/n",20-j,u[j],v[j]);  else   printf("x[%d]=%.12lf/n",20-j,u[j]); }}

     

     

     

     


    最新回复(0)