数值分析计算实习及遇到的问题分析

以上是计算实习的题目。

说说我遇到的问题吧。

1、初始向量的选择对迭代结果的影响问题

我在刚开始做第一题时就遇到了一个大问题：得到的结果虽然收敛，但是却与正确答案不相符。刚开始我设置的初始迭代向量 为(1,0,0,……,0)，用幂法计算λ得到的结果是-0.208098108534e+01，迭代次数为280次。这个错误几乎令我想要放弃，因为我已经通过在代码的各个地方设置断点不断缩小可能出错的地方，但是直到第二天我仍没有找到错在哪。直到我上网看到了幂法的收敛与否不仅与特征值分布有关，还与初始迭代向量有关，于是我修改了迭代向量 为(1,1,1,……,1)，得到的答案果然完全正确！之后的题目就可以说是势如破竹，在很短的时间内就完成了。

2、关于预编译时发生堆栈溢出的问题

第一次调试，系统就提示Stack overflow，随后编译停留在了汇编指令test dword ptr[ecx],eax上。于是我重新编译，采取单步执行，随后报错：please enter the file for :CHKSTK.ASM。上网查找这两个错误后，得知原因是声明的函数大量使用了传值，要想解决问题，必须将函数的传参方式尽可能改为传递指针（引用）。我发现在第一次编写代码时，例如矩阵与向量相乘的函数，我对函数参数的定义就是直接将矩阵里所有的值传递到函数中，若改为将指针传进函数，将解决预编译栈溢出的问题。

3、向量可重复利用的问题

修改这个问题之前，我声明了一个大小为MAXSTEP=10000的向量数组μ，每个向量中又有501个double型的分量，占用了大量的空间，导致调试时发生了卡顿现象，因为内存被分配了太多。

幂法中用到的向量μ虽然确实有迭代步数k次，但经过一步后，前一步用过的向量将不再会用到，所以只需要声明两个向量 进行交替存储即可。

以下是最终输出的答案：

lamda1的值为-0.107001136151E+2。

lamda501的值为0.972463409986E+1。

lamda_s的值为-0.555791079422E-2。

lamda_i1的值为-0.101829340331E+2。

lamda_i2的值为-0.958570742507E+1。

lamda_i3的值为-0.917267242393E+1。

lamda_i4的值为-0.865228400790E+1。

lamda_i5的值为-0.809348380868E+1。

lamda_i6的值为-0.765940540769E+1。

lamda_i7的值为-0.711968464869E+1。

lamda_i8的值为-0.661176433940E+1。

lamda_i9的值为-0.606610322660E+1。

lamda_i10的值为-0.558510105263E+1。

lamda_i11的值为-0.511408352981E+1。

lamda_i12的值为-0.457887217687E+1。

lamda_i13的值为-0.409647092626E+1。

lamda_i14的值为-0.355421121575E+1。

lamda_i15的值为-0.304109001813E+1。

lamda_i16的值为-0.253397031113E+1。

lamda_i17的值为-0.200323076956E+1。

lamda_i18的值为-0.150355761123E+1。

lamda_i19的值为-0.993558606008E-0。

lamda_i20的值为-0.487042673885E-0。

lamda_i21的值为0.223173624957E-1。

lamda_i22的值为0.532417474207E-0。

lamda_i23的值为0.105289896269E+1。

lamda_i24的值为0.158944588188E+1。

lamda_i25的值为0.206033046027E+1。

lamda_i26的值为0.255807559707E+1。

lamda_i27的值为0.308024050931E+1。

lamda_i28的值为0.361362086769E+1。

lamda_i29的值为0.409137851045E+1。

lamda_i30的值为0.460303537828E+1。

lamda_i31的值为0.513292428390E+1。

lamda_i32的值为0.559490634808E+1。

lamda_i33的值为0.608093385703E+1。

lamda_i34的值为0.668035409211E+1。

lamda_i35的值为0.729387744813E+1。

lamda_i36的值为0.771711171424E+1。

lamda_i37的值为0.822522001405E+1。

lamda_i38的值为0.864866606519E+1。

lamda_i39的值为0.925420034458E+1。

CondA2的值为0.192520427393E+4。

detA的值为0.277278614175E+119。

以下是实现代码：

#include "c1.h"
#include "struct.h"

double powerMethod(double lamda);//幂法求A-λI的主特征值
double reversePowerMethod(double lamda);//反幂法求A-μk按模最小的特征值
double det();//使用LU分解求矩阵A的特征值

