北航数值分析计算实习报告一

航空航天大学《数值分析》计算实习报告第一大题学院：自动化科学与电气工程学院专业：控制科学与工程学生姓名：学号：教师：电话：完成日期： 2015年11月6日航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。

矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1.求1λ，501λ和s λ的值。

2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。

3.求A 的(谱数)条件数2)A (cond 和行列式detA 。

说明：1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。

2.选择算法时，应使矩阵A 的所有零元素都不储存。

3.打印以下容： （1）全部源程序；（2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。

4.采用e 型输出实型数，并且至少显示12位有效数字。

一、算法设计方案1、求1λ，501λ和s λ的值。

由于||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤，可知绝对值最大特征值必为1λ和501λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。

将矩阵A 进行一下平移：I -A A'λ= （1）对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。

s λ为按模最小特征值，||min ||5011i i s λλ≤≤=，可对A 使用反幂法求得。

2、求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。

计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。

因此对A 进行平移变换：)39,,2,1k -A A k k ==（I μ （2）对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。

3、求A 的(谱数)条件数2)(A cond 和行列式detA 。

由矩阵A 为非奇异对称矩阵可得：||)(min max2λλ=A cond （3）其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。

在进行反幂法求解时，要对A 进行LU 分解得到。

因L 为单位下三角阵，行列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。

二、 算法实现1、矩阵存储原矩阵A 为一个上、下半带宽都为2的501×501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501×501的空间保存矩阵的话会浪费很多空间。

因此，为了节省存储量，A 的带外元素不给存储，值存储带元素，如下C 矩阵所示：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=000000501500499321cccc b b b b b a a a a a a b bb b bc c c c C（4） C 是一个5×501的矩阵，相比A 大大节省了存储空间，在数组C 中检索矩阵A 的带元素ij a 的方法是：j 1,s j -i c a A ++=中的元素的带内元素C ij （5）2、幂法幂法迭代公式如下：⎪⎪⎪⎩⎪⎪⎪⎨⎧====∈-------k T k kk k k k k u y Ay u u y u u 111k 11211k n0/R βηη任取非零向量 （6）其中λβ=∞→k k lim ，不断迭代当εβββ≤--||/||1k k k 时即可认为其满足精度要求，令k βλ=。

在程序中计算1u -=k k Ay 时，根据A 矩阵的特点，简化如下：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⋅⋅⋅=+++++-+-=++=+++=+++=++=)499,,3()2()1()()()1()2()()501()500()499()501()501()500()499()498()500()4()3()2()2(C )1()2()3()2()1(C )1(]][3[]501][3[]500][3[]2][3[]1][3[i i cy i by i y i C i by i cy i u y C by cy u by y C by cy u cy by y by u cy by y u i （7）3、反幂法反幂法迭代公式如下：⎪⎪⎪⎩⎪⎪⎪⎨⎧====∈-------kT k k k k k k k u y y Au u y u u 111k 11211k n0/R βηη任取非零向量 （8）当k 足够大时，kβλ1s ≈。

在求解1-=k k y Au 时，可先对A 进行Doolittle 分解，由于A 是带状结构，所以分解出的L 、U 也是带状结构，利用C 矩阵进行Doolittle 分解并求向量k u 的算法如下： （1）作分解A=LU对于n ,,2,1k =执行：)],min(,,2,1[/)(:)],min(,,1,[:,11),,1max(,1,1,1,11),,1max(,1,1,1,1n r k k k i c ccc c n s k k k j ccc c ks k s k r i t k s k t t s t i k s k i k s k i k s j r k t js j t t s t k j s j k j s j k +++=-=++=-=+---=++-++-++-++----=++-++-++-++-∑∑ （9）由于C 语言中数组下标是从0开始的，所以在程序中矩阵元素c 的下标都减1。

（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y,在程序中b 和y 都保存在数组y[501]中。

）)1,,2,1(/)(/),,3,2(,1),min(1,1,11),1max(,1 --=-===-=+++=++-+--=++-∑∑n n i c x cb xc b x n i b cb b i s tn s i i t t s t i i i ns n n i r i t tt s t i i i （10）求出k u 后，其他部分与幂法求解相同。

三、结果分析实验表明，本程序中，初始向量[]Tu u u u 501211501,,, =⨯对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。

（实验结果截图见附录）表1 不同初始向量对应的1λ和501λ及其迭代次数由表1可以得到如下结论：1. 不同的初始向量对本程序的1λ影响大，对501λ没有影响，都能保证收敛到正确值。

2. 初始向量中必须保证501462~u u 中至少有一个为1才能保证1λ收敛到正确值。

3. 初始向量非零值的多少和大小对迭代次数并没有明显影响。

4.为解决初始向量对程序的影响，可以先对A做平移变换再求1四、实验程序#include <stdio.h>#include <math.h>#include <stdlib.h>static double b=0.16,c=-0.064;#define Precision 1e-12void copy(double b[501],double y[501]);double dianji(double x[],double y[]);//计算两个向量积void InitMatrix(double *p);//Init Matrix Adouble NeiJi(double a[],double b[]);//get 2数void get_y(double *y,double *u);//get yvoid get_u(double *u,double *y,double *a);//get uvoid Initu(double *p);//初始化初始向量udouble Get_Fabs_Eigenvalue(double *a,double *u,int *iterations);//循环迭代得到绝对值最大特征值void A_sub_minI(double *a,double min);//A-min*Ivoid InitC(double C[5][501]);//初始化数组Cvoid DoolittleC(double C[5][501],int n,int s,int r);//进行Doolittle分解int min(int a,int b);//取返回a,b最小值int max(int a,int b);//求最大值void Doolittle_getx(double C[5][501],double y[501],double u[501],int n,int s,int r);double Get_min_Eigenvalue(double C[5][501],double *u,int n,int s,int r,int *iterations);double detA(double C[5][501]);//求A的行列式struct FinalValue{double min;//特征值最小值double max;//特征值最大值double abs_min;//模最小特征值double abs_max;//模最大特征值double detA;//A的行列式double cond2;//A的条件数int min_iterations;//最小值迭代次数int max_iterations;//最大值迭代次数int abs_min_iterations;//求绝对值最小的迭代次数};int main(){FinalValue main_num= {0,0,0,0};double temp;//两值交换中间变量int temp1;double u[501]={0};//,y[501];//为u0赋初值double Norm_u=0;//数double C[5][501] = {0};InitC(C);//初始化Cint i=0,*iterations;Initu(u);iterations = &main_num.min_iterations;main_num.min = Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.abs_max = main_num.min;A_sub_minI(C[2],main_num.min);for(i=0;i<501;i++)u[i]=i+1;iterations = &main_num.max_iterations;main_num.max = Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.max += main_num.min;if(main_num.min>main_num.max){temp = main_num.min;main_num.min = main_num.max;main_num.max = temp;temp1 = main_num.min_iterations;main_num.min_iterations = main_num.max_iterations;main_num.max_iterations = temp1;}printf("最小特征值为：%.12e,迭代次数为：%d\n",main_num.min,main_num.min_iterations);printf("最大特征值为：%.12e,迭代次数为：%d\n",main_num.max,main_num.max_iterations);/**********************************//*以下利用反幂法求解模最小的特征值*//**********************************/InitC(C);//初始化CInitu(u);DoolittleC(C,501,2,2);main_num.detA = detA(C);iterations = &main_num.abs_min_iterations;main_num.abs_min = Get_min_Eigenvalue(C,u,501,2,2,iterations);main_num.cond2 = fabs(main_num.abs_max/main_num.abs_min);printf("绝对值最小特征值为：%.12e,迭代次数为：%d\n",main_num.abs_min,main_num.abs_min_iterations);printf("A的行列式为：%.12e；",main_num.detA);printf("A的条件数cond(A)2为：%.12e\n",main_num.cond2);/**********************************************//*以下利用反幂法求与数列Uk中元素最相近的特征值*//**********************************************/double Uk[39];//保存Uk的值并保存与Uk[i]最接近的λdouble B[39];int diedai;iterations = &diedai;for( i=1;i<=39;i++){Uk[i-1] = main_num.min + i*(main_num.max - main_num.min)/40;InitC(C);Initu(u);for(int j =0;j<501;j++)C[2][j] -= Uk[i-1];DoolittleC(C,501,2,2);B[i-1] = Get_min_Eigenvalue(C,u,501,2,2,iterations);B[i-1] +=Uk[i-1];printf("μ%-2d=%-3.12e,与μ%-2d最相近的特征值λ=%-3.12e\n",i,Uk[i-1],i,B[i-1]);}return 0;}double detA(double C[5][501]){int i =0;double e =1;for(i =0;i<501;i++)e*=C[2][i];return e;}void InitMatrix(double *p){for(int i=1;i<=501;i++){*p=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);p++;}}void Initu(double *p){for(int i=1;i<=501;i++){if(i<=500)*p=0;else *p=1;p++;}// *p=1;}void A_sub_minI(double *a,double min) {for(int i=1;i<=501;i++){*a=*a-min;a++;}}double NeiJi(double *a,double *b) {double e=0.0;for(int i=0;i<501;i++){e=e+(*a)*(*b);a++;b++;}return e;}void get_y(double *y,double *u){ double Norm_u=sqrt(NeiJi(u,u));for(int i=0;i<501;i++){*y = *u/Norm_u;y++;u++;}}void get_u(double *u,double *y,double *a){u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3];for(int i=2;i<499;i++)u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];u[500]=c*y[498]+b*y[499]+a[500]*y[500];}void InitC(double C[5][501]){int i;for(i = 2;i < 501;i++){C[0][i] = -0.064;C[4][i-2] = -0.064;}for(i = 1;i < 501;i++){C[1][i] = 0.16;C[3][i-1] = 0.16;}for(i = 1;i <= 501;i++){C[2][i-1] = (1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);}}double Get_Fabs_Eigenvalue(double *a,double *u,int *iterations) {double y[501],B_k0=0,B_k1=0;double wucha;int i=0;while(1){++i;get_y(y,u);get_u(u,y,a);B_k1=NeiJi(y,u);//是否判断B_K1是否为0？wucha=(fabs(B_k1-B_k0))/(fabs(B_k1));//get wuchaif(wucha<=Precision)break;else if(i>10000){printf("迭代次数超长，请更改初始向量\n");break;}else B_k0 = B_k1;}*iterations = i;return B_k1;}double Get_min_Eigenvalue(double C[5][501],double *u,int n,int s,int r,int *iterations){double y[501] = {0},B_k0=0,B_k1=0;double b[501] = {0};//储存y值的中间向量double wucha;int i=0;while(1){++i;get_y(y,u);copy(b,y);//保护y向量Doolittle_getx(C,y,u,501,2,2);copy(y,b);B_k1=NeiJi(y,u);//是否判断B_K1是否为0？wucha=(fabs(1/B_k1-1/B_k0))/(1/fabs(B_k1));//get wucha// wucha=(fabs(B_k1-B_k0))/(fabs(B_k1));//get wuchaif(wucha<=Precision)break;else if(i>10000){printf("迭代次数超长，请更改初始向量\n");break;}else B_k0 = B_k1;}*iterations = i;return (1/B_k1);}void Doolittle_getx(double C[5][501],double y[501],double u[501],int n,int s,int r){int i,t;double e;for(i = 2;i <= n;i++){e = 0;for(t = max(1,i-r);t <= i-1;t++){e+=C[i-t+s+1-1][t-1]*y[t-1];}y[i-1] -= e;}u[n-1] = y[n-1]/C[s+1-1][n-1];for(i = n-1;i>=1;i--){e = 0;for(t = i+1;t<=min(i+s,n);t++)e += C[i-t+s+1-1][t-1]*u[t-1];u[i-1] = (y[i-1] - e)/C[s+1-1][i-1];}}void copy(double b[501],double y[501]){for(int i=0;i<501;i++){b[i] = y[i];}}void DoolittleC(double C[5][501],int n,int s,int r){int k,j,i,t;double e;for(k = 1;k <= n;k++){for(j = k;j <= min(k+s,n);j++){ e = 0;for(t = max(max(1,k-r),j-s);t <= k-1;t++){e = e + C[k-t+s+1-1][t-1]*C[t-j+s+1-1][j-1];}C[k-j+s+1-1][j-1] = C[k-j+s+1-1][j-1] - e;}if(k == n)break;for(i = k+1;i <= min(k+r,n);i++){ e = 0;for(t = max(max(1,i-r),k-s);t <= k-1;t++){e = e + C[i-t+s+1-1][t-1]*C[t-k+s+1-1][k-1];}C[i-k+s+1-1][k-1] = (C[i-k+s+1-1][k-1] - e)/C[s+1-1][k-1];}}}int min(int a,int b){if(a<b)return a;else return b;}int max(int a,int b){if(a>=b)return a;elsereturn b;}附录：部分实验程序截图1、[]Tu 1,,1,15011 =⨯2、[]Tu 501,,2,11501 =⨯3、[]Tu 0,,0,11501 =⨯4、[]Tu 1,0,,01501 =⨯5、[]Tu 0,0,1,11501 ，=⨯6、146121====u u u7、140021====u u u8、另14621===u u ，其余都为09、1462 u ，其余为010、1501462===u u ，其余为011、1481 u ，其余为0。