北航数值分析1_Jacobi法计算矩阵特征值
北航研究生数值分析作业第一题

北航研究⽣数值分析作业第⼀题北航研究⽣数值分析作业第⼀题:⼀、算法设计⽅案1.要求计算矩阵的最⼤最⼩特征值,通过幂法求得模最⼤的特征值,进⾏⼀定判断即得所求结果;2.求解与给定数值接近的特征值,可以该数做漂移量,新数组特征值倒数的绝对值满⾜反幂法的要求,故通过反幂法即可求得;3.反幂法计算时需要⽅程求解中间过渡向量,需设计Doolite分解求解;4.|A|=|B||C|,故要求解矩阵的秩,只需将Doolite分解后的U矩阵的对⾓线相乘即为矩阵的Det。
算法编译环境:vlsual c++6.0需要编译函数:幂法,反幂法,Doolite分解及⽅程的求解⼆、源程序如下:#include#include#include#includeint Max(int value1,int value2);int Min(int value1,int value2);void Transform(double A[5][501]);double mifa(double A[5][501]);void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]); double Det(double A[5][501]);/***定义2个判断⼤⼩的函数,便于以后调⽤***/int Max(int value1,int value2){return((value1>value2)?value1:value2);}int Min(int value1,int value2){return ((value1}/*****************************************//***将矩阵值转存在⼀个数组⾥,节省空间***/void Transform(double A[5][501],double b,double c){int i=0,j=0;A[i][j]=0,A[i][j+1]=0;for(j=2;j<=500;j++)A[i][j]=c;i++;j=0;A[i][j]=0;for(j=1;j<=500;j++)A[i][j]=b;i++;for(j=0;j<=500;j++)A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++;for(j=0;j<=499;j++)A[i][j]=b;A[i][j]=0;i++;for(j=0;j<=498;j++)A[i][j]=c;A[i][j]=0,A[i][j+1]=0;}/***转存结束***///⽤于求解模最⼤的特征值,幂法double mifa(double A[5][501]){int s=2,r=2,m=0,i,j;double b2,b1=0,sum,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+A[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=exp(-12));return b2;}//带状DOOLITE分解,并且求解出⽅程组的解void daizhuangdoolite(double A[5][501],double x[501],double b[501]) { int i,j,k,t,s=2,r=2;double B[5][501],c[501];for(i=0;i<=4;i++){for(j=0;j<=500;j++)B[i][j]=A[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; }for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//⽤于求解模最⼤的特征值,反幂法double fanmifa(double A[5][501]){int s=2,r=2,m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);daizhuangdoolite(A,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)>=fabs(b1)*exp(-12));return 1/b2;}//⾏列式的LU分解,U的主线乘积即位矩阵的DET double Det(double A[5][501]) {int i,j,k,t,s=2,r=2;for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k];A[i-k+s][k]=A[i-k+s][k]/A[s][k];}}double det=1;for(i=0;i<=500;i++)det*=A[s][i];return det;}void main(){double b=0.16,c=-0.064,p,q;int i,j;double A[5][501];Transform(A,b,c); //进⾏A的赋值cout.precision(12); //定义输出精度double lamda1,lamda501,lamdas;double k=mifa(A);if(k>0) //判断求得最⼤以及最⼩的特征值.如果K>0,则它为最⼤特征值值,//并以它为偏移量再⽤⼀次幂法求得新矩阵最⼤特征值,即为最⼤ //与最⼩的特征值的差{lamda501=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda1=mifa(A)+lamda501;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}else //如果K<=0,则它为最⼩特征值值,并以它为偏移量再⽤⼀次幂法//求得新矩阵最⼤特征值,即为最⼤与最⼩的特征值的差{lamda1=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda501=mifa(A)+lamda1;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}lamdas=fanmifa(A);FILE *fp=fopen("result.txt","w");fprintf(fp,"λ1=%.12e\n",lamda1);fprintf(fp,"λ501=%.12e\n",lamda501);fprintf(fp,"λs=%.12e\n\n",lamdas);fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n");for(i=1;i<=39;i++) //反幂法求得与给定值接近的特征值{p=lamda1+(i+1)*(lamda501-lamda1)/40;for(j=0;j<=500;j++)A[2][j]=A[2][j]-p;q=fanmifa(A)+p;for(j=0;j<=500;j++)A[2][j]=A[2][j]+p;fprintf(fp,"µ%d: %.12e λi%d: %.12e\n",i,p,i,q);}double cond=fabs(mifa(A)/fanmifa(A));double det=Det(A);fprintf(fp,"\ncond(A)=%.12e\n",cond);fprintf(fp,"\ndetA=%.12e\n",det);}三、程序运⾏结果λ1=-1.069936345952e+001λ501=9.722283648681e+000λs=-5.557989086521e-003要求接近的值实际求得的特征值µ1: -9.678281104107e+000 λi1: -9.585702058251e+000µ2: -9.167739926402e+000 λi2: -9.172672423948e+000µ3: -8.657198748697e+000 λi3: -8.652284007885e+000µ4: -8.146657570993e+000 λi4: -8.0934********e+000µ5: -7.636116393288e+000 λi5: -7.659405420574e+000µ6: -7.125575215583e+000 λi6: -7.119684646576e+000µ7: -6.615034037878e+000 λi7: -6.611764337314e+000µ8: -6.104492860173e+000 λi8: -6.0661********e+000µ9: -5.593951682468e+000 λi9: -5.585101045269e+000µ10: -5.0834********e+000 λi10: -5.114083539196e+000µ11: -4.572869327058e+000 λi11: -4.578872177367e+000µ12: -4.062328149353e+000 λi12: -4.096473385708e+000µ13: -3.551786971648e+000 λi13: -3.554211216942e+000µ14: -3.0412********e+000 λi14: -3.0410********e+000µ15: -2.530704616238e+000 λi15: -2.533970334136e+000µ16: -2.020*********e+000 λi16: -2.003230401311e+000µ17: -1.509622260828e+000 λi17: -1.503557606947e+000µ18: -9.990810831232e-001 λi18: -9.935585987809e-001µ19: -4.885399054182e-001 λi19: -4.870426734583e-001µ20: 2.200127228676e-002 λi20: 2.231736249587e-002µ21: 5.325424499917e-001 λi21: 5.324174742068e-001µ22: 1.043083627697e+000 λi22: 1.052898964020e+000µ23: 1.553624805402e+000 λi23: 1.589445977158e+000µ24: 2.064165983107e+000 λi24: 2.060330427561e+000µ25: 2.574707160812e+000 λi25: 2.558075576223e+000µ26: 3.0852********e+000 λi26: 3.080240508465e+000µ27: 3.595789516221e+000 λi27: 3.613620874136e+000µ28: 4.106330693926e+000 λi28: 4.0913********e+000µ29: 4.616871871631e+000 λi29: 4.603035354280e+000µ30: 5.127413049336e+000 λi30: 5.132924284378e+000µ31: 5.637954227041e+000 λi31: 5.594906275501e+000µ32: 6.148495404746e+000 λi32: 6.080933498348e+000µ33: 6.659036582451e+000 λi33: 6.680354121496e+000µ34: 7.169577760156e+000 λi34: 7.293878467852e+000µ35: 7.680118937861e+000 λi35: 7.717111851857e+000µ36: 8.190660115566e+000 λi36: 8.225220016407e+000µ37: 8.701201293271e+000 λi37: 8.648665837870e+000µ38: 9.211742470976e+000 λi38: 9.254200347303e+000µ39: 9.722283648681e+000 λi39: 9.724634099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、分析如果初始向量选择不当,将导致迭代中X1的系数等于零.但是,由于舍⼊误差的影响,经若⼲步迭代后,.按照基向量展开时,x1的系数可能不等于零。
《实验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_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
雅克比法解特征值

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
北航数值分析计算实习第一题编程