void main()
{
freopen("output.txt","w",stdout);
eTypeRealNumber Answer1_1,Answer1_2,Answer1_3,Answer2_1[39],Answer3_1,Answer3_2;
double lamda=powerMethod(0);
double lamda_=powerMethod(lamda);
double lamda1=lamda<0?lamda:lamda_;
double lamda501=lamda>0?lamda:lamda_;
double lamda_s=reversePowerMethod(0);
double lamda_ik[39];
for(int i=0;i<39;i++)
{
lamda_ik[i]=reversePowerMethod(UK);
}
double CondA2=dabs(lamda)/dabs(lamda_s);
double detA=det();
//转换
Answer1_1=doubleToEtype(lamda1);
Answer1_2=doubleToEtype(lamda501);
Answer1_3=doubleToEtype(lamda_s);
for(i=0;i<39;i++)
{
Answer2_1[i]=doubleToEtype(lamda_ik[i]);
}
Answer3_1=doubleToEtype(CondA2);
Answer3_2=doubleToEtype(detA);
//输出
printf("lamda1的值为%.12lfE%c%d\n",Answer1_1.bottom,Answer1_1.power>0?'+':'-',abs(Answer1_1.power));
printf("lamda501的值为%.12lfE%c%d\n",Answer1_2.bottom,Answer1_2.power>0?'+':'-',abs(Answer1_2.power));
printf("lamda_s的值为%.12lfE%c%d\n",Answer1_3.bottom,Answer1_3.power>0?'+':'-',abs(Answer1_3.power));
for(i=0;i<39;i++)
{
printf("lamda_i%d的值为%.12lfE%c%d\n",i+1,Answer2_1[i].bottom,Answer2_1[i].power>0?'+':'-',abs(Answer2_1[i].power));
}
printf("CondA2的值为%.12lfE%c%d\n",Answer3_1.bottom,Answer3_1.power>0?'+':'-',abs(Answer3_1.power));
printf("detA的值为%.12lfE%c%d\n",Answer3_2.bottom,Answer3_2.power>0?'+':'-',abs(Answer3_2.power));
}

double powerMethod(double lamda)//幂法求A-λI的主特征值
{
//
//主对角线的所有ai的定义如下
//ai=(1.64-0.024i)sin(0.2i)-0.64e^{0.1/i};(i=1 to 501);
//
RealSymmetricMatricePoint A=(RealSymmetricMatricePoint)malloc(sizeof(RealSymmetricMatrice));
int size=MAXSIZE;//以下代码中所有矩阵的阶数和向量的维数均为size;
for(int i=1;i<=size;i++)//初始化主对角线ai的值
{
A->mainDiagonal[i-1]=AI-lamda;
}
A->b=B;
A->c=C;

/*接下来对A使用幂法求得λ*/
Vector *u1,*u2;//向量u
Vector *y1;//向量y
u1=(Vector*)malloc(sizeof(Vector));u1->partNumber=(double*)malloc(size*sizeof(double));
y1=(Vector*)malloc(sizeof(Vector));y1->partNumber=(double*)malloc(size*sizeof(double));
u2=(Vector*)malloc(sizeof(Vector));u2->partNumber=(double*)malloc(size*sizeof(double));
for(i=0;i<size;i++)//初始化u0为(1,1,....,1)
{
u1->partNumber[i]=1;
}

double beta[MAXSTEP];//标量β
double tempAnswer1;
double hrk2=u1->partNumber[maxIndexinVector(u1,size)];//向量u1中最大的分量
double hrk1=0;
for(i=1,beta[0]=0;i<MAXSTEP;i++)//对A用幂法迭代的第二种格式
{
*y1=vectorMultipleConst(u1,(double)(1)/dabs(hrk2),size);
*u2=RSMattriceMultipleVector(A,y1,size);
hrk1=hrk2;
hrk2=u2->partNumber[maxIndexinVector(u1,size)];
beta[i]=(hrk1>0?1:-1)*hrk2;
if(i>2&&(dabs(beta[i]-beta[i-1])
4000
/dabs(beta[i]))<=KSI)
{
break;
}
u1=u2;
}
tempAnswer1=beta[i-1]+lamda;
return tempAnswer1;
}

