矩阵古典雅克比迭代算法

古典雅克比算法一、实验目的1、学习和掌握古典的jacobi 算法。

2、 运用古典jacobi 算法进行计算。

二、方法原理在正交相似变换下，矩阵元素的平方和保持不变，寻找相似变换，是对称矩阵A 经过变换之后所得矩阵的非对角线元素的平方和减少，对角线元素的平方和增大，且保持对称性不变，不断的施行这种正交变换，最终是非对角元素的平方和任意的接近与零，对角线元素平方的取极大值。

三、算法描述： step1:1)1(---=k qq k pp a a y 1)1()1(2*)sgn(----=kpq k qq k pp a a a x step 2:)/1(2/1cos 222y x y ++=θ θθcos /1*/*2/1sin 22y x x += step 3: 由于k A 得对称性，实际计算时只需对k A 的上三角元素进行即可。

()()(1)(1)cos sin k k k k ip pi ip iq a a a a θθ--==+()()(1)(1)sin cos k k k k iq qi ip iq a a a a θθ--==-+()()(1)(1)(1)22()sin cos (cos sin )k k k k k pq qp qq pp pq a a a a a θθθθ---==-+-()(1)2(1)(1)2cos 2sin cos sin k k k k pp pp pq qq a a a a θθθθ---=++()(1)2(1)(1)2sin 2sin cos cos k k k k qq pp pq qq a a a a θθθθ---=-+四、程序原代码如下：头文件：#include<iostream.h>#include<math.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();Matrix& operator =(Matrix &A);void max(int a[]);void get_cos_sin(double b[],int a[]);void Jacobi(Matrix &X,double b[],int a[]);bool BeZero(int a[]);double SUM();};CPP文件：#include<iostream.h>#include<malloc.h>#include<iomanip.h>#include"Jacobi.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);}}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;}Matrix& Matrix::operator =(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;return *this;}void Matrix::print()////////////////////////打印函数{for(int i=0;i<M;++i){for(int j=0;j<N;++j){cout<<setw(15)<<setprecision(5)<<ar[i][j];}cout<<endl;}}void Matrix::Init()/////////////////////////初始化函数{for(int i=0;i<M;++i){cout<<"input the data of "<<i+1<<" line"<<endl;for(int j=0;j<N;++j){cin>>ar[i][j];}}}void Matrix::max(int a[])////////////////按模最大{int i=0,j=0;Datatype max=0;for(i=0;i<M;++i){for(j=0;j<N;++j){if(i==j){continue;}if(max<ar[i][j]){a[0]=i;a[1]=j;max=ar[i][j];}}}}bool Matrix::BeZero(int a[])/////////////若分母为零{if(ar[a[0]][a[0]]==ar[a[1]][a[1]]){return true;//取值45度}else{return false;//按公式计算}}void Matrix::get_cos_sin(double b[],int a[])/////计算cos sin {double y=fabs(ar[a[0]][a[0]]-ar[a[1]][a[1]]);double x=0;if(ar[a[0]][a[0]]-ar[a[1]][a[1]]>0){x=2*ar[a[0]][a[1]];}else{x=-2*ar[a[0]][a[1]];}b[0]=sqrt((y/sqrt(pow(x,2)+pow(y,2))+1)/2);b[1]=x/(2*sqrt(pow(x,2)+pow(y,2))*b[0]);}void Matrix::Jacobi(Matrix &X,double b[],int a[])/////////进行矩阵迭代得到新矩阵X{X=*(this);for(int i=0;i<M;++i){if(i==a[0]||i==a[1]){continue;}else{X.ar[a[0]][i]=X.ar[i][a[0]]=ar[i][a[0]]*b[0]+ar[i][a[1]]*b[1];X.ar[a[1]][i]=X.ar[i][a[1]]=ar[i][a[1]]*b[0]-ar[i][a[0]]*b[1];}}X.ar[a[0]][a[0]]=ar[a[0]][a[0]]*pow(b[0],2)+2*ar[a[0]][a[1]]*b[0]*b[1]+ar[a[1]][ a[1]]*pow(b[1],2);X.ar[a[1]][a[1]]=ar[a[0]][a[0]]*pow(b[1],2)-2*ar[a[0]][a[1]]*b[0]*b[1]+ar[a[1]][ a[1]]*pow(b[0],2);X.ar[a[0]][a[1]]=X.ar[a[1]][a[0]]=(ar[a[1]][a[1]]-ar[a[0]][a[0]])*b[0]*b[1]+ar[a[0 ]][a[1]]*(pow(b[0],2)-pow(b[1],2));}double Matrix::SUM()/////////////////迭代结束判断{double sum=0;for(int i=0;i<M;++i){for(int j=0;j<N;++j){if(i==j){continue;}else{sum+=pow(ar[i][j],2);}}}return sum;}主程序：#include"Jacobi.h"#define E pow(10,-3)void main(){int row;cout<<"请输入阶数"<<endl;cin>>row;int a[2]={0};double b[2]={0};Matrix A(row,row),X(row,row);cout<<"请输入系数矩阵"<<endl;A.Init();//系数int i=1;while(1){A.max(a);if(A.BeZero(a)){b[0]=b[1]=sqrt(2)/2;}else{A.get_cos_sin(b,a);}A.Jacobi(X,b,a);cout<<"第"<<i<<"次迭代结束后矩阵为："<<endl;X.print();cout<<endl;if(A.SUM()<E){break;}A=X;++i;} } 五、数值计算:102021211A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦六、参考文献：[1]刑志栋，矩阵数值分析， 陕西： 陕西科学技术出版社， 2005。

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

[3]翁惠玉，C 语言程序设计思想与方法，北京：人民邮电出版社，2008.。