北航数值分析1_Jacobi法计算矩阵特征值
《实验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)任是实对称
雅克比法求矩阵特征值特征向量
C语言课程设计报告课程名称：计算机综合课程设计学院：土木工程学院设计题目：矩阵特征值分解级别：B学生姓名：学号：同组学生：无学号：无指导教师：2012年9月5日C语言课程设计任务书（以下要求需写入设计报告书）学生选题说明:以所发课程设计要求为准，请同学们仔细阅读；本任务书提供的设计案例仅供选题参考；也可自选，但难易程度需难度相当；鼓励结合本专业（土木工程、力学）知识进行选题，编制程序解决专业实际问题。
限2人选的题目可由1-2人完成（A级）；限1人选的题目只能由1人单独完成（B级）；设计总体要求:采用模块化程序设计；鼓励可视化编程；源程序中应有足够的注释；学生可自行增加新功能模块（视情况可另外加分）；必须上机调试通过；注重算法运用，优化存储效率与运算效率；需提交源程序(含有注释)及相关文件(数据或数据库文件)；(cpp文件、txt或dat文件等) 提交设计报告书，具体要求见以下说明。
设计报告格式：目录1．课程设计任务书（功能简介、课程设计要求）；2．系统设计（包括总体结构、模块、功能等，辅以程序设计组成框图、流程图解释）；3．模块设计（主要模块功能、源代码、注释（如函数功能、入口及出口参数说明，函数调用关系描述等）；4．调试及测试：（调试方法，测试结果的分析与讨论，截屏、正确性分析）；5．设计总结：（编程中遇到的问题及解决方法）；6．心得体会及致谢；参考文献1.课程设计任务书功能简介：a）输入一个对称正方矩阵A，从文本文件读入；b）对矩阵A进行特征值分解，将分解结果：即U矩阵、S矩阵输出至文本文件；c）将最小特征值及对应的特征向量输出至文本文件；d）验证其分解结果是否正确。
提示：A=USUT,具体算法可参考相关文献。
功能说明：矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。
注：以三阶对称矩阵为例2.系统设计总体结构main函数tezheng函数从文本文找矩阵A递推求矩递推求矩向屏幕和向屏幕和求最小特阵S件中读入阵U中非对角txt文件输txt文件输征值及其数组A入矩阵U入矩阵S元素中的对应特征最大值，向量，并并记下其输出到屏位置幕和txt文件中3.模块设计#include<stdio.h>#include<stdlib.h>#include<math.h>intmain(){FILE*fp;inttezheng(double*a,intn,double*s,double*u,doubleeps,intitmax);//函数调用声明inti,j,p,itmax=1000;//itmax为最大循环次数doubleeps=1e-7,s[3][3],u[3][3];//eps为元素精度，s为对角矩阵S，u为矩阵U doublea[9];//a为待分解矩阵Ai=tezheng(a,3,s,u,eps,1000);WORD格式if(i>0)//i对应函数中的返回值it{if((fp=fopen("juzhen.txt","w"))==NULL)//打开待输入txt文件{printf("无法打开文件.\n");return;}printf("U矩阵为:\n");//下几句分别向屏幕和txt文件输入矩阵Ufprintf(fp,"U矩阵为:\n");for(i=0;i<3;i++){for(j=0;j<3;j++){printf("%10.6f",u[i][j]);fprintf(fp,"%10.6f",u[i][j]);}printf("\n");fprintf(fp,"\n");}printf("S对角矩阵为:\n");//下几句分别向屏幕和txt文件输入矩阵Sfprintf(fp,"S对角矩阵为:\n");for(i=0;i<3;i++){for(j=0;j<3;j++){printf("%10.6f",s[i][j]);fprintf(fp,"%10.6f",s[i][j]);}printf("\n");fprintf(fp,"\n");}p=0;for(i=0;i<3;i++)//下面几句为求最小特征值及其对应特征向量，并输出到屏幕和txt文件中if(s[i][i]<s[0][0])p=i;printf("最小特征值为:%10.6f\n",s[p][p]);fprintf(fp,"最小特征值为:%10.6f\n",s[p][p]);printf("对应特征向量为:\n");fprintf(fp,"对应特征向量为:\n");for(i=0;i<3;i++){printf("%10.6f\n",u[i][p]);fprintf(fp,"%10.6f\n",u[i][p]);}}}inttezheng(double*a,intn,doubles[3][3],doubleu[3][3],doubleeps,intitmax){FILE*A;charline[1000];//存放从文件juzhenA.txt读入的数据char*k="";inti,j,p,q,it;doublesint,cost,sin2t,cos2t,d,tmp,t1,t2,t3,fm;it=0;if((A=fopen("juzhenA.txt","r"))==NULL){printf("无法打开文件\n");return;}fgets(line,1000,A);//从文件指针A指向的文件，即juzhenA.txt中读入字符数据到数组line中strtok(line,k);//原型char*strtok(char*s,constchar*delim);说明：strtok()用来将字符串分割成一个个片段。
北航数值分析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.6630353#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;double cosz;double sinz;double MAX;int kk;//#define PTS pts#ifdef PTSvoid PTS(double *m){printf("-----------------------------------------------------------------------\n");printf(" 迭代第%d次\n",kk);for(int i = 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(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 1 ; i <= PFrols+2 ; i++){printf(" ... ");}putchar(10);}#elsevoid PTS(double *m){}#endifvoid recounti(int i , 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;}}}void refreshMetrix(double *m){int ipr,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(int i = 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);}//void calCosSin(double *m){double app = m[(p+1)*p/2];double aqq = m[(q+1)*q/2];double apq = m[(p-1)*p/2+q];cosz = cos(atan(2 * apq / (app - aqq))/2);sinz = sin(atan(2 * apq / (app - aqq))/2); }//void find_pq(double *m){double max = 0.0;int pp = 0;int qq = 0;for(int i = 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;}void init(double *m){for(int i = 1 ; i <= N ;i++)m[(i+1)*i/2] = a(i);for(int i = 2 ; i <= N ; i++)m[(i-1)*i/2+i-1] = b;for(int i = 3 ; i <= N ; i++)m[(i-1)*i/2+i-2] = c;PTS(m);}void calFinal(double *m){printf("---------------------------------------------------------------------------------------------------\n");printf("结果输出:\n\n");double conda;double deta = 1.0;double minlumda = pow((double)10.0,12);double maxlumda = pow((double)10.0,-12);double absminlumda = pow((double)10.0,12);for(int i = 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(int i = 1 ; i <= FK ; i++){double muk = minlumda + i * (maxlumda - minlumda) / 40;double lumdak = 0.0;double tempabsmin = 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(int argc, _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(int i = 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(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 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;}运行结果如下：中间运行状态如下：结果分析➢有效性分析1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2.代码中有定义预编译宏//#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G存，64位windows10操作系统）。
用共轭梯度法解方程,用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
矩阵特征值问题的数值计算
矩阵特征值问题的计算方法特征值问题：A V=λV¾直接计算：A的阶数较小,且特征值分离得较好 特征值：det(λI-A)=0,特征向量：(λI-A)V=0¾迭代法：幂法与反幂法¾变换法：雅可比方法与QR方法内容：一、 特征值的估计及其误差问题二、 幂法与反幂法三、 雅可比方法四、 QR方法一、 特征值的估计及其误差问题 （一）特征值的估计结论 1.1：n 阶矩阵()ij n n A a ×=的任何一个特征值必属于复平面上的n 个圆盘：1,||||,1,2,ni ii ij j j i D z z a a i n =≠⎧⎫⎪⎪=−≤=⎨⎬⎪⎪⎩⎭∑"（10.1） 的并集。
结论1.2：若（10.1）中的m个圆盘形成一个连通区域D，且D与其余的n-m个圆盘不相连，则D中恰有A的m个特征值。
（二）特征值的误差问题结论1.3：对于n 阶矩阵()ij n n A a ×=，若存在n 阶非奇异矩阵H ，使得11(,,)n H AH diag λλ−=Λ="， （10.2）则11min ||||||||||||||i p p p i nH H A λλ−≤≤−≤∆ （10.3）其中λ是A A +∆的一个特征值，而(1,,)i i n λ="是A 的特征值，1,2,p =∞。
结论1.4：若n 阶矩阵A 是实对称的，则1min ||||||i p i nA λλ≤≤−≤∆。
（10.4）注：（10.4）表明，当A 是实对称时，由矩阵的微小误差所引起的特征值摄动也是微小的。
但是对于非对称矩阵而言，特别是对条件数很大的矩阵，情况未必如此。
二、 幂法与反幂法（一） 幂法：求实矩阵按模最大的特征值与特征向量假设n 阶实矩阵A 具有n 个线性无关的特征向量,1,iV i n ="，则对于任意的0nX R ∈，有 01ni ii X a V ==∑，从而有01111112((/))n nk k k i i i i ii i nk k i i i i A X a A V a V a V a V λλλλ======+∑∑∑.若A 的特征值分布如下：123||||||||n λλλλ>≥≥≥"，则有01111()k kk A X a V λλ→∞⎯⎯⎯→为对应的特征向量须注意的是，若1||1λ<，则10kλ→，出现“下溢”，若1||1λ>，则1kλ→∞，出现“上溢”，为避免这些现象的发生，须对0kA X 进行规范化。
八、矩阵特征值的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。
北航数值分析-实习作业1(C语言详细注释)
《数值分析》计算实习作业《一》北航第一题 设有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，再重复上述步骤即可求得两个特征值。
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
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
jacobian矩阵的特征值
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。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。
1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。
2)使用反幂法求λs ，其中需要解线性方程组。
因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。
（2）与140k λλμλ-5011=+k 最接近的特征值λik 。
通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。
（3）2cond(A)和det A 。
1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。
2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。
由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。
2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。
Jacobi__矩阵特征值和特征向量
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
用Jacobi方法求正交矩阵的特征值与特征向量
用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迭代法是一种经典的求解复矩阵特征值和特征向量的方法。
在数值分析领域，特征值和特征向量的求解是一个十分重要且常见的问题。
它不仅在理论上有重要意义，还在实际应用中有着广泛的应用，比如在物理、工程、金融等领域。
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迭代法以其简洁、直观的特点，在复矩阵特征值和特征向量的求解中得到了广泛的应用。
雅克比法求矩阵特征值特征向量
C语言课程设计报告课程名称：计算机综合课程设计学院：土木工程学院设计题目：矩阵特征值分解级别： B学生姓名：学号：同组学生：无学号：无指导教师：2012年 9 月 5 日C语言课程设计任务书（以下要求需写入设计报告书）学生选题说明:➢以所发课程设计要求为准，请同学们仔细阅读；➢本任务书提供的设计案例仅供选题参考；也可自选，但难易程度需难度相当；➢鼓励结合本专业（土木工程、力学）知识进行选题，编制程序解决专业实际问题。
➢限2人选的题目可由1-2人完成（A级）；限1人选的题目只能由1人单独完成（B级）；设计总体要求:➢采用模块化程序设计；➢鼓励可视化编程；➢源程序中应有足够的注释；➢学生可自行增加新功能模块（视情况可另外加分）；➢必须上机调试通过；➢注重算法运用，优化存储效率与运算效率；➢需提交源程序(含有注释)及相关文件(数据或数据库文件)；(cpp文件、txt或dat文件等) ➢提交设计报告书，具体要求见以下说明。
设计报告格式：目录1．课程设计任务书（功能简介、课程设计要求）；2．系统设计（包括总体结构、模块、功能等，辅以程序设计组成框图、流程图解释）；3．模块设计（主要模块功能、源代码、注释（如函数功能、入口及出口参数说明，函数调用关系描述等）；4．调试及测试：（调试方法，测试结果的分析与讨论，截屏、正确性分析）；5．设计总结：（编程中遇到的问题及解决方法）；6．心得体会及致谢；参考文献1.课程设计任务书功能简介：a）输入一个对称正方矩阵A，从文本文件读入；b）对矩阵A进行特征值分解，将分解结果：即U矩阵、S矩阵输出至文本文件；c）将最小特征值及对应的特征向量输出至文本文件；d）验证其分解结果是否正确。
提示：A=USU T,具体算法可参考相关文献。
功能说明：矩阵特征值分解被广泛运用于土木工程问题的数值计算中，如可用于计算结构自振频率与自振周期、结构特征屈曲问题等。
注：以三阶对称矩阵为例2.系统设计3.模块设计#include<stdio.h>#include<stdlib.h>#include<math.h>int main(){FILE *fp;int tezheng(double *a,int n,double *s,double *u,double eps,int itmax); //函数调用声明int i,j,p,itmax=1000; //itmax为最大循环次数double eps=1e-7,s[3][3],u[3][3]; //eps为元素精度，s为对角矩阵S，u为矩阵Udouble a[9];//a为待分解矩阵Ai=tezheng(a,3,s,u,eps,1000);if(i>0) //i对应函数中的返回值it{if((fp=fopen("juzhen.txt","w"))==NULL) //打开待输入txt文件{printf("无法打开文件.\n");return;}printf("U矩阵为:\n"); //下几句分别向屏幕和txt文件输入矩阵Ufprintf(fp,"U矩阵为:\n");for(i=0;i<3;i++){for(j=0;j<3;j++){printf("%10.6f",u[i][j]);fprintf(fp,"%10.6f",u[i][j]);}printf("\n");fprintf(fp,"\n");}printf("S对角矩阵为:\n"); //下几句分别向屏幕和txt文件输入矩阵Sfprintf(fp,"S对角矩阵为:\n");for(i=0;i<3;i++){for(j=0;j<3;j++){printf("%10.6f",s[i][j]);fprintf(fp,"%10.6f",s[i][j]);}printf("\n");fprintf(fp,"\n");}p=0;for(i=0;i<3;i++)//下面几句为求最小特征值及其对应特征向量，并输出到屏幕和txt 文件中if(s[i][i]<s[0][0])p=i;printf("最小特征值为:%10.6f\n",s[p][p]);fprintf(fp,"最小特征值为:%10.6f\n",s[p][p]);printf("对应特征向量为:\n");fprintf(fp,"对应特征向量为:\n");for(i=0;i<3;i++){printf("%10.6f\n",u[i][p]);fprintf(fp,"%10.6f\n",u[i][p]);}}}int tezheng(double *a,int n,double s[3][3],double u[3][3],double eps,int itmax){FILE *A;char line[1000]; //存放从文件juzhen A.txt读入的数据char *k=" ";int i,j,p,q,it;double sint,cost,sin2t,cos2t,d,tmp,t1,t2,t3,fm;it=0;if((A=fopen("juzhen A.txt","r"))==NULL){printf("无法打开文件\n");return;}fgets(line,1000,A); //从文件指针A指向的文件，即juzhen A.txt中读入字符数据到数组line 中strtok(line,k); //原型char *strtok(char *s, const char *delim);说明：strtok()用来将字符串分割成一个个片段。
jacobi迭代法求复矩阵特征值和特征向量
题目：深入探究jacobi迭代法求复矩阵特征值和特征向量上线性代数的学习过程中，我们经常会遇到求解复矩阵的特征值和特征向量的问题。
而jacobi迭代法则是一种被广泛应用的方法之一。
本文将深入探讨jacobi迭代法的原理、应用以及个人观点和理解。
### 1. jacobi迭代法的原理和概念jacobi迭代法是一种通过不断相似变换将矩阵对角化的方法，它可以被用于求解实对称矩阵的特征值和特征向量，而在这篇文章中，我们将着重讨论其在求解复矩阵时的应用。
### 2. jacobi迭代法的算法步骤在使用jacobi迭代法求解复矩阵特征值和特征向量时，我们需要经历一系列的算法步骤。
我们可以通过对角线元素的绝对值大小来判断矩阵是否已经对角化，然后进行迭代，直到满足精度要求为止。
### 3. jacobi迭代法的实际应用在实际应用中，jacobi迭代法除了可以求解复矩阵的特征值和特征向量外，还可以在解决其他涉及特征值和特征向量的问题时发挥重要作用。
通过简单的算法步骤和迭代过程，我们可以有效地得到复矩阵的特征值和特征向量，为进一步的分析和计算提供便利。
### 4. 个人观点和理解从个人的角度来看，jacobi迭代法在求解复矩阵特征值和特征向量时具有一定的优势，尤其在算法实现的过程中，我们可以通过简单的迭代步骤快速得到结果。
然而，对于大规模复矩阵的计算，可能还需要考虑其他更高效的方法或并行计算的应用。
### 结论通过本文的深入探讨，我们对jacobi迭代法求解复矩阵特征值和特征向量有了更深入的了解。
在实际应用中，我们需要灵活运用不同的方法和算法，以便更好地解决实际问题。
总结来说，jacobi迭代法是一种常用的求解复矩阵特征值和特征向量的方法，它通过简单的算法步骤和迭代过程，能够快速有效地得到结果。
然而，在实际应用中，我们还需要综合考虑不同的因素，以便获得更好的计算效果。
通过本文的阐述，希望读者能够更加深入地理解jacobi迭代法以及其在求解复矩阵特征值和特征向量中的应用，为进一步的学习和研究打下良好的基础。
雅可比迭代法与矩阵的特征值
实验五矩阵的lu分解法，雅可比迭代法班级：学号：姓名：实验五 矩阵的LU 分解法，雅可比迭代一、目的与要求：熟悉求解线性方程组的有关理论和方法；会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序； 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
二、实验内容：会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。
三、程序与实例列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示 1) 消元过程 对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。
③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i kj k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算kk ik ik a a l /=对j=l+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。
②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算ii ni j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下：#include <iostream>#include <cmath>using namespace std;void disp(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void disp(double* q,int n){cout<<"====================================="<<endl; for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMax(double** p,int start,int end){int max=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[max][start]))max=i;}return max;}void swapRow(double** p,int one,int other,int col){double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMax(p,i,row); //找列主元行号if(p[flag][i]==0) return false;swapRow(p,i,flag,col); //交换行for(int j=i+1;j<row;j++){double elem=p[j][i]/p[i][i]; //消元因子 p[j][i]=0;for(int k=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}return true;}double sumRow(double** p,double* q,int row,int col){ double sum=0;for(int i=0;i<col-1;i++){if(i==row)continue;sum+=(q[i]*p[row][i]);}return sum;}void back(double** p,int row,int col,double* q){for(int i=row-1;i>=0;i--){q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; }}int main(){cout<<"Input n:";int n; //方阵的大小cin>>n;double **p=new double* [n];for(int i=0;i<n;i++){p[i]=new double [n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){ cout<<"奇异"<<endl; return 0;}double* q=new double[n]; for(int i=0;i<n;i++) q[i]=0;back(p,n,n+1,q);disp(q,n); delete[] q;for(int i=0;i<n;i++) delete[] p[i]; delete[] p; }1. 用列主元消去法解方程⎪⎩⎪⎨⎧=++-=++-=++035.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.1833.555x 2.304x 0.101x 321321321例2 解方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧=++++=++++=++++=++++=++++-12.041.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E 2.92D 1.46C3.53B -20.071.0F 1.10E 4.48D 1.21C 4.93B -32.041.0F 1.55E 5.66D 2.40C 8.77B 计算结果如下B=-1.461954 C= 1.458125 D=-6.004824E=-2.209018F= 14.719421矩阵直接三角分解法算法：将方程组A x=b 中的A 分解为A =LU ，其中L 为单位下三角矩阵，U 为上三角矩阵，则方程组A x=b 化为解2个方程组Ly =b ，Ux =y ，具体算法如下： ①对j=1,2,3,…,n 计算j j a u 11=对i=2,3,…,n 计算1111/a a l i i =②对k=1,2,3,…,n:a. 对j=k,k+1,…,n 计算∑-=-=11k q qj kq kj kj u l a ub. 对i=k+1,k+2,…,n 计算kk k q qk iq ik ik u u l a l /)(11∑-=-=③11b y =，对k=2,3,…,n 计算∑-=-=11k q q kq k k y l b y④nn n n u y x /=,对k=n-1,n-2,…,2,1计算kk nk q q kqk k u x uy x /)(1∑+=-=注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n n n na a a a a a a a a a a a施行算法②，③，此时U 的第n+1列元素即为y 。
(完整版)北京航空航天大学数值分析课程知识点总结.docx
1.2 误差知识与算法知识1.2.2 绝对误差、相对误差与有效数字设 a 是准确值 x 的一个近似值，记 ex a ，称 e 为近似值 a 的绝对误差，简称误差。
如果 |e |的一个上界已知，记为 ，即 | e |，则称 为近似值 a 的绝对误差限或绝对误差界，简称误差限或误差界。
记 e re x a，称 e r 为近似值 a 的相对误差。
由于 x 未知，实际上总把e作为 a 的xxae x ae 的上界，即 r相对误差，并且也记为 e r，相对误差一般用百分比表示。
aar| a |称为近似值 a 的相对误差限或相对误差界。
定义 设数 a 是数 x 的近似值。
如果 a 的绝对误差限时它的某一位的半个单位，并且从该位 到它的第一位非零数字共有 n 位，则称用 a 近似 x 时具有 n 位有效数字。
1.2.3 函数求值的误差估计~设 uf (x) 存在足够高阶的导数， a 是 x 的近似值， 则 uf (a) 是 u f (x) 的近似值。
~若 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) (a) 0 ，且比值~f( k)(a)ke(u)k! e( a)大，则有误差估计。
f ( k) (a)~k(u)(a)k !~nf (a 1, a 2,..., a n )e(a )e(u)i 1 x i i对于 n 元函数，有误差估计~nf ( a 1 ,a 2 ,..., a n )(u)(a i )i 1x if (k 1) (a) / f (k ) (a) 不很；若一阶偏导全为零或很小，则要使用高阶项。
1.2.4 算法及其计算复杂性（ 1）要有数值稳定性，即能控制舍入误差的传播。
- 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.6630353#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;double cosz;double sinz;double MAX;int kk;//#define PTS pts#ifdef PTSvoid PTS(double *m){printf("-----------------------------------------------------------------------\n");printf(" 迭代第%d次\n",kk);for(int i = 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(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 1 ; i <= PFrols+2 ; i++){printf(" ... ");}putchar(10);}#elsevoid PTS(double *m){}#endifvoid recounti(int i , 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;}}}void refreshMetrix(double *m){int ipr,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(int i = 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);}//void calCosSin(double *m){double app = m[(p+1)*p/2];double aqq = m[(q+1)*q/2];double apq = m[(p-1)*p/2+q];cosz = cos(atan(2 * apq / (app - aqq))/2);sinz = sin(atan(2 * apq / (app - aqq))/2); }//void find_pq(double *m){double max = 0.0;int pp = 0;int qq = 0;for(int i = 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;}void init(double *m){for(int i = 1 ; i <= N ;i++)m[(i+1)*i/2] = a(i);for(int i = 2 ; i <= N ; i++)m[(i-1)*i/2+i-1] = b;for(int i = 3 ; i <= N ; i++)m[(i-1)*i/2+i-2] = c;PTS(m);}void calFinal(double *m){printf("---------------------------------------------------------------------------------------------------\n");printf("结果输出:\n\n");double conda;double deta = 1.0;double minlumda = pow((double)10.0,12);double maxlumda = pow((double)10.0,-12);double absminlumda = pow((double)10.0,12);for(int i = 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(int i = 1 ; i <= FK ; i++){double muk = minlumda + i * (maxlumda - minlumda) / 40;double lumdak = 0.0;double tempabsmin = 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(int argc, _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(int i = 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(int i = 1 ; i <= PFrols+1 ; i++){printf(" ... ");}putchar(10);printf(" . .\n");printf(" . .\n");printf(" . .\n");for(int i = 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;}运行结果如下：中间运行状态如下：结果分析➢有效性分析1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；2.代码中有定义预编译宏//#define PTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G存，64位windows10操作系统）。