《数值分析A》计算实习题目第一题
一．算法设计方案：
1.矩阵A的存储与检索
将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] .
由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：
A的带内元素a ij=C中的元素c i-j+2,j
2.求解λ1，λ501，λs
①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs；
如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。
②使用带原点平移的幂法（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和行列式d etA。
①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有
对角线上元素的乘积。
二．源程序(VS2010环境下，C++语言)
#include
#include
#include
#include
#include
#include
#include
#define E 1.0e-12 /*定义全局变量相对误差限*/
int max2(int a,int b) /*求两个整型数最大值的子程序*/
{
if(a>b)
return a;
else
return b;
}
int min2(int a,int b) /*求两个整型数最小值的子程序*/
{
if(a>b)
return b;
else
return a;
}
int max3(int a,int b,int c) /*求三整型数最大值的子程序*/
{ int t;
if(a>b)
t=a;
else t=b;
if(treturn(t);
}
void assignment(double array[5][501]) /*将矩阵A转存为数组C[5][501]*/ {
int i,j,k;
//所有元素归零
for(i=0;i<=4;)
{
for(j=0;j<=500;)
{
array[i][j]=0;
j++;
}
i++;
}
//第0,4行赋值
for(j=2;j<=500;)
{
k=500-j;
array[0][j]=-0.064;
array[4][k]=-0.064;
j++;
}
//第1,3行赋值
for(j=1;j<=500;)
{
k=500-j;
array[1][j]=0.16;
array[3][k]=0.16;
j++;
}
//第2行赋值
for(j=0;j<=500;)
{ k=j;
j++;
array[2][k]=(1.64-0.024*j)*sin((double)(0.2*j))-0.64*exp((double)(0.1/j));
}
}
double mifa(double u[501],double array[5][501],double p) /*带原点平移的幂法*/ {
int i,j; /* u[501]为初始迭代向量*/
double a,b,c=0; /* array[5][501]为矩阵A的转存矩阵*/
double y[501]; /*p为平移量*/
for(;;)
{
a=0;
b=0;
/*选用第一种迭代格式*/ //求ηk-1
for(i=0;i<=500;i++)
{
a=a+u[i]*u[i];
}
a=sqrt(a);
//求y k-1
for(i=0;i<=500;i++)
{
y[i]=u[i]/a;
}
//求u k
for(i=0;i<=500;i++)
{
u[i]=0;
for(j=max2(i-2,0);j<=min2(i+2,500);j++)
{
u[i]+=array[i-j+2][j]*y[j];
}
u[i]=u[i]-p*y[i]; /*引入平移量*/
}
//求βk
for(i=0;i<=500;i++)
{
b+=y[i]*u[i];
}
if(fabs((b-c)/b)<=E) /*达到精度水平，迭代终止*/ break;
c=b;
}
return (b+p); /*直接返回A的特征值*/
}
void chuzhi(double a[]) /*用随机数为初始迭代向量赋值*/ {
int i;
srand((int)time(0));
for(i=0;i<=500;i++)
{
a[i]=(10.0*rand()/RAND_MAX); /*生成0~10的随机数*/ }