i − t + s +1,t t − k + s +1, k t = max(1,i − r ,k − s )
∑c
c
) / cs +1, k
[i = k + 1, k + 2,⋯ , min( k + r , n); k < n]
(2) 求解 Ly = b,Ux = y (数组 b 先是存放原方程右端向量,后来存放中间向量 y)
0 b a2
b c
c b a3 b c
⋯ ⋯ ⋯ ⋯ ⋯
c b a499 b c
c b a500 b 0
c ⎤ b ⎥ ⎥ a501 ⎥ ⎥ 0 ⎥ 0 ⎥ ⎦
在数组 C 中检索矩阵 A 的带内元素 aij 的方法是: A 的带内元素 aij =C 中的元素 ci − j + s +1, j
2
数值分析计算实习题目一
i −1
bi := bi −
பைடு நூலகம்
i − t + s +1,t t t = max(1,i − r )
∑c
b
(i = 2,3,⋯ , n)
xn := bn / cs +1, n
min( i + s )
xi := (bi −
t = i +1
∑c
i −t + s +1,t t
x ) / cs +1,i
(i = n − 1, n − 2,⋯ ,1)
3、Doolittle 分解求解 n 元带状线性方程组(doolittle()函数)
按照上述对带状矩阵 A 的存储方法和元素 aij 的检索方法,并且把三角分解的结果 ukj 和 lik 分 别存放在 akj 和 aik 原先的存储单元内,那么用 Doolittle 分解法求解 n 元带状线性方程组的算法 可重新表述如下(其中“:=”表示赋值) : (1) 作分解 A = LU 。 对于 k=1,2, ……,n 执行
北京航空航天大学数值分析课程知识点总结

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方法的原理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。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值

