北航数值分析实习题目第一题

北航数值分析实习题目第一题

《数值分析B》大作业一

ZY1515105 樊雪松

一． 算法设计方案：

1.矩阵A的存储与检索

将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] 。在数组MatrixC[5][501]中检索A的带内元素aij的方法是：A的带内元素aij=C中的元素ci-j+2,j 。

2.求解λ1，λ501，λs

1、首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs；如果λ max>0,则λ501=λmax；如果λmax<0,则λ1=λmax。

2、使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出 对应的按摸最大的特征值λ’max， 如果λ max>0,则λ1=λ’max+p；如果λmax<0,则λ501=λ’max+p。

3、求解A的与数μk=λ1+k（λ501-λ1）/40 的最接近的特征值λik （k=1,2，…，39）。

使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λ ik。

4、求解A的（谱范数）条件数cond（A）2和行列式detA。

cond（A）其中λ1和λn分别是矩阵 A的模最大和最 小2=|λ1/λn|，

特征值。求解矩阵A的行列式，可先对矩阵A进行LU分解后，detA

北航数值分析实习题目第一题

等于U所有对角线上元素的乘积。

二． 源程序

#include<stdio.h>

#include<math.h>

#include<conio.h>

//定义A中元素

double C[5][501];

double a[501];

double b;

double c;

//声明所有函数

void YaSuoJZ(double C[5][501],double a[501],double b,double c) ;//压缩矩阵函数

double mifa(double C[5][501]); //幂法函数

void daizhuangLU(double A[5][501]); //带状矩阵的LU分解

double fanmifa(double C[5][501]);//反幂法函数

//最值函数

int max2(int x,int y);

int max3(int x,int y,int z);

int min(int x,int y);

//最值函数

int max2(int x,int y) //求2个数的最大值

{

int z;

z=x>y?x:y;

return(z);

}

int max3(int x,int y,int z) //求3个数的最大值

{

int w;

w = z > max2(x,y)? z:max2(x,y);

return(w);

}

int min(int x,int y) //求2个数的最小值

{

int z;

z=x>y?y:x;

北航数值分析实习题目第一题

return(z);

}

//将矩阵A压缩存储在矩阵C中

void YaSuoJZ(double C[5][501],double a[501],double b,double c)

{

int i;

for(i=0;i<=500;i++)

{

if(i>=2) C[0][i]=c;

else C[0][i]=0;

if(i>=1) C[1][i]=b;

else C[1][i]=0;

if(i<=499) C[3][i]=b;

else C[3][i]=0;

if(i<=498) C[4][i]=c;

else C[4][i]=0;

C[2][i]=a[i];

}

}

//幂法函数：用幂法求矩阵模最大的特征值

double mifa(double C[5][501])

{

double u[501];

double y[501]={0},η=0;

double β,βk=0;

double ε=1; // ε为精度

double sumu=0,sumAY=0;

int i,j,k=1; //k为循环次数

for (i=0;i<=500;i++) //取任一非零向量u0

u[i] = 1.0;

while(ε>=1e-12)

{

for(i=0;i<=500;i++) //求u(k-1)的2范数η

sumu=sumu+u[i]*u[i];

η=sqrt(sumu);

sumu=0;

for(i=0;i<=500;i++) //求y(k-1)

y[i]=u[i]/η;

for(i=0;i<=500;i++) //求u(k)的各分量u[i]

{

for(j=max2(0,i-2);j<=min(i+2,500);j++)

北航数值分析实习题目第一题

sumAY=sumAY+C[i-j+2][j]*y[j];

u[i]=sumAY;

sumAY=0;

}

//求幂法中的βk

β=βk; //将β(k-1)放在β中

βk=0;

for(i=0;i<=500;i++) //求βk

βk=βk+y[i]*u[i];

if(k>=2)

ε=fabs(βk-β)/fabs(βk);

k++;

}

return(βk);

}

//带状矩阵的LU分解

void daizhuangLU(double A[5][501])

