Jacobi迭代法求解非齐次线性方程组Ax=b
消元法这种直接解法适合低维方程组求解, 一般情况高维方程组的求解用迭代法
Jacobi迭代法求解Ax=b中
Jacobi.h
#include <iostream>
/////////////////////////
class JACOBI {
public :
JACOBI();
public :
void red_matA(int n); ////读取矩阵mat_A和mat_b, n为方阵的阶数
void creatLUD(void); //创建LUD
void creatD_iv(void); //求mat_D的逆矩阵
void cal(double); //求解方程组
private :
int N; //方阵的阶数
double **mat_A, *mat_b; //初始系数矩阵和右值矩阵
double **mat_D; //mat_A的下三角矩阵mat_L,上三角矩阵mat_U,和对角矩阵mat_D
double **mat_D_iv; //mat_D的逆矩阵
double *result; //结果矩阵
double eps; //迭代求解的精度
};
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((2*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<2*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");
for(i=0; i<N; i++) {
*(temp+i) = *(result+i);
*(result+i) = 0;
}
flag=fabs(sum_pre-sum_next);
}while(flag > eps);
/*---迭代结束---*/
fclose(fp);
std::cout << "write all the vectors" << '\n';
}
Jacobi.cpp
#include <iostream>
#include "JacobFile.h"
////////////////////////////////
int main() {
JACOBI jcb;
std::cout << "matAb.txt: 3\n";
int n;
std::cout << "input dimention of matA\n";
std::cin >> n;
jcb.red_matA(n);
jcb.creatLUD();
jcb.creatD_iv();
double esp;
std::cout << "input esp\n";
std::cin >> esp;
jcb.cal(esp);
return 0;
}
验证：
文本文件
结果：
浙公网安备 33010602011771号