北航数值分析1-Jacobi法计算矩阵特征值
数值分析—矩阵特征值问题的数值计算
( vk +1 )i vk = α1 x1 , lim = λ1 . k k →∞ λ k →∞ ( v ) k i 1
lim
可见，当 k 充分大时， vk 近似于主特征向量（相差一个常数倍） ， vk +1 与 vk 的对应非零分 量的比值近似于主特征值。 在实际计算中，需要对计算结果进行规范化。因为当 λ1 <1 时， vk 趋于零；当 λ1 >1 时 ， vk 的 非 零 分 量 趋 于 无 穷 ， 从 而 计 算 时 会 出 现 下 溢 或 上 溢 。 为 此 ， 对 向 量
λ1 ≥ λ2 ≥ L ≥ λn ,
∑x y
i =1 i
n
i
为向量 x 和 y 的内积。
定理 8.3 设 A 为 n 阶实对称矩阵，其特征值都为实数，排列为 对应的特征向量 x1 , x2 ,L , xn 组成正交向量组，则有 1) 对任何非零向量 x ∈ R n ，有 λn ≤ R( x) ≤ λ1 ， 2) λ1 = max R ( x) = R( x1 ) ，
Ζ = ( z1 , z2 ,L, zn )T ∈ R n ， 记 max( Ζ) = zi ，其中 zi = Ζ ∞ ，这样， 我们有 如 下 乘幂 法 的实 用
的计算公式： 任取 v0 = u0 ≠ 0 ，对于 k = 1, 2,L 分别计算 vk = Auk −1 , uk = vk / max(vk ). 求出对应矩阵的主特征向量和特征值的近似值，有下面的定理。 定理 8.4
m1 0 M = 0 0 0
称为质量矩阵，而
0 m2 0 0 0
0 0 m3 0 0
0 0 0 m4 0
0 0 0 0 m5 0 0 −k4 k 4 + k5 − k5 0 0 0 − k5 k5
《实验8：Jacobi法求实对称矩阵的特征值及特征向量》
实验名称实验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。
矩阵特征值与特征向量的计算-Jacobi方法
所在平面旋转了一个角度 , 其它坐标保持不 变, 称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)任是实对称
北京航空航天大学数值分析课程知识点总结
北京航空航天大学数值分析课程知识点总结1.2 误差知识与算法知识1.2.2 绝对误差、相对误差与有效数字设a 是准确值x 的一个近似值，记e x a =-，称e 为近似值a 的绝对误差，简称误差。
如果||e 的一个上界已知，记为ε，即||e ε≤，则称ε为近似值a 的绝对误差限或绝对误差界，简称误差限或误差界。
记r e x ae x x-==，称r e 为近似值a 的相对误差。
由于x 未知，实际上总把e a 作为a 的相对误差，并且也记为r e x a e a a -==，相对误差一般用百分比表示。
r e 的上界，即||r a εε=称为近似值a 的相对误差限或相对误差界。
定义设数a 是数x 的近似值。
如果a 的绝对误差限时它的某一位的半个单位，并且从该位到它的第一位非零数字共有n 位，则称用a 近似x 时具有n 位有效数字。
1.2.3 函数求值的误差估计设()u f x =存在足够高阶的导数，a 是x 的近似值，则~()u f a =是()u f x =的近似值。
若'()0f a ≠且|''()|/|'()|f a f a 不很大，则有误差估计~~()'()()()'()()e uf a e a u f a a εε≈≈。
若(1)()'()''()...()0,()0k k f a f a fa f a -====≠，且比值(1)()()/()k k f a fa +不很大，则有误差估计[][]()~()~()()()!()()()!k kk k f a e u e a k f a u a k εε≈≈。
对于n 元函数，有误差估计~121~121(,,...,)()()(,,...,)()()nn i i i nn i i if a a a e u e a x f a a a u a x εε==?≈??≈?∑∑；若一阶偏导全为零或很小，则要使用高阶项。
用共轭梯度法解方程,用Jacobi方法求矩阵的全部特征值和特征向量
end
A(1,1)=4;
A(n,n)=4;
A(n,n-1)=1;
X=zeros(n,1);
fori=1:n
X(i,1)=1;
end
%用共轭梯度法求解方程
fprintf('方程的精确解\n');
X
fprintf('用共轭法求解方程\n');
x=cg(A,b)
%用方法求解方程的特征值和特征向量
-0.1859 -0.0000 0.0000 0.0000 5.6180 -0.0000 -0.0000
-0.2236 0.0000 -0.0000 0.0000 0.0000 5.4142 0.0000
-0.2558 0.0000 -0.0000 0.0000 0.0000 0.0000 5.1756
0.1859 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000
0.1436 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000
-0.0977 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000
0.6634
1.1794
0.6188
1.3453
用方法Jacobi求解矩阵的全部特征值及特征向量
D =
Columns 1 through 7
4.0000 0 0 0 0 0 0
0.1382 5.9021 0.0000 0.0000 -0.0000 0.0000 0.0000
-0.2629 0.0000 5.6180 0.0000 0.0000 -0.0000 -0.0000
八、矩阵特征值的jacobi方法
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方法的原理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。
jacobi方法求取实对称矩阵的特征值和特征向量
#ifnde f ksf loat#defin e ksf loat doub le#e ndif#ifn def i nt16#defin e int16 i nt#e ndif#ifn def u int16#defi ne ui nt16unsig ned i nt#e ndif#ifn def k snew#defi ne ks new(s,t) ((t*)m alloc((s)*sizeo f(t)))#en dif#ifnd ef ks free#defin e ksf ree(p) if(p!=0){free(p);p=0;}#endi f/*ksJac obiGe tReal SymMa trixF eatur eValu e用jac obi方法求取*实对称*矩阵的特征值和特征向量 aMa trix[in]: n 阶实对称阵 sid e[in]: 矩阵阶大小 eps[in]:控制精度 ma xTime s[in]: 最大迭代次数返回值：前sid e*sid e个单元存储特征向量，一列为一特征向量；后side个存储单元为特征值，也就是总大小为 siz e = s ide*(side+1) 存储单元*/ksfl oat *ksJac obiGe tReal SymMa trixF eatur eValu e( ks float aMat rix[] , ui nt16side,ksf loateps , uint16 ma xTime s ){ ksfl oat *t, *s, *s1, *ve ctor, *q1; floa t a1, c , s2 , t1 , maxE ps =0 , d1 , d2, r= 1;u int16 i ,j , p , q, m , time s = 0 , le n = 0; uin t16 s tartI ndex= 0 , endI ndex; endI ndex= sid e + s tartI ndex;side+= st artIn dex;len = side*side; t = ksne w(len , ks float); s = ksn ew(le n , k sfloa t); s1 = ks new(l en ,ksflo at);q1 = k snew(len , ksfl oat);vecto r = k snew(len + side , ks float); mem set ( vect or ,0 , (len+s ide)*sizeo f(ksf loat)); for( i = star tInde x ; i < en dInde x ; i++) { ve ctor[i*sid e+i]= 1;//将初始Q[i][j]置为单位矩阵 }w hile(r >=eps && tim es <maxTi mes){p = q= sta rtInd ex; m axEps = 0;f or(i= sta rtInd ex; i < en dInde x; i++) {for(j = st artIn dex;j < e ndInd ex; j++) {i f( (j != i) &&fabs(aMatr ix[i*side+j]) > maxE ps) { ma xEps= (fl oat)f abs(a Matri x[i*s ide+j]);//找非对角元素绝对值最大的元p = i; q =j;} } }r = ma xEps;//重置当前误差 a1= (fl oat)((aMat rix[q*side+q]-a Matri x[p*s ide+p])/(2*aMat rix[p*side+q]));if(a1 >= 0){t1 = (float)(1/(fabs(a1)+s qrt(1+a1*a1)));}e lse { t1 = (flo at)(-1/(fa bs(a1)+sqr t(1+a1*a1))); } c = (flo at)(1/sqrt(1+t1*t1));s2 =t1*c;memse t ( s , 0, len*size of(ks float)); for( i = star tInde x ; i < en dInde x ; i++) { s[i*side+i] =1;//将初始Q[i][j]置为单位矩阵}s[p*s ide+p] = c ;s[p*sid e+q]=s2; s[q*side+p] = -s2;s[q*s ide+q]=c; fo r(i = star tInde x; i< end Index; i++){f or(j= sta rtInd ex; j < en dInde x; j++){ s1[i*si de+j] = s[j*sid e+i];//将矩阵s1[i][j]化为s[i][j]的转置 } }f or(i= sta rtInd ex; i < en dInde x; i++) {for(j = st artIn dex;j < e ndInd ex; j++) {d1 = 0; fo r(m = star tInde x; m< end Index; m++) {d1 += (flo at)(s1[i*s ide+m]*aMa trix[m*sid e+j]);//计算s1[i][j]*a Matri x[i][j] } t[i*side+j] = d1; } } fo r(i = star tInde x; i< end Index; i++){f or(j= sta rtInd ex; j < en dInde x; j++){ d1 = d2 = 0; for(m =start Index; m < endI ndex; m++) {d1 +=(floa t)(t[i*sid e+m]*s[m*s ide+j]);//计算t[i][j]*s[i][j] d2 += (float)(vec tor[i*side+m]*s[m*si de+j]);//计算Q[i][j]*s[i][j] } aMat rix[i*side+j] = d1; q1[i*side+j] = d2; } } me mcpy( vec tor , q1 , len*sizeo f(ksf loat)); time s +=1; }// 对特征值与特征向量进行整合for(i = st artIn dex;i < e ndInd ex; i++) { ve ctor[len+i] = a Matri x[i*s ide+i]; }k sfree ( t); ksf ree ( s );ksfre e ( s1 );k sfree ( q1 );i f(tim es >maxTi mes){k sfree ( ve ctor); retu rn NU LL; }retur n vec tor;}。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
准备工作算法设计矩阵特征值的求法有幂法、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操作系统）。