36阶的矩阵和向量，形式为Ax=b
原始的数据：
这是左端的系数矩阵A
View Code
4 -2 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 4 0 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 0 4 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -2 4
右端的系数矩阵b
View Code
-0.57548 0.85828 -0.15069 -0.9514 -0.43732 1.9378 2.4544 2.076 -0.36449 -2.3013 -1.0578 2.4544 -5.3396 -4.5164 0.79296 5.0068 2.3013 -5.3396 0.84572 0.71536 -0.12559 -0.79296 -0.36449 0.84572 4.8168 4.0744 -0.71536 -4.5164 -2.076 4.8168 0.57548 -0.85828 0.15069 0.9514 0.43732 -1.9378
此时的矩阵的条件数为 cond(A) = 1.2630e+016 非常的大
如果将矩阵的第一行和第一列去掉，那么我们的数据变为
左端系数矩阵A
View Code
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 4 0 0 0 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 4 -2 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 4 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -2 4 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 0 4 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2 0 0 0 0 -2 4
此时矩阵的条件数变为 370.4672 （就改变这么一点，发生如此大的改变，太奇妙了）
右端的向量b暂时保持不变（仅仅为了验证是否可解，严格意义上b也相应的改变，待续）
通过选全选主元高斯消去法，参考书籍：徐士良 编著 《常用算法程序集(C语言描述)》(第三版) 清华大学出版社
主程序（6GAUS00.C）为
View Code
#include "stdio.h" #include "6gaus.c" main() { FILE *fptr1, *fptr2,*fptr3; int i,j; double a[36][36]; double b[36]; fptr1 = fopen("input1.txt", "r"); for (i=0; i<36; i++) { for (j=0; j<36; j++) { fscanf(fptr1,"%lf", &a[i][j]); } fscanf(fptr1,"\n"); } fclose(fptr1); fptr2 = fopen("input2.txt", "r"); for (i=0; i<36; i++) { fscanf(fptr2,"%lf", &b[i]); } fclose(fptr2); fptr3 = fopen("output.txt", "w"); if (gaus(a,b,36)!=0) for (i=0;i<=35;i++) { printf("x(%d)=%e\n",i,b[i]); fprintf(fptr3,"%18.5e",b[i]); fprintf(fptr3,"\n"); } fclose(fptr3); }
函数接口（6GAUS.C）为
View Code
#include "stdlib.h" #include "math.h" #include "stdio.h" int gaus(a,b,n) int n; double a[],b[]; { int *js,l,k,i,j,is,p,q; double d,t; js=malloc(n*sizeof(int)); 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) { free(js); printf("fail\n"); 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) { free(js); printf("fail\n"); 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;} free(js); return(1); }
求解结果为
View Code
-5.75480e-001 1.02260e+000 7.62579e-001 4.33064e-001 8.31612e-001 1.63594e+000 1.09184e+000 1.23477e+000 8.72671e-001 5.44733e-001 8.47381e-001 1.47137e+000 -5.56563e-001 -1.24047e-001 1.31310e+000 2.32712e+000 1.59960e+000 1.00393e-001 2.26960e+000 2.02891e+000 1.38369e+000 8.44233e-001 8.22222e-001 1.07059e+000 4.73143e+000 3.87104e+000 1.47411e+000 -3.63137e-001 1.38953e-001 1.69180e+000 4.09724e+000 3.17531e+000 1.72019e+000 6.06561e-001 4.80925e-001 6.01914e-001
上述问题的编译环境全部是在VC6.0下测试过的 。
如果要在vs2008下通过 需要改动一下
创建一个项目，命名为dongtaishuru
在源文件下 添加 6GAUS00.C 6GAUS.C
在资源文件下 添加 input1.txt input2.txt
其组成如图所示：
其6GAUS00.C文件 代码如下
View Code
#include <stdio.h> //#include <fstream> // #include "6gaus.c" int main() { FILE *fptr1, *fptr2,*fptr3; int i,j; double a[36][36]; double b[36]; //ifstream inAMatrix; //inAMatrix.open("F:\\fluid\\test_poissonEqution\\AMatrix.txt"); fptr1 = fopen("input1.txt", "r"); for (i=0; i<36; i++) { for (j=0; j<36; j++) { fscanf(fptr1,"%lf", &a[i][j]); } fscanf(fptr1,"\n"); } fclose(fptr1); fptr2 = fopen("input2.txt", "r"); for (i=0; i<36; i++) { fscanf(fptr2,"%lf", &b[i]); } fclose(fptr2); fptr3 = fopen("output.txt", "w"); if (gaus(a,b,36)!=0) for (i=0;i<=35;i++) { printf("x(%d)=%e\n",i,b[i]); fprintf(fptr3,"%18.5e",b[i]); fprintf(fptr3,"\n"); } fclose(fptr3); //inAMatrix.close(); return 0; }
其他不变