jacobi旋轉(zhuǎn)法,Jacobi過關(guān)法

1、實(shí)驗(yàn)名稱：項(xiàng)目一Jacobi旋轉(zhuǎn)法，Jacobi過關(guān)法實(shí)驗(yàn)題目：用Jacobi旋轉(zhuǎn)法，Jacobi過關(guān)法計(jì)算下面矩陣A的全部特征值以及特征向量2-100-12-10A=0-12-1_00-121實(shí)驗(yàn)過程：1、 Jacobi旋轉(zhuǎn)法程序如下:functionEigValMat,EigVecMat=JocobiRot(A,eps)%本程序是根據(jù)jacobi旋轉(zhuǎn)法求實(shí)對稱矩陣的全部特征值和特征向量%輸入變量:A為對稱實(shí)矩陣，eps為允許誤差.%輸出變量：EigValMat為A的特征值,EigVecMat為特征向量n=size(A);n=n(1);P=eye(n);trc=1;num=0;whileab

2、s(trc)>=epsMaxMes=FinMax(A);l=MaxMes(1);s=MaxMes(2);trc=MaxMes(3);%n為矩陣A的階數(shù)%P為旋轉(zhuǎn)矩陣，賦初值為單位陣%trc為矩陣A的非對角元素的平方和，賦初值為1;%num設(shè)置為累加器，記錄迭代的次數(shù)%進(jìn)行正交變換A=PAP'將A的兩個絕對值最大的非對角元素化零，直到所有非對角元素的平方和小于給定的eps則結(jié)束循環(huán).%尋找絕對值最大的非對角元素的位置及所有非對角元素的平方和%l為絕對值最大的非對角元素的行號%s為絕對值最大的非對角元素的列號%trc為矩陣A的非對角元素的平方和RotMat=ComRotMat(A,l

3、,s);%計(jì)算此次旋轉(zhuǎn)的平面旋轉(zhuǎn)矩陣RotMatA=RotMat*A*RotMat'%對當(dāng)前A進(jìn)行一次旋轉(zhuǎn)變換將A的兩個絕對值最大的非對角元素化零，并仍記為AP=RotMat*P;%記錄到目前為止所有旋轉(zhuǎn)矩陣的乘積num=num+1;%記錄已經(jīng)進(jìn)行旋轉(zhuǎn)的次數(shù)endEigValMat=A;EigVecMat=P;numfunctionMaxMes=FinMax(A)%對上三角進(jìn)行掃描，尋找矩陣A的絕對值最大的非對角元素的位置及所有非對角元素的平方和%輸入量：實(shí)對稱矩陣A%輸出量:MaxMes記錄矩陣A的絕對值最大的非對角元素的行號列號及所有非對角元素的平方和n=size(A);n=n(1

4、);trc=0;MaxVal=0;%MaxVal記錄矩陣A的絕對值最大的非對角元素賦初值為0fori=1:n-1forj=i+1:ntrc=trc+2*A(i,j)A2;ifabs(A(i,j)>MaxValMaxVal=abs(A(i,j);l=i;%l為絕對值最大的非對角元素的行號s=j;%s為絕對值最大的非對角元素的列號endendendMaxMes=l,s,trc;functionRotMat=ComRotMat(A,l,s)%計(jì)算平面旋轉(zhuǎn)矩陣%輸入量：實(shí)對稱矩陣A及矩陣A的絕對值最大的非對角元素的行號l列號s%輸出量：旋轉(zhuǎn)的平面旋轉(zhuǎn)矩陣RotMatn=size(A);n=n(1

5、);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)A2);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過關(guān)法程序如下:functionEigValMat,EigVecMat=JocobiGG(A,eps)%本程序是根據(jù)jacobi過關(guān)法求實(shí)對稱矩陣的

6、全部特征值和特征向量%輸入變量:A為對稱實(shí)矩陣，eps為允許誤差.%輸出變量：EigValMat為A的特征值，EigVecMat為特征向量n=size(A);n=n(1);P=eye(n);trc=0;num=0;fori=1:n-1forj=i+1:n%n為矩陣A的階數(shù)%P為旋轉(zhuǎn)矩陣，賦初值為單位陣%trc為矩陣A的非對角元素的平方和，賦初值為0;%num設(shè)置為累加器，記錄迭代的次數(shù)trc=trc+2*A(i,j)A2;end%計(jì)算矩陣A的非對角元素平方和endr0=trcA(1/2);r1=r0/n;whiler1>=epsfori=1:n-1forj=i+1:nifabs(A(i,

7、j)>=r1%對abs(A(i,j)>=r1的元素進(jìn)行旋轉(zhuǎn)變換，將A(i,j)化為0,其余元素視為“過關(guān)”，不作相應(yīng)的變換l=i;s=j;RotMat=ComRotMat(A,l,s);A=RotMat*A*RotMat'P=RotMat*P;num=num+1;end計(jì)算此次旋轉(zhuǎn)的平面旋轉(zhuǎn)矩陣RotMat記錄到目前為止所有旋轉(zhuǎn)矩陣的乘積記錄已經(jīng)進(jìn)行旋轉(zhuǎn)的次數(shù)endendr1=r1/n;endEigValMat=A;EigVecMat=P;num%當(dāng)所有非對角元素絕對值都小于r1后，縮小閥值functionRotMat=ComRotMat(A,l,s)%計(jì)算平面旋轉(zhuǎn)矩陣%輸

8、入量：實(shí)對稱矩陣A及矩陣A的絕對值最大的非對角元素的行號l列號s%輸出量：旋轉(zhuǎn)的平面旋轉(zhuǎn)矩陣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)A2);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;實(shí)驗(yàn)結(jié)果:1、Jacobi旋轉(zhuǎn)法

9、運(yùn)行結(jié)果如下:>>A=2-100;-12-10;0-12-1;00-12;eps=0.0001;EigValMat,EigVecMat=JocobiRot(A,eps)num=EigValMat=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過關(guān)法運(yùn)行結(jié)果如下：>>EigValMat,

10、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.37170.60150.6015-0.6015-0.37170.37170.6015-0.3717-0.3718實(shí)驗(yàn)名稱：項(xiàng)目二Romberg方法求數(shù)值積分實(shí)驗(yàn)題目：用Romberg方法計(jì)算下列積分：3(

11、1) Le'sinxdx,要求準(zhǔn)確到101 2_.-(2)1°e"dx,要求準(zhǔn)確到10實(shí)驗(yàn)過程：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=2A(m-1);h=hI=Ifori=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;whileM>1T(m+1,k+1)=(4Ak*T(m+1,k)-T(m,k)/(4Ak-1);M=M/2;k=k+1;endifabs(T(k,k)-T(k,k-1)<epsbreak;endm=m+1;endI=T(k,k);實(shí)驗(yàn)結(jié)果：>>fu