矩阵雅克比迭代算法

雅克比迭代
实验目的：
1.学习和掌握线性代数方程组的jacobi 迭代法。

2.运用jacobi 迭代法进行计算。

方法原理：
设方程组Ax=b 的系数矩阵A 非奇异而且),...,2,1(0n i a ii =≠，将A 分裂为 A=D+L+U,可以使计算简便。

其中),,...,,(2211nn a a a diag D =
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=
0...
............0...00 (002)
1
21n n a a a L ，⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡=0...
...............00...02112n
n a a a U A=D+L+U ，其中),,...,,(2211nn a a a diag D =
⎥⎥⎥⎥⎦⎤⎢⎢⎢
⎢⎣⎡=0...
............0...00 (002)
1
21
n n a a a L ，⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡=0...
...............00...02112n
n a a a U 将方程组
n ,...,2,1i ,b x a i n
1j j ij ==∑
=乘以
ii
a 1，得到等价的方程组
⎪⎪⎪
⎭
⎫ ⎝
⎛-
=∑
≠=n
i
j 1
j j ij i ii i x a b a 1x ，i=1，2，…n ，简记为x Bx f
=+。

其中 11()B I D A D L U --=-=-+, 1f D b -=. 我们称
x Bx f ϕ=+为迭代函数。

任取初始向量(0)
x x
=，按照
(1)
(
)
k k x
Bx f
+=+形成迭代格式，称这种迭代方法为Jacobi 迭代法。

算法描述：
Step1：给定一组x ，即初值。

Step2：用for 循环计算：
x[k+1]=(b[i]-∑∑+=-=-
n
1
i j 1
i 1
j ]j [x ]j ][i [a ]j [x ]j ][i [a )/a[i][i].
Step3:当fabs(x[k+1]-x[k])<eps时停止。

程序代码：
头文件：
#include<iostream.h>
typedef double Datatype;
class Matrix
{
private:
Datatype **ar;
int M;
int N;
public:
Matrix(int a=0,int b=0);
Matrix(Matrix &A);
~Matrix();
void print();
void Init();
void Jacobi(Matrix b,Matrix x,Matrix &xx);
bool compare(Matrix X);
void change(Matrix xx);
};
CPP文件：
#include<iostream.h>
#include<malloc.h>
#include<iomanip.h>
#include"Jacobi.h"
#define E 0.0000001
#include<math.h>
Matrix::Matrix(int a,int b)///////////////////构造函数
{
M=a;
N=b;
ar=(Datatype **)malloc(sizeof(Datatype *)*M);
for(int i=0;i<M;++i)
{
ar[i]=(Datatype *)malloc(sizeof(Datatype)*N);
for(int j=0;j<N;++j)
{
ar[i][j]=0;
}
}
}
Matrix::~Matrix()////////////////////////////析构函数
{
for(int i=0;i<M;++i)
{
free(ar[i]);
}
free(ar);
}
Matrix::Matrix(Matrix &A)
{
M=A.M;
N=A.N;
Datatype **p=(Datatype **)malloc(sizeof(Datatype *)*M);
for(int i=0;i<M;++i)
{
p[i]=(Datatype *)malloc(sizeof(Datatype)*N);
for(int j=0;j<N;++j)
{
p[i][j]=A.ar[i][j];
}
}
ar=p;
}
void Matrix::print()////////////////////////打印函数
{
for(int i=0;i<M;++i)
{
for(int j=0;j<N;++j)
{
cout<<setw(5)<<ar[i][j];
}
cout<<endl;
}
}
void Matrix::Init()/////////////////////////初始化函数
{
for(int i=0;i<M;++i)
{
cout<<"请输入第"<<i+1<<"行数据"<<endl;
for(int j=0;j<N;++j)
{
cin>>ar[i][j];
}
}
}
void Matrix::Jacobi(Matrix b,Matrix x,Matrix &xx) {
for(int i=0;i<N;++i)
{
Datatype sum=0;
for(int j=0;j<N;++j)
{
if(i==j)
continue;
sum+=ar[i][j]*x.ar[j][0];
}
xx.ar[i][0]=(b.ar[i][0]-sum)/ar[i][i];
}
}
bool Matrix::compare(Matrix X)
{
int flag=0;
for(int i=0;i<M;++i)
{
for(int j=0;j<N;++j)
{
if(fabs(ar[i][j]-X.ar[i][j])>E)
{
flag=1;
break;
}
}
break;
}
if(flag==1)
return false;//程序需继续迭代
else
return true;//迭代结束
}
void Matrix::change(Matrix xx)
{
for(int i=0;i<M;++i)
{ for(int j=0;j<N;++j) { ar[i][j]=xx.ar[i][j]; } } }
主程序：
#include"Jacobi.h"
void main() { int row; cout<<"*****"<<endl; cin>>row; Matrix A(row,row),b(row,1),x(row,1),xx(row,1); cout<<"请输入系数矩阵"<<endl; A.Init();//系数 cout<<"请输入矩阵b"<<endl; b.Init();//b cout<<"请输入迭代初始值"<<endl; x.Init();//迭代初始值 while(1) { A.Jacobi(b,x,xx); if(pare(xx)) { xx.print(); break; } x.change(xx); } }
测试数据：
1231231
232213225
x x x x x x x x x +-=⎧⎪
++=⎨⎪++=⎩ 测试结果：
参考文献：
[1]刑志栋，矩阵数值分析，陕西：陕西科学技术出版社，2005。

[2]谭浩强，C语言程序设计，北京：清华大学出版社，2005。