您所在位置：主题：急啊！大作业要交了，可程序编译没问题，运行就出错，哪位仁兄给看一下啊
dong1982628
[专家分：0] 发布于 2006-11-11 22:19:00
原题是这样的，设有501*501的矩阵
A= a1 b c
b a2 b c
c b a3 b c
.......
c b a499 b c
c b a500 b
c b a501
其中
ai=（1.64-0.024i）sin（0.2i）-0.64exp(0.1/i)
(i=1,2,...,501);b=0.16,c=-0.064.矩阵A的特征值lamda（i）（i=1，2,...，501）满足
lamda1<lamda2<...<lamda501,|lamdas|=mim|lamdai|试求1.lamda1,lamda501,lamdas
2.A的与数miu（k）=lamda1+k(lamda501-lamda1)/40最接近的特征值lamda（ik）（k1,2,...,39)
3.A的（谱范数）条件数cond（A)2和行列式detA
应该用幂法与反幂法求最大最小特征值和模最小特征值，以及与某数距离最近的特征值，要求39个，我编的程序如下，可运行时有问题,哪位仁兄给看一下啊
#include"stdio.h"
#include"math.h"
void main()
{double MI(double Array[][501]);
double FMI(double A[][501]);
double C[5][501],B[5][501],E[39],L,s,cond,lamda,lamda0,lamda1,lamdas,lamda501;int i,j,k;
for(i=0;i<5;i=i+4) //稀疏矩阵的存储
for(j=0;j<501;j++)
C[i][j]=-0.064;
for(i=1;i<5;i=i+2)
for(j=0;j<501;j++)
C[i][j]=0.16;
for(j=1;j<502;j++)
C[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);
lamda0=MI(C); //求lamda1，lamda501
for(i=0;i<5;i++)
for(j=0;j<501;j++)
B[i][j]=C[i][j];
for(j=0;j<501;j++)
B[2][j]-=lamda0;
lamda=MI(B);
if(lamda>0) lamda501=lamda0;lamda1=lamda0+lamda;
if(lamda<0) lamda1=lamda0;lamda501=lamda0+lamda;
for(k=1;k<40;k++) //求lamda(ik);
{ L=lamda1+k*(lamda501-lamda1)/40;
for(j=0;j<501;j++)
B[2][j]=C[2][j]-L;
E[k-1]=FMI(B);E[k-1]=E[k-1]+L;
}
lamdas=FMI(C);lamdas=fabs(lamdas); //求lamdas
cond=fabs(lamda0/lamdas); //求cond(A)
for(i=0;i<501;i++) //求detA
{s=1; s=s*C[2][i];}
printf("%lf,%lf,%lf,%lf,%lf",lamda1,lamda501,lamdas,cond,s);
for(i=0;i<40;i++)
printf("%lf\n",E[i]);
}
double MI(double Array[][501]) //幂法子函数
{ int i,j;double beta0=0.0,beta=0.0,yita=0.0,u[501], y[501],s;
for(i=0;i<501;i++) u[i]=1; //u0;
while(1)
{ s=0.0;
for(i=0;i<501;i++)
s+=u[i]*u[i];
yita=sqrt(s);
for(i=0;i<501;i++)
y[i]=u[i]/yita;
for(j=0;j<3;j++) //uk=A*y(k-1)
u[0]+=Array[2-j][j]*y[j];
for(j=0;j<4;j++)
u[1]+=Array[3-j][j]*y[j];
for(i=2;i<499;i++)
for(j=i-2;j<3+i;j++)
u[i]+=Array[2+i-j][j]*y[j];
for(j=497;j<501;j++)
u[499]+=Array[501-j][j]*y[j];
for(j=498;j<501;j++)
u[500]+=Array[502-j][j]*y[j];
for(i=0;i<501;i++)
beta+=y[i]*u[i];
if (fabs(beta-beta0)/fabs(beta)<=(1.0e-12))
{break;return(beta);}
beta0=beta;beta=0.0;
}
}
double FMI(double A[][501]) //反幂法子函数
{int i,j,k,t;double beta0=0.0,beta=0.0,yita=0.0,s,x;double u[501], y[501],yy[501];
for(k=0;k<501;k++) //LU分解
{if(k+1<500) //L矩阵
{for(j=k;j<=k+2;j++)
if(j==1) {for(t=0;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
else {for(t=j-2;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
}
else{ for(j=k;j<=500;j++)
for(t=j-2;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
if((k+2<=500)&&(k<500)) //U矩阵
{for(i=k+1;i<=k+2;i++)
{for(t=i-2;t<=k-1;t++)
A[i-k+2][k]-=A[i-t+2][t]*A[t-k+2][k];
A[i-k+2][k]/=A[2][k];}
}
else{for(t=i-1;t<=k;t++)
A[i-k+2][k]-=A[i-t+2][t]*A[t-k+2][k];
A[i-k+2][k]/=A[2][k];
}
}
for(i=0;i<501;i++) u[i]=1; //u0;
while(1)
{s=0.0;
for(i=0;i<501;i++)
s+=u[i]*u[i];
yita=sqrt(s);
for(i=0;i<501;i++)
{y[i]=u[i]/yita;
yy[i]=y[i];}
y[1]=y[1]-A[3][0]*y[0];
for(i=2;i<=500;i++)
for(t=i-2;t<=i-1;i++)
y[i]-=A[i-t+2][t]*y[t];
u[500]=y[500]/A[2][500];
for(i=499;i>=0;i--)
{if(i+2<500)
for(t=i+1;t<=i+2;t++)
{u[i]=y[i]-A[i-t+2][t]*u[t];
u[i]/=A[2][i];
}
else for(t=i+1;t<=500;t++)
{u[i]=y[i]-A[i-t+2][t]*u[t];
u[i]/=A[2][i];
}
}
for(i=0;i<501;i++)
beta+=yy[i]*u[i];
x=fabs(beta-beta0)/fabs(beta);
if(x<=(1.0e-12))
{return(beta);break;}
beta0=beta;beta=0.0;
}
}
运行结果：Loaded 'ntdll.dll', no matching symbolic information found.
Loaded 'C:\WINDOWS\system32\kernel32.dll', no matching symbolic information found.
The thread 0xE38 has exited with code -1073741510 (0xC000013A).
The thread 0xC7C has exited with code -1073741510 (0xC000013A).
The program 'F:\Debug\shu.exe' has exited with code -1073741510 (0xC000013A).
A= a1 b c
b a2 b c
c b a3 b c
.......
c b a499 b c
c b a500 b
c b a501
其中
ai=（1.64-0.024i）sin（0.2i）-0.64exp(0.1/i)
(i=1,2,...,501);b=0.16,c=-0.064.矩阵A的特征值lamda（i）（i=1，2,...，501）满足
lamda1<lamda2<...<lamda501,|lamdas|=mim|lamdai|试求1.lamda1,lamda501,lamdas
2.A的与数miu（k）=lamda1+k(lamda501-lamda1)/40最接近的特征值lamda（ik）（k1,2,...,39)
3.A的（谱范数）条件数cond（A)2和行列式detA
应该用幂法与反幂法求最大最小特征值和模最小特征值，以及与某数距离最近的特征值，要求39个，我编的程序如下，可运行时有问题,哪位仁兄给看一下啊
#include"stdio.h"
#include"math.h"
void main()
{double MI(double Array[][501]);
double FMI(double A[][501]);
double C[5][501],B[5][501],E[39],L,s,cond,lamda,lamda0,lamda1,lamdas,lamda501;int i,j,k;
for(i=0;i<5;i=i+4) //稀疏矩阵的存储
for(j=0;j<501;j++)
C[i][j]=-0.064;
for(i=1;i<5;i=i+2)
for(j=0;j<501;j++)
C[i][j]=0.16;
for(j=1;j<502;j++)
C[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);
lamda0=MI(C); //求lamda1，lamda501
for(i=0;i<5;i++)
for(j=0;j<501;j++)
B[i][j]=C[i][j];
for(j=0;j<501;j++)
B[2][j]-=lamda0;
lamda=MI(B);
if(lamda>0) lamda501=lamda0;lamda1=lamda0+lamda;
if(lamda<0) lamda1=lamda0;lamda501=lamda0+lamda;
for(k=1;k<40;k++) //求lamda(ik);
{ L=lamda1+k*(lamda501-lamda1)/40;
for(j=0;j<501;j++)
B[2][j]=C[2][j]-L;
E[k-1]=FMI(B);E[k-1]=E[k-1]+L;
}
lamdas=FMI(C);lamdas=fabs(lamdas); //求lamdas
cond=fabs(lamda0/lamdas); //求cond(A)
for(i=0;i<501;i++) //求detA
{s=1; s=s*C[2][i];}
printf("%lf,%lf,%lf,%lf,%lf",lamda1,lamda501,lamdas,cond,s);
for(i=0;i<40;i++)
printf("%lf\n",E[i]);
}
double MI(double Array[][501]) //幂法子函数
{ int i,j;double beta0=0.0,beta=0.0,yita=0.0,u[501], y[501],s;
for(i=0;i<501;i++) u[i]=1; //u0;
while(1)
{ s=0.0;
for(i=0;i<501;i++)
s+=u[i]*u[i];
yita=sqrt(s);
for(i=0;i<501;i++)
y[i]=u[i]/yita;
for(j=0;j<3;j++) //uk=A*y(k-1)
u[0]+=Array[2-j][j]*y[j];
for(j=0;j<4;j++)
u[1]+=Array[3-j][j]*y[j];
for(i=2;i<499;i++)
for(j=i-2;j<3+i;j++)
u[i]+=Array[2+i-j][j]*y[j];
for(j=497;j<501;j++)
u[499]+=Array[501-j][j]*y[j];
for(j=498;j<501;j++)
u[500]+=Array[502-j][j]*y[j];
for(i=0;i<501;i++)
beta+=y[i]*u[i];
if (fabs(beta-beta0)/fabs(beta)<=(1.0e-12))
{break;return(beta);}
beta0=beta;beta=0.0;
}
}
double FMI(double A[][501]) //反幂法子函数
{int i,j,k,t;double beta0=0.0,beta=0.0,yita=0.0,s,x;double u[501], y[501],yy[501];
for(k=0;k<501;k++) //LU分解
{if(k+1<500) //L矩阵
{for(j=k;j<=k+2;j++)
if(j==1) {for(t=0;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
else {for(t=j-2;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
}
else{ for(j=k;j<=500;j++)
for(t=j-2;t<=k-1;t++)
A[k-j+2][j]-=A[k-t+2][t]*A[t-j+2][j];}
if((k+2<=500)&&(k<500)) //U矩阵
{for(i=k+1;i<=k+2;i++)
{for(t=i-2;t<=k-1;t++)
A[i-k+2][k]-=A[i-t+2][t]*A[t-k+2][k];
A[i-k+2][k]/=A[2][k];}
}
else{for(t=i-1;t<=k;t++)
A[i-k+2][k]-=A[i-t+2][t]*A[t-k+2][k];
A[i-k+2][k]/=A[2][k];
}
}
for(i=0;i<501;i++) u[i]=1; //u0;
while(1)
{s=0.0;
for(i=0;i<501;i++)
s+=u[i]*u[i];
yita=sqrt(s);
for(i=0;i<501;i++)
{y[i]=u[i]/yita;
yy[i]=y[i];}
y[1]=y[1]-A[3][0]*y[0];
for(i=2;i<=500;i++)
for(t=i-2;t<=i-1;i++)
y[i]-=A[i-t+2][t]*y[t];
u[500]=y[500]/A[2][500];
for(i=499;i>=0;i--)
{if(i+2<500)
for(t=i+1;t<=i+2;t++)
{u[i]=y[i]-A[i-t+2][t]*u[t];
u[i]/=A[2][i];
}
else for(t=i+1;t<=500;t++)
{u[i]=y[i]-A[i-t+2][t]*u[t];
u[i]/=A[2][i];
}
}
for(i=0;i<501;i++)
beta+=yy[i]*u[i];
x=fabs(beta-beta0)/fabs(beta);
if(x<=(1.0e-12))
{return(beta);break;}
beta0=beta;beta=0.0;
}
}
运行结果：Loaded 'ntdll.dll', no matching symbolic information found.
Loaded 'C:\WINDOWS\system32\kernel32.dll', no matching symbolic information found.
The thread 0xE38 has exited with code -1073741510 (0xC000013A).
The thread 0xC7C has exited with code -1073741510 (0xC000013A).
The program 'F:\Debug\shu.exe' has exited with code -1073741510 (0xC000013A).