数值分析计算实习作业一

数值分析计算实习题一

学号：

：

院系：

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的矩阵C

A 的带元素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

#include

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

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;