数值作业1.docx

《数值作业1.docx》由会员分享，可在线阅读，更多相关《数值作业1.docx（16页珍藏版）》请在冰豆网上搜索。

数值作业1

1.第一题　2

算法设计说明：

一.在计算最小特征值和最大特征值是的算法设计如下：

1.对A用幂法计算出矩阵的按模为最大的特征值：

λm；

2.对矩阵B＝A-λｍI用幂法计算出矩阵的按模为最大的特征值：

λn；

3.λp＝λm＋λn，则有：

λ1＝min（λm，λp），λ501＝max（λm，λp）.

二.幂法迭代公式：

　　当｜βk－βk-1｜／｜βk｜≤1ｅ－12时，迭代中止，以当前的βk作为按模最大特征值的近似值．

三．反幂法迭代公式：

每迭代一次都要求解一次线性方程组Auk=y（k-1）．ｋ足够大时，当｜βk－βk-1｜／｜βk｜≤1ｅ－12时，迭代中止，此时可取λｎ≈1／βk．

四．在对给定矩阵存储的时，由于A是上半带宽度为ｓ＝2，下半带宽度ｒ＝2的带状矩阵，A的零元素较多，为了节省存储量，A的带外元素不给存储，仅存A的带内元素，将A以及进行平移后矩阵存在一个二维数组C（5,501），数组C的第ｊ列存放A的第ｊ列带内元素，并使A的主对角线元素存放C的第ｓ＋1＝3行中．在数组C中检索矩阵A的带内元素ａｉｊ的方法是：

　　　A的带内元素ａｉｊ＝C的元素ci－j＋s＋1，j

五．在求解一次线性方程组Auk=y（k-1）时，需要对A进行Doolittle分解来求解新的迭代向量uk，算法可分解两步：

（1）作分解A＝LU

对于k＝1，2，…ｎ（501）执行

　　　　ck-j+s+1,j:

=ck-j+s+1,j-

ck-t+s+1,jct-j+s+1,j

[j=k,k+1,…min（k+s,n）]

Ci-k+s+1,k:

=（ci-k+s+1,k-

ci-t+s+1,tct-k+s+1,k）/cs+1,k

