北航数值分析编程大作业一

《数值分析b》大作业一。

4系 zy1104219 许柱山。

一．算法设计方案：

1 .首先应该创建原矩阵，但由于题目所给的矩阵是501*501的大型矩阵，计算机在计算时可能会运行不起来。而我们可以发现原矩阵a是一个带状矩阵，因此可以利用这一特点将其压缩存储。

2．由于λ1λ2…λ501，所以在以所有特征值建立的数轴上，λ1、λ501位于数轴的两端，两者之一必为按模最大。利用幂法，可以求出来按模最大的特征值，但该值可能为λ1也可能为λ501；

3．上步求出按模最大的特征值λm后，将原矩阵平移λm,再利用幂法求一次平移后矩阵的按模最大的特征值，即是数轴上另一端点值λm′。

4．比较λm与λm+λm′的大小，大的为矩阵a的最大特征值，小的为a矩阵的最小的特征值；

5．利用反幂法，求矩阵a的按模最小的特征值。但是反幂法中要用到线性方程组的求解，而原矩阵a又是带状矩阵，采用lu分解。所以在这之前要定义一个lu分解子程序，将a矩阵分解为单位下三角矩阵l和上三角矩阵u的乘积。

6．先利用循环求出k从1到39变化的uk的值。当循环次数为39时，在每次循环时都将压缩后的矩阵a的第三行减去相应的，然后调用lu分解的子程序，利用反幂法求出与uk最接近的特征值，该特征值等于利用反幂法求出的值与uk的和。

7．a的谱范数条件数等于按模最大的特征值的绝对值与按模最小的特征值的绝对值之比，按模最大的特征值与按模最小的特征值已分别在前面求出。

6．a的行列式的值就是矩阵a进行lu分解后u的对角线元素的乘积。先把带状线性矩阵a[501][501]转存为一个矩阵c[5][501].

二．源程序。

#include<>

全局变量、函数声明。

#define n 501

#define epsi 1.0e-12

#define r 2

#define s 2

double c[5][na非零元素的压缩存储矩阵*/

double fuzhi对a赋值函数*/

void ludet利用lu分解求解矩阵a的行列式*/

int int_max2(int a,int b求两个数字中最大值的*/

int int_min2(int a,int b求两个数字中最小值的*/

int int_max3(int a,int b,int c);/求三个数字中最大值的*/

double mifa幂法计算矩阵a按模最大的特征值*/

double fmifa反幂法求矩阵a按模最大的特征值*/

/*主程序*//

void main()

int i;

/*利用幂法计算矩阵a的最大特征值和最小特征值*/

double a1,a2;

fuzhi();

a1=mifa();

if (a1<0)

cout<<"矩阵a最小的特征值lambda1："=0)

cout<<"矩阵a最大的特征值lambda501："

a2=mifa()+a1;

if (a2<0)

cout<<"矩阵a最小的特征值lambda1："=0)

cout<<"矩阵a最大的特征值lambda501："

double a3;

fuzhi();

a3=fmifa();

cout<<"矩阵a按模最小的特征值lambdas:"

fuzhi();

double u[39]=;

cout<<"与数uk最接近的特征值："

/*计算矩阵a的条件数，取2范数*/

double cond_a;

if(a1>0)

cond_a=fabs(a1/a3);

else if(a1<0)

cond_a=fabs(a2/a3);

cout<<"矩阵a的条件数："

double det=1.0;

fuzhi();

ludet对矩阵a作lu分解

for(i=0;i det*=c[2][i];

cout<<"矩阵a的行列式值："

int int_max2(int a,int b)

return(a>b?a:b); 求两个数字中最大值的子程序*/

int int_min2(int a,int b)

return(a}

int int_max3(int a,int b,int c)

int t;

if(a>b)

t=a求三个数字中最大值的子程序*/

else t=b;

if(t return(t);

*幂法求最大特征值*/

double mifa()

double u[n],unitu[n];

double sum,lengthu,beta,temp;

int i,j;

beta=0.0;

for(i=0;i u[i]=0.5迭代初始向量*/

do while(fabs(beta-temp)>=epsi);

return(beta);

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

double fmifa()

double u[n],unitu[n],unitu_u[n];

double sum,lengthu,beta,temp;

int i,j;

beta=0.0;

for(i=0;i u[i]=0.5迭代初始向量*/

do ludet();

for(i=2;i<=n;i++)

u[n-1]=unitu_u[n-1]/c[2][n-1];

for(i=n-1;i>=1;i--)

sum=0;

for(i=0;isum=sum+unitu[i]*u[i

beta=sum;

}while(fabs(beta-temp)>=epsi);

beta=1/beta;

return(beta);

double fuzhi()

int i,j;

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

for(j=0;j

return(c[i][j]);

*lu三角分解*/

void ludet()

double sum;

int k,i,j;

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

for(j=k;j<=int_min2(k+s,n);j++)

sum=0;

for(i=int_max3(1,k-r,j-s);i<=k-1;i++)

北航数值分析编程大作业一

数值分析编程大作业工业工程

北航数值分析大作业第二次

北航数值分析第二次大作业

其他用户还读了