北航数值分析计算实习一

  • 格式:pdf
  • 大小:327.74 KB
  • 文档页数:22

下载文档原格式

  / 22
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

λ1=λm1 λ501=λm2
结束
图 1 求最大和最小特征值算法流程图
7. 程序中用到的子程序 (1) 幂法迭代 设置一维数组 uk[dim]和 yk[dim],分别存放每一步迭代产生的向量和其单位化的向 量。选取非零初始向量 u0,进行幂法迭代uk = A ∗ yk,每次迭代产生一个特征值的近似 值 b1 和其对应的单位特征向量 yk,直到特征值近似值 b1 满足精度要求,迭代结束,输 出 b1。 流程图如图 2 所示。 (2) 反幂法迭代 设置一维数组 uk[dim]和 yk[dim],分别存放每一步迭代产生的向量和其单位化的向 量。选取非零初始向量 u0,进行反幂法迭代uk = ������������−1 ∗ yk(用 Doolittle 分解法求解A ∗ uk = yk得到 uk) ,每次迭代产生一个特征值的近似值 b1 和其对应的单位特征向量 yk, 直到特征值近似值 b1 满足精度要求,迭代结束,输出 b1。 流程图如图 3 所示。
double mifadd2(double x[dim],double aa[r+s+1][dim],double e0) { b0=0,b1=0,e=1; long k=0; long i,j=0; for(j=0;j<dim;j++) //初始向量->u0 { uk[j]=x[j]; } //printf("第%d 次\n",k); //display(uk);
似值,通过原点平移,先用反幂法求 A-μI 的按模最小特征值λ‘������������ ,进而求出 A 的特征值 λ������������ = λ‘������������ + µ。 带原点平移的反幂法迭代流程图如图 4 所示。
开始
读取矩阵A参数 设置平移量p
对A进行原点平移A-pI A主对角线元素减p,其余不变
//幂法迭代
do { k++; unit(uk); //单位化 uk->yk //printf("特征值 %12.12e\n",b1); //printf("第%d 次\n",k); mtx_vtr2(aa,yk); //迭代 A*yk->u(k+1) b1=vtr_vtr(yk,uk); //求特征向量 b(k+1)=(yk,u(k+1)) if(k>1) { e=(b1-b0)/b0; e=fabs(e); } b0=b1; }while(e>e0); //printf("迭代次数:%d\n",k-1); //printf("\n\n%12.12e\n",b1); return b1; } double mifadd2_py(double x[dim],double aa[r+s+1][dim],double p,double e0) 平移的幂法迭代 { b0=0,b1=0,e=1; double apy[r+s+1][dim]={{0},{0}}; double b2=0; long k=0; long i,j=0; for(j=0;j<dim;j++) //初始向量->u0 { uk[j]=x[j]; } //printf("第%d 次\n",k); //display(uk); do { k++; unit(uk); //单位化 uk->yk //printf("特征值 %12.12e\n",b1); //带原点 //判断是否满足精度要求 //计算误差
最接近的特征值λ������������������������ (k=1,2,…,39)
cond(������������)2 =
�����������1 �����������
开始
幂法求A 按模最大特征值λm1
原点平移A+λm1I
反幂法求A 特征值λm2
λm1>λm2?
Y N
λ1=λm2 λ501=λm1
输出b1
结束
图 3 反幂法迭代流程图
(3) 原点平移 若 λ 为矩阵 A 的特征值,则 λ-p 是矩阵 A-pI 的特征值;反之,若 λ-p 是矩阵 A-pI 的 特征值,则 λ 为矩阵 A 的特征值,特征向量相同。 实际中常用带原点平移的反幂法求 A 的某个特征向量。 设 μ 为 A 的某个特征值的近
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;
数值分析计算实习(一)
姓名:张时毓 学号: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 = � ������������������������������������
开始
初始化 读取A的压缩矩阵、 相对误差精度e0 选取非零初始向量u0, 并存入uk 迭代次数k=0
k++
uk单位化—>yk
迭代产生下一uk A*yk—>uk 求特征值 b1=ykTuk
计算误差 e=(b1-b0)/b0
判断是否 满足精度要求 e<e0? Y N
取当前b1作为特征值Baidu Nhomakorabeam的近似值 yk作为λm对应的特征向量
输出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
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 //维数 //上半带宽 //下半带宽 //允许误差 //压缩的系数矩阵 //迭代初始向量 //迭代过程中特征向量序列 //迭代过程中特征向量序列(单位化) //迭代过程中特征值 //相对误差 //向量输出函数
������������=1 501
5. 求 A 的(谱范数)条件数cond(������������)2 6. 求 A 的与数������������������������ = λ1 + ������������
λ501 −λ1 40
作原点平移 A-μkI,求出 A-μkI 按模最小特征值λ = λ������������������������ − ������������������������ , λ������������������������ = λ + ������������������������
//计算向量的模
//向量单位化
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 向量相乘