北航数值分析作业第一题

数值分析作业第一题

一、 算法设计方案

利用带状Dollittle分解，将A[501][501]转存到数组C[5][501]，以节省存储空间

1、计算λ1和λ501

首先使用幂法求出矩阵的按模最大的特征值λ0：如果λ0>0，则其必为按模最大值，因此λ501=λ0，然后采用原点平移法，平移量为λ501，使用幂法迭代求出矩阵A-λ501I的按模最大的特征值，其特征值按从小到大排列应为λ1-λ501、λ2-λ501、……、0。因此A-λ501I的按模最大的特征值应为λ1-λ501。又因为λ501的值已求得，由此可直接求出λ1。

2、计算λS

λS为矩阵A按模最小的特征值，可以通过反幂法直接求出。

3、计算λik

λik是对矩阵A进行λik平移后，再用反幂法求出按模最小的特征值λλik=λik+λmin。

4、计算矩阵A的条件数计算cond（A）2和行列式det(A)

矩阵A的条件数为cond（A）2 min，1，其中λn1和λn分别是矩阵A的模最大

和最小特征值，直接利用上面求得的结果直接计算。

矩阵A的行列式可先对矩阵A进行LU分解后，det(A)等于U所有对角线上元素的乘积。

二、源程序：

#include<math.h>

#include<stdio.h>

#include<stdlib.h>

#include<iostream.h>

#define s 2

#define r 2

int Max(int v1,int v2);

int Min(int v1,int v2);

int maxt(int v1,int v2,int v3);

void storage(double C[5][501],double b,double c);

double mifa(double C[5][501]);

void LU(double C[5][501]);

double fmifa(double C[5][501]);

int Max(int v1,int v2) //求两个数的最大值

{ return((v1>v2)?v1:v2);

}

int Min(int v1,int v2) //求两个数最小值

{

}

int maxt(int v1,int v2,int v3) //求三个数最大值

{ int t;

if(v1>v2) t=v1;

}

/***将矩阵值转存在一个数组里，以节省存储空间***/

void storage(double C[5][501],double b,double c)

{

int i=0,j=0; C[i][j]=0,C[i][j+1]=0; for(j=2;j<=500;j++) i++; j=0; C[i][j]=0; for(j=1;j<=500;j++) i++; for(j=0;j<=500;j++) i++; for(j=0;j<=499;j++) i++; for(j=0;j<=498;j++) C[i][j]=c; C[i][j]=0,C[i][j+1]=0; C[i][j]=b; C[i][j]=0; C[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); C[i][j]=b; C[i][j]=c; else t=v2; if(t<v3) t=v3; return(t); return ((v1<v2)?v1:v2);

}

//用于求解最大的特征值，幂法

double mifa(double C[5][501])

