数值分析大作业
一、算法设计方案
1、矩阵初始化
矩阵A??aij?501?501的下半带宽r=2,上半带宽s=2，设置矩阵C?r?s?1??501?，在矩阵C中检索矩阵A中的带内元素aij的方法是：aij?ci?j?s?1,j。这样所需要的存储单元数大大减少，从而极大提高了运算效率。 2、利用幂法求出?1，?501 幂法迭代格式：
?非零向量u0?Rn?T??k?1?uk?1uk?1??yk?1?uk?1?k?1 ?u?Ayk?1?kT??k?yk?1uk?当?k??/?k?10?12时，迭代终止。
首先对于矩阵A利用幂法迭代求出一个?，然后求出矩阵B，其中
，对矩阵B进行幂法迭代，求出??，之后令????????，B?A??I(I为单位矩阵）
大小，大者为?501，小者为?1。 比较?与???的3、利用反幂法求出?s,?ik反幂法迭代格式：
?非零向量u0?Rn???k?1?uTk?1u?k1? ?yk?1?uk?1??k1
?Au?yk?1?kT??k?yk?1uk?当?k??/?k?10?12时，迭代终止，?s?1?k。
每迭代一次都要求解一次线性方程组Auk?yk?1，求解过程为： （1）作分解A?LU 对于k?1,2,...,n执行
1
k?1ck?j?s?1,j:?ck?j?s?1,j?k?t?s?1,tt?j?s?1,jt?max(1,k?r,j?s)?ck?1c?j?k,k?1,...,min(k?s,n)?ci?k?s?1,k:?(ci?k?s?1,k?i?t?s?1,tt?k?s?1,kt?max(1,i?r,k?s)?cc)/cs?1,k
?i?k?1,k?2,...,min(k?r,n);k?s?（2）求解Ly?b,Ux?y（数组b先是存放原方程组右端向量，后来存放中间向量y)
bi:?bi?i?t?s?1,ttt?max(1,i?r)i?1?cb(i?2,3,...,n)xn:?bn/cs?1.nmin(i?s,n)
i?t?s?1.ttxi:?(bi?t?i?1?cx)/cs?1,i(i?n?1,n?2,...,1) 使用反幂法，直接可以求得矩阵按模最小的特征值?s。 求与数?k??1?k?501??140(k?1,2,...,39)最接近的特征值?ik，对矩阵A??kI实行反幂法，即可求出对应的?k?1/?k,?ik??k??k。 4、求出A的条件数和行列式 根据cond(A)2??max，其中分子分母分别对应按模最大和最小的特征值。 ?sdet(A)的计算：由于A?LU,其中L为下三角矩阵，且对角线元素为1，故det(L)?1，所以有A?LU?U，又U为上三角矩阵，故det(U)为对其对角线
上各元素的乘积，最后可得det(A)?det(U)。
2
二、程序源代码
（1）定义所需要的函数：
#include
int min(int a,int b); // 求最小值
int max(int a,int b,int c); // 求最大值
double Fan_two(double x[N]);//计算二范数
void FenjieLU(double (*C)[N]);//解线性方程组的LU分解过程
void Solve(double (*C)[N], double *b,double *x);//解线性方程组的求解过程
double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//幂法
double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D);//反幂法 };
（2）程序的主函数，Main.cpp代码如下： void main() {
double C[R+S+1][N]; double u[N]; double y[N]; double miu[39];
double C1[R+S+1][N]; double bta = 1.0;
double Namda1,Namda501,NamdaS; double Namda[39]; double CondA2; double detA = 1.0; double D = 1.0e-12; int i, j, k; FILE * fp;
fp = fopen(\ //对数组进行初始化// int i, j;
for (i = 0; i < N; i++) { u[i] = 1;
3
}
for (i = 0;i< R + S + 1;i++) { for (j = 0;j< N;j++) { if (i==0||i==4) { C[i][j]=-0.064; } else if (i==1||i==3) { C[i][j]=0.16; } else if (i==2) { C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1)) -0.64*exp(0.1/(j+1)); } } }
//幂法求Namda1//
Namda1 = PowerMethod(C, u, y, bta, D);
printf(\printf(\
printf(\//幂法求Namda501// bta = 1.0;
for (i = 0; i < R + S + 1; i++) { for (j = 0; j < N; j++) { if (i == 2) C1[i][j] = C[i][j] -Namda1; else C1[i][j] = C[i][j]; } }
Namda501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1;
printf(\printf(\
printf(\//反幂法求NamdaS// bta = 1.0;
4
NamdaS = InversePowerMethod(C, u, y, bta, D);
printf(\ printf(\
printf(\ //反幂法求Namda[k]//
printf(\ for (k = 0; k < 39; k++) { miu[k] = Namda1 + (k + 1) * (Namda501 - Namda1) / 40.0; bta = 1.0; for (i = 0; i < R + S + 1; i++) { for (j = 0; j < N; j++) { if (i == 2) C1[i][j] = C[i][j] - miu[k]; else C1[i][j] = C[i][j]; } } Namda[k] = InversePowerMethod(C1, u, y, bta, D) + miu[k]; fprintf(fp,\与.11e最接近的特征值为:.11e\\n\ }
printf(\求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.txt中\
printf(\ //求A的谱范数//
printf(\ printf(\的谱范数为:.11e\
printf(\ //求A的条件数//
CondA2 = fabs( Namda1 / NamdaS);
printf(\ printf(\的谱范数的条件数Cond(A)2为:.11e\
printf(\ //求det(A)2的值// for (j = 0; j < N; j++) detA *= C[2][j];
printf(\ printf(\行列式A的值为:.11e\
printf(\ fclose(fp); _getch(); return;
5