Physics 版 (精华区)

发信人: zjliu (秋天的萝卜), 信区: Physics
标  题: 一元3次方程 zz
发信站: 哈工大紫丁香 (Thu Mar 11 13:13:08 2004), 站内信件

发信站: 吉林大学牡丹园站


用高思消元法可解非奇异的线性方程组
#include<iostream.h>
#include<math.h>
int agaus(double a[],double b[],int n);
int main()
{
  int i;
  static double a[4][4]=
  {{0.2368,0.2471,0.2567,1.2571}  //系数距阵
  ,{0.1968,0.2071,1.2168,0.2271}
  ,{0.1581,1.1675,0.1768,0.1871}
  ,{1.1161,0.1254,0.1397,0.1490}
  };
  static double b[4]={1.8471,1.7471,1.6471,1.5471};//存解向量
  if (agaus(&a[0][0],b,4)!=0)
   for (i=0;i<=3;i++)
      cout<<"x"<<"("<<i<<")"<<"= "<<b[i]<<endl;
   return 0;
}
int agaus(double *a,double *b,int n)
{
    int *js,l,k,i,j,is,p,q;
    double d,t;

    js =new int[n];
    l=1;
    for (k=0;k<=n-2;k++)
      { d=0.0;
        for (i=k;i<=n-1;i++)
          for (j=k;j<=n-1;j++)
            { t=fabs(a[i*n+j]);
              if (t>d) { d=t; js[k]=j; is=i;}
            }
        if (d+1.0==1.0) l=0;
        else
          { if (js[k]!=k)
              for (i=0;i<=n-1;i++)
                { p=i*n+k; q=i*n+js[k];
                  t=a[p]; a[p]=a[q]; a[q]=t;
                }
            if (is!=k)
              { for (j=k;j<=n-1;j++)
                  { p=k*n+j; q=is*n+j;
                    t=a[p]; a[p]=a[q]; a[q]=t;
                  }
                t=b[k]; b[k]=b[is]; b[is]=t;
              }
          }
        if (l==0)
          {
            delete []js;
            return(0);
          }
        d=a[k*n+k];
        for (j=k+1;j<=n-1;j++)
          { p=k*n+j; a[p]=a[p]/d;}
        b[k]=b[k]/d;
        for (i=k+1;i<=n-1;i++)
          { for (j=k+1;j<=n-1;j++)
              { p=i*n+j;
                a[p]=a[p]-a[i*n+k]*a[k*n+j];
              }
            b[i]=b[i]-a[i*n+k]*b[k];
          }
      }
    d=a[(n-1)*n+n-1];
    if (fabs(d)+1.0==1.0)
      {
        delete []js;
        return(0);
      }
    b[n-1]=b[n-1]/d;
    for (i=n-2;i>=0;i--)
      { t=0.0;
        for (j=i+1;j<=n-1;j++)
          t=t+a[i*n+j]*b[j];
        b[i]=b[i]-t;
      }
    js[n-1]=n-1;
    for (k=n-1;k>=0;k--)
      if (js[k]!=k)
        { t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}

      delete []js;
    return(1);
 }
//用牛顿二分法求某些方程某个区间的一个根。
//2x3-4x2+3x-6=0
#include<iostream.h>
#include<math.h>

double newton(double & ,double &);

int main()
{
  double  low,up;
  cout<<"请输入求根的范围"<<endl;
  cin>>low>>up;
  cout<<newton(low,up)<<endl;;
  return 0;
}

double newton(double &x1, double &x2)
{
 double x3,y1,y2,y3;
 while(fabs(x1-x2)>1e-6)
 {
   y1=2*x1*x1*x1-4*x1*x1+3*x1-6;
   y2=2*x2*x2*x2-4*x2*x2+3*x2-6;
   if(y1*y2<0)
   {
   x3=(x1+x2)/2;
   y3=2*x3*x3*x3-4*x3*x3+3*x3-6;
   if(y1*y3<0)
       x2=x3;
   if(y2*y3<0)
       x1=x3;
   }
 }
 return x3;
}
--
╔═══════════════════╗
║★★★★★友谊第一  比赛第二★★★★★║
╚═══════════════════╝

※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 202.118.229.162]