《数值分析》计算实习题目第一题: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)作为初始向量进行迭代,可得出以上结果。
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。
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方法

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.
§3 计算实对称矩阵特征值的Jacobi方法

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方法求正交矩阵的特征值与特征向量

用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。
北航数值分析第一次大作业

一、算法的设计方案1、求矩阵最大特征值,最小特征值与按模最小特征值的方法首先用幂法求出矩阵A 的一个特征值λ,则其必为最大特征值与最小特征值二者其一,之后对矩阵A 进行一次移项,即A-λI ,然后再次用幂法求出另一个按模最大特征值,再比较这两个值的大小,则较大的为矩阵A 的最大特征值,较小的为矩阵A 的最小特征值。
用反幂法可以求得矩阵的按模最小特征值λs 2、求矩阵A 与k μ最接近的特征值k i λ可以先对矩阵A 进行移项,即A-k μI ,对这个移项后的矩阵用反幂法求出按模最小的特征值,然后再加上k μ,就求出所要求的k i λ。
3、求矩阵A 的条件数cond(A)2和行列式detA由于矩阵A 是非奇异的实对称矩阵,所以可以用以下公式方便地求出矩阵A 的条件数cond(A)2=sλλ501对于矩阵A 行列式的求法也比较简单。
由于在用反幂法的过程中对A 进行了Doolittle LU 分解,所以detA=detL*detU ,而detL=1,detU 可以用对角线元素相乘方便地算出,所以detA 就是U 阵对角线元素的乘积。
4、几点说明由于A 中的零元素都不存储,所以在存储矩阵的时候采用书上26页的压缩存储方式。
在反幂法中采用LU 分解求解带状线性方程组的算法来求解每一次迭代的方程组,由于每一次方程左边的系数都相同,所以只要进行一次LU 分解即可。
因为幂法,反幂法,LU 分解,求最大值与最小值在程序编写的过程中多次用到,所以这几项作为子函数单独进行编写。
二、源程序如下:#include "stdio.h"#include "math.h"# define s 2# define r 2# define N 501double c[5][N]={0};double lameda[40];double max(double x,double y);double min(double x,double y);double mifa();double fanmifa();void LUfenjie();void main(){int i=0,j=0;/*============对数组进行赋值==============*/ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;/*========幂法求最大和最小特征值==============*/ double a=mifa();for(j=1;j<=N;j++)c[2][j-1]-=a;double b=mifa()+a;double lameda501=max(a,b);double lameda1=min(a,b);printf("矩阵A最大的特征值=%13.11e\n",lameda501);printf("矩阵A最小的特征值=%13.11e\n",lameda1);/*========反幂法求绝对值最小特征值===========*/for(j=1;j<=N;j++)c[2][j-1]+=a;double lamedas=fanmifa();printf("矩阵A按模最小的特征值=%13.11e\n",lamedas);/*========求条件数和行列式的值===========*/double detA=1;for(j=1;j<=N;j++)detA*=c[2][j-1];printf("矩阵A的行列式=%13.11e\n",detA);double condA=fabs(lameda501/lamedas);printf("矩阵A的条件数=%13.11e\n",condA);/*========反幂法求与uk最接近的特征值========*/for(int k=1;k<40;k++){ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;for(j=1;j<=N;j++)c[2][j-1]-=(lameda1+k*(lameda501-lameda1)/40);lameda[k]=fanmifa()+(lameda1+k*(lameda501-lameda1)/40);printf("矩阵A最接近u%d的特征值=%13.11e\n",k,lameda[k]);}}double max(double x,double y) //求两数中的最大值{double z;z=x>y ? x:y;return(z);}double min(double x,double y) //求两数中的最小值{double z;z=x<y ? x:y;return(z);}double mifa() //幂法求按模最大特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b1=0,b2=0;int i=0,j=0;for(i=0;i<N;i++)u[i]=1;do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<N;i++){ u[i]=0;for(j=max(0,i-2);j<=min(i+2,N-1);j++){u[i]+=c[i-j+s][j]*y[j];}}b2=b1;zero=0;for(i=0;i<N;i++)zero+=y[i]*u[i];b1=zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}double fanmifa() //反幂法求按模最小特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b[N]={0};double b1=0,b2=0;int i=0,j=0,t=0;for(i=0;i<N;i++)u[i]=1;LUfenjie();do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++){b[i]=u[i]/sqrt(sum);y[i]=b[i];}for(i=1;i<N;i++)for(t=max(0,i-r);t<i;t++)y[i]-=c[i-t+s][t]*y[t];u[N-1]=y[N-1]/c[s][N-1];for(i=N-2;i>=0;i--){ for(t=i+1;t<=min(i+s,N-1);t++)y[i]-=c[i-t+s][t]*u[t];u[i]=y[i]/c[s][i];}b2=b1;zero=0;for(i=0;i<N;i++)zero+=b[i]*u[i];b1=1/zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}void LUfenjie() //对矩阵做LU分解{ int k=0;int j=0;int i=0;int t=0;for(k=0;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++)for(t=max(max(0,k-r),j-s);t<k;t++)c[k-j+s][j]-=c[k-t+s][t]*c[t-j+s][j];for(i=k+1;i<=min(k+r,N-1);i++){ for(t=max(max(0,i-r),k-s);t<k;t++)c[i-k+s][k]-=c[i-t+s][t]*c[t-k+s][k];c[i-k+s][k]/=c[s][k];}}}三、运行结果如下四、初始向量对计算结果的影响在本程序的编写中,为了方便起见,所以迭代向量的初值选为u=[1 1 1 ...1]。
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迭代法以其简洁、直观的特点,在复矩阵特征值和特征向量的求解中得到了广泛的应用。
北航数值分析-实习作业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,再重复上述步骤即可求得两个特征值。
- 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操作系统)。