数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。

这样所需要的存储单元数大大减少，从而极大提高了运算效率。

2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。

3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。

求与数)39,...,2,1(4015011=-+=k k k λλλμ最接近的特征值ik λ，对矩阵IA k μ-实行反幂法，即可求出对应的k k ik k k μλλβλ+==,/1。

4、求出A 的条件数和行列式根据max 2()scond A λλ=，其中分子分母分别对应按模最大和最小的特征值。

det()A 的计算：由于A LU =,其中L 为下三角矩阵，且对角线元素为1，故det()1L =，所以有A LU U ==，又U 为上三角矩阵，故det()U 为对其对角线上各元素的乘积，最后可得det()det()A U =。

二、程序源代码（1）定义所需要的函数：#include <stdio.h>#include <conio.h>#include <math.h>#define N 501#define R 2#define S 2int 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("Namda.txt","w");//对数组进行初始化//int i, j;for (i = 0; i < N; i++){u[i] = 1;}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("\n================================================\n"); printf("Namda1 = %12.11e", Namda1);printf("\n================================================\n"); //幂法求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;elseC1[i][j] = C[i][j];}}Namda501 = algorism.PowerMethod(C1, u, y, bta, D) +Namda1;printf("\n================================================\n"); printf("Namda501 = %12.11e", Namda501);printf("\n================================================\n"); //反幂法求NamdaS//bta = 1.0;NamdaS = InversePowerMethod(C, u, y, bta, D);printf("\n================================================\n");printf("NamdaS = %12.11e", NamdaS);printf("\n================================================\n");//反幂法求Namda[k]//printf("\n================================================\n");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];elseC1[i][j] = C[i][j];}}Namda[k] = InversePowerMethod(C1, u, y, bta, D) + miu[k];fprintf(fp,"与%12.11e最接近的特征值为:%12.11e\n",miu[k],Namda[k]);}printf("求与miu[k]最接近的Namda[k]的计算结果已经输出到文件Namda.txt 中");printf("\n================================================\n");//求A的谱范数//printf("\n================================================\n");printf("A的谱范数为:%12.11e", sqrt(Namda501));printf("\n================================================\n");//求A的条件数//CondA2 = fabs( Namda1 / NamdaS);printf("\n================================================\n");printf("A的谱范数的条件数Cond(A)2为:%12.11e",CondA2);printf("\n================================================\n");//求det(A)2的值//for (j = 0; j < N; j++)detA *= C[2][j];printf("\n================================================\n");printf("行列式A的值为:%12.11e",detA);printf("\n================================================\n");fclose(fp);_getch();return;}（3）成员函数的实现int min(int a,int b){return a < b ? a : b;}int max(int a,int b,int c){int temp;temp = a > b ? a : b;return temp > c ? temp : c;}double Fan_two(double x[N]){double sum = 0.0;int i;for (i = 0; i < N; i++){sum += pow(x[i],2);}return sqrt(sum);}void FenjieLU(double (*C)[N]){double sum = 0;int i, j, k,t;for (k = 0; k < N; k++){j = k;i = k + 1;while (1){if (j == min(k + S + 1, N))break;for (t = max(0, k - R, j - S); t <= k - 1; t++){sum += C[k-t+S][t] * C[t-j+S][j];}C[k-j+S][j] = C[k-j+S][j] - sum;sum = 0.0;j++;if (k == N-1)break;if (i == min(k + R + 1, N))break;for (t = max(0, i - R,k - S); t <= k - 1; t++){sum += C[i-t+S][t] * C[t-k+S][k];}C[i-k+S][k] = (C[i-k+S][k] - sum) / C[S][k];sum = 0;i++;}}}void Solve(double (*C)[N], double *b,double *x){double sum = 0;int i, t;sum = 0;for (i = 1; i < N; i++){for (t = max(0, i - R); t <= i - 1; t++){sum += C[i-t+S][t] * b[t];}b[i] = b[i] - sum;sum = 0;}x[N-1] = b[N-1] / C[S][N-1];for (i = N - 2; i >= 0; i--){for (t = i+1; t <= min(i + S, N - 1); t++){sum += C[i-t+S][t] * x[t];}x[i] = (b[i] - sum) / C[S][i];sum = 0;}}double PowerMethod(double C[][N],double u[N],double y[N],double bta,double D) {double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;while (fabs(bta - temp) / fabs(bta) > D){temp = bta;ita = Fan_two(u);for (i = 0; i < N; i++){y[i] = u[i] / ita;}for (i = 0; i < N; i++){for (j = max(0,i - R); j < min(i + S + 1,N); j++){sum += C[i - j + S][j] * y[j];}u[i] = sum;sum = 0;}for (i = 0; i < N; i++){sum += y[i] * u[i];}bta = sum;sum = 0;k++;}return bta;}double InversePowerMethod(double C[][N],double u[N],double y[N],double bta,double D){double TC[R+S+1][N];double ty[N];double ita;double sum = 0;double temp = 0.0;int i,j,k = 0;FenjieLU(C);while (abs(1/bta - 1/temp) / abs(1/bta) > D){temp = bta;ita = Fan_two(u);for (i = 0; i < N; i++){y[i] = u[i] / ita;}//用到临时存储数组TC[][]和ty[][]是因为函数Solve执行过程中会改变A[][]和y[][]for (i = 0; i < R + S + 1; i++){for (j = 0; j < N; j++)TC[i][j] = C[i][j];}for (i = 0; i < N; i++)ty[i] = y[i];Solve(C, y, u);for (i = 0; i < R+S+1; i++){for (j = 0; j < N; j++)C[i][j] = TC[i][j];}for (i = 0; i < N; i++)y[i] = ty[i];for (i = 0; i < N; i++){sum += y[i] * u[i];}bta = sum;sum = 0;k++;}bta = 1.0 / bta;return bta;}三、程序运行结果下图为主程序运行结果的结果输出在Namda.txt文件中，结果如下：其中ik四、分析迭代初始向量对计算结果的影响选择不同的初始向量[]u N可能会得到不同的特征值。