输入线性方程组的维数,然后随即生成一定有解的线性方程组的增广矩阵,求出解,然后输出时间和方程组的解,以及和标准解的误差(其实就是和标准解的方差)。
1 #include <cstdio> 2 #include <algorithm> 3 #include <cstdlib> 4 #include <iostream> 5 #include <cstring> 6 #include <cmath> 7 #include <ctime> 8 using namespace std; 9 const double eps=1e-9; 10 const int MAX=20000; 11 double ans[MAX]; 12 int n; 13 double **inputdata; double *result; 14 /* 15 输出标准解 16 */ 17 void printresult() { 18 for(int i=0;i<n;i++) 19 printf("%.4lf ",result[i]); 20 printf("\n"); 21 } 22 /* 23 打印出生成的的增广矩阵 24 */ 25 void printb() { 26 for(int i=0;i<n;i++) 27 printf("%.4lf ",inputdata[i][n]); 28 printf("\n"); 29 } 30 /* 31 随机生成增广矩阵 32 */ 33 void pro_inputdata() 34 { 35 int i,j,k; 36 double res=0; 37 printf("please input the number of the element:\n"); 38 scanf("%d",&n); 39 srand(time(0)); 40 result=(double*)malloc(n*sizeof(double)); 41 inputdata=(double**)malloc(n*sizeof(double)); 42 for(i=0;i<n;i++) 43 result[i]=(double)(rand()%10000)/1000.0; 44 for(i=0;i<n;i++) { 45 inputdata[i]=(double*)malloc((n+1)*sizeof(double)); 46 res=0; 47 for(j=0;j<n;j++) { 48 k=(rand()%2==0)?1:(-1); 49 inputdata[i][j]=k*(double)(rand()%10000)/1000.0; 50 res=res+inputdata[i][j]*result[j]; 51 } 52 inputdata[i][n]=res; 53 } 54 } 55 56 /* 57 高斯消元函数 58 */ 59 inline void solve(double **a, double ans[], int n) { 60 for (int i = 0; i < n; ++i) { 61 for (int j = i; j < n; ++j) 62 if (fabs(a[j][i]) > eps) { 63 for (int k = i; k <= n; ++k) swap(a[j][k], a[i][k]); 64 break; 65 } 66 for (int j = i+1; j < n; ++j) { 67 if (fabs(a[j][i]) > eps) { 68 double tmp = a[j][i]/a[i][i]; 69 for (int k = i; k <= n; ++k) a[j][k] -= tmp * a[i][k]; 70 } 71 } 72 } 73 for (int i = n - 1; i >= 0; --i) { 74 double sum = 0.0; 75 for (int j = n - 1; j > i; --j) 76 sum += ans[j] * a[i][j]; 77 ans[i] = (a[i][n] - sum) / a[i][i]; 78 } 79 for (int i = 0; i < n; ++i) printf("%.4lf ", ans[i]); 80 printf("\n"); 81 return; 82 } 83 int main(void) { 84 freopen("in.txt", "r", stdin); 85 clock_t begin, end; 86 begin = clock(); 87 pro_inputdata(); 88 solve(inputdata, ans, n); 89 end = clock(); 90 printf("Time is: %lf seconds\n",(double)(end-begin)/CLOCKS_PER_SEC); 91 double sum=0; 92 for (int i = 0; i < n; ++i) { 93 sum += (ans[i]-result[i])*(ans[i]-result[i]); 94 } 95 printf("Error is %.4lf\n",sqrt(sum)); 96 97 return 0; 98 }
囧,等下在贴代码。==
刚才系统剪切板里面的内容不能粘贴到浏览器里面了,注销了一下就可以了=_=