数值第一次大作业

《数值分析B》计算实习题第一题姓名：乔迈学号：SY1004319学院：能源与动力工程学院一、算法的设计方案：1、首先创建原矩阵，对矩阵A进行幂法运算，得出按模最大的特征值。

2、由于λ1‹λ2‹…‹λ501，所以在以所有特征值建立的数轴上，λ1、λ501位于数轴的两端，两者之一必为按模最大。

利用幂法，可以求出来按模最大的特征值，但该值可能为λ1也可能为λ501；判断该特征值的符号，若为正，则为λ501，否则，为λ1。

3、A平移，再用幂法，即可得另一个最大或最小的特征值。

4、对矩阵A进行反幂法运算，得出按模最小的特征值λs。

5、A平移相应量μk，进行反幂法运算，即可得相应的特征值λik，但是反幂法中要用到线性方程组的求解，而原矩阵A又是带状矩阵，采用LU分解。

6、A的谱范数条件数cond（A）为按模最大的特征值与按模最小的特征值之比按模最大的特征值与按模最小的特征值已分别在前面求出。

7、detA为LU分解后U的对角线元素的乘积。

二、源程序：#include<stdio.h>#include<math.h>#include<iostream.h>#define N 501/*定义全局变量*/double a[5][N];double u[N],y[N];/*定义函数*/void input_a() //输入矩阵A{int i;for(i=2;i<N;i++)a[0][i]=-0.064;for(i=1;i<N;i++)a[1][i]=0.16;for(i=1;i<=N;i++)a[2][i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);for(i=0;i<N-1;i++)a[3][i]=0.16;for(i=0;i<N-2;i++)a[4][i]=-0.064;}void input_u() //输入初始向量u{for(int i=0;i<N;i++)u[i]=1;}int max(int a,int b) //求最大值{if(a>b) return a;else return b;}int min(int a,int b) //求最小值{if(a<b) return a;else return b;}void matrix_pingyi(double fun_a[5][N],double b) //矩阵的主对角线元素平移{for(int i=0;i<N;i++)a[2][i]=a[2][i]-b;}void doolittle(double fun_a[N][N],double fun_b[N]) //Doolittle分解{int k,i,j,t;int m,l;/*分解*/for(k=0;k<N;k++){m=min(k+2,N-1);for(j=k;j<=m;j++){double sum=0;l=max(max(0,k-2),j-2);for(t=l;t<=k-1;t++)sum+=fun_a[k-t+2][t]*fun_a[t-j+2][j];fun_a[k-j+2][j]=fun_a[k-j+2][j]-sum;}if(k<N-1){m=min(k+2,N-1);for(i=k+1;i<=m;i++){double sum=0;l=max(max(0,i-2),k-2);for(t=l;t<=k-1;t++)sum+=fun_a[i-t+2][t]*fun_a[t-k+2][k];fun_a[i-k+2][k]=(fun_a[i-k+2][k]-sum)/fun_a[2][k];}}}/*求解*/double copy_b[N];for(i=0;i<N;i++)copy_b[i]=fun_b[i];for(i=1;i<N;i++){double sum=0;m=max(0,i-2);for(t=m;t<=i-1;t++)sum+=fun_a[i-t+2][t]*copy_b[t];copy_b[i]=copy_b[i]-sum;}u[N-1]=copy_b[N-1]/fun_a[2][N-1];for(i=N-2;i>=0;i--){double sum=0;l=min(i+2,N-1);for(t=i+1;t<=l;t++)sum+=fun_a[i-t+2][t]*u[t];u[i]=(copy_b[i]-sum)/fun_a[2][i];}}double vector_chengji(double fun_a[N],double fun_b[N]) //计算向量乘积{double sum=0;for(int i=0;i<N;i++)sum=sum+fun_a[i]*fun_b[i];return sum;}void vector_bianxing(double fun_a[N],double b) //向量变形{for(int i=0;i<N;i++)y[i]=fun_a[i]/b;}void matix_vector(double fun_a[5][N],double fun_b[N]) //计算矩阵与向量的乘积{int i,j;double m;for(i=0;i<N;i++)u[i]=0;for(i=0;i<N;i++){m=max(0,i-2);for(j=min(i+2,N-1);j>=m;j--)u[i]+=a[i+2-j][j]*y[j];}/*for(i=0;i<2;i++){for(j=i+2;j>=0;j--)u[i]=u[i]+a[i+2-j][j]*y[j];}for(i=2;i<N-2;i++){for(j=i+2;j>=i-2;j--)u[i]+=a[i+2-j][j]*y[j];}for(i=N-2;i<N;i++){for(j=N-1;j>=i-2;j--)u[i]=u[i]+a[i+2-j][j]*y[j];}*/}double mifa(double a[5][N]) //幂法{double v0=0,v1;double v_middle;double xiangdui;input_u();for(int k=0;;k++){v_middle=vector_chengji(u,u);v_middle=sqrt(v_middle);vector_bianxing(u,v_middle);matix_vector(a,y);v1=vector_chengji(y,u);xiangdui=fabs(v1-v0)/fabs(v1);v0=v1;if(xiangdui<=1.0e-12)break;}return v1;}double fanmifa(double a[5][N],double b) //反幂法{double v0=0,v1;double v_middle;double xiangdui;input_u();for(int k=0;;k++){v_middle=vector_chengji(u,u);v_middle=sqrt(v_middle);vector_bianxing(u,v_middle);input_a();matrix_pingyi(a,b);doolittle(a,y);v1=vector_chengji(y,u);v1=1/v1;xiangdui=fabs(v1-v0)/fabs(v1);v0=v1;if(xiangdui<=1.0e-12)break;}return v1;}double det(double fun_a[5][N]) //求行列式的值{int k,i,j,t;int m,l;for(k=0;k<N;k++){m=min(k+2,N-1);for(j=k;j<=m;j++){double sum=0;l=max(max(0,k-2),j-2);for(t=l;t<=k-1;t++)sum+=fun_a[k-t+2][t]*fun_a[t-j+2][j];fun_a[k-j+2][j]=fun_a[k-j+2][j]-sum;}if(k<N-1){m=min(k+2,N-1);for(i=k+1;i<=m;i++){double sum=0;l=max(max(0,i-2),k-2);for(t=l;t<=k-1;t++)sum+=fun_a[i-t+2][t]*fun_a[t-k+2][k];fun_a[i-k+2][k]=(fun_a[i-k+2][k]-sum)/fun_a[2][k];}}}double sum=1;for(j=0;j<N;j++)sum=sum*a[2][j];return sum;}/*定义主函数*/int main(int argc, char* argv[]){printf("\n");printf("第1题答案是：\n");printf("\n");double t_max,t_min,t_s;input_a();t_max=mifa(a);matrix_pingyi(a,t_max);t_min=mifa(a);if(t_max>=0){t_min=t_min+t_max;}else{double t;t=t_max;t_max=t_min+t_max;t_min=t;}printf("矩阵的最小特征值λ1=%.11e\n",t_min);printf("矩阵的最大特征值λ501=%.11e\n",t_max);t_s=fanmifa(a,0);printf("矩阵的按模最小特征值λs=%.11e\n",t_s);printf("\n");printf("第2题答案是：\n");printf("\n");double miu[39];for(int k=1;k<=39;k++){double l;l=t_min+k*(t_max-t_min)/40;miu[k-1]=fanmifa(a,l)+l;printf("与%.11e 最接近的特征值λ%2d=%.11e\n",l,k,miu[k-1]);}printf("\n");printf("第3题答案是：\n");printf("\n");double cond;if(fabs(t_min)>fabs(t_max))cond=fabs(t_min);elsecond=fabs(t_max);cond=cond/fabs(t_s);printf("矩阵A的条件数cond(A)=%.11e\n",cond);input_a();printf("矩阵的行列式值detA=%.11e\n",det(a));printf("\n");return 0;}三、计算结果：四．讨论迭代初始向量的选取对于计算结果的影响：1.影响迭代速度。