高斯方程解线性方程组函数,c++版的是正确的, 可是改成delphi语言同一个方程解得结果就变了这是为什么???
int gauss(double* a,double *b,double *x,int n)//  aij*xi=bi  n是方程的维数
{
int *js,l,k,i,j,is,p,q;
    float d,t;
    js=new int[40];                // 为js分配n个整形空间;
    l=1;
    for(k=0;k<=n-2;k++)       /* 进行第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;}   /* js[k]代表主元所在
                                     列,is代表主元所在行 */
}
        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; cout<<"fail"<<endl;
            return(0);
}
        d=a[k*n+k];
        for(j=k+1;j<=n-1;j++)       /*第k行除以系数akk*/
{
p=k*n+j;a[p]=a[p]/d;}
            b[k]=b[k]/d;
            for(i=k+1;i<=n-1;i++)        /*计算aij=aij-aik*akj*/
{
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];
}
}//end of k for!
    d=a[(n-1)*n+n-1];
    if(fabs(d)+1.0==1.0)            /*判断方程组是否有奇异解*/
{
delete js;  cout<<"fail"<<endl;
        return(0);
}
    x[n-1]=b[n-1]/d;
    for(i=n-2; i>=0; i--)                   /*回代过程,求出x[k] */
{
t=0.0;
        for (j=i+1; j<=n-1; j++)
        t=t+a[i*n+j]*x[j];
        x[i]=b[i]-t;
}    js[n-1]=n-1;
    for (k=n-1; k>=0; k--)       /* 如果交换了列,就必须换回来 */
       if (js[k]!=k)
       {
   t=x[k]; x[k]=x[js[k]]; x[js[k]]=t;
   }
    delete js;
    return(1);}
delphi版的, 一模一样, 就只是把表达方式变了,解就不同了……   
Gauss(var a,b,x:Array of real;Const n:integer);    //高斯消元法解方程组     问题问题问题?????
var                                                                 //a[]  左值            ???、/、
i,j,l,k,iss,p,q:integer;                                            //b[]右值             ???/
js: array of integer;                                               //x[]解
d,t:Real;
begin
  Setlength(js,n);
  l:=1;
  for i:=0 to n-1 do   js[i]:=0;   //初试化数据
  for k:=0 to n-2 do     // 进行第K次消元
    begin
      d:=0.0;
      for i:=k to n-1 do
        for j:=k to n-1 do
          begin
            t:=abs(a[i*n+j]);
            if t>d then
              begin
                d:=t;js[k]:=j;iss:=i;
              end;   // js[k]代表主元所在
          end;                         //列,iss代表主元所在行
      if(d+1.0=1.0)  then  l:=0
        else
          begin
            if js[k]<>k  then
              for i:=0  to n-1 do        //交换列
                begin
                  p:=i*n+k;q:=i*n+js[k];
                  t:=a[p];a[p]:=a[q];a[q]:=t;
                end;
            if iss<>k  then
              begin
                for j:=k to n-1 do         //交换行
          begin
                    p:=k*n+j;q:=iss*n+j;
                    t:=a[p];a[p]:=a[q];a[q]:=t;
                  end;
                t:=b[k];b[k]:=b[iss];b[iss]:=t;
              end ;
          end; {else}
     if  l=0  then
       begin
         showmessage('方程无解');
         exit;
       end;
     d:=a[k*n+k];
     for j:=k+1 to n-1 do       //第k行除以系数akk
       begin
         p:=k*n+j;a[p]:=a[p]/d;
       end;
     b[k]:=b[k]/d;
     for i:=k+1 to n-1 do        //计算aij=aij-aik*akj
       begin
         for j:=k+1 to n-1 do
           begin
             p:=i*n+j;
             a[p]:=a[p]-a[i*n+k]*a[k*n+j];
           end;
         b[i]:=b[i]-a[i*n+k]*b[k];
       end;
   end;//end of k for!
   d:=a[(n-1)*n+n-1];
   if (abs(d)+1.0)=1.0  then           //判断方程组是否有奇异解
     begin
       showmessage('方程无解');
       exit;
     end;
   x[n-1]:=b[n-1]/d;
   for i:=n-2  downto 0 do                   //回代过程,求出x[k]
     begin
       t:=0.0;
       for j:=i+1 to n-1 do  t:=t+a[i*n+j]*x[j];
       x[i]:=b[i]-t;
     end;
    js[n-1]:=n-1;
    for  k:=n-1 downto 0 do       // 如果交换了列,就必须换回来 
      if js[k]<>k  then
        begin
          t:=x[k]; x[k]:=x[js[k]]; x[js[k]]:=t;
        end;
   Finalize(js);
end;

解决方案 »

  1. 不对啊,我把Real换成了double类型结果也不一样啊, 我查了,两个编译器的double型都是15~16有效位,1.7×10308
      

  2. 呵呵, 也不是哦, 刚才我把实型数据全部替换成了double两个解就是不同, 郁闷啊!    看来要单独写一个小程序来验证了。
      

  3. 可能是这两句:
      if (d+1.0=1.0) then 
      if (abs(d)+1.0)=1.0 then
    这里的浮点数判断出现了问题。这两处都改成这样试试: if abs(d) < 1E-9 then
      

aliyun

类似问题 »