北航数值分析1-Jacobi法计算矩阵特征值

准备工作算法设计矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。

分析一：由题目中所给条件λ1≤λ2≤…≤λn，可得出λ1、λn按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现|λs|<λ1|<|λn |或|λs|<λn|<|λ1 |的情况，导致按幂法和反幂法无法求解λ1或λn二者中的一者；分析二：题目要求求解与数μk =λ1+k(λn-λ1)/40最接近的特征值λik(k=1,2,3…39)，这个问题其实可以转换为求A-μk 按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得λ1和λn，所以第二个问题暂先搁浅；分析三：cond(A)2= ||A|| * ||A-1|| =|λ|max *|λ|min,这可以用幂法和反幂法求得，det(A) =λ1 *λ2 * … *λn，这需要求得矩阵A的所有特征值。

由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。

所以该题可以用Jacobi法求解。

模块设计由数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。

但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。

基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。

完整代码如下（编译环境windows10 + visual studio2010）：完整代码// math.cpp : 定义控制台应用程序的入口点。

//#include "stdafx.h"#include<stdio.h>#include<math.h>#include<time.h>#define N 501#define V (N+1)*N/2+1#define e 2.718281828459045235360287471352662497757247093699959574966967627724076630353#define a(i) (1.64 - 0.024 * (i)) * sin(0.2 * (i)) - 0.64 * pow(e , 0.1 / (i))#define b 0.16#define c -0.064#define eps pow((double)10.0,-12)#define PFbits "%10.5f "#define PFrols 5#define PFe %.11e#define FK 39int p;int q;doublecosz;doublesinz;double MAX;intkk;//#define PTS pts#ifdef PTSvoid PTS(double *m){printf("-----------------------------------------------------------------------\n");printf(" 迭代第%d次\n",kk);for(inti = 1 ; i<= PFrols ; i++){for(int j = (i-1)*i/2+1 ; j <= (i+1)*i/2 ; j++){printf(PFbits,m[j]);}putchar(10);}for(inti = 1 ; i<= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(inti = 1 ; i<= PFrols+2 ; i++){printf(" ... ");}putchar(10);}#elsevoid PTS(double *m){}#endifvoidrecounti(inti , int *pp, int *qq){for(int j = 0 ; j <= N-1 ; j++){if( (i - (j+1)*j/2) <= j+1){*pp = j+1;*qq = i - (j+1)*j/2;break;}}}voidrefreshMetrix(double *m){intipr,ipc,iqr,iqc;m[(p+1)*p/2] = m[(p+1)*p/2] * pow(cosz,2) + m[(q+1)*q/2] * pow(sinz,2) + 2 * m[(p-1)*p/2+q] * cosz * sinz;m[(q+1)*q/2] = m[(p+1)*p/2] * pow(sinz,2) + m[(q+1)*q/2] * pow(cosz,2) - 2 * m[(p-1)*p/2+q] * cosz * sinz;for(inti = 1; i<= N ;i++){if(i != p &&i != q){if(i> p){ipr = i;ipc = p;}else{ipr = p;ipc = i;}if(i> q){iqr = i;iqc = q;}else{iqr = q;iqc = i;}m[(ipr-1)*ipr/2+ipc] = m[(ipr-1)*ipr/2+ipc] * cosz + m[(iqr-1)*iqr/2+iqc] * sinz;m[(iqr-1)*iqr/2+iqc] = -m[(ipr-1)*ipr/2+ipc] * sinz + m[(iqr-1)*iqr/2+iqc] * cosz;}}m[(p-1)*p/2+q] = 0;PTS(m);}//voidcalCosSin(double *m){double app = m[(p+1)*p/2];doubleaqq = m[(q+1)*q/2];doubleapq = m[(p-1)*p/2+q];cosz = cos(atan(2 * apq / (app - aqq))/2);sinz = sin(atan(2 * apq / (app - aqq))/2); }//voidfind_pq(double *m){double max = 0.0;int pp = 0;intqq = 0;for(inti = 1 ; i<= V ; i++){if(fabs(m[i]) > max){recounti(i,&pp,&qq);if(pp != qq){max = fabs(m[i]);p = pp;q = qq;}}}MAX = max;}voidinit(double *m){for(inti = 1 ; i<= N ;i++)m[(i+1)*i/2] = a(i);for(inti = 2 ; i<= N ; i++)m[(i-1)*i/2+i-1] = b;for(inti = 3 ; i<= N ; i++)m[(i-1)*i/2+i-2] = c;PTS(m);}voidcalFinal(double *m){printf("---------------------------------------------------------------------------------------------------\n");printf("结果输出:\n\n");doubleconda;doubledeta = 1.0;doubleminlumda = pow((double)10.0,12);doublemaxlumda = pow((double)10.0,-12);doubleabsminlumda = pow((double)10.0,12);for(inti = 1 ; i<=N ;i++){if(m[(i+1)*i/2] >maxlumda)maxlumda = m[(i+1)*i/2];if(m[(i+1)*i/2] <minlumda)minlumda = m[(i+1)*i/2];if(fabs(m[(i+1)*i/2]) <absminlumda)absminlumda = fabs(m[(i+1)*i/2]);deta *= m[(i+1)*i/2];}if(fabs(minlumda) <fabs(maxlumda))conda = fabs(maxlumda) / absminlumda;elseconda = fabs(minlumda) / absminlumda;printf(" Lumda(1)=%.11e Lumda(%d)=%.11e Lumda(s)=%.11e\n",minlumda,N,maxlumda,absminlumda);printf(" Cond(A)=%.11e\n",conda);printf(" Det(A)=%.11e\n\n",deta);for(inti = 1 ; i<= FK ; i++){doublemuk = minlumda + i * (maxlumda - minlumda) / 40;doublelumdak = 0.0;doubletempabsmin = pow((double)10.0,12);for(int j = 1 ; j <= N ;j++){if(fabs(muk - m[(j+1)*j/2]) <tempabsmin){lumdak = m[(j+1)*j/2];tempabsmin = fabs(muk - m[(j+1)*j/2]);}}printf(" Lumda(i%d)=%.11e ",i,lumdak);if(i%3==0)putchar(10);}putchar(10);printf("------------------------------------------------------------------------------------------------------\n");putchar(10);putchar(10);}int _tmain(intargc, _TCHAR* argv[]){double m[(N+1)*N/2+1] = {0.0};kk=0;MAX=1.0;time_t t0,t1;t0 = time(&t0);init(m);#ifndef PTSprintf("正在计算...\n\n");#endifwhile(true){kk++;find_pq(m);if(MAX<eps)break;#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n"); #endifcalCosSin(m);refreshMetrix(m);}#ifdef PTSprintf(" p=%d q=%d |max|=%e\n",p,q,MAX);printf("-----------------------------------------------------------------------\n\n");#endifprintf("矩阵最终形态...\n");for(inti = 1 ; i<= PFrols ; i++){for(int j = (i-1)*i/2+1 ; j <= (i+1)*i/2 ; j++){printf(PFbits,m[j]);}putchar(10);}for(inti = 1 ; i<= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(inti = 1 ; i<= PFrols+2 ; i++){printf(" ... ");}putchar(10);t1 = time(&t1);#ifdef PTSprintf("计算并输出用时%.2f秒\n\n",difftime(t1,t0));#elseprintf("迭代次数%d,计算用时%.2f秒\n\n",kk,difftime(t1,t0));#endifcalFinal(m);return 0; }运行结果如下：中间运行状态如下：结果分析数值分析2015/11/10有效性分析1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2.代码中有定义预编译宏//#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G内存，64位windows10操作系统）。

