病态方程求根

来源：互联网发布：明道软件官网编辑：程序博客网时间：2024/04/30 01:43

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