jacobi旋转法-Jacobi过关法

实验名称：项目一Jacobi旋转法，Jacobi过关法实验题目：用Jacobi旋转法，Jacobi过关法计算下面矩阵A的全部特征值以及特征向量2100121001210012A实验过程：1、Jacobi旋转法程序如下：function[EigValMat,EigVecMat]=JocobiRot(A,eps)%本程序是根据jacobi旋转法求实对称矩阵的全部特征值和特征向量%输入变量:A为对称实矩阵,eps为允许误差.%输出变量:EigValMat为A的特征值,EigVecMat为特征向量n=size(A);n=n(1);%n为矩阵A的阶数P=eye(n);%P为旋转矩阵,赋初值为单位阵trc=1;%trc为矩阵A的非对角元素的平方和,赋初值为1;num=0;%num设置为累加器,记录迭代的次数whileabs(trc)=eps%进行正交变换A=PAP'将A的两个绝对值最大的非对角元素化零,直到所有非对角元素的平方和小于给定的eps,则结束循环.MaxMes=FinMax(A);%寻找绝对值最大的非对角元素的位置及所有非对角元素的平方和l=MaxMes(1);%l为绝对值最大的非对角元素的行号s=MaxMes(2);%s为绝对值最大的非对角元素的列号trc=MaxMes(3);%trc为矩阵A的非对角元素的平方和RotMat=ComRotMat(A,l,s);%计算此次旋转的平面旋转矩阵RotMatA=RotMat*A*RotMat';%对当前A进行一次旋转变换将A的两个绝对值最大的非对角元素化零,并仍记为AP=RotMat*P;%记录到目前为止所有旋转矩阵的乘积num=num+1;%记录已经进行旋转的次数endEigValMat=A;EigVecMat=P;numfunctionMaxMes=FinMax(A)%对上三角进行扫描，寻找矩阵A的绝对值最大的非对角元素的位置及所有非对角元素的平方和%输入量:实对称矩阵A%输出量:MaxMes记录矩阵A的绝对值最大的非对角元素的行号列号及所有非对角元素的平方和n=size(A);n=n(1);trc=0;MaxVal=0;%MaxVal记录矩阵A的绝对值最大的非对角元素赋初值为0fori=1:n-1forj=i+1:ntrc=trc+2*A(i,j)^2;ifabs(A(i,j))MaxValMaxVal=abs(A(i,j));l=i;%l为绝对值最大的非对角元素的行号s=j;%s为绝对值最大的非对角元素的列号endendendMaxMes=[l,s,trc];functionRotMat=ComRotMat(A,l,s)%计算平面旋转矩阵%输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s%输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);ifA(l,l)==A(s,s)c1=1/sqrt(2);s1=1/sqrt(2);elsetan2=2*A(l,s)/(A(l,l)-A(s,s));c2=1/sqrt(1+(tan2)^2);s2=tan2*c2;c1=sqrt((1+c2)/2);s1=s2/(2*c1);endRotMat=eye(n);RotMat(l,l)=c1;RotMat(s,s)=c1;RotMat(l,s)=s1;RotMat(s,l)=-s1;2、Jacobi过关法程序如下：function[EigValMat,EigVecMat]=JocobiGG(A,eps)%本程序是根据jacobi过关法求实对称矩阵的全部特征值和特征向量%输入变量:A为对称实矩阵,eps为允许误差.%输出变量:EigValMat为A的特征值,EigVecMat为特征向量n=size(A);n=n(1);%n为矩阵A的阶数P=eye(n);%P为旋转矩阵,赋初值为单位阵trc=0;%trc为矩阵A的非对角元素的平方和,赋初值为0;num=0;%num设置为累加器,记录迭代的次数fori=1:n-1forj=i+1:ntrc=trc+2*A(i,j)^2;%计算矩阵A的非对角元素平方和endendr0=trc^(1/2);r1=r0/n;whiler1=epsfori=1:n-1forj=i+1:nifabs(A(i,j))=r1%对abs(A(i,j))=r1的元素进行旋转变换，将A(i,j)化为0，其余元素视为“过关”，不作相应的变换l=i;s=j;RotMat=ComRotMat(A,l,s);%计算此次旋转的平面旋转矩阵RotMatA=RotMat*A*RotMat';P=RotMat*P;%记录到目前为止所有旋转矩阵的乘积num=num+1;%记录已经进行旋转的次数endendendr1=r1/n;%当所有非对角元素绝对值都小于r1后，缩小阀值endEigValMat=A;EigVecMat=P;numfunctionRotMat=ComRotMat(A,l,s)%计算平面旋转矩阵%输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s%输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);ifA(l,l)==A(s,s)c1=1/sqrt(2);s1=1/sqrt(2);elsetan2=2*A(l,s)/(A(l,l)-A(s,s));c2=1/sqrt(1+(tan2)^2);s2=tan2*c2;c1=sqrt((1+c2)/2);s1=s2/(2*c1);endRotMat=eye(n);RotMat(l,l)=c1;RotMat(s,s)=c1;RotMat(l,s)=s1;RotMat(s,l)=-s1;实验结果：1、Jacobi旋转法运行结果如下：A=[2-100;-12-10;0-12-1;00-12];eps=0.0001;[EigValMat,EigVecMat]=JocobiRot(A,eps)num=7EigValMat=0.382000-0.000003.6180000-0.00001.382000002.6180EigVecMat=0.37170.60150.60150.3717-0.37170.6015-0.60150.3717-0.6015-0.37170.37170.60150.6015-0.3717-0.37170.60152、Jacobi过关法运行结果如下：[EigValMat,EigVecMat]=JocobiGG(A,eps)num=16EigValMat=0.38200.0000-0.0000-0.00000.00003.6180-0.00000.0000-0.0000-0.00001.38200.00000.00000.00000.00002.6180EigVecMat=0.37170.60150.60150.3717-0.37180.6015-0.60150.3717-0.6015-0.37170.37170.60150.6015-0.3717-0.37180.6015实验名称：项目二Romberg方法求数值积分实验题目：用Romberg方法计算下列积分：（1）30sinxexdx，要求准确到610（2）210xedx，要求准确到610实验过程：functionI=romberg(fun,a,b,eps)m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while1N=2^(m-1);h=h/2;I=I/2;fori=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;whileM1T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;endifabs(T(k,k)-T(k,k-1))epsbreak;endm=m+1;endI=T(k,k);实验结果：fun=inline('exp(x)*sin(x)');a=0;b=3;eps=0.000001;I=romberg(fun,a,b,eps)I=11.8595fun=inline('exp(-x^2)');a=0;b=1;eps=0.000001;I=romberg(fun,a,b,eps)I=0.7468