double reversePowerMethod(double miu)//反幂法求A-μk按模最小的特征值
{
RealSymmetricMatricePoint A=(RealSymmetricMatricePoint)malloc(sizeof(RealSymmetricMatrice));
int size=MAXSIZE;
for(int i=1;i<=size;i++)//初始化主对角线ai的值
{
A->mainDiagonal[i-1]=AI-miu;
}
A->b=B;
A->c=C;
double L[3][501]={0};
double U[501][3]={0};
LU(A->mainDiagonal,A->b,A->c,L,U);

Vector *u1,*u2;//向量u
Vector *y1;//向量y
u1=(Vector*)malloc(sizeof(Vector));u1->partNumber=(double*)malloc(size*sizeof(double));
y1=(Vector*)malloc(sizeof(Vector));y1->partNumber=(double*)malloc(size*sizeof(double));
u2=(Vector*)malloc(sizeof(Vector));u2->partNumber=(double*)malloc(size*sizeof(double));

double beta[MAXSTEP];//标量β
double tempAnswer1;

for(i=0;i<size;i++)//初始化u0为(1,1,....,1)
{
u1->partNumber[i]=1;
}
for(i=1,beta[0]=0;i<MAXSTEP;i++)
{
*y1=vectorMultipleConst(u1,(double)(1)/sqrt(vectorMultipleVector(u1,u1,size)),size);
findUk(L,U,y1,u2);
beta[i]=vectorMultipleVector(y1,u2,size);
if(i>2&&(dabs(beta[i]-beta[i-1])/dabs(beta[i]))<=KSI)
{
break;
}
u1=u2;
}
tempAnswer1=(double)(1.0)/beta[i-1]+miu;
return tempAnswer1;
}
double det()//使用LU分解求矩阵A的特征值
{
double ret=1;
RealSymmetricMatricePoint A=(RealSymmetricMatricePoint)malloc(sizeof(RealSymmetricMatrice));
int size=MAXSIZE;
for(int i=1;i<=size;i++)//初始化主对角线ai的值
{
A->mainDiagonal[i-1]=AI;
}
A->b=B;
A->c=C;
double L[3][MAXSIZE]={0};
double U[MAXSIZE][3]={0};
LU(A->mainDiagonal,A->b,A->c,L,U);//LU分解

for(i=0;i<MAXSIZE;i++)
{
ret*=U[i][0];
}
return ret;
}

以下是常用函数文件struct.h

#define MAXSIZE 501//最大阶（维）数
#define MAXSTEP 10000//最大迭代次数
#define KSI 0.0000000000001//允许误差ε
#define AI (1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/(double)(i))//ai
#define UK lamda1+(i+1)*(lamda501-lamda1)/40//μk
#define B 0.16//b
#define C -0.064//c
/*存储类型*/
typedef struct RealSymmetricMatrice//用于本题的特殊实对称矩阵存储结构
{
double mainDiagonal[MAXSIZE];//主对角线上的元素，待初始化
double b;
double c;
}*RealSymmetricMatricePoint;
typedef struct Vector//定义向量
{
double *partNumber;//分量
}*VectorPoint;
typedef struct eTypeRealNumber//定义要求输出的形式
{
double bottom;//底数
int power;//幂
}eTypeRealNumber;
/*函数*/
double dabs(double a)//自定义double型求绝对值函数
{
if(a<0)
return -a;
else
return a;
}
eTypeRealNumber doubleToEtype(double number)//将double型转换为要求输出的形式
{
double bottom=dabs(number);
int power=0;
eTypeRealNumber ret;
while(bottom<0.1||bottom>=1)
{
if(bottom<0.1)
{
power--;
bottom*=10;
}
else
{
power++;
bottom/=10;
}
}
ret.bottom=(number>0?1:-1)*bottom;
ret.power=power;
return ret;
}
double vectorMultipleVector(Vector* a,Vector* b,int size)//向量乘向量，返回标量
{
double ret=0;
for(int i=0;i<size;i++)
{
ret+=a->partNumber[i]*b->partNumber[i];
}
return ret;
}
Vector vectorMultipleConst(Vector* a,double b,int size)//向量乘以常数，返回向量
{
Vector ret;
ret.partNumber=(double*)malloc(size*sizeof(double));
for(int i=0;i<size;i++)
{
ret.partNumber[i]=a->partNumber[i]*b;
}
return ret;
}
Vector RSMattriceMultipleVector(RealSymmetricMatrice* A,Vector* b,int size)//特殊矩阵乘以向量，返回向量
{
Vector ret;
ret.partNumber=(double*)malloc(size*sizeof(double));
for(int i=0;i<size;i++)
{
double sum=0;
for(int j=-2;j<=2;j++)
{
if(i+j<0||i+j>=size)
{
continue;
}
switch(j)
{
case -2:sum+=A->c*b->partNumber[i+j];break;
case -1:sum+=A->b*b->partNumber[i+j];break;
case 0:sum+=A->mainDiagonal[i]*b->partNumber[i];break;
case 1:sum+=A->b*b->partNumber[i+j];break;
case 2:sum+=A->c*b->partNumber[i+j];break;
}
}
ret.partNumber[i]=sum;
}
return ret;
}

