Jacobi迭代法求解方程组

最新推荐文章于 2021-06-07 09:36:06 发布
shiguangrenran1
最新推荐文章于 2021-06-07 09:36:06 发布
阅读量1.1k
点赞数 1
分类专栏： C/C++
本文链接：https://blog.csdn.net/shiguangrenran1/article/details/8159056
版权
C/C++ 专栏收录该内容
8 篇文章 0 订阅
订阅专栏
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';	                                                                                                                                                                                                                                                                                                                                                                    
}