[i=k+1,k+2,…min（k+r,n）;k

（2）求解Lx＝y，Uuk＝x（数组ｙ先是存放原方程组右端的向量，后来存放中间变量，最后求解新的迭代向量uk）

　　yi:

=yi-

ci-t+s+1,tyt

[i=2,3,…,n]

un:

=yn/cs+1,n

un:

=（yi-

ci-t+s+1,tut）/cs+1,i

[i=n-1,n-2,…1]

六．算法说明：

　　1.选择不同的初始向量会得到不同的特征值，这是因为收敛速度的不同引起的．通过选择几组不同的初始向量比较后，发现当初始向量按照程序中的取值时比较合理，得到的特征值较为精确．

　　　2.利用矩阵U所有主元素乘积求detA，cond（A）

可有公式cond（A）

＝

求得。

，

分别利用幂法和反幂法求的。

#include

#include

#include

intMax（inta,intb,intc）;

intMin（inta,intb）;

doublea[501],C[5][501];/*声明全局变量*/

/******************************************/

/*函数：

power*/

/*功能：

幂法按模最大求特征值*/

/*说明：

p是矩阵运算是需要平移的量*/

/******************************************/

doublepower（doublep）

{

intm,k,t,n;

doubledPreVal,dVal;

doubley[501],u[501],s=0,b=0.16,c=-0.064,sum;

dPreVal=0.0;

dVal=0.0;

t=0;

/*********初始化u[n]*******/

for（n=0;n

u[n]=1;

while（t=1e-12）/*迭代循环*/

{sum=0;

dPreVal=dVal;

dVal=0.0;

for（k=0;k

{

sum+=u[k]*u[k];

}

/**********计算迭代向量Y*********/

for（k=0;k

{

y[k]=u[k]/sqrt（sum）;

}

/*******矩阵与迭代向量相乘***********/

/*******求出新的迭代向量***********/

/*说明：

p是矩阵运算是需要平移的量*/

for（m=0;m

{

if（m==0）

{

s=（a[m]-p）*y[0]+b*y[1]+c*y[2];

}

elseif（m==1）

{

s=b*y[0]+（a[m]-p）*y[1]+b*y[2]+c*y[3];

}

elseif（m==499）

{

s=c*y[497]+b*y[498]+（a[m]-p）*y[499]+b*y[500];

}

elseif（m==500）

{

s=c*y[498]+b*y[499]+（a[m]-p）*y[500];

}

elseif（m!

=0&&m!

=1&&m!

=499&&m!

=500）

{

s=c*y[m-2]+b*y[m-1]+（a[m]-p）*y[m]+b*y[m+1]+c*y[m+2];

}

u[m]=s;

dVal=dVal+u[m]*y[m];

}

t++;

}

returndVal;

}

/******************************************/

/*函数：

Max*/

/*功能：

求三个数的最大值*/

/******************************************/

intMax（inta,intb,intc）

{

intmax;

max=a;

if（b>max）max=b;

if（c>max）max=c;

returnmax;

}

/******************************************/

/*函数：

Min*/

/*功能：

求两个数的最小值*/

/******************************************/

intMin（inta,intb）

{

if（a

elsereturnb;

}

/******************************************/

/*函数：

revpower*/

/*功能：

反幂法按模最小求特征值*/

/******************************************/

doublerevpower（）

{

Intk,t,n,m,x;

doubledPreVal,dVal;

doubley[501],u[501],v[501],s,sum;

dPreVal=0.0;

dVal=0.0;

t=0;

/********初始化迭代向量u[n]**********/

for（n=0;n

u[n]=1;

/*************迭代循环（根据精度要求选择适当的迭代次数）********/

while（t=1e-12）

{sum=0;

dPreVal=dVal;

dVal=0.0;

for（k=0;k

{

sum+=u[k]*u[k];

}

/***********计算迭代向量Y**************/

for（k=0;k

{

y[k]=u[k]/sqrt（sum）;

v[k]=y[k];

}

/**********根据Doolittle分解来求解方程组得到新的迭代向量*********/

for（k=2;k

{

s=0;

m=Max（1,1,k-2）;

for（x=m;x

{

s=s+C[k-x+2][x-1]*y[x-1];

}

y[k-1]=y[k-1]-s;

}

u[500]=y[500]/C[2][500];

for（k=500;k>0;k--）

{

s=0;

m=Min（k+2,501）;

for（x=k+1;x

{

s=s+C[k-x+2][x-1]*u[x-1];

}

u[k-1]=（y[k-1]-s）/C[2][k-1];

}

for（k=0;k

{

dVal=dVal+v[k]*u[k];

}

t++;

}

return1.0/dVal;

}

/******************************************/

/*函数：

SAVE*/

/*功能：

存储矩阵A带内元素，不存非零元素*/

/*说明：

p是矩阵运算是需要平移的量*/

/******************************************/

voidSave（doublep）

{

inti,j;

doubleb=0.16,c=-0.064;

for（i=0;i

{

for（j=0;j

{

if（i==0||i==4）

{

if（（i==0&&j=498））C[i][j]=0;

elseC[i][j]=c;

}

elseif（i==1||i==3）

{

if（（i==1&&j==0）||（i==3&&j==500））C[i][j]=0;

elseC[i][j]=b;

}

elseif（i==2）

{

C[i][j]=a[j]-p;

}

}

}

}

/******************************************/

/*函数：

Doolittle*/

/*功能：

Doolittle分解矩阵A带内元素*/

/******************************************/

voidDoolittle（）

{

ints,t,x,k,j,i;

doublesum;

for（k=1;k

{

s=Min（k+2,501）;

for（j=k;j

{

t=Max（1,k-2,j-2）;

sum=0.0;