数值分析计算实习作业一

数值分析计算实习题一学号：：院系：2015年11月5日一、分析1.1算法分析题目要求求出：1）特征值从小到大排列的最小特征值1λ和最大特征值501λ。

2）特征值中模最小的特征值s λ。

3）靠近一组数k μ的一组特征值k i λ。

4）矩阵A 的条件数cond(A)2。

5）行列式detA 。

解决方法：1）若将所有行列式按模的大小排列则模最大的特征值一定是1λ和501λ中的一个，因此利用幂法求出模最大的特征值1m λ。

然后利用带原点平移的幂法，将系数矩阵变为1m A I λ-即将所有特征值都减去1m λ，则特征值按大小顺序排列的次序不变，模最大的特征值依然在整个排列的两端，再用一次幂法得到模最大的特征值21=m m λλλ-，其中λ为带原点平移的幂法求出的特征值，最后两个特征值1m λ、2m λ比较大小，大的为501λ，小的为1λ。

2）因为s λ为按模最小的特征值，因此用反幂法可求的其特征值。

3）因为k i λ靠近数k μ，因此k i k λμ-一定是所有的k λμ-中模最小的，因此可利用带原点平移的反幂法求出特征值k i λ，此时的系数矩阵变为k A I μ-。

4）条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值，因此利用1和2中求出的1m λ和s λ可解出条件数。

5）可对矩阵A 进行LU 分解，即A LU =则det()det()det()A L U =⨯，又因为矩阵L 对角线元素为1，则det()L =1，所以det()det()A U =，U 为上三角阵，行列式为对角线元素的乘积，因此可得A 的行列式。

1.2程序分析1.2.1 因为A 为拟三角阵，储存时零元素不储存，因此将矩阵A 压缩为5*501的矩阵CA 的带元素ij a =C 中的元素1,i j s j c -++ 程序中A[5][501]即为压缩后的矩阵。

1.1.2 程序中的B[5][501]为过渡矩阵，在幂法迭代、反幂法迭代以及LU 分解中均用矩阵B 来计算，计算之间对B 进行适当的赋值。

先令B=A ，调用幂法函数求1m λ；再令B=A-1m λI ，调用幂法函数求出2m λ；再令B=A ，调用反幂法函数求s λ；在令B=k A I μ-，调用反幂法函数求k i λ，最后令B=A ，将B 进行LU 分解，计算A 的行列式。

1.2.3幂法求解过程： 1）取初始向量u 0=[1,1,…1]; 2）进行k 次迭代，k=1,2,3…3）对每一次迭代，计算上一次迭代中的(1)k u -的模（2数）1k η-=(1)k u -单位化后赋值给1k y -，即1(1)1/k k k y u η---=，计算矩阵B 与向量1k y -的乘积，1k k u By -=，计算1k y -与k u 的积1T k k k y u β-=。

4）如果两次相邻两次迭代的β满足：1/k k k βββε--≤（允许误差），则结束迭代，k β的值可认为是矩阵B 对应的模最大的特征值。

如果不满足误差条件，则重复3）的迭代直到达到误差允许值。

1.2.4 反幂法求解过程：1）取初始向量u 0=[1,1,…1];对矩阵B 进行LU 分解。

2）进行k 次迭代，k=1,2,3…3）对每一次迭代，计算上一次迭代中的(1)k u -的模（2数）1k η-=(1)k u -单位化后赋值给1k y -，即1(1)1/k k k y u η---=，利用LU 分解后的矩阵B ，求线性方程组1k k Bu y -=，计算1k y -与k u 的积1T k k k y u β-=。

4）如果两次相邻两次迭代的β满足：1/k k k βββε--≤（允许误差），则结束迭的值可认为是矩阵B对应的模最大的特征值。

如果不满足误差条件，则代，k重复3）的迭代直到达到误差允许值。

1.2.5 拟三角矩阵的LU分解过程可参照数值分析书26页的算法得到。

