今天很高兴因为自己调试出了一个程序,调试的好辛苦。但完成时很有成就感。
今天上机实现Crout算法,结合参考书很快就完成了。没有遇到什么大的错误。
#include<iostream.h>
#include<math.h>
#define MAX_N 20
int main()
{
int n;
int i,j,k;
static double a[MAX_N][MAX_N], b[MAX_N], x[MAX_N], y[MAX_N];
static double l[MAX_N][MAX_N], u[MAX_N][MAX_N];
cout<< "Input n value(dim of ax=b):"<<endl;
cin>> n;
if (n>MAX_N)
cout<< "The input n is large then MAX_N,please redefine the MAX_N.";
if (n<=0)
cout << "Please input a number between 1 and"<<MAX_N;
cout << "Now input the matrix a(i,j), i,j=0,...,"<<n-1<<endl;
for (i=0; i<n; i++)
for (j=0; j<n; j++)
cin>> a[i][j]; //输入ax=b的a矩阵
cout<< "Now input the matrix b(i), i=0,...,"<<n-1<<endl;
for (i=0; i<n; i++)
cin>> b[i]; //输入不b矩阵
for (i=0; i<n; i++)
u[i][i]=1;
for (k=0; k<n; k++) //计算L矩阵
{
for (i=k; i<n; i++)
{
l[i][k]=a[i][k];
for (j=0; j<=k-1; j++)
l[i][k]-=(l[i][j]*u[j][k]);
}
for (j=k+1; j<n; j++) //计算U矩阵
{
u[k][j]=a[k][j];
for (i=0; i<=k-1; i++)
u[k][j]-=(l[k][i]*u[i][j]);
u[k][j]/=l[k][k];
}
}
for(i=0; i<n; i++)
{
y[i]=b[i];
for (j=0; j<=i-1; j++)
y[i]-=(l[i][j]*y[j]);
y[i]/=l[i][i];
}
for (i=n-1; i>=0; i--)
{
x[i]=y[i];
for (j=i+1; j<n; j++)
x[i]-=(u[i][j]*x[j]);
}
cout.setf(ios::showpoint );
cout << "Solve... x_i="<<endl;
for (i=0; i<n; i++)
cout <<x[i]<<endl;
return 0;
}