北航 数值分析B 作业题一研究生.docx
《北航 数值分析B 作业题一研究生.docx》由会员分享，可在线阅读，更多相关《北航 数值分析B 作业题一研究生.docx（29页珍藏版）》请在冰点文库上搜索。
北航数值分析B作业题一研究生
《数值分析》计算实习题目
目录
一算法的设计方案1
二全部源代码2
2.1“daizhuang.cpp”中源代码2
2.2“matrix-tool.cpp”中源代码6
2.3”convert_mufa.cpp”中源代码8
2.4”doolittle.cpp”中源代码10
2.5”mufa.cpp”中源代码12
2.6“main.cpp”的源代码14
三：
特征值以及
，行列式的值16
四：
结果分析与讨论17
一算法的设计方案
根据题目要求一，可以使用幂法求出矩阵A的最小特征值
，而通过原点平移原理可以求得矩阵A的最大特征值
，而
则可使用反幂法求得。
同理，通过原点平移原理，我们也可以得到与数
最接近的特征值
。
这样就满足了要求二。
而根据矩阵条件数的定义
，只要知道了模值为最大和最小的特征值便可以求得。
而将A通过QR分解，便可以得到A所有的特征值（位于上三角矩阵的主对角线上）。
由题目要求，矩阵A的（s=2,r=2）所以可以选用矩阵
来存储矩阵A中的非零元素。
在算法实现过程中，我使用了5个头文件和相应的源文件来完成不同的矩阵操作，这5个源文件和头文件分别为：
”convert_mufa.h”.”convert_mufa.c””doolittle.h”.”doolittle.c”.“matrix-tool.h”.“matrix-tool.c””mufa.h”.”mufa.c”“daizhuang.h”.”daizhuang.c”
其中：
”mufa.h”，”mufa.c”中的函数为：
voiditeration（double*error1）;
这个函数完成幂法的一次迭代过程
doublemufa_root（）;
这个函数通过幂法求得模最大的特征值
”convert_mufa.h”.”convert_mufa.c”中的函数为：
voidConvert_Iteration（double*error）;
这个函数完成反幂法的迭代操作。
doubleConvert_root（）;
这个函数完成反幂法的特征根求解。
”doolittle.h”.”doolittle.c”.中的函数为：
intinput（intn,double*a,double*b）;
这个函数完成原始系统矩阵a（n*n）的数据、维数的输入，方程右边向量b（n）的输入。
intDLLUdecompostion（double*a,intn）;
这个函数是进行doolittle的LU分解。
intROOT（intn,double*a,double*b）;
这个函数进行求根运算
voiddisplay（intn,double*a,double*b）;
这个函数用于显示处理前和处理后的矩阵。
“daizhuang.h”.”daizhuang.c”中的函数为：
voidINITIALDAIZHUANG（）;
这个函数主要用于对A[501][501]的存储矩阵C[（1+r+s）][n]的初始化和按题目要求赋值。
intTRI_DLLUdecompostion（）;
这个函数完成对拟三角矩阵的DOLITTILE分解。
intTRI_ROOT（double*xb,double*B）;
这个函数完成对拟上三角矩阵方程组的求根运算。
voidTRI_display（）;
这个函数完成对拟上三角矩阵的显示，以方便调试。
voidTRI_ROOT_display（double*xb）;
这个函数完成对拟上三角矩阵方程组的根显示。
“matrix-tool.h”.“matrix-tool.c”.中的函数主要实现矩阵和向量的一般操作，方便其他函数调用。
函数的定义和作用如下所示：
voidMatrix_M_Vector（double*A,double*b,double*y,intn）;
函数实现N*N的矩阵乘以N*1的向量，得到N*1向量Y
voidMatrix_Converted_M_Vector（double*A,double*b,double*y,intn）;
函数实现转置后N*N矩阵乘以N*1向量,得到N*1向量Y
voidColumn_M_Row（double*a,double*b,double*A,intn）;
函数实现列乘以行得到一个矩阵A
voidRow_M_Column（double*a,double*b,double*y,intn）;
函数实现行乘以列得到一个数值y
voidPrint_Matrix（double*a,intn）;
函数实现将矩阵a显示出来
voidPrint_Vector（double*b,intn）;
函数实现将向量b显示出来
二全部源代码
2.1“daizhuang.cpp”中源代码
#include"daizhuang.h"
#include"iostream"
#include"math.h"
#include"matrix-tool.h"
usingnamespacestd;
//带状线性方程组中相关变量的定义
chars=2;
charr=2;
intn=501;//后来可以修改成为501
doubleformula=0;
double*c,*d;
doubleb[501];
doublexa[501];
voidINITIALDAIZHUANG（）
{
//把A[N][N]中非零数据放入C[M][N]中进行存储，以下存储方式为题目特定方式
c=newdouble[（1+r+s）*n];
d=newdouble[（1+r+s）*n];//得把c的地址赋给d
inti,j;
//把c中元素全部初始化
for（i=0;i{
for（j=0;j{
if（i==j）
{
formul=（（1.64-0.024*（i+1））*（double）sin（0.2*（i+1））-0.64*（double）exp（0.1/（i+1）））;
（*（c+（i-j+s）*n+j））=formula;
formula=0;
}
elseif（abs（i-j）==1）
{
（*（c+（i-j+s）*n+j））=0.16;
}
elseif（abs（i-j）==2）
{
（*（c+（i-j+s）*n+j））=-0.064;
}
}
}
//把ｃ中其他元素赋值０
for（i=0;i<（1+r+s）;i++）
{
for（j=0;j{
if（（（j-s+i）<0）||（（j-s+i）>（n-1）））
（*（c+i*n+j））=0;
}
}
//把c中所有的元素都赋给d
for（i=0;i<（1+r+s）;i++）
{
for（j=0;j{
（*（d+i*n+j））=（*（c+i*n+j））;
}
}
}
intTRI_DLLUdecompostion（）//进行doolittle的LU分解
{
intt0=0;
intj0=0;
inti0=0;
intk=0;
intt=0;
intj=0;
inti=0;
doublesum=0;
doubletemp=0;
for（k=0;k{
j0=（k+s）;
if（j0>（n-1））j0=（n-1）;
//赋值循环一
for（j=k;j<=j0;j++）
{
//使t0取得1,（k-r）,（j-s）三个数中的最大值
t0=0;
if（t0<（k-r））t0=（k-r）;
if（t0<（j-s））t0=（j-s）;
for（t=t0;t<（k）;t++）
{
sum+=（*（d+（k-t+s）*n+t））*（*（d+（t-j+s）*n+j））;//d[][]中元素的数组下标表示有可能要减一
}
（*（d+（k-j+s）*n+j））-=sum;
d["<sum=0;
}
//赋值循环二
//使i0获得k+r和n两者当中的最小值
i0=（k+r）;
if（i0>（n-1））i0=（n-1）;
for（i=（k+1）;i<=i0;i++）
{
//使t0获得1,i-r,k-s三者间的最大值
t0=0;
if（t0<（i-r））t0=（i-r）;
if（t0<（k-s））t0=（k-s）;
for（t=t0;t<（k）;t++）
{
sum+=（*（d+（i-t+s）*n+t））*（*（d+（t-k+s）*n+k））;//d[][]中元素的数组下标表示减一
}
*（d+（i-k+s）*n+k）=（（*（d+（i-k+s）*n+k））-sum）/（*（d+s*n+k））;
sum=0;
}
}
return0;
}
intTRI_ROOT（double*xb,double*B）//进行对角线形线性方程组求根运算L*y=B;U*x=y;
{
inti=0;
intt0=0;
intt=0;
doublesum=0;
for（i=1;i{
//使t0取得1和i-r两者的最大值
t0=0;
if（t0<（i-r））t0=（i-r）;
for（t=t0;t<（i）;t++）
{
cout<<"a["<<（i-1）<<"]"<<"["<sum+=（*（d+（i-t+s）*n+t））*（*（B+t））;//c[i-t+s+1][t]没有写成（*（c+（i-t+s+1）*n+t））的形式
}
（*（B+t））-=sum;
sum=0;
}
（*（xb+n-1））=（*（B+n-1））/（*（d+（s）*n+n-1））;//c[s+1][t]没有写成（*（c+（s+1）*n+t））的形式
for（i=（n-2）;i>=0;i--）
{
t0=（i+s）;
if（t0>n）t0=（n-1）;
for（t=（i+1）;t<=t0;t++）
{
sum+=（*（d+（i-t+s）*n+t））*（*（xb+t））;
}
（*（xb+i））=（（*（B+i））-sum）/（*（d+（s）*n+i））;
sum=0;
}
return0;
}
voidTRI_display（）
{
inti,j;
//cout<<"LU变换后的A矩阵为"<for（i=0;i{
for（j=0;j{
if（（i-j+s）>=0）
{
cout<<（*（d+（abs（i-j+s））*n+j））<<"";
}
}
cout<}
return;
}
voidTRI_ROOT_display（double*xb）
{
inti;
for（i=0;i{
cout<<"["<}
cout<return;
}
2.2“matrix-tool.cpp”中源代码
#include"matrix-tool.h"
#include
#include
usingnamespacestd;
voidMatrix_M_Vector（double*A,double*b,double*y,intn）
{
//此时需要使用全局变量n也需要给函数三个指针型的形参赋值
inti,j;
doublesum=0;
for（i=0;i{
for（j=0;j{
sum+=（*（A+i*n+j））*（*（b+j））;
}
（*（y+i））=sum;
sum=0;
}
}
voidMatrix_Converted_M_Vector（double*A,double*b,double*y,intn）
{
inti,j;
doublesum=0;
for（i=0;i{
for（j=0;j{
sum+=（*（A+j*n+i））*（*（b+j））;
}
（*（y+i））=sum;
sum=0;
}
}
voidColumn_M_Row（double*a,double*b,double*A,intn）//列乘以行得到一个矩阵A
{
inti,j;
for（i=0;i{
for（j=0;j{
（*（A+i*n+j））=（*（a+i））*（*（b+j））;
}
}
}
voidRow_M_Column（double*a,double*b,double*y,intn）//行乘以列得到一个数值y
{
inti;
doublesum=0;
for（i=0;i{
sum+=（*（a+i））*（*（b+i））;
}
cout<<"Row_M_Columnsum"<y=∑
sum=0;
}
voidPrint_Matrix（double*a,intn）
{
inti,j;
for（i=0;i{
for（j=0;j{
cout<<（*（a+i*n+j））<<"";
}
cout<}
}
voidPrint_Vector（double*b,intn）
{
inti;
for（i=0;icout<<（*（b+i））<<""<}
2.3”convert_mufa.cpp”中源代码
#include"math.h"
#include"doolittle.h"
#include"daizhuang.h"
#include
usingnamespacestd;
//变量初始化
intflag2=0;
doubleyb[501]={0};
doubleub[501]={0};
doublemyb[501]={0};
doubleerror_b[2]={1,100};
externintn;
externdouble*d;
voidConvert_Iteration（double*error）
{
intk;
doubleeta=0;
doublebeta=0;
doublesum3=0;
doublesum4=0;
for（k=0;k{
sum3+=（ub[k]*ub[k]）;
}
eta=sqrt（sum3）;
sum3=0;
for（k=0;k{
yb[k]=ub[k]/eta;
myb[k]=yb[k];
}
//线性方程组求解子程序段
TRI_ROOT（ub,yb）;
for（k=0;k{
sum4=（*（myb+k））*（*（ub+k））;
beta+=sum4;
}
（*error）=（*（error+1））;
（*（error+1））=（1/beta）;
beta=0;
sum4=0;
flag2++;
}
doubleConvert_root（）
{
doublesubtract=10;
while（subtract/fabs（error_b[1]）>=（1E-12））
{
Convert_Iteration（error_b）;
subtract=fabs（error_b[0]-error_b[1]）;
}
returnerror_b[1];
}
2.4”doolittle.cpp”中源代码
#include
#include"doolittle.h"
usingnamespacestd;
intinput（intn,double*a,double*b）//原始系统矩阵a（n*n）的数据、维数的输入，方程右边向量b（n）的输入
{
inti,j;
cout<<"PleaseentervalueforeachelementfortheMatrix"<for（i=0;i{
for（j=0;j{
cout<<"thea["<cin>>*（a+i*n+j）;
}
}
for（i=0;i{
cout<<"theb["<cin>>（*（b+i））;
}
return0;
}
intDLLUdecompostion（double*a,intn）//进行doolittle的LU分解
{
intk,i,j,t;
doublesum1=0;
doublesum2=0;
for（k=0;kfor（j=k;j{
for（t=0;t<（k）;t++）
{
sum1+=（*（a+k*n+t））*（*（a+t*n+j））;
}
*（a+k*n+j）=（（*（a+k*n+j））-sum1）;//求的是U[k][j]
sum1=0;
};
for（i=（k+1）;i{
for（t=0;t<（k）;t++）
{
sum2+=（*（a+i*n+t））*（*（a+t*n+k））;
}
*（a+i*n+k）=（（*（a+i*n+k））-sum2）/（*（a+k*n+k））;//求的是L[i][k]
sum2=0;
}
}
return0;
}
intROOT（intn,double*a,double*b）//进行求根运算
{
inti,t;
doublesum3=0;
doublesum4=0;
//求取Y[n]
for（i=1;i{
for（t=0;t<（i）;t++）
{
sum3+=（*（a+i*n+t））*（*（b+t））;
}
（*（b+i））-=sum3;
sum3=0;
}
//求取X[n]
（*（b+n-1））/=（*（a+n*n-1））;
for（i=（n-2）;i>=0;i--）
{
for（t=（i+1）;t{
sum4+=（*（a+i*n+t））*（*（b+t））;
}
（*（b+i））=（*（b+i）-sum4）/（*（a+n*i+i））;//实为求取X[I];
sum4=0;
}
return0;
}
voiddisplay（intn,double*a,double*b）
{
inti,j;
for（i=0;i{
for（j=0;j{
cout<<（*（a+i*n+j））<<"";
}
cout<}
for（i=0;i{
cout<<"b["<}
cout<}
2.5”mufa.cpp”中源代码
#include"mufa.h"
#include"math.h"
#include
usingnamespacestd;
//幂法中相关变量的定义
doublesum1=0;
doublesum2=0;
doublemed=0;
doublespec=0;
doubleya[501]={0};
doubleua[501];
//externdoubleerror[2];
externdouble*c;
externchars;
externcharr;
externintn;
intflag1=0;
voiditeration（double*error）
{
inti,j;
//实现U[K-1]的标准化
for（i=0;i{
su
升级会员 