{

int i,j,k,m,t;

double sumukj=0,sumlik=0;

for(k=0;k<=500;k++)

{

for(j=k;j<=min(k+2,500);j++) //求ukj并存在A[k-j+2][j]中 {

for(t=max3(0,k-2,j-2);t<=k-1;t++)

sumukj=sumukj+A[k-t+2][t]*A[t-j+2][j];

A[k-j+2][j]=A[k-j+2][j]-sumukj;

sumukj=0;

}

if(k<500)

for(i=k+1;i<=min(k+2,500);i++) //求lik并存在A[i-k+2][k]中 {

for(m=max3(0,i-2,k-2);m<=k-1;m++)

sumlik=sumlik+A[i-m+2][m]*A[m-k+2][k];

A[i-k+2][k]=(A[i-k+2][k]-sumlik)/A[2][k];

sumlik=0;

}

}

}

//反幂法函数：用反幂法求矩阵的模最小的特征值

double fanmifa(double M[5][501])

{

北航数值分析实习题目第一题

double u[501];

double y[501]={0},x[501],η=0;

double fβ,fβk=0;

double ε=1;

double fsumu=0,sumLX=0,sumUu=0;

int i,t,m,k=1;

for(i=0;i<=500;i++) //任取一非零向量u0

u[i]=1;

daizhuangLU(M); //对A进行LU分解A=LU，Au(k)=y(k-1)等价于Uu(k)=x和Lx=y(k-1) while(ε>=1e-12)

{

for(i=0;i<=500;i++) //求u(k-1)的2范数η

fsumu=fsumu+u[i]*u[i];

η=sqrt(fsumu);

fsumu=0;

for(i=0;i<=500;i++) //求y(k-1)

y[i]=u[i]/η;

for(i=0;i<=500;i++) //求中间向量x

x[i]=y[i];

for(i=1;i<=500;i++)

{

for(t=max2(0,i-2);t<=i-1;t++)

sumLX=sumLX+M[i-t+2][t]*x[t];

x[i]=x[i]-sumLX;

sumLX=0;

}

u[500]=x[500]/C[2][500]; //求u(k)的各分量u[i]

for(i=499;i>=0;i--)

{

for(m=i+1;m<=min(i+2,500);m++)

sumUu=sumUu+M[i-m+2][m]*u[m];

u[i]=(x[i]-sumUu)/M[2][i];

sumUu=0;

}

//求反幂法中的βk

fβ=fβk; //将fβ(k-1)放在fβ中

fβk=0;

for(i=0;i<=500;i++) //求fβk

fβk=fβk+y[i]*u[i];

if(k>=2)

ε=fabs(1/fβk-1/fβ)/fabs(1/fβk);

k++;

}

return(1/fβk);

北航数值分析实习题目第一题

}

//主函数

void main()

{

int i,j,k;

double λ1,λ501,λm,λm1,λm2,λs,λ,p;

double cond,detA=1;

for(i=1;i<=501;i++)

a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);

b=0.16;

c=-0.064;

YaSuoJZ(C,a,b, c); //将矩阵A中元素压缩存储在C中

λm1=mifa(C); //对A用幂法求出模最大的特征值λm1

λs=fanmifa(C); //对A用反幂法求出模最小的特征值λs

YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中

for(j=0;j<=500;j++) //对A进行平移，平移量为λm1，平移后矩阵元素压缩存储在C中 C[2][j]=C[2][j]-λ?m1;

λm=mifa(C);

λm2=λm1+λm; //λm1与λm2是矩阵的最大最小特征值

if(λm1>λm2) //判断A最大最小特征值

{

λ501=λm1;λ1=λm2;

}

else

{

λ501=λm2;λ1=λm1;

}

printf("数值分析计算实习第一题\n\n ZY1515105 樊雪松\n\n （1）A的最大最小以及模最小的特征值\n");

printf("A的最小特征值λ1=%.13e\n",λ1);

printf("A的最大特征值λ501=%.13e\n",λ501);

printf("A的模最小特征值λs=%.13e\n",λs);

printf("\n（2）与数μk最接近的特征值\n");

printf("\t要求接近的值\t\t\t实际求得的特征值\n");

YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中

for(k=1;k<=39;k++)

{

p=λ1+k*(λ501-λ1)/40;

for(j=0;j<=501;j++)

C[2][j]=C[2][j]-p;

λ=fanmifa(C)+p;

printf("μ%d=%.13e λ%d=%.13e\n",k,p,k,λ);

YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中

北航数值分析实习题目第一题

}

printf("\n（3）计算A的条件数cond(A)和行列式detA\n");

cond=λm1/λs;

daizhuangLU(C);

for(j=0;j<=500;j++)

detA=detA*C[2][j];

printf("A的条件数cond(A)=%.13e\n",cond);