{

int m=0,i,j; double b2,b1=0,sum; double u[501],y[501]; for (i=0;i<501;i++)

{ u[i] = 1.0;

}

do { sum=0;

}

/*****行列式LU分解*****/

void LU(double C[5][501])

{

}

/***带状DOOLITE分解，并且求解出方程组的解***/

void solve(double C[5][501],double x[501],double b[501])

{

int i,j,k,t; double B[5][501],c[501]; for(i=0;i<=4;i++) { for(j=0;j<=500;j++) B[i][j]=C[i][j]; double sum; int k,i,j; for(k=1;k<=501;k++) { for(j=k;j<=Min(k+s,501);j++) } { sum=0; for(i=maxt(1,k-r,j-s);i<=k-1;i++) } { sum=0; for(i=maxt(1,j-r,k-s);i<=k-1;i++) } sum+=C[j-i+s][i-1]*C[i-k+s][k-1]; C[j-k+s][k-1]=(C[j-k+s][k-1]-sum)/C[s][k-1]; sum+=C[k-i+s][i-1]*C[i-j+s][j-1]; C[k-j+s][j-1]-=sum; } while(fabs(b2-b1)/fabs(b2)>=1.0e-12); return b2; m++; for(i=0;i<=500;i++) sum+=u[i]*u[i]; y[i]=u[i]/sqrt(sum); for(i=0;i<=500;i++) for(i=0;i<=500;i++) { u[i]=0; } b2=0; for(i=0;i<=500;i++) b2=b2+y[i]*u[i]; for(j=Max(i-r,0);j<=Min(i+s,500);j++) u[i]=u[i]+C[i-j+s][j]*y[j]; for(j=k+1;j<=Min(k+r,501);j++)

}

//用于求解模最大的特征值，反幂法

double fmifa(double C[5][501])

{

int m=0,i; double b2,b1=0,sum=0,u[501],y[501]; for (i=0;i<=500;i++) { [i] = 1.0; } do { if(m!=0)b1=b2; m++; sum=0; for(i=0;i<=500;i++) sum+=u[i]*u[i]; y[i]=u[i]/sqrt(sum); for(i=0;i<=500;i++) solve(C,u,y); b2=0; for(i=0;i<=500;i++) b2+=y[i]*u[i]; for(i=0;i<=500;i++) { } for(i=1;i<=500;i++) for(t=Max(0,i-r);t<=i-1;t++) c[i]=c[i]-B[i-t+s][t]*c[t]; x[500]=c[500]/B[s][500]; for(i=499;i>=0;i--) { } x[i]=c[i]; for(t=i+1;t<=Min(i+s,500);t++) x[i]=x[i]-B[i-t+s][t]*x[t]; x[i]=x[i]/B[s][i]; c[i]=b[i]; for(j=k;j<=Min(k+s,500);j++) { } for(i=k+1;i<=Min(k+r,500);i++) { for(t=Max(0,Max(i-r,k-s));t<=k-1;t++) } B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k]; for(t=Max(0,Max(k-r,j-s));t<=k-1;t++) B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; for(k=0;k<=500;k++)

}

/***主程序***/

void main()

{

} double b=0.16,c=-0.064,det=1.0; int i; double C[5][501],cond; storage(C,b,c); //进行C的赋值 cout.precision(12); //定义输出精度 double k1=mifa(C); //利用幂法计算矩阵的最大特征值和最小特征值 if(k1<0) else if(k1>=0) for(i=0;i<501;i++) C[2][i]=C[2][i]-k1; double k2=mifa(C)+k1; if(k2<0) printf("λ1=%.12e\n",k2); printf("λ501=%.12e\n",k2); else if(k2>=0) storage(C,b,c); double k3=fmifa(C); //利用反幂法计算矩阵A的按模最小特征值 storage(C,b,c); //计算最接近特征值 double u[39]={0}; for(i=0;i<39;i++) { u[i]=k1+(i+1)*(k2-k1)/40; } if(k1>0) //计算矩阵A的条件数，取2范数 cond=fabs(k1/k3); C[2][i]=C[2][i]-u[i]; u[i]=fmifa(C)+u[i]; return 1/b2; printf("λ1=%.12e\n",k1); printf("λ501=%.12e\n",k1); printf("λs=%.12e\n",k3); printf("与数u%d 最接近的特征值λ%d: %.12e\n",i+1,i+1,u[i]); else if(k1<0) cond=fabs(k2/k3); storage(C,b,c); LU(C); //利用LU分解计算矩阵A的行列式 for(i=0;i<501;i++) det*=C[2][i]; printf("\ncond(A)=%.12e\n",cond); printf("\ndet(A)=%.12e\n",det);

三、 计算结果：

四、结果分析

迭代初始向量的选择对果有一定的影响，选择不同的初始向量可能会得到不同阶的特征值。

若改初始向量u[1]=1，u[i]=0，(i=2,3,...,501)。则得出结果

此结果与正确结果相差较多。

…… 此处隐藏：2590字，全部文档内容请下载后查看。喜欢就下载吧 ……