实验名称实验8实验地点6A-XXX 实验类型设计实验学时 2 实验日期20 /X/X ★撰写注意：版面格式已设置好（不得更改），填入内容即可。

一、实验目的
1.Jacobi法求实对称矩阵的特征值及特征向量
二、实验内容
1.实验任务
1.Jacobi法求实对称矩阵的特征值及特征向量
2.程序设计
1）数据输入（输入哪些数据、个数、类型、来源、输入方式）
double a[N][N], int n
2）数据存储（输入数据在内存中的存储）
函数
void Jacobi(double a[N][N], int n)
3）数据处理（说明处理步骤。

若不是非常简单，需要绘制流程图）
1.输入要处理的数据进入变量中
2.进行函数处理
3.输出函数处理结果
4）数据输出（贴图：程序运行结果截图。

图幅大小适当，不能太大）
三、实验环境
1.操作系统：WINDOWS 7及以上
2.开发工具：VS 2015
3.实验设备：PC。

所在平面旋转了一个角度 , 其它坐标保持不 变, 称Upq为平面旋转矩阵.
数值分析
基于Upq的相似变换
用Upq对A作正交相似变换得到矩阵 A(1) ，即：
U
T pq
AU
pq
=
A(1)
aaq((p11qp))
= =
a pp a pp
cos2 + aqq sin2 sin2 + aqq cos2
i , j=1
则有 lki→m k = 0成立.
i j
数值分析
UTAU=D 其中 D＝diag[λ1, λ2, …, λn]，即D的对角元素为 A 的 特征值，对应的U的列向量即为相应的特征向量。
Jacobi方法的思路：通过一系列的旋转变换（正交 变换）把A中非对角线上的非零元变为零 。
数值分析
定义下面的 n 阶正交矩阵：
1
cos
U ( p, q, ) =
素apq； Step 2. 根据 cot 2 = app − a qq ，求出相应的旋转矩
2a pq
阵Upq；
Step 3. 根据公式计算矩阵A(1)的元素；
Step 4. 若
，则停机；否则，令A=A(1)
并返回Step 1.
数值分析
1. 令k＝1，R（1）＝I，给定矩阵A（＝A（1）），收敛条件ε
=
a(1) qp
=
1 2 (aqq
−
app )sin 2
+
a pq
cos 2
矩阵A(a1i()j1的) =第a(j1ip) 行= a，ij ,i第, j p列p, q和第q行，第q列的元素
发生了变化，其余元素不变。 且A(1)任是实对称

A1 仍然是实对称阵，且 A1与 A 的特征值相同。

把式(3.12)的右端乘开，并与左端比较，得到A1 的元素计算公式：
A
(1)
a a (p11) a (1) q1 a (1) n1
(1) 11
a a a a
(1) 1p

4
sgn(a pq ).
2 tan 1 2 1 tan c

4
， 故 tan 1,
故可取 tan
sgn(c ) c c2 1
t

t sin t cos 1 t 2 则 1 1 cos 1 tan 2 1 t 2

i2 ( B) B .
2 i 1
F
n

对于 n 维列向量 x， Upq x 相当于把坐标轴Oxp和 Oxq于所在平面内旋转角度 。 记矩阵 A = [aij]n×n ，对A作一次正交相似变换， 得到矩阵A1。


记
A1= UpqT A Upq = [aij(1)]n×n
(3.12)
故可选择角θ，使
c12 c21 (a22 a11 ) sin cos a12 (cos sin ) 0
2 2
c12 c21 (a22 a11 ) sin cos a12 (cos2 sin 2 ) 0
即 tan 2 2a12 ， a11 a22
1 2 a aij n(n 1) i j
2 pq
从而
a
i j
(1) 2 ij
2 2 1 aij n(n 1) i j

1、用jacobi方法计算对称矩阵A的特征值和对应的特征向量。