二、程序整个数值分析在c语言中的程序如下：#include <stdio.h>#include <math.h>double mifa();double fanmifa();void fenjie();//声明三个函数，mifa是用幂法求模最大的特征值，fanmifa 是用反幂法求模最小的特征值，fenjie是对系数矩阵做LU分解double B[5][501]={0};double A[5][501]={0};//定义两个系数矩阵为全局变量,B作为求解过程中的系数矩阵使用double sigma=1e-12;//定义误差允许值/********主程序********/void main(){int i,j,t;double lamk,lamm,min,max,miu,lam,x,lams;double cond2;double det=1;for(i=0;i<=500;i++)//输入矩阵A{j=i+1;x=0.1/j;A[2][i]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(x);A[3][i]=0.16;A[4][i]=-0.064;}for(i=1;i<=500;i++){A[1][i]=0.16;}for(i=2;i<=500;i++){A[0][i]=-0.064;}A[3][500]=0;A[4][500]=0;A[4][499]=0;for(i=0;i<=4;i++)//先使B=A{for(j=0;j<=500;j++){B[i][j]=A[i][j];}}lamk=mifa();//幂法求模最大的特征值for(i=0;i<=4;i++)//赋值计算B=（A-lamk*I）{for(j=0;j<=500;j++){if(i==2) B[i][j]=A[i][j]-lamk;else B[i][j]=A[i][j];}}lamm=mifa()+lamk;//带原点位移的幂法求位移后模最大的特征值if(lamk<lamm) {min=lamk;max=lamm;}else {min=lamm;max=lamk;}printf("lam1=%13.11le\n",min);//显示12位有效数字printf("lam501=%13.11le\n",max);for(i=0;i<=4;i++)//赋值B=A{for(j=0;j<=500;j++){B[i][j]=A[i][j];}}lams=fanmifa();printf("lams=%13.11le\n",lams);//反幂法求模最小的特征值for(t=1;t<=39;t++)//循环求与miu接近的特征值{miu=min+t*(max-min)/40;for(i=0;i<=4;i++)//赋值B=(A-miu*I){for(j=0;j<=500;j++){if(i==2) B[i][j]=A[i][j]-miu;else B[i][j]=A[i][j];}}lam=fanmifa()+miu;//用带原点的反幂法求与miu接近的特征值printf("lami%d=%13.11le\n",t,lam); }cond2=fabs(lams/lamk);//计算A的条件数printf("condA=%13.11le\n",cond2);for(i=0;i<=4;i++)//赋值B=A{for(j=0;j<=500;j++){B[i][j]=A[i][j];}}fenjie();//对B进行LU分解，其第三行的元素即为矩阵U的对角线元素，即乘积为A的行列式的值for(i=0;i<=500;i++){det=det*B[2][i];}printf("det(A)=%13.11le\n",det);//计算输出A的行列式}/********幂法程序********/double mifa()//幂法的函数{int i,j,p,s;double u[501],y[501]={0};double ita,a,beta1=0,beta2=0,b,c,lam;for(i=0;i<=500;i++) u[i]=1;while(1)//开始迭代{a=0;for(i=0;i<=500;i++){a=a+u[i]*u[i];}ita=sqrt(a);//计算2数for(i=0;i<=500;i++)//计算向量y(k-1){y[i]=u[i]/ita;}for(i=0;i<=500;i++){u[i]=0;}for(i=2;i<=502;i++)//计算u(k){if(i<4) p=i;else p=4;if(i<500) s=i;else s=500;for(j=i-p;j<=s;j++){u[i-2]=u[i-2]+B[i-j][j]*y[j];}}for(i=0;i<=500;i++) beta2=beta2+y[i]*u[i];//计算beta(T) b=fabs(beta2-beta1);c=fabs(beta2);if((b/c)<=sigma)//如果达到精度则结束并返回值{lam=beta2;goto finish;}else{beta1=beta2;beta2=0;}}finish:return(lam);}/********反幂法程序********/double fanmifa()//反幂法函数{int i,k,t,p;double u[501],y[501]={0},d[501]={0};double ita,a,beta1=0,beta2=0,b,c,lam;for(i=0;i<=500;i++) u[i]=1;fenjie();//对系数矩阵进行LU分解/*for(i=0;i<=500;i++)printf("%12le ",B[0][i]);*/for(k=1;;k++)//开始迭代{a=0;for(i=0;i<=500;i++){a=a+u[i]*u[i];}ita=sqrt(a);//计算2数for(i=0;i<=500;i++)//计算向量y（k-1）{y[i]=u[i]/ita;}for(i=0;i<=500;i++){u[i]=0;}d[0]=y[0];//用LU分解解拟三角矩阵的方法求解u(k) for(i=1;i<=500;i++){t=1;if(i-1>t) t=i-1;else t=t;double temp1=0;for(;t<=i;t++){temp1+=B[i-t+3][t-1]*d[t-1];}d[i]=y[i]-temp1;}u[500]=d[500]/B[2][500];for(i=499;i>=0;i--){if(i+2>500) p=500;else p=i+2;double temp2=0;for(t=i+1;t<=p;t++){temp2+=B[i-t+2][t]*u[t];}u[i]=(d[i]-temp2)/B[2][i];}for(i=0;i<=500;i++) beta2=beta2+y[i]*u[i];b=fabs(beta2-beta1);c=fabs(beta2);if((b/c)<=sigma)//如果达到精度则结束并返回值{lam=1/beta2;goto finish;}else{beta1=beta2;beta2=0;}}finish:return(lam);}/********矩阵分解程序********/void fenjie()//LU分解函数{int k,t,j,i,p;for(i=3;i<=4;i++)//对于k=1时对应的矩阵元素先赋值{B[i][0]=B[i][0]/B[2][0];}for(k=2;k<=501;k++)//从k=2开始将U中的元素存到系数矩阵中 {if(k+2<501) p=k+2;else p=501;for(j=k;j<=p;j++){t=1;if(k-2>t) t=k-2;else t=t;if(j-2>t) t=j-2;else t=t;double temp3=0;for(;t<=k-1;t++){temp3+=B[k-t+2][t-1]*B[t-j+2][j-1];}B[k-j+2][j-1]=B[k-j+2][j-1]-temp3;}if(k<501){if(k+2<501) p=k+2;else p=501;for(i=k+1;i<=p;i++){t=1;if(i-2>t) t=i-2;else t=t;if(k-2>t) t=k-2;else t=t;double temp4=0;for(;t<=k-1;t++){temp4+=B[i-t+2][t-1]*B[t-k+2][k-1];}B[i-k+2][k-1]=B[i-k+2][k-1]-temp4;B[i-k+2][k-1]=B[i-k+2][k-1]/B[2][k-1];}}}}三、计算结果计算结果如下图所示四、编程中遇到的问题在幂法与反幂法的求解过程中初始值选取不同得到的结果也不同，比如取u0=[1,0,…0]和取u0=[1,1,….1]得到的结果不同，书上说可以任意取初始向量，但这里却不可以，原因未知。