北航数值分析计算实习一

开始
初始化 读取A的压缩矩阵、 相对误差精度e0 选取非零初始向量u0， 并存入uk 迭代次数k=0
k++
uk单位化—>yk
迭代产生下一uk A*yk—>uk 求特征值 b1=ykTuk
计算误差 e=(b1-b0)/b0
判断是否 满足精度要求 e<e0? Y N
取当前b1作为特征值λm的近似值 yk作为λm对应的特征向量
//计算向量的模
//向量单位化
void mtx_vtr2(double aa[r+s+1][dim],double x[dim]) { double mv[dim]={0}; long i=0,j=0;
//幂法迭代过程中 A*y->u
mv[0]=aa[s][0]*x[0]+b*x[1]+c*x[2]; mv[1]=aa[s][1]*x[1]+b*(x[0]+x[2])+c*x[3]; for(i=2;i<dim-2;i++) //计算 3-499 { mv[i]=aa[s][i]*x[i]+b*(x[i-1]+x[i+1])+c*(x[i-2]+x[i+2]); } mv[dim-2]=aa[s][dim-2]*x[1]+b*(x[dim-3]+x[dim-1])+c*x[dim-4]; mv[dim-1]=aa[s][dim-1]*x[dim-1]+b*x[dim-2]+c*x[dim-3]; for(i=0;i<dim;i++) { uk[i]=mv[i]; } } double vtr_vtr(double x1[dim],double x2[dim]) { double vv=0; long i=0; for(i=0;i<dim;i++) { vv=vv+x1[i]*x2[i]; } return vv; } //1*n 向量与 n*1 向量相乘
double acc[s+r+1][dim]={{0},{0}}; double u0[dim]={0}; double b=0.16,c=-0.064; double yk[dim]={0}; double uk[dim]={0}; double b0=0,b1=0; double e=1; void display(double x[dim]) { long i=0;
long j=0; for(i=0;i<dim;i++) { printf("%12.12e\t",x[i]); } printf("\n"); } void display_pq(double aa[r+s+1][dim]) { long i=0; long j=0; for(i=0;i<r+s+1;i++) { for(j=0;j<dim;j++) { printf("%f\t",aa[i][j]); } printf("\n"); } printf("\n"); } double length(double x[dim]) { double leng=0; long i=0; for(i=0;i<dim;i++) { leng=leng+x[i]*x[i]; } leng=sqrt(leng); return leng; } void unit(double x[dim]) { long i=0; double ll=length(x); for(i=0;i<dim;i++) { yk[i]=x[i]/ll; } } //矩阵输出函数
进行反幂法迭代 求A-pI的近似特征值λ’s
求A的近似特征值 λs=λ’s+p
结束
图 4 带原点平移的反幂法迭代流程图
二、源代码 #include<stdio.h> #include<math.h> #define dim 501 #define s 2 #define r 2 #define er0 1E-12 //维数 //上半带宽 //下半带宽 //允许误差 //压缩的系数矩阵 //迭代初始向量 //迭代过程中特征向量序列 //迭代过程中特征向量序列（单位化） //迭代过程中特征值 //相对误差 //向量输出函数
输出b1
结束
图 2 幂法迭代流程图
开始
初始化 读取A的压缩矩阵、 相对误差精度e0 选取非零初始向量u0， 并存入uk 迭代次数k=0
k++
uk单位化—>yk
用Doolittle分解法求解A*uk=yk 迭代产生下一uk
求特征值 b1=1/(ykTuk)
计算误差 e=(b1-b0)/b0
判断是否 满足精度要求 e<e0? Y 取当前b1作为特征值λs的近似值 yk作为λs对应的特征向量 N
输出b1
结束
图 3 反幂法迭代流程图
（3） 原点平移 若 λ 为矩阵 A 的特征值，则 λ-p 是矩阵 A-pI 的特征值；反之，若 λ-p 是矩阵 A-pI 的 特征值，则 λ 为矩阵 A 的特征值，特征向量相同。 实际中常用带原点平移的反幂法求 A 的某个特征向量。 设 μ 为 A 的某个特征值的近
数值分析计算实习（一）
姓名：张时毓 学号：SY1403531
一、设计方案 1． A 为大型带状矩阵，维数 dim=501，上半带宽 s=2，下半带宽 r=2。 将 A 压缩储存，设置二维数组 acc[r+s+1][dim]储存 A 的带内元素，acc 的第 j 列存放 A 的 第 j 列带内元素，并使 A 的主对角线元素存放在 acc 的第 s+1 行中，多出的单元取零。 2． 求最小特征值 λ1 和最大特征值 λ501 先用幂法迭代求出按模最大特征值 λm1，则 λm1 为最小特征值 λ1 或最大特征值 λ501。 （1） 若λ������������1 = λ1 ，则λ������������1 < 0，因为|λ501 | < |λ������������1 |，所以 λ1 + λ������������1 ≤ λ2 + λ������������1 ≤ ⋯ ≤ λ501 + λ������������1 < 0 所以 |λ1 + λ������������1 | ≥ |λ2 + λ������������1 | ≥ ⋯ ≥ |λ501 + λ������������1 | 用原点平移的反幂法求 A+λm1I 按模最小的特征值λ’������������2 = λ501 + λ������������1 则λ������������2 = λ501 = λ’������������2 −λ������������1 则λ1 = λ������������1 ，λ501= λ������������2 （2） 若λ������������1 = λ501 ，则λ������������1 > 0，因为|λ1 | < |λ������������1 |，所以 0 < λ1 + λ������������1 ≤ λ2 + λ������������1 ≤ ⋯ ≤ λ501 + λ������������1 所以 |λ1 + λ������������1 | ≤ |λ2 + λ������������1 | ≤ ⋯ ≤ |λ501 + λ������������1 | 用原点平移的反幂法求 A+λm1I 按模最小的特征值λ’������������2 = λ1 + λ������������1 则λ������������2 = λ1 = λ’������������2 −λ������������1 则λ1 = λ������������2 ，λ501= λ������������1 综上所述，先用幂法迭代求出按模最大特征值 λm1，再通过原点平移 A+λm1I 用反幂法求 特征值 λm2，比较 λm1 和 λm2 的大小，大的为最大特征值 λ501，小的为最小特征值 λ1。 流程图如图 1 所示。 同理可用 A-λm1I 作原点平移，再用幂法求取 λm2。 3． 求按模最小特征值 用反幂法，求取 A 的按模最小特征值 λs。 4． 求 A 的行列式 detA 将 A 做 Doolittle 分解， det A = detL ∗ detU 即 U 对角线元素的乘积 det A = � ������������������������������������
最接近的特征值λ������������������������ （k=1，2，…，39）
cond(������������)2 =
�����������1 �����������
开始
幂法求A 按模最大特征值λm1
原点平移A+λm1I
反幂法求A 特征值λm2
λm1>λm2?
Y N
λ1=λm2 λ501=λm1
������������=1 501
5． 求 A 的（谱பைடு நூலகம்数）条件数cond(������������)2 6． 求 A 的与数������������������������ = λ1 + ������������
λ501 −λ1 40
作原点平移 A-μkI，求出 A-μkI 按模最小特征值λ = λ������������������������ − ������������������������ ， λ������������������������ = λ + ������������������������