printf("A的行列式detA=%.13e\n",detA);

getch();

}

三、运行结果

数值分析计算实习第一题

ZY1515105 樊雪松

（1）A的最大最小以及模最小的特征值

A的最小特征值λ1=-1.0700113615018e+001

A的最大特征值λ501=9.7246340987773e+000

A的模最小特征值λs=-5.5579107942295e-003

（2）与数μk最接近的特征值

要求接近的值 实际求得的特征值

μ1=-1.0189494922173e+001 λ1=-1.0182934033146e+001

μ2=-9.6788762293280e+000 λ2=-9.5857074250676e+000

μ3=-9.1682575364831e+000 λ3=-9.1726724239280e+000

μ4=-8.6576388436383e+000 λ4=-8.6522840078976e+000

μ5=-8.1470201507934e+000 λ5=-8.0934838086753e+000

μ6=-7.6364014579485e+000 λ6=-7.6594054076924e+000

μ7=-7.1257827651036e+000 λ7=-7.1196846486912e+000

μ8=-6.6151640722588e+000 λ8=-6.6117643393973e+000

μ9=-6.1045453794139e+000 λ9=-6.0661032265951e+000

μ10=-5.5939266865690e+000 λ10=-5.5851010526284e+000

μ11=-5.0833079937241e+000 λ11=-5.1140835298122e+000

μ12=-4.5726893008792e+000 λ12=-4.5788721768651e+000

μ13=-4.0620706080344e+000 λ13=-4.0964709262599e+000

μ14=-3.5514519151895e+000 λ14=-3.5542112157508e+000

μ15=-3.0408332223446e+000 λ15=-3.0410900181333e+000

μ16=-2.5302145294997e+000 λ16=-2.5339703111304e+000

μ17=-2.0195958366549e+000 λ17=-2.0032307695635e+000

北航数值分析实习题目第一题

μ18=-1.5089771438100e+000 λ18=-1.5035576112274e+000

μ19=-9.9835845096511e-001 λ19=-9.9355860600754e-001

μ20=-4.8773975812023e-001 λ20=-4.8704267388496e-001

μ21=2.2878934724645e-002 λ21=2.2317362495748e-002

μ22=5.3349762756952e-001 λ22=5.3241747420686e-001

μ23=1.0441163204144e+000 λ23=1.0528989626935e+000

μ24=1.5547350132593e+000 λ24=1.5894458818809e+000

μ25=2.0653537061042e+000 λ25=2.0603304602743e+000

μ26=2.5759723989490e+000 λ26=2.5580755970728e+000

μ27=3.0865910917939e+000 λ27=3.0802405093071e+000

μ28=3.5972097846388e+000 λ28=3.6136208676923e+000

μ29=4.1078284774837e+000 λ29=4.0913785104506e+000

μ30=4.6184471703285e+000 λ30=4.6030353782791e+000

μ31=5.1290658631734e+000 λ31=5.1329242838984e+000

μ32=5.6396845560183e+000 λ32=5.5949063480833e+000

μ33=6.1503032488632e+000 λ33=6.0809338570269e+000

μ34=6.6609219417080e+000 λ34=6.6803540921116e+000

μ35=7.1715406345529e+000 λ35=7.2938774481266e+000

μ36=7.6821593273978e+000 λ36=7.7171117142356e+000

μ37=8.1927780202427e+000 λ37=8.2252200140502e+000

μ38=8.7033967130876e+000 λ38=8.6486660651935e+000

μ39=9.2140154059324e+000 λ39=9.2542003445750e+000

（3）计算A的条件数cond(A)和行列式detA

A的条件数cond(A)=1.9252042739022e+003

A的行列式detA=2.7727861417521e+118

四、结果分析

设A的n个线性无关的特征向量为x，x2,…,1xn，其相对应的特征值满足的关系为 1 2 3 n。

1在幂法计算过程中，如果初始向量选择不当,将导致迭代中在x方向的分量等于零，但

是由于舍入误差的影响,必然会在迭代的某一步生成的迭代向量，按照基向量展开时,在x1方向的分量不等于零，把这一向量看作初始向量，重新开始迭代，用幂法继续求向量序列,仍然会得出预期的结果。但是迭代次数会变大，收敛速度较慢，所以初始向量的选取会影响到迭代的次数。

本文来源：https://www.bwwdw.com/article/tzl4.html