北航数值分析作业第一题题解

北航数值分析作业第一题：

一、算法设计方案

1.要求计算矩阵的最大最小特征值，通过幂法求得模最大的特征值，进行一定

判断即得所求结果；

2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝

对值满足反幂法的要求，故通过反幂法即可求得；

3.反幂法计算时需要方程求解中间过渡向量，需设计Doolite分解求解；

4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对角线相

乘即为矩阵的Det。

算法编译环境：vlsual c++6.0

需要编译函数：幂法，反幂法，Doolite分解及方程的求解

二、源程序如下：

#include

#include

#include

#include

int Max(int value1,int value2);

int Min(int value1,int value2);

void Transform(double A[5][501]);

double mifa(double A[5][501]);

void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]);

double Det(double A[5][501]);

/***定义2个判断大小的函数，便于以后调用***/

int Max(int value1,int value2)

{

return((value1>value2)?value1:value2);

}

int Min(int value1,int value2)

{

return ((value1

}

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

/***将矩阵值转存在一个数组里，节省空间***/

void Transform(double A[5][501],double b,double c)

{

int i=0,j=0;

A[i][j]=0,A[i][j+1]=0;

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

A[i][j]=c;

i++;

j=0;

A[i][j]=0;

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

A[i][j]=b;

i++;

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

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

i++;

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

A[i][j]=b;

A[i][j]=0;

i++;

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

A[i][j]=c;

A[i][j]=0,A[i][j+1]=0;

}

/***转存结束***/

//用于求解模最大的特征值，幂法

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

{

int s=2,r=2,m=0,i,j;

double b2,b1=0,sum,u[501],y[501];

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

{

u[i] = 1.0;

}

do

{

sum=0;

if(m!=0)b1=b2;

m++;

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

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

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

y[i]=u[i]/sqrt(sum);

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

{

u[i]=0;

for(j=Max(i-r,0);j<=Min(i+s,500);j++)

u[i]=u[i]+A[i-j+s][j]*y[j];

}

b2=0;

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

b2=b2+y[i]*u[i];

}

while(fabs(b2-b1)/fabs(b2)>=exp(-12));

return b2;

}

//带状DOOLITE分解，并且求解出方程组的解

void daizhuangdoolite(double A[5][501],double x[501],double b[501]) {

int i,j,k,t,s=2,r=2;

double B[5][501],c[501];

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

{

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

B[i][j]=A[i][j];

}

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

c[i]=b[i];

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

{

for(j=k;j<=Min(k+s,500);j++)

{

for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)

B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j];

}

for(i=k+1;i<=Min(k+r,500);i++)

{

for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)

B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k];

B[i-k+s][k]=B[i-k+s][k]/B[s][k];

}

}

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

for(t=Max(0,i-r);t<=i-1;t++)

c[i]=c[i]-B[i-t+s][t]*c[t];

x[500]=c[500]/B[s][500];

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

{

x[i]=c[i];

for(t=i+1;t<=Min(i+s,500);t++)