北航数值分析大作业一
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
《数值分析B》大作业一
SY1103120 朱舜杰
一.算法设计方案:
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所有对角线上元素的乘积。
二.源程序
#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(t return(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的随机数*/