雅可比迭代法与矩阵的特征值.docx
《雅可比迭代法与矩阵的特征值.docx》由会员分享,可在线阅读,更多相关《雅可比迭代法与矩阵的特征值.docx(33页珍藏版)》请在冰点文库上搜索。
雅可比迭代法与矩阵的特征值
实验五
矩阵的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;j
cout<
cout<}
}
voiddisp(double*q,intn){
cout<<"====================================="<for(inti=0;icout<<"X["<
cout<<"====================================="<}
voidinput(double**p,introw,intcol){
for(inti=0;icout<<"输入第"<
";
for(intj=0;j
cin>>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;i
temp=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;k
p[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;j
cin>>p[i][j];
}
}
voiddispMatrix(double**p,introw,intcol){
for(inti=0;ifor(intj=0;j
cout<
cout<}
}
voiddispVector(double*q,intn){
cout<<"====================================="<for(inti=0;icout<<"X["<
cout<<"====================================="<}
voidinitMatrix(double**p,introw,intcol){
for(inti=0;ifor(intj=0;j
p[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;i
U[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;j
if(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;i
sum+=(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;j
cin>>p[i][j];
}
}
voiddispMatrix(double**p,introw,intcol){
for(inti=0;ifor(intj=0;j
cout<
cout<}
}
voiddispVector(double*q,intn){
for(inti=0;icout<<"X["<
cout<}
doublesumRow(double**A,double*x1,introw,intcol){
doublesum=0;
for(inti=0;i
if(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)
{
floata[10][10
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|