int maxIndexinVector(Vector *a,int size)//寻找向量中所有分量求绝对值后最大分量的编号
{
int maxIndex=0;
double max=0;
for(int i=0;i<size;i++)
{
if(dabs(a->partNumber[i])>max)
{
max=dabs(a->partNumber[i]);
maxIndex=i;
}
}
return maxIndex;
}

void LU(double a[],double b, double c, double (*L)[MAXSIZE] , double (*U)[3])//计算特殊矩阵A的矩阵的LU分解
{
U[0][0]=a[0];
U[0][1]=b;
U[0][2]=c;
L[1][0]=b/U[0][0];
L[2][0]=c/U[0][0];
U[1][0]=a[1]-L[1][0]*U[0][1];

U[1][1]=b-L[1][0]*U[0][2];
U[1][2]=c;
L[1][1]=(b-L[2][0]*U[0][1])/U[1][0];
L[2][1]=c/U[1][0];

for (int k = 2; k < (MAXSIZE-2); ++k)
{
U[k][0]=a[k]-L[2][k-2]*U[k-2][2]-L[1][k-1]*U[k-1][1];
U[k][1]=b-L[1][k-1]*U[k-1][2];
U[k][2]=c;
L[1][k]=(b-L[2][k-1]*U[k-1][1])/U[k][0];
L[2][k]=c/U[k][0];
}
k=(MAXSIZE-2);
U[k][0]=a[k]-L[2][k-2]*U[k-2][2]-L[1][k-1]*U[k-1][1];
U[k][1]=b-L[1][k-1]*U[k-1][2];
L[1][k]=(b-L[2][k-1]*U[k-1][1])/U[k][0];
k=(MAXSIZE-1);
U[k][0]=a[k]-L[2][k-2]*U[k-2][2]-L[1][k-1]*U[k-1][1];
}
void findUk(double (*L)[MAXSIZE], double (*U)[3], VectorPoint b, VectorPoint u2)//由LU分解所得矩阵数组L,U,根据yk-1反解uk
{
double y[MAXSIZE];
for (int i = 0; i < MAXSIZE; ++i)//解方程Ly=b
{
y[i]=b->partNumber[i];
for (int t = (0<(i-2)?(i-2):0); t <= i-1; ++t)
{
y[i]=y[i]-L[i-t][t]*y[t];
}
}
u2->partNumber[MAXSIZE-1]=y[MAXSIZE-1]/U[MAXSIZE-1][0];//解方程Ux=y
for (i = (MAXSIZE-2); i >=0; --i)
{
u2->partNumber[i]=y[i];
int z=((i+2)<MAXSIZE-1?(i+2):MAXSIZE-1);
for (int t = i+1; t <= ((i+2)<(MAXSIZE-1)?(i+2):MAXSIZE-1); ++t)
{
u2->partNumber[i]=u2->partNumber[i]-U[i][t-i]*u2->partNumber[t];
}
u2->partNumber[i]=u2->partNumber[i]/U[i][0];
}
}

c1.h /* c1.h (程序名) */
#include<string.h>
#include<ctype.h>
#include<malloc.h> /* malloc()等 */
#include<limits.h> /* INT_MAX等 */
#include<stdio.h> /* EOF(=^Z或F6),NULL */
#include<stdlib.h> /* atoi() */
#include<io.h> /* eof() */
#include<math.h> /* floor(),ceil(),abs() */
#include<process.h> /* exit() */
/* 函数结果状态代码 */
#define TRUE 1
#define FALSE 0
#define OK 1
#define ERROR 0
#define INFEASIBLE -1
/* #define OVERFLOW -2 因为在math.h中已定义OVERFLOW的值为3,故去掉此行 */
typedef int Status; /* Status是函数的类型,其值是函数结果状态代码，如OK等 */
typedef int Boolean; /* Boolean是布尔类型,其值是TRUE或FALSE */