雅可比迭代法与矩阵的特征值DOC.docx
《雅可比迭代法与矩阵的特征值DOC.docx》由会员分享，可在线阅读，更多相关《雅可比迭代法与矩阵的特征值DOC.docx（28页珍藏版）》请在冰豆网上搜索。
雅可比迭代法与矩阵的特征值DOC
实验五
矩阵的lu分解法，雅可比迭代法
班级：
学号：
姓名：
实验五矩阵的LU分解法，雅可比迭代
一、目的与要求：
Ø熟悉求解线性方程组的有关理论和方法；
Ø会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序；
Ø通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
二、实验内容：
Ø会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。
三、程序与实例
Ø列主元高斯消去法
算法：
将方程用增广矩阵[A∣b]=（
表示
1）消元过程
对k=1,2,…,n-1
①选主元，找
使得
=
②如果
，则矩阵A奇异，程序结束；否则执行③。
③如果
，则交换第k行与第
行对应元素位置，
j=k,┅,n+1
④消元，对i=k+1,┅,n计算
对j=l+1,┅,n+1计算
2）回代过程
①若
，则矩阵A奇异，程序结束；否则执行②。
②
；对i=n-1,┅,2,1,计算
程序与实例
程序设计如下：
#include
#include
usingnamespacestd;
voiddisp（double**p,introw,intcol）{
for（inti=0;ifor（intj=0;jcout<cout<}
}
voiddisp（double*q,intn）{
cout<<"====================================="<for（inti=0;icout<<"X["<cout<<"====================================="<}
voidinput（double**p,introw,intcol）{
for（inti=0;icout<<"输入第"<";
for（intj=0;jcin>>p[i][j];
}
}
intfindMax（double**p,intstart,intend）{
intmax=start;
for（inti=start;iif（abs（p[i][start]）>abs（p[max][start]））
max=i;
}
returnmax;
}
voidswapRow（double**p,intone,intother,intcol）{
doubletemp=0;
for（inti=0;itemp=p[one][i];
p[one][i]=p[other][i];
p[other][i]=temp;
}
}
booldispel（double**p,introw,intcol）{
for（inti=0;iintflag=findMax（p,i,row）;//找列主元行号
if（p[flag][i]==0）returnfalse;
swapRow（p,i,flag,col）;//交换行
for（intj=i+1;jdoubleelem=p[j][i]/p[i][i];//消元因子
p[j][i]=0;
for（intk=i+1;kp[j][k]-=（elem*p[i][k]）;
}
}
}
returntrue;
}
doublesumRow（double**p,double*q,introw,intcol）{
doublesum=0;
for（inti=0;iif（i==row）
continue;
sum+=（q[i]*p[row][i]）;
}
returnsum;
}
voidback（double**p,introw,intcol,double*q）{
for（inti=row-1;i>=0;i--）{
q[i]=（p[i][col-1]-sumRow（p,q,i,col））/p[i][i];
}
}
intmain（）
{
cout<<"Inputn:
";
intn;//方阵的大小
cin>>n;
double**p=newdouble*[n];
for（inti=0;ip[i]=newdouble[n+1];
}
input（p,n,n+1）;
if（!
dispel（p,n,n+1））{
cout<<"奇异"<return0;
}
double*q=newdouble[n];
for（inti=0;iq[i]=0;
back（p,n,n+1,q）;
disp（q,n）;
delete[]q;
for（inti=0;idelete[]p[i];
delete[]p;
}
1.用列主元消去法解方程
例2解方程组
计算结果如下
B=-1.461954
C=1.458125
D=-6.004824
E=-2.209018
F=14.719421
Ø矩阵直接三角分解法
算法：
将方程组Ax=b中的A分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组Ax=b化为解2个方程组Ly=b，Ux=y，具体算法如下：
①对j=1,2,3,…,n计算
对i=2,3,…,n计算
②对k=1,2,3,…,n:
a.对j=k,k+1,…,n计算
b.对i=k+1,k+2,…,n计算
③
，对k=2,3,…,n计算
④
对k=n-1,n-2,…,2,1计算
注：
由于计算u的公式于计算y的公式形式上一样，故可直接对增广矩阵
[A∣b]=
施行算法②，③，此时U的第n+1列元素即为y。
程序与实例
例3求解方程组Ax=b
A=
b=
结果为
X[0]=3.000001
X[1]=-2.000001
X[2]=1.000000
X[3]=5.000000
#include
usingnamespacestd;
double**newMatrix（introw,intcol）{
double**p=newdouble*[row];//行
for（inti=0;ip[i]=newdouble[col];
returnp;
}
voiddelMatrix（double**p,introw）{
for（inti=0;idelete[]p[i];
delete[]p;
}
voidinputMatrix（double**p,introw,intcol）{
for（inti=0;icout<<"输入第"<";
for（intj=0;jcin>>p[i][j];
}
}
voiddispMatrix（double**p,introw,intcol）{
for（inti=0;ifor（intj=0;jcout<cout<}
}
voiddispVector（double*q,intn）{
cout<<"====================================="<for（inti=0;icout<<"X["<cout<<"====================================="<}
voidinitMatrix（double**p,introw,intcol）{
for（inti=0;ifor（intj=0;jp[i][j]=0;
}
doublesum（double**L,double**U,introw,intcol）{
intmin=（row>=col?
col:
row）;
doublesum=0;
for（inti=0;isum+=（U[i][col]*L[row][i]）;
returnsum;
}
voidresolve（double**A,double**L,double**U,introw,intcol）{
for（inti=0;iU[0][i]=A[0][i];//把A的第一行给U的第一行
L[0][0]=A[0][0];
for（inti=1;iL[i][0]=A[i][0]/A[0][0];
for（inti=1;ifor（intj=1;jif（i<=j）
U[i][j]=A[i][j]-sum（L,U,i,j）;
else
L[i][j]=（A[i][j]-sum（L,U,i,j））/U[j][j];
}
for（inti=0;iL[i][i]=1;
}
doublesumRowY（double**L,double*y,introw）{
doublesum=0;
for（inti=0;isum+=（L[row][i]*y[i]）;
returnsum;
}
voidbackY（double**L,double*b,double*y,introw）{
for（inti=0;iy[i]=0;//初始化y向量全为零
for（inti=0;iy[i]=b[i]-sumRowY（L,y,i）;
}
doublesumRowX（double**U,double*x,introw,intcol）{
doublesum=0;
for（inti=row+1;isum+=（U[row][i]*x[i]）;
returnsum;
}
voidbackX（double**U,double*y,double*x,introw）{
for（inti=0;ix[i]=0;//初始化x向量全为零
for（inti=row-1;i>=0;i--）
x[i]=（y[i]-sumRowX（U,x,i,row））/U[i][i];
}
intmain（）
{
cout<<"输入方阵大小n:
";
intn=0;
cin>>n;
cout<<"====================================="<cout<<"开始录入方阵数据....."<double**A=newMatrix（n,n）;//开辟矩阵A
inputMatrix（A,n,n）;//录入数据到A中
cout<<"录入方阵数据完毕......"<cout<<"====================================="<cout<<"开始录入列阵....."<cout<<"输入列阵:
";
double*b=newdouble[n];//开辟向量b
for（inti=0;icin>>b[i];//录入数据到b中
cout<<"录入列阵数据完毕....."<double**L=newMatrix（n,n）;//开辟方阵L
initMatrix（L,n,n）;//初始化L全为零
double**U=newMatrix（n,n）;//开辟方阵U
initMatrix（U,n,n）;//初始化U
resolve（A,L,U,n,n）;//分解A为L与U
double*y=newdouble[n];//开辟向量y
backY（L,b,y,n）;//回代求出y
double*x=newdouble[n];//开辟向量X
backX（U,y,x,n）;//回代求出X
dispVector（x,n）;
//释放空间
delMatrix（A,n）;
delMatrix（L,n）;
delMatrix（U,n）;
delete[]b;
delete[]y;
delete[]x;
}
Ø迭代法
雅可比迭代法
算法：
设方程组Ax=b系数矩阵的对角线元素
M为迭代次数容许的最大值,ε为容许误差。
①取初始向量x=
，令k=0。
②对i=1,2,…,n计算
③如果
，则输出
，结束；否则执行④。
④如果k≥M,则不收敛，终止程序；否则
转②。
程序与实例
例4用雅可比迭代法解方程组
结果为
迭代次数为20
X[0]=1.000000
X[1]=2.000000
X[2]=-1.000000
#include
#include
usingnamespacestd;
double**newMatrix（introw,intcol）{
double**p=newdouble*[row];//行
for（inti=0;ip[i]=newdouble[col];
returnp;
}
voiddelMatrix（double**p,introw）{
for（inti=0;idelete[]p[i];
delete[]p;
}
voidinputMatrix（double**p,introw,intcol）{
for（inti=0;icout<<"输入第"<";
for（intj=0;jcin>>p[i][j];
}
}
voiddispMatrix（double**p,introw,intcol）{
for（inti=0;ifor（intj=0;jcout<cout<}
}
voiddispVector（double*q,intn）{
for（inti=0;icout<<"X["<cout<}
doublesumRow（double**A,double*x1,introw,intcol）{
doublesum=0;
for（inti=0;iif（row==i）
continue;
sum+=（A[row][i]*x1[i]）;
}
returnsum;
}
voiditerat（double**A,double*b,double*x1,double*x2,introw）{
for（inti=0;ix2[i]=（b[i]-sumRow（A,x1,i,row））/A[i][i];
}
boolblow_error（double*x1,double*x2,introw,doublee）{
doublesum=0;
for（inti=0;isum+=abs（x2[i]-x1[i]）;
}
if（sumelsereturnfalse;
}
intmain（）
{
cout<<"输入方阵大小n:
";
intn=0;
cin>>n;
cout<<"====================================="<cout<<"开始录入方阵数据....."<double**A=newMatrix（n,n）;//开辟矩阵A
inputMatrix（A,n,n）;//录入数据到A中
cout<<"录入方阵数据完毕......"<cout<<"====================================="<cout<<"开始录入列阵....."<cout<<"输入列阵:
";
double*b=newdouble[n];//开辟向量b
for（inti=0;icin>>b[i];//录入数据到b中
cout<<"录入列阵数据完毕....."<cout<<"====================================="<cout<<"输入迭代次数容许的最大值:
";
intM=0;
cin>>M;
cout<<"输入容许误差:
";
doublee=0;
cin>>e;
cout<<"====================================="<double*x1=newdouble[n];//开辟x1向量
double*x2=newdouble[n];//开辟x2向量
cout<<"输入初始向量:
";
for（inti=0;icin>>x1[i];
cout<<"====================================="<intk=0;//迭代计数器
for（;;）{
iterat（A,b,x1,x2,n）;//迭代一次
if（blow_error（x1,x2,n,e））{
dispVector（x2,n）;
cout<<"迭代次数为:
"<return0;
}else{
if（k>=M）{
cout<<"不收敛"<return0;
}else{
k++;
for（inti=0;ix1[i]=x2[i];
}
}
}
//释放空间
delMatrix（A,n）;
delete[]b;
delete[]x1;
delete[]x2;
}
Ø高斯-塞尔德迭代法
算法:
设方程组Ax=b的系数矩阵的对角线元素,
M为迭代次数容许的最大值,ε为容许误差
①取初始向量
令k=0。
②对i=1,2,…,n计算
③如果
，则输出
结束；否则执行④。
④如果
则不收敛，终止程序；否则
转②。
程序与实例
例5用高斯-塞尔德迭代法解方程组
结果为
X[0]=3.000000
X[1]=2.000000
X[2]=1.000000
#include
#include
#include
#include
#defineN100
main（）
{
inti;
float*x;
floatc[12]={8.0,-3.0,2.0,20.0,
4.0,11.0,-1.0,33.0,
6.0,3.0,12.0,36.0};
float*GauseSeidel（float*,int）;
x=GauseSeidel（c,3）;
for（i=0;i<=2;i++）printf（"x[%d]=%f\n",i,x[i]）;
getch（）;
}
float*GauseSeidel（float*a,intn）
{
inti,j,nu=0;
float*x,dx,d;
x=（float*）malloc（n*sizeof（float））;
for（i=0;i<=n-1;i++）x[i]=0.0;
do
{
for（i=0;i<=n-1;i++）
{
d=0.0;
for（j=0;j<=n-1;j++）
d+=*（a+i*（n+1）+j）*x[j];
dx=（*（a+i*（n+1）+n）-d）/（*（a+i*（n+1）+i））;
x[i]+=dx;
}
if（nu>=N）
{
printf（"foldnumber\n"）;
}
nu++;
}
while（fabs（dx）>1e-6）;
returnx;
}
例6用雅可比迭代法解方程组
迭代4次得解
若用高斯-塞尔德迭代法则发散。
#include
voidmain（void）
{
intk,n;
doublex[3]={7,2,5};
for（k=0;k<5;k++）
{
doublea,b;
a=x[0];
b=x[1];
x[0]=（7-2.0*x[1]+2*x[2]）/1;
x[1]=（2-a-x[2]）/1;
x[2]=（5-2*a-2*b）/1;
}
for（n=0;n<3;n++）
{
printf（"x[%d]=%8.6f\n",n,x[n]）;
}
}
用高斯-塞尔德迭代法解方程组
迭代84次得解
，若用雅克比迭代法则发散。
#include
#include
voidLOOP（floata[10][10],floatb[10],floatx[10],int）;
voidmain（void）
升级会员 