Jacobi迭代法求解方程組
阿新 • • 發佈：2019-02-18
JacobiFile.h
#include <iostream> #include <stdlib.h> #include <malloc.h> #include <math.h> #include "Jacobi.h" //////////////////////////////////////// char matAb[20]; //存放矩陣mat_A和mat_b的檔案 char matRes[20]; //存放每次迭代結果 /////////////////////////////////////////////// JACOBI::JACOBI() { N=0; eps=0.0; mat_A=NULL; mat_D=NULL; mat_D_iv=NULL; result=NULL; } ///////////////////////////////////////////////////// void JACOBI::red_matA(int n) { //讀取矩陣mat_A和mat_b N = n; std::cout << "方陣的階數為 " << N << '\n'; std::cout << "read the mat_A" << '\n'; mat_A = (double **)malloc(N*sizeof(double)); for(int i=0; i<N; i++) *(mat_A+i) = (double *)malloc(N*sizeof(double)); mat_b = (double *)malloc(N*sizeof(double)); std::cout << "input filrname of matAb\n"; std::cin >> matAb; FILE *fp; char ch; if((fp=fopen(matAb,"rt"))==NULL) { //covariance檔案中每一行代表一個樣本向量,列數為維數 printf("Cannot open file strike any key exit!"); exit(0); } ch=fgetc(fp); std::string str=""; //這裡str是類string的一個物件 while(ch!=EOF) { for(i=0; i<N; i++) { while(ch==' ' || ch=='v' || ch=='\n') ch=fgetc(fp); //讀取mat_A for(int j=0; j<N; j++) { while(ch==' ') ch=fgetc(fp); while(ch!=' ') { str+=ch; ch=fgetc(fp); } double aa=atof(str.c_str()); //將字串轉換成double型 *(*(mat_A+i)+j) = aa; //讀取*(*(mat_A+i)+j) str=""; } //讀取mat_A while(ch==' ') ch=fgetc(fp); while(ch!=' ' && ch!='\n') { str+=ch; ch=fgetc(fp); } double aa=atof(str.c_str()); //將字串轉換成double型 *(mat_b+i) = aa; //讀取mat_b+i str=""; //ch=fgetc(fp); } ch=fgetc(fp); } std::cout << "read all the vectors" << '\n'; fclose(fp); } ///////////////////////////////////////////////////////////////////// void JACOBI::creatLUD() { //建立D和L+U mat_D = (double **)malloc(N*sizeof(double)); for(int i=0; i<N; i++) *(mat_D+i) = (double *)malloc(N*sizeof(double)); /*----建立D-----*/ for(i=0; i<N; i++) { for(int j=0; j<N; j++) if(j==i) *(*(mat_D+i)+j) = *(*(mat_A+i)+j); else *(*(mat_D+i)+j) = 0; } /*----建立D-----*/ /*----建立L+U-----*/ for(i=0; i<N; i++) *(*(mat_A+i)+i) =0.0; /*----建立L+U-----*/ //L U建立結束 std::cout << "output mat_D \n"; for(i=0; i<N; i++) { for(int j=0; j<N; j++) std::cout << *(*(mat_D+i)+j) << " "; std::cout << '\n'; } std::cout << "output L+U \n"; for(i=0; i<N; i++) { for(int j=0; j<N; j++) std::cout << *(*(mat_A+i)+j) << " "; std::cout << '\n'; } } ////////////////////////////////////////////////////// void JACOBI::creatD_iv() { //建立mat_D的逆矩陣 /*-----初始化擴充套件矩陣mat_DI(mat_D和單位矩陣I組成)----*/ double **mat_DI = (double **)malloc(N*sizeof(double)); for(int i=0; i<N; i++) *(mat_DI+i) = (double *)malloc(N*sizeof(double)); std::cout << "擴充套件矩陣為" << '\n'; for(i=0; i<N; i++) { for(int j=0; j<N; j++) *(*(mat_DI+i)+j) = *(*(mat_D+i)+j); for(j=N; j<2*N; j++) { if((j-N)==i) *(*(mat_DI+i)+j) = 1; else *(*(mat_DI+i)+j) = 0; } } for(i=0; i<N; i++) { for(int j=0; j<2*N; j++) std::cout << *(*(mat_DI+i)+j) << " "; std::cout << '\n'; } /*-----初始化擴充套件矩陣mat_DI(mat_D和單位矩陣I組成)----*/ /*****************求逆模組***********************/ for(i=0;i<N;i++) { if(*(*(mat_DI+i)+i)==0) { for(int k=i;k<N;k++) { if(*(*(mat_DI+k)+k)!=0) { for(int j=0;j<2*N;j++) { double temp; temp=*(*(mat_DI+i)+j); *(*(mat_DI+i)+j)=*(*(mat_DI+k)+j); *(*(mat_DI+k)+j)=temp; } break; } } if(k==N) { std::cout<<"該矩陣不可逆！"<< '\n'; } } for(int j=2*N-1;j>=i;j--) { *(*(mat_DI+i)+j)/=*(*(mat_DI+i)+i); } for(int k=0;k<N;k++) { if(k!=i) { double temp=*(*(mat_DI+k)+i); for(int j=0;j<N;j++) { *(*(mat_DI+k)+j)-=temp*(*(*(mat_DI+i)+j)); } } } } /*****************求逆模組***********************/ //儲存逆矩陣 mat_D_iv = (double **)malloc(N*sizeof(double)); for(i=0; i<N; i++) *(mat_D_iv+i) = (double *)malloc(N*sizeof(double)); for(i=0;i<N;i++) for(int j=N;j<2*N;j++) *(*(mat_D_iv+i)+j-N) = *(*(mat_DI+i)+j); std::cout << "mat_D的逆矩陣為" << '\n'; for(i=0; i<N; i++) { for(int j=0; j<N; j++) std::cout << *(*(mat_D_iv+i)+j) << " "; std::cout << '\n'; } } //////////////////////////////////////////////////////// void JACOBI::cal(double _eps) { //求解方程組 eps=_eps; double *temp = (double *)malloc(N*sizeof(double)); //迭代求解時用於儲存臨時值的陣列 for(int i=0; i<N; i++) *(temp+i) = 0; //初始化temp作為迭代起始值 //迭代公式是x^(k+1)=(-mat_D_iv*(mat_L+mat_U))x^k+mat_D_iv*mat_b /*---求(-mat_D_iv*(mat_L+mat_U))---*/ double **mat_LU = (double **)malloc(N*sizeof(double)); for(i=0; i<N; i++) *(mat_LU+i) = (double *)malloc(N*sizeof(double)); for(i=0; i<N; i++) for(int j=0; j<N; j++) { double sum=0.0; for(int k=0; k<N; k++) sum += *(*(mat_D_iv+i)+k) * (*(*(mat_A+k)+j)); *(*(mat_LU+i)+j) = -1*sum; } std::cout << "\nthe matrix D^-1*(L+U)\n"; for(i=0; i<N; i++) { for(int j=0; j<N; j++) std::cout << *(*(mat_LU+i)+j) << " "; std::cout << '\n'; } /*---求(-mat_D_iv*(mat_L+mat_U))---*/ /*---求mat_D_iv*mat_b---*/ double *mat_Db = (double *)malloc(N*sizeof(double)); for(i=0; i<N; i++) { double sum=0.0; for(int j=0; j<N; j++) sum += *(*(mat_D_iv+i)+j) * (*(mat_b+j)); *(mat_Db+i) = sum; } std::cout << "\nthe matrix D^-1*b\n"; for(i=0; i<N; i++) { std::cout << *(mat_Db+i) << '\n'; std::cout << '\n'; } /*---求mat_D_iv*mat_b---*/ /*---迭代開始---*/ std::cout << "output the result\n"; std::cout << "input the fileName to save the results\n"; std::cin >> matRes; FILE *fp; if((fp=fopen(matRes,"w"))==NULL) { printf("Cannot build file strike any key exit!"); exit(0); } result = (double *)malloc(N*sizeof(double)); double sum_pre = 0.0; //計算x^k的和 double sum_next = 0.0; //計算x^(k+1)的和 double flag = 0; do { for(int i=0; i<N; i++) { double sum=0.0; for(int j=0; j<N; j++) sum += *(*(mat_LU+i)+j) * (*(temp+j)); sum += *(mat_Db+i); *(result+i) = sum; } sum_pre = 0.0; //計算x^k的和 sum_next = 0.0; for(i=0; i<N; i++) { sum_pre += *(temp+i); sum_next += *(result+i); } //將迭代結果輸入文字檔案 fprintf(fp,"%s"," v "); for(i=0; i<N; i++) fprintf(fp,"%.4f%s",*(temp+i)," "); fprintf(fp,"%s","\n"); fprintf(fp,"%s"," v "); for(i=0; i<N; i++) fprintf(fp,"%.4f%s",*(result+i)," "); fprintf(fp,"%s","\n"); for(i=0; i<N; i++) { *(temp+i) = *(result+i); *(result+i) = 0; } flag=fabs(sum_pre-sum_next); }while(flag > eps); /*---迭代結束---*/ std::cout << "output the result \n"; for(i=0; i<N; i++) std::cout << *(temp+i)<< '\n'; fclose(fp); std::cout << "write all the vectors" << '\n'; }