北航数值分析计算实习报告一————————————————————————————————作者:————————————————————————————————日期:北京航空航天大学《数值分析》计算实习报告第一大题学院:自动化科学与电气工程学院专业: 控制科学与工程学生姓名:学号:教师:电话:完成日期: 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 和行列式d etA 。
说明: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 和行列式det A。
由矩阵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 矩阵所示:⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=000000501500499321cc cc b b b b b a a a a a a b b b b b c c c cC (4) C 是一个5×501的矩阵,相比A 大大节省了存储空间,在数组C 中检索矩阵A 的带内元素ij a 的方法是:j 1,s j -i c a A ++=中的元素的带内元素C ij (5)2、幂法幂法迭代公式如下:⎪⎪⎪⎩⎪⎪⎪⎨⎧====∈-------kT k k k 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 进行D ool it tl e分解,由于A 是带状结构,所以分解出的L 、U也是带状结构,利用C矩阵进行Dool it tle分解并求向量k u 的算法如下: (1)作分解A=L U ﻩ 对于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是不同初始向量的情况下的得到的部分结果。
(实验结果截图见附录)1501⨯u1λ迭代次数501λ迭代次数 11=u000370809810853.2-+e159 000787246340993.9+e381 121==u u000e 380809810853.2-+160 000787246340993.9+e381 14001===u u000962085531594.9-+e535 000087246340989.9+e57314611===u u000659539789789.9-+e350 000407246340988.9+e592 14621===u u001509.0700113611-+e674 000727246340987.9+e6111462=u 001502.0700113611-+e311 000727246340987.9+e611 1481=u001501.0700113611-+e304 000727246340987.9+e6111501462===u u001502.0700113611-+e343 000727246340987.9+e611 15011===u u001502.0700113611-+e343 000727246340987.9+e611 []T 501,,2,1001502.0700113611-+e343 000727246340987.9+e611表1 不同初始向量对应的1λ和501λ及其迭代次数由表1可以得到如下结论:1. 不同的初始向量对本程序的1λ影响大,对501λ没有影响,都能保证收敛到正确值。
2. 初始向量中必须保证501462~u u 中至少有一个为1才能保证1λ收敛到正确值。
3. 初始向量非零值的多少和大小对迭代次数并没有明显影响。
4. 为解决初始向量对程序的影响,可以先对A 做平移变换再求1λ。
四、实验程序#in clude <st dio.h>#i nclude <mat h.h> #i ncl ude <stdl ib.h>static double b=0.16,c=-0.064; #define P recision 1e-12void copy (doub le b[501],dou bl e y[501]);doubl e dianji (dou ble x[],doubl e y[]);//计算两个向量内积 voi d InitMatri x(double *p );//Init Matrix A doubl e NeiJi(dou ble a[],doubl e b[]);//ge t 2范数 voi d g et_y(d ou ble *y,doubl e *u);//g et y void get _u(dou ble *u ,do uble *y,double *a);//get u void I nitu(doubl e *p);//初始化初始向量udoubl e Ge t_Fa bs_Eigen value(dou bl e *a,double *u ,in t *it era ti on s);//循环迭代得到绝对值最大特征值void A_sub_minI(doubl e *a,doub le min );//A-mi n*I void InitC (do ub le C[5][501]);//初始化数组Cvoid Doolit tleC(double C[5][501],in t n,int s ,int r);//进行Doolit tle分解 i nt m in(int a ,int b );//取返回a,b 最小值 in t m ax(int a,i nt b);//求最大值voi d Dool ittle_getx(double C[5][501],dou ble y[501],double u[501],int n ,in t s,int r);dou ble G et_min_Eigenva lue(dou ble C[5][501],dou ble *u ,int n ,int s,i nt r,int*iterations);double detA(double C[5][501]);//求A的行列式structFinalValue{doublemin;//特征值最小值ﻩdouble max;//特征值最大值doubleabs_min;//模最小特征值ﻩdouble abs_max;//模最大特征值ﻩdouble detA;//A的行列式ﻩdouble cond2;//A的条件数ﻩint min_iterations;//最小值迭代次数int max_iterations;//最大值迭代次数int abs_min_iterations;//求绝对值最小的迭代次数};intmain(){ﻩﻩFinalValue main_num= {0,0,0,0};double temp;//两值交换中间变量ﻩinttemp1;ﻩdoubleu[501]={0};//,y[501];//为u0赋初值doubleNorm_u=0;//范数ﻩdoubleC[5][501]= {0};InitC(C);//初始化Cﻩint 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.m ax_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];ﻩintdiedai;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;}voidA_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(inti=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++;}}voidget_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(doubleC[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 wuchaﻩif(wucha<=Precision)break;ﻩﻩelseif(i>10000)ﻩﻩ{ﻩﻩprintf("迭代次数超长,请更改初始向量\n");ﻩﻩbreak;ﻩ}ﻩelse B_k0 =B_k1;}ﻩ*iterations =i;ﻩreturn B_k1;}double Get_min_Eigenvalue(doubleC[5][501],double *u,int n,int s,int r,int *iterations){ﻩdouble y[501] ={0},B_k0=0,B_k1=0;doubleb[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 wuchaﻩﻩif(wucha<=Precision)break;ﻩﻩelse if(i>10000)ﻩ{ﻩprintf("迭代次数超长,请更改初始向量\n");ﻩbreak;ﻩﻩ}else B_k0 =B_k1;ﻩ}ﻩ*iterations =i;return(1/B_k1);}void Doolittle_getx(doubleC[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],doubley[501]){for(int i=0;i<501;i++)ﻩ{ﻩb[i] = y[i];ﻩ}}voidDoolittleC(doubleC[5][501],int n,int s,int r){intk,j,i,t;doublee;ﻩ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];}ﻩ}}intmin(inta,int b){if(a<b)ﻩreturn a;ﻩelsereturn b;}int max(int a,intb){ﻩif(a>=b)ﻩﻩreturn a;elsereturn b;}附录:部分实验程序截图1、[]T u 1,,1,15011 =⨯2、[]Tu 501,,2,11501 =⨯3、[]T u 0,,0,11501 =⨯4、[]T u 1,0,,01501 =⨯5、ﻬ[]T u 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。