高斯方程解线性方程组函数,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;
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;
if (d+1.0=1.0) then
if (abs(d)+1.0)=1.0 then
这里的浮点数判断出现了问题。这两处都改成这样试试: if abs(d) < 1E-9 then