function [k,Bk,V,D,Wc]=jacobite(A,jd,max1[n,n]=size(A;P0=eye(n;Vk=eye(n;Bk=A;k=1;state=1;while ((k<=max1&(state==1aij=abs(Bk-diag(diag(Bk;[m1 i]=max(abs(aij;[m2 j]=max(m1;i=i(j;Aij=(Bk-diag(diag(Bk;mk=m2*sign(Aij(i,j;Wc=m2;Dk=diag(diag(Bk;Pk=P0;c=(Bk(j,j-Bk(i,i/(2*Bk(i,j;t=sign(c/(abs(c+sqrt(1+c^2;pii=1/(sqrt(1+t^2;pij=t/(sqrt(1+t^2;Pk(i,i=pii;Pk(i,j=pij;Pk(j,j=pii;Pk(j,i=-pij;B1=Pk'*Bk;B2=B1*Pk;Vk=Vk*Pk;Bk=B2;state=0;if (Wc>jdstate=1;endk,k=k+1;Pk,Vk,Bk=B2,Wc,endif (k>max1disp('迭代次数k已经达到最大迭代次数ma1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下:'elsedisp('迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下:'endk=k-1;Bk=B2;V=Vk;D=diag(diag(Bk;Wc;[V1,D1]=eig(A,'nobalance'>> A=[12 -56 3 -1;-56 7 2 0;3 2 5 1;-1 0 1 12];>> [k,Bk,V,D,Wc]=jacobite(A,0.0001,100k =1Pk =0.7227 0.6912 0 0-0.6912 0.7227 0 00 0 1.0000 00 0 0 1.0000 Vk =0.7227 0.6912 0 0-0.6912 0.7227 0 00 0 1.0000 00 0 0 1.0000 Bk =65.5558 0 0.7858 -0.7227-0.0000 -46.5558 3.5189 -0.69120.7858 3.5189 5.0000 1.0000-0.7227 -0.6912 1.0000 12.0000 Wc = 56k =2Pk =1.0000 0 0 00 0.9977 0.0678 00 -0.0678 0.9977 00 0 0 1.0000 Vk =0.7227 0.6896 0.0468 0-0.6912 0.7210 0.0490 00 -0.0678 0.9977 00 0 0 1.0000 Bk =65.5558 -0.0533 0.7840 -0.7227-0.0533 -46.7948 0 -0.75740.7840 0.0000 5.2391 0.9509-0.7227 -0.7574 0.9509 12.0000 Wc = 3.5189k =3Pk =1.0000 0 0 00 1.0000 0 00 0 0.9906 0.13670 0 -0.1367 0.9906 Vk =0.7227 0.6896 0.0464 0.0064-0.6912 0.7210 0.0485 0.00670 -0.0678 0.9883 0.13640 0 -0.1367 0.9906Bk =65.5558 -0.0533 0.8754 -0.6088-0.0533 -46.7948 0.1035 -0.75020.8754 0.1035 5.1079 -0.0000-0.6088 -0.7502 -0.0000 12.1312 Wc = 0.9509k =4Pk =0.9999 0 -0.0145 00 1.0000 0 00.0145 0 0.9999 00 0 0 1.0000 Vk =0.7233 0.6896 0.0359 0.0064-0.6904 0.7210 0.0585 0.00670.0143 -0.0678 0.9882 0.1364-0.0020 0 -0.1367 0.9906 Bk =65.5685 -0.0518 -0.0000 -0.6087-0.0518 -46.7948 0.1043 -0.7502-0.0000 0.1043 5.0952 0.0088-0.6087 -0.7502 0.0088 12.1312 Wc = 0.8754k =5Pk =1.0000 0 0 00 0.9999 0 -0.01270 0 1.0000 00 0.0127 0 0.9999 Vk =0.7233 0.6896 0.0359 -0.0024-0.6904 0.7211 0.0585 -0.00250.0143 -0.0660 0.9882 0.1372-0.0020 0.0126 -0.1367 0.9905 Bk = 65.5685 -0.0595 -0.0000 -0.6080-0.0595 -46.8044 0.1044 0.0000-0.0000 0.1044 5.0952 0.0075-0.6080 0.0000 0.0075 12.1407 Wc = 0.7502k =6Pk =0.9999 0 0 0.01140 1.0000 0 00 0 1.0000 0-0.0114 0 0 0.9999 Vk =0.7233 0.6896 0.0359 0.0059-0.6903 0.7211 0.0585 -0.01030.0127 -0.0660 0.9882 0.1374-0.0132 0.0126 -0.1367 0.9905 Bk = 65.5754 -0.0595 -0.0001 -0.0000-0.0595 -46.8044 0.1044 -0.0007-0.0001 0.1044 5.0952 0.0075-0.0000 -0.0007 0.0075 12.1338 Wc =0.6080k =7Pk =1.0000 0 0 00 1.0000 0.0020 00 -0.0020 1.0000 00 0 0 1.0000 Vk =0.7233 0.6895 0.0373 0.0059-0.6903 0.7209 0.0600 -0.01030.0127 -0.0680 0.9881 0.1374-0.0132 0.0129 -0.1366 0.9905 Bk = 65.5754 -0.0595 -0.0002 -0.0000-0.0595 -46.8046 -0.0000 -0.0007-0.0002 -0.0000 5.0954 0.0075-0.0000 -0.0007 0.0075 12.1338 Wc = 0.1044k =8Pk =1.0000 0.0005 0 0-0.0005 1.0000 0 00 0 1.0000 00 0 0 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9881 0.1374-0.0133 0.0129 -0.1366 0.9905 Bk = 65.5754 0.0000 -0.0002 0.0000-0.0000 -46.8046 -0.0000 -0.0007-0.0002 -0.0000 5.0954 0.00750.0000 -0.0007 0.0075 12.1338 Wc = 0.0595k =9Pk =1.0000 0 0 00 1.0000 0 00 0 1.0000 0.00110 0 -0.0011 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9880 0.1384-0.0133 0.0129 -0.1377 0.9903 Bk = 65.5754 0.0000 -0.0002 0.0000-0.0000 -46.8046 0.0000 -0.0007-0.0002 0.0000 5.0954 0.00000.0000 -0.0007 -0.0000 12.1338 Wc = 0.0075k =10Pk =1.0000 0 0 00 1.0000 0 -0.00000 0 1.0000 00 0.0000 0 1.0000 Vk =0.7229 0.6899 0.0373 0.0059-0.6907 0.7206 0.0600 -0.01030.0128 -0.0680 0.9880 0.1384-0.0133 0.0129 Bk = 65.5754 0.0000 0.0000 -46.8046 -0.0002 0.0000 0.0000 -0.0000 Wc = 6.9206e-004 k= 11 Pk = 1.0000 0 0 1.0000 -0.0000 0 0 0 Vk = 0.72290.6899 -0.6907 0.7206 0.0128 -0.0680 -0.0133 0.0129 Bk = 65.5754 -0.0000 -0.0000 -46.8046 -0.0000 0.0000 0.0000 -0.0000 Wc = 2.0482e-004 k= 12 Pk = 1.0000 0 0 1.0000 0 -0.0000 0 0 Vk = 0.7229 0.6899 -0.6907 0.7206 0.0128 -0.0680 -0.0133 0.0129 Bk = 65.5754 -0.0000 -0.0000 -46.8046 -0.0000 0.0000 0.0000 -0.0000 -0.1377 -0.00020.0000 5.0954 -0.0000 0.9903 0.0000 0.0000 -0.0000 12.1338 0.0000 0 1.0000 0 0.0373 0.0600 0.9880 -0.1377 0.0000 0.0000 5.0954 -0.0000 0 0 0 1.0000 0.0059 -0.0103 0.1384 0.9903 0.0000 0.0000 -0.0000 12.1338 0 0.0000 1.0000 0 0.0373 0.0600 0.9880 -0.1377 0.0000 -0.0000 5.0954 -0.0000 0 0 0 1.0000 0.0059 -0.0103 0.1384 0.9903 0.0000 0.0000 -0.0000 12.1338Wc = 6.2740e-007 迭代次数 k，对称矩阵 Bk，以特征向量为列向量的矩阵 V，特征值为对角元的对角矩阵 D 如下： V1 = 0.6899 -0.0373 0.0059 -0.7229 0.7206 -0.0600 -0.0103 0.6907 -0.0680 -0.9880 0.1384 -0.0128 0.0129 0.1377 0.9903 0.0133 D1 = -46.8046 0 0 0 0 5.0954 0 0 0 0 12.1338 0 0 0 0 65.5754 k= 12 Bk = 65.5754 -0.0000 0.0000 0.0000 -0.0000 -46.8046 -0.0000 0.0000 -0.0000 0.0000 5.0954 -0.0000 0.0000 -0.0000 -0.0000 12.1338 V= 0.7229 0.6899 0.0373 0.0059 -0.6907 0.7206 0.0600 -0.0103 0.0128 -0.0680 0.9880 0.1384 -0.0133 0.0129 -0.1377 0.9903 D= 65.5754 0 0 0 0 -46.8046 0 0 0 0 5.0954 0 0 0 0 12.1338 Wc = 6.2740e-007。

实验名称：项目一 Jacobi 旋转法，Jacobi 过关法实验题目：用Jacobi 旋转法，Jacobi 过关法计算下面矩阵A 的全部特征值以及特征向量2100121001210012A -⎡⎤⎢⎥--⎢⎥=⎢⎥--⎢⎥-⎣⎦ 实验过程：1、 Jacobi 旋转法程序如下：function [EigV alMat,EigV ecMat]=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 设置为累加器,记录迭代的次数while abs(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;numfunction MaxMes=FinMax(A)% 对上三角进行扫描，寻找矩阵A的绝对值最大的非对角元素的位置及所有非对角元素的平方和% 输入量:实对称矩阵A% 输出量:MaxMes记录矩阵A的绝对值最大的非对角元素的行号列号及所有非对角元素的平方和n=size(A);n=n(1);trc=0;MaxVal=0; % MaxVal记录矩阵A的绝对值最大的非对角元素赋初值为0for i=1:n-1for j=i+1:ntrc=trc+2*A(i,j)^2;if abs(A(i,j))>MaxValMaxVal=abs(A(i,j));l=i; % l为绝对值最大的非对角元素的行号s=j; % s为绝对值最大的非对角元素的列号endendendMaxMes=[l,s,trc];function RotMat=ComRotMat(A,l,s)% 计算平面旋转矩阵% 输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s% 输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);if A(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 [EigV alMat,EigV ecMat]=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设置为累加器,记录迭代的次数for i=1:n-1for j=i+1:ntrc=trc+2*A(i,j)^2; % 计算矩阵A的非对角元素平方和endendr0=trc^(1/2);r1=r0/n;while r1>=epsfor i=1:n-1for j=i+1:nif abs(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;numfunction RotMat=ComRotMat(A,l,s)% 计算平面旋转矩阵% 输入量:实对称矩阵A及矩阵A的绝对值最大的非对角元素的行号l列号s% 输出量:旋转的平面旋转矩阵RotMatn=size(A);n=n(1);if A(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 -1 0 0;-1 2 -1 0;0 -1 2 -1;0 0 -1 2];eps=0.0001;[EigValMat,EigV ecMat]=JocobiRot(A,eps) num =7EigValMat =0.3820 0 0 -0.00000 3.6180 0 00 -0.0000 1.3820 00 0 0 2.6180EigVecMat =0.3717 0.6015 0.6015 0.3717-0.3717 0.6015 -0.6015 0.3717-0.6015 -0.3717 0.3717 0.60150.6015 -0.3717 -0.3717 0.60152、Jacobi过关法运行结果如下：>> [EigValMat,EigVecMat]=JocobiGG(A,eps)num =16EigValMat =0.3820 0.0000 -0.0000 -0.00000.0000 3.6180 -0.0000 0.0000-0.0000 -0.0000 1.3820 0.00000.0000 0.0000 0.0000 2.6180EigVecMat =0.3717 0.6015 0.6015 0.3717-0.3718 0.6015 -0.6015 0.3717-0.6015 -0.3717 0.3717 0.60150.6015 -0.3717 -0.3718 0.6015实验名称：项目二 Romberg 方法求数值积分 实验题目：用Romberg 方法计算下列积分：（1）30sin x e xdx ⎰，要求准确到610-（2）210x e dx -⎰，要求准确到610-实验过程：function I=romberg(fun,a,b,eps)m=1;h=b-a;I=h/2*(feval(fun,a)+feval(fun,b));T(1,1)=I;while 1N=2^(m-1);h=h/2;I=I/2;for i=1:NI=I+h*feval(fun,a+(2*i-1)*h);endT(m+1,1)=I;M=2*N;k=1;while M>1T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1);M=M/2;k=k+1;endif abs(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.8595>> fun=inline('exp(-x^2)');a=0;b=1;eps=0.000001;I=romberg(fun,a,b,eps)I =0.7468。

1 ，其中 1 和 n 分别是矩阵 A 的 n
2.4 迭代法
2.4.1 迭代法的一般形式及其收敛性
x ( k 1) Gx ( k ) d (k 0,1,...)
定义 设 n n 矩阵 G 的特征值是 1 , 2 ,..., n ，称 (G ) max | i | 为矩阵 G 的谱半径。
n T
x 1 xi
i 1
n
x2 x
则 1 ， 2 和 都是向量范数。 定理 1.2 设


x
i 1 1 i n
n
2 i
max xi
和

是 R 上的任意两种向量范数，则存在与向量 x 无关的常数 m 和
n
M(0<m<M)，使下列关系式成立
m x
1.3.2 矩阵范数
~
若 f '(a ) 0 且 | f ''( a ) | / | f '( a ) | 不很大，则有误差估计
e(u ) f '(a )e(a )
~
(u ) f '(a) (a)
~
。
若 f '(a ) f ''(a ) ... f
( k 1)
(a ) 0, f
(k )
... ... ... ... ln ,n 1
为节省空间，用 C(m,n)存储 A 的带内元素，其中 m=r+s+1，并且 aij ci j s 1, j 。 2.2.5 拟三对角线性方程组的求解方法
a1 d 2 A cn p1 d 2 r1
e xa e ，称 er 为近似值 a 的相对误差。由于 x 未知，实际上总把 作为 a 的 x x a e xa ， 相对误差一般用百分比表示。er 的上界， 即r a a |a|

一、引言Jacobi方法是一种用于计算矩阵特征值和特征向量的迭代数值方法。

它是数值线性代数中的重要算法之一，广泛应用于科学计算、工程技术和金融领域。

本文将通过一个例题来介绍Jacobi方法的原理和求解过程，并分析其在实际问题中的应用。

二、Jacobi方法的原理Jacobi方法是一种通过迭代对矩阵进行相似变换，使得原矩阵逐步转化为对角矩阵的方法。

通过数值迭代，可以逐步逼近矩阵的特征值和对应的特征向量。

其基本原理如下：1. 对称矩阵特征值问题：对于对称矩阵A，存在一个正交矩阵P，使得P^T * A * P = D，其中D为对角矩阵，其对角线上的元素为A的特征值。

所以我们可以通过迭代找到P，使得P逼近正交矩阵，从而逼近A的特征值和特征向量。

2. Jacobi迭代：Jacobi方法的基本思想是通过正交相似变换，逐步将矩阵对角化。

具体来说，对于矩阵A，找到一个旋转矩阵G，使得A' = G^T * A * G为对角矩阵，然后递归地对A'进行相似变换，直到达到精度要求。

三、Jacobi方法求解特征值和特征向量的例题考虑以下矩阵A：A = [[4, -2, 2],[-2, 5, -1],[2, -1, 3]]我们将通过Jacobi方法来计算矩阵A的特征值和特征向量。

1. 对称化矩阵我们需要对矩阵A进行对称化处理。

对称化的思路是找到正交矩阵P，使得P^T * A * P = D，其中D为对角矩阵。

我们可以通过迭代找到逼近P的矩阵序列，直到达到一定的精度。

2. Jacobi迭代在Jacobi迭代的过程中，我们需要找到一个旋转矩阵G，使得A' =G^T * A * G为对角矩阵。

具体的迭代过程是：找到矩阵A中绝对值最大的非对角元素a[i][j]，然后构造一个旋转矩阵G，将a[i][j]置零。

通过迭代地对A'进行相似变换，最终使得A'的非对角元素逼近零，即达到对角化的目的。

3. 计算特征值和特征向量经过一定次数的Jacobi迭代后，得到了对称矩阵A的对角化矩阵D和正交矩阵P。

#ifnd‎e f ks‎f loat‎‎#defi‎n e ks‎f loat‎ dou‎b le#‎e ndif‎#if‎n def ‎i nt16‎‎#defi‎n e in‎t16 ‎i nt#‎e ndif‎#if‎n def ‎u int1‎6‎#def‎i ne u‎i nt16‎unsi‎g ned ‎i nt#‎e ndif‎#if‎n def ‎k snew‎‎#def‎i ne k‎s new(‎s,t) ‎((t*)‎m allo‎c((s)‎*size‎o f(t)‎))#e‎n dif‎#ifn‎d ef k‎s free‎‎#defi‎n e ks‎f ree(‎p) if‎(p!=0‎){fre‎e(p);‎p=0;}‎#end‎i f/‎*‎ksJa‎c obiG‎e tRea‎l SymM‎a trix‎F eatu‎r eVal‎u e‎用ja‎c obi方‎法求取*实‎对称*矩阵‎的特征值和‎特征向量‎‎ aM‎a trix‎[in]:‎ n 阶‎实对称阵‎‎ si‎d e[in‎]: 矩阵‎阶大小‎‎ eps‎[in]:‎控制精度‎‎ m‎a xTim‎e s[in‎]: 最大‎迭代次数‎返‎回值：‎‎前si‎d e*si‎d e个单元‎存储特征向‎量，一列为‎一特征向量‎；后sid‎e个存储单‎元为特征值‎，‎‎也就是总大‎小为 si‎z e = ‎s ide*‎(side‎+1) 存‎储单元*‎/ksf‎l oat ‎*ksJa‎c obiG‎e tRea‎l SymM‎a trix‎F eatu‎r eVal‎u e( k‎s floa‎t aMa‎t rix[‎] , u‎i nt16‎side‎,ks‎f loat‎eps ‎, uin‎t16 m‎a xTim‎e s )‎{‎ ksf‎l oat ‎*t, *‎s, *s‎1, *v‎e ctor‎, *q1‎;‎ flo‎a t a1‎, c ‎, s2 ‎, t1 ‎, max‎E ps =‎0 , ‎d1 , ‎d2, r‎= 1;‎‎u int1‎6 i ,‎j , ‎p , q‎, m ‎, tim‎e s = ‎0 , l‎e n = ‎0;‎ ui‎n t16 ‎s tart‎I ndex‎= 0 ‎, end‎I ndex‎;‎ end‎I ndex‎= si‎d e + ‎s tart‎I ndex‎;‎side‎+= s‎t artI‎n dex;‎‎len ‎= sid‎e*sid‎e;‎ t ‎= ksn‎e w(le‎n , k‎s floa‎t);‎ s ‎= ks‎n ew(l‎e n , ‎k sflo‎a t);‎ s‎1 = k‎s new(‎l en ,‎ksfl‎o at);‎‎q1 = ‎k snew‎(len ‎, ksf‎l oat)‎;‎vect‎o r = ‎k snew‎(len ‎+ sid‎e , k‎s floa‎t);‎ me‎m set ‎( vec‎t or ,‎0 , ‎(len+‎s ide)‎*size‎o f(ks‎f loat‎));‎ fo‎r( i ‎= sta‎r tInd‎e x ; ‎i < e‎n dInd‎e x ; ‎i++)‎ {‎‎ v‎e ctor‎[i*si‎d e+i]‎= 1;‎//将初始‎Q[i][‎j]置为单‎位矩阵‎ }‎‎w hile‎(r >=‎eps ‎&& ti‎m es <‎maxT‎i mes)‎‎{‎‎p = q‎= st‎a rtIn‎d ex; ‎m axEp‎s = 0‎;‎‎f or(i‎= st‎a rtIn‎d ex; ‎i < e‎n dInd‎e x; i‎++)‎‎ {‎‎‎for(‎j = s‎t artI‎n dex;‎j < ‎e ndIn‎d ex; ‎j++)‎‎‎ {‎‎‎‎i f( (‎j != ‎i) &&‎fabs‎(aMat‎r ix[i‎*side‎+j]) ‎> max‎E ps) ‎‎‎‎ {‎‎‎‎ m‎a xEps‎= (f‎l oat)‎f abs(‎a Matr‎i x[i*‎s ide+‎j]);/‎/找非对角‎元素绝对值‎最大的元‎‎‎‎‎p = i‎; q =‎j;‎‎‎‎}‎‎ }‎‎ }‎‎‎r = m‎a xEps‎;//重置‎当前误差‎‎ a1‎= (f‎l oat)‎((aMa‎t rix[‎q*sid‎e+q]-‎a Matr‎i x[p*‎s ide+‎p])/(‎2*aMa‎t rix[‎p*sid‎e+q])‎);‎‎if(a‎1 >= ‎0)‎‎{‎‎‎t1 = ‎(floa‎t)(1/‎(fabs‎(a1)+‎s qrt(‎1+a1*‎a1)))‎;‎‎}‎‎e lse‎‎ {‎‎‎ t1 ‎= (fl‎o at)(‎-1/(f‎a bs(a‎1)+sq‎r t(1+‎a1*a1‎)));‎‎ }‎‎ c ‎= (fl‎o at)(‎1/sqr‎t(1+t‎1*t1)‎);‎‎s2 =‎t1*c‎;‎‎mems‎e t ( ‎s , 0‎, le‎n*siz‎e of(k‎s floa‎t));‎‎ fo‎r( i ‎= sta‎r tInd‎e x ; ‎i < e‎n dInd‎e x ; ‎i++)‎‎ {‎‎‎ s[i‎*side‎+i] =‎1;//‎将初始Q[‎i][j]‎置为单位矩‎阵‎‎}‎‎s[p*‎s ide+‎p] = ‎c ;s‎[p*si‎d e+q]‎=s2;‎‎ s[‎q*sid‎e+p] ‎= -s2‎;s[q*‎s ide+‎q]=c;‎‎ f‎o r(i ‎= sta‎r tInd‎e x; i‎< en‎d Inde‎x; i+‎+)‎‎{‎‎‎f or(j‎= st‎a rtIn‎d ex; ‎j < e‎n dInd‎e x; j‎++)‎‎‎{‎‎‎ s‎1[i*s‎i de+j‎] = s‎[j*si‎d e+i]‎;//将矩‎阵s1[i‎][j]化‎为s[i]‎[j]的转‎置‎‎ }‎‎ }‎‎‎f or(i‎= st‎a rtIn‎d ex; ‎i < e‎n dInd‎e x; i‎++)‎‎ {‎‎‎for(‎j = s‎t artI‎n dex;‎j < ‎e ndIn‎d ex; ‎j++)‎‎‎ {‎‎‎‎d1 = ‎0;‎‎‎ f‎o r(m ‎= sta‎r tInd‎e x; m‎< en‎d Inde‎x; m+‎+)‎‎‎ {‎‎‎‎‎d1 +‎= (fl‎o at)(‎s1[i*‎s ide+‎m]*aM‎a trix‎[m*si‎d e+j]‎);//计‎算s1[i‎][j]*‎a Matr‎i x[i]‎[j] ‎‎‎‎}‎‎‎ t[‎i*sid‎e+j] ‎= d1;‎‎‎ }‎‎ }‎‎ f‎o r(i ‎= sta‎r tInd‎e x; i‎< en‎d Inde‎x; i+‎+)‎‎{‎‎‎f or(j‎= st‎a rtIn‎d ex; ‎j < e‎n dInd‎e x; j‎++)‎‎‎{‎‎‎ d‎1 = d‎2 = 0‎;‎‎‎ fo‎r(m =‎star‎t Inde‎x; m ‎< end‎I ndex‎; m++‎)‎‎‎ {‎‎‎‎‎d1 +=‎(flo‎a t)(t‎[i*si‎d e+m]‎*s[m*‎s ide+‎j]);/‎/计算t[‎i][j]‎*s[i]‎[j]‎‎‎‎ d‎2 += ‎(floa‎t)(ve‎c tor[‎i*sid‎e+m]*‎s[m*s‎i de+j‎]);//‎计算Q[i‎][j]*‎s[i][‎j]‎‎‎ }‎‎‎‎ aMa‎t rix[‎i*sid‎e+j] ‎= d1;‎‎‎‎ q1[‎i*sid‎e+j] ‎= d2;‎‎‎ }‎‎ }‎‎ m‎e mcpy‎( ve‎c tor ‎, q1 ‎, len‎*size‎o f(ks‎f loat‎));‎‎ tim‎e s +=‎1;‎ }‎‎// 对特‎征值与特征‎向量进行整‎合‎for(‎i = s‎t artI‎n dex;‎i < ‎e ndIn‎d ex; ‎i++)‎ {‎‎ v‎e ctor‎[len+‎i] = ‎a Matr‎i x[i*‎s ide+‎i];‎ }‎‎k sfre‎e ( t‎);‎ ks‎f ree ‎( s )‎;‎ksfr‎e e ( ‎s1 );‎‎k sfre‎e ( q‎1 );‎‎i f(ti‎m es >‎maxT‎i mes)‎‎{‎‎k sfre‎e ( v‎e ctor‎); ‎‎ ret‎u rn N‎U LL;‎ }‎‎retu‎r n ve‎c tor;‎}‎。

jacobian矩阵的特征值摘要：1.引言2.Jacobian 矩阵的定义和性质3.特征值和特征向量的概念4.Jacobian 矩阵的特征值求解方法5.应用实例6.结论正文：1.引言在机器学习和人工智能领域，Jacobian 矩阵被广泛应用于优化算法、反向传播算法以及链式法则等。

了解Jacobian 矩阵的特征值对于分析算法的收敛性和稳定性具有重要意义。

本文将介绍Jacobian 矩阵的特征值及其求解方法。

2.Jacobian 矩阵的定义和性质Jacobian 矩阵是用于描述多元函数在某点偏导数的矩阵。

设f(x) 是一个多元函数，其对变量x_i 的偏导数为f_x_i(x)，则Jacobian 矩阵J_f(x) 可以表示为：J_f(x) = | f_x_1(x) f_x_2(x)...f_x_n(x) |其中，n 为函数f(x) 的维数。

Jacobian 矩阵具有以下性质：(1) Jacobian 矩阵是方阵，其行数等于列数，即n×n。

(2) Jacobian 矩阵的元素是函数f(x) 的偏导数。

(3) Jacobian 矩阵的行列式表示函数f(x) 在点x 处的方向导数。

3.特征值和特征向量的概念特征值和特征向量是线性代数中的基本概念。

对于给定的矩阵A，如果存在非零向量x 和标量λ，使得Ax = λx，则λ称为矩阵A 的特征值，x 称为对应于特征值λ的特征向量。

4.Jacobian 矩阵的特征值求解方法求解Jacobian 矩阵的特征值和特征向量，可以采用以下步骤：(1) 计算Jacobian 矩阵。

(2) 计算Jacobian 矩阵的行列式。

(3) 求解行列式为0 的方程，得到特征值。

(4) 构造特征值方程，求解特征向量。

5.应用实例考虑函数f(x) = x^2 + 3x + 2，求其在点x=1 处的Jacobian 矩阵特征值。

首先，计算函数的偏导数：f_x_1(x) = 2x + 3f_x_2(x) = 2然后，计算Jacobian 矩阵：J_f(1) = | 2*1 + 3 2 |= | 5 2 |接下来，计算行列式：det(J_f(1)) = 5 * 2 - 2 * 5 = 0由此得到特征值λ=0。

a pp aqq 1 sin 2 aqq a pp 0 2 2
为保证单值，限定 | | / 4 。
Jacobi 算法
2a pq tan 2 aqq a pp
(1) (2) (3) (4)
在 A 的非主对角线元素中，找到最大元 apq . 用式(3.14)计算 tan2，求 cos, sin 及矩阵Upq . 用公式 (3.13) 求 A1 。 若 max | bij | ,停止计算. 否则, 令 A = A1 , 重复执行 (1) ~ (4).

平面旋转矩阵 对于 p ≠ q,下面定义的n 阶矩阵 Upq 是平面旋转矩阵。
第p列 第q列 1 0 0 1 cos sin 第p 行 1 1 第q 行 sin cos 1 0 0
A1 仍然是实对称阵，因为，UpqT =Upq-1 ，知 A1与 A 的 特征值相同.


下面，以 4 阶矩阵为例，来计算 (3.12)
取 p 2, q 3, A1 U T AU pq . pq
A= aij ,
0 0 1 0 cos sin U 23 0 sin cos 0 0 0
2 2 2 bpi bqi a 2 aqi , i p, q (2 ) pi
由矩阵 A1，A 均是对称阵，得
2 2 2 2 bip biq aip aiq , i p, q (3 )
正交阵
再注意，相似阵的迹不变，我们得到下面的关系式：
bpq bpp bpp bpq cos cos sin a2pq cos sin cos sin sin a pp 2 B2 (4) a A2 sin cos (5) bqq cos sin cos cos bqp bqp bqq sin sin qp aqq

xCni 1且x 0
x x
x H Ax x H Ex max max H xCni 1且x 0 x x xCni 1且x 0 x H x
x Ax max i 由定理9.2.3，xCni1且x 0 x H x
又由定理9.2.2，对任意x≠0，有
... a1 a0
它是关于参数λ的n次多项式，称为矩阵A的特 征多项式， 其中a0=(-1)n｜A｜.
显然，当λ是A的一个特征值时，它必然
是 f ( ) 0的根. 反之，如果λ是 f ( ) 的根， 0
那么齐次方程组(2)有非零解向量x，使(1)式
成立. 从而，λ是A的一个特征值.
k min R( x)或k
xCk 且x 0
xCnk 1且x 0
min
R( x )
9.2.2 极值定理
定理9.2.4(极值定理)n 设Hermite矩阵的n个 1 2 ... 特征值为 , u ,..., u ，其相应的标准酉 u1 2 n 交特征向量为 . 用Ck表示酉空 间Cn中任意的k维子空间，那么
London, England: Millennium ('Wobbly') Bridge (1998-2002, Norman Foster and Partners and Arup Associates)
I decide that I have to write something today, otherwise I would not know how to speak English here. • This is a very quick story about a bridge. • London launched three major construction projects to celebrate the arrival of the Millennium. After all, Greenwich (pronounced green-ich) is supposed to be (supposed to be?!) where the prime meridian lies, and the place where the Millennium officially starts in the world. The three projects are the Millennium Dome in North Greenwich, so far the largest single roofed structure in the world, London Eye right across Westminster, which becomes so far the largest observation wheel in the world, and the Millennium Bridge that links Southeast London with St. Paul’s Cathedral, which is currently…well...not swinging any more, it is said.

定消元容限.常用的一种方法是在第一轮取容限 1 1 1 n 2 2 1 E0 F [2 aij ] , n n i , j 1 i j 第二轮取容限
1 2 E1 n
F
1 n (s) 2 2 [2 (aij ) ] , n i , j 1
i j
1
(s) s 其中 aij 是做完第一轮消元后得到的 As 的(i,j)元素，
(3.3)
为使bpq=0,必须
(aqq a pp ) sin 2 2a pq cos 2 0,
即旋转角 应满足关系式: tg 2 2a pq /(a pp aqq ).
, ] 上,若 a pp aqq ,则可取 sign(a pq ) . 4 4 4 称A的(p,q)位置元素 a pq 为旋转主元.
前面介绍计算矩阵的部分特征值和特征向量的一些方法.现介绍的
Jacobi方法是计算一个实对称矩阵的全部特征值和特征向量的方法. 设 A [aij ] 是任一n阶实对称矩阵,则必存在一个直交相似变换矩 阵U将它化为一个对角阵,即
UAU T D diag(1 , 2 ,, n ),
其中 1 , 2 ,, n 就是A的n个特征值,而U的第j列向量是与 j 相应 的特征值向量.但是这种直交相似变换矩阵U必须用一个无限迭代过程 求得.也就是说,必须通过一系列的直交相似变换 U1 ,U 2 ,,U k , 把矩阵A化为一个对角阵:
常将 限制在区间 [

(3.4) (3.5)
(0) A [aij ] ,即令 A [ a 用 0 ij ]表示原来给定的n阶实对称矩阵
A0 A
(0) 其中 aij aij , i, j 1,, n .古典的Jacobi方法首先选取矩阵 A0 主对角

用Jacobi方法求正交矩阵的特征值与特征向量PROGRAM TEZHENG_Jacobi* 用Jacobi方法求正交矩阵的特征值与特征向量* 就是用平面旋转矩阵U不断对矩阵A作正交相似变换把A化为对角矩阵, * 从而求出A的特征值与特征向量** Made by ZhaoZunsheng,2005.12.20*IMPLICIT NONEINTEGER NMAX,N,I,J,P,QPARAMETER(NMAX=50)REAL A(NMAX,NMAX),AMAX,TEMP,ZEMP,COO,SII,CO,SI,APP,AQQ,APQ,API,AQI REAL R(NMAX,NMAX),RIP,RIQCHARACTER NAME*12,NAMEO*12,CHR*1* 从文件中读入实对称矩阵AWRITE(*,*)'输入实对称矩阵维数n(n<51):'READ(*,*) NWRITE(*,*)'输入矩阵文件:'READ(*,*) NAMEOPEN(6,FILE=NAME)DO I=1,NREAD(6,*) (A(I,J),J=1,I)DO J=1,IA(J,I)=A(I,J)ENDDOENDDOCLOSE(6)* R矩阵存放正交变换矩阵U,在这先初始化,即单位矩阵DO I=1,NDO J=1,NR(I,J)=0ENDDOR(I,I)=1ENDDO* 在矩阵A的非主对角线元素中,找出按模最大的元素Apq 100 AMAX=ABS(A(2,1))P=2Q=1DO I=2,NDO J=1,I-1IF(ABS(A(I,J)).GT.AMAX) THENAMAX=ABS(A(I,J))P=IQ=JENDIFENDDOENDDO* do i=1,n* WRITE(*,*)(A(I,J),J=1,N)* enddo* 当非主对角线元素化为0,即小于给定精度时,输出特征值与特征向量IF(AMAX.LE.1.0E-7) THENWRITE(*,*) 'A的特征值为:'WRITE(*,*) (A(I,I),I=1,N)WRITE(*,*) 'A的特征向量为:'WRITE(*,*) ' X1 X2 X3 ...:'DO I=1,NWRITE(*,*)(R(I,J),J=1,N)ENDDOWRITE(*,*) '是否将结果存入文件(Y/N)?'READ(*,*) CHRIF(CHR.EQ.'Y'.OR.CHR.EQ.'y') THENWRITE(*,*) '输入文件名(小于12字符):'READ(*,*) NAMEOOPEN(8,FILE=NAMEO)WRITE(8,*) 'A的特征值为:'WRITE(8,*) (A(I,I),I=1,N)WRITE(8,*) 'A的特征向量为:'WRITE(8,*) ' X1 X2 X3 ...:'DO I=1,NWRITE(8,*)(R(I,J),J=1,N)ENDDOCLOSE(8)ENDIFSTOPENDIF* 开始准备计算平面旋转矩阵UTEMP=2*A(P,Q)/(A(P,P)-A(Q,Q)+1.0e-30) ZEMP=(A(P,P)-A(Q,Q))/(2*A(P,Q))IF(ABS(TEMP).LT.1.0) THENCOO=(1+TEMP**2)**(-0.5)SII=TEMP*(1+TEMP**2)**(-0.5)ELSECOO=ABS(ZEMP)*(1+ZEMP**2)**(-0.5)SII=SIGN(1.0,ZEMP)*(1+ZEMP**2)**(-0.5) ENDIFCO=SQRT(0.5*(1+COO))SI=SII/(2*CO)* 计算平面旋转矩阵UDO I=1,NRIP=R(I,P)*CO+R(I,Q)*SIRIQ=-R(I,P)*SI+R(I,Q)*COR(I,P)=RIPR(I,Q)=RIQENDDO* 对A进行变换APP=A(P,P)*CO**2+A(Q,Q)*SI**2+2*A(P,Q)*CO*SI AQQ=A(P,P)*SI**2+A(Q,Q)*CO**2-2*A(P,Q)*CO*SI APQ=0.5*(A(Q,Q)-A(P,P))*SII+A(P,Q)*COOA(P,P)=APPA(Q,Q)=AQQA(P,Q)=APQA(Q,P)=A(P,Q)DO I=1,NIF(I.EQ.P.OR.I.EQ.Q) THENELSEAPI=A(P,I)*CO+A(Q,I)*SIAQI=-A(P,I)*SI+A(Q,I)*COA(P,I)=APIA(Q,I)=AQIA(I,P)=A(P,I)A(I,Q)=A(Q,I)ENDIFENDDOGOTO 100END。

jacobi迭代法求复矩阵特征值和特征向量Jacobi迭代法是一种经典的求解复矩阵特征值和特征向量的方法。

在数值分析领域，特征值和特征向量的求解是一个十分重要且常见的问题。

它不仅在理论上有重要意义，还在实际应用中有着广泛的应用，比如在物理、工程、金融等领域。

Jacobi迭代法的提出，极大地简化了这个复杂问题的求解过程，为研究人员和工程师提供了一个高效、可靠的数值计算工具。

我们需要了解什么是特征值和特征向量。

对于一个n阶方阵A，如果存在数λ和一个非零向量x，使得Ax=λx成立，则称λ是A的特征值，x是对应于特征值λ的特征向量。

特征值和特征向量的求解十分重要，因为它们包含了矩阵A的重要特性和信息，对于矩阵的对角化、矩阵的稳定性、矩阵的特征分解等问题有着重要的作用。

接下来，让我们来介绍Jacobi迭代法的基本思想和步骤。

Jacobi迭代法的核心思想是通过一系列相似变换，将原始矩阵对角化，从而得到其特征值和特征向量。

具体步骤如下：1. 我们选择一个n阶方阵A，将其初始化为对角矩阵D，将初始的特征向量矩阵初始化为单位矩阵I。

2. 我们选择两个不同的下标i和j(1≤i，j≤n，i≠j)，使得矩阵A的元素aij为非零元素，即aij≠0。

这两个下标表示我们要进行的相似变换的维度。

3. 我们构造一个旋转矩阵P，使得通过P的相似变换，可以将aij对应的元素变为0。

这一步是Jacobi迭代法的核心步骤，旋转矩阵P的构造涉及到对称双射矩阵的变换和特征值的迭代计算。

4. 我们通过P的相似变换，更新矩阵A和特征向量矩阵I，得到新的对角矩阵D和新的特征向量矩阵。

5. 我们检查新得到的对角矩阵D的非对角线元素是否足够小，如果满足要求，则停止迭代，否则继续进行第2步的操作。

通过这样一系列的迭代操作，我们可以逐步地将矩阵A对角化，并得到其特征值和特征向量。

Jacobi迭代法以其简洁、直观的特点，在复矩阵特征值和特征向量的求解中得到了广泛的应用。

题目：深入探究jacobi迭代法求复矩阵特征值和特征向量上线性代数的学习过程中，我们经常会遇到求解复矩阵的特征值和特征向量的问题。

而jacobi迭代法则是一种被广泛应用的方法之一。

本文将深入探讨jacobi迭代法的原理、应用以及个人观点和理解。

### 1. jacobi迭代法的原理和概念jacobi迭代法是一种通过不断相似变换将矩阵对角化的方法，它可以被用于求解实对称矩阵的特征值和特征向量，而在这篇文章中，我们将着重讨论其在求解复矩阵时的应用。

### 2. jacobi迭代法的算法步骤在使用jacobi迭代法求解复矩阵特征值和特征向量时，我们需要经历一系列的算法步骤。

我们可以通过对角线元素的绝对值大小来判断矩阵是否已经对角化，然后进行迭代，直到满足精度要求为止。

### 3. jacobi迭代法的实际应用在实际应用中，jacobi迭代法除了可以求解复矩阵的特征值和特征向量外，还可以在解决其他涉及特征值和特征向量的问题时发挥重要作用。

通过简单的算法步骤和迭代过程，我们可以有效地得到复矩阵的特征值和特征向量，为进一步的分析和计算提供便利。

### 4. 个人观点和理解从个人的角度来看，jacobi迭代法在求解复矩阵特征值和特征向量时具有一定的优势，尤其在算法实现的过程中，我们可以通过简单的迭代步骤快速得到结果。

然而，对于大规模复矩阵的计算，可能还需要考虑其他更高效的方法或并行计算的应用。

### 结论通过本文的深入探讨，我们对jacobi迭代法求解复矩阵特征值和特征向量有了更深入的了解。

在实际应用中，我们需要灵活运用不同的方法和算法，以便更好地解决实际问题。

总结来说，jacobi迭代法是一种常用的求解复矩阵特征值和特征向量的方法，它通过简单的算法步骤和迭代过程，能够快速有效地得到结果。

然而，在实际应用中，我们还需要综合考虑不同的因素，以便获得更好的计算效果。

通过本文的阐述，希望读者能够更加深入地理解jacobi迭代法以及其在求解复矩阵特征值和特征向量中的应用，为进一步的学习和研究打下良好的基础。

《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。

变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值 4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。

对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。

如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。