实验五矩阵的LU分解法雅可比迭代.docx
《实验五矩阵的LU分解法雅可比迭代.docx》由会员分享，可在线阅读，更多相关《实验五矩阵的LU分解法雅可比迭代.docx（16页珍藏版）》请在冰豆网上搜索。
实验五矩阵的LU分解法雅可比迭代
实验五矩阵的LU分解法，雅可比迭代
实
验
报
告
学院：
计算机科学与软件学院
班级：
116班
姓名：
薛捷星
学号：
112547
一、目的与要求：
Ø熟悉求解线性方程组的有关理论和方法；
Ø会编制列主元消去法、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,计算
程序与实例
例1解方程组
输出结果如下：
X[0]=-0.398234
X[1]=0.013795
X[2]=0.335144
程序如下：
#include
#include
main（）
{
inti,j,p,o,l,q;
doublea[3][4]={{0.101,2.304,3.555,1.183},{-1.347,3.712,4.623,2.137},{-2.835,1.072,5.643,3.035}};
doublex[3],z[4];
printf（"列主元消去法\n"）;
for（j=0;j<2;j++）
{
for（i=j+1;i<3;i++）
{
if（fabs（a[j][j]）{
for（p=0;p<4;p++）
{
z[p]=a[j][p];
a[j][p]=a[i][p];
a[i][p]=z[p];
}/*交换得最大主元*/
}
}
for（l=j+1;l<3;l++）
{
for（q=3;q>=j;q--）
{
a[l][q]=（a[l][q]-（a[l][j]/a[j][j]）*a[j][q]）;
}
}
printf（"进行消去：
\n"）;
for（o=0;o<3;o++）
{
for（p=0;p<4;p++）
{
printf（"%12.6f",a[o][p]）;
}
printf（"\n"）;
}
}
x[2]=a[2][3]/a[2][2];
x[1]=（a[1][3]-x[2]*a[1][2]）/a[1][1];
x[0]=（a[0][3]-x[2]*a[0][2]-x[1]*a[0][1]）/a[0][0];
printf（"最后的解:
\n"）;
for（i=0;i<3;i++）
{
printf（"x[%d]=%12.6f\n",i,x[i]）;
}
}
结果如下：
例2解方程组
计算结果如下
B=-1.161954
C=1.458125
D=-6.004824
E=-2.209018
F=14.719421
程序如下：
#include
#include
voidmain（void）
{
inti,j,p,o,l,q;
doublea[5][6]={{8.77,2.40,5.66,1.55,1.0,-32.04},{4.93,1.21,4.48,1.10,1.0,-20.07},{3.53,1.46,2.92,1.21,1.0,-8.53},{5.05,4.04,2.51,2.01,1.0,-6.30},{3.54,1.04,3.47,1.02,1.0,-12.04}};
doublex[5],z[6];
printf（"列主元消去法求五元一次方程组：
\n"）;
for（j=0;j<4;j++）
{
for（i=j+1;i<5;i++）
{
if（fabs（a[j][j]）{
for（p=0;p<6;p++）
{
z[p]=a[j][p];
a[j][p]=a[i][p];
a[i][p]=z[p];
}/*交换得最大主元*/
}
}
for（l=j+1;l<5;l++）
{
for（q=5;q>=j;q--）
{
a[l][q]=（a[l][q]-（a[l][j]/a[j][j]）*a[j][q]）;
}
}
printf（"消去一列：
\n"）;
for（o=0;o<5;o++）
{
for（p=0;p<6;p++）
{
printf（"%12.6f",a[o][p]）;
}
printf（"\n"）;
}
}
x[4]=a[4][5]/a[4][4];
x[3]=（a[3][5]-x[4]*a[3][4]）/a[3][3];
x[2]=（a[2][5]-x[4]*a[2][4]-x[3]*a[2][3]）/a[2][2];
x[1]=（a[1][5]-x[4]*a[1][4]-x[3]*a[1][3]-x[2]*a[1][2]）/a[1][1];
x[0]=（a[0][5]-x[4]*a[0][4]-x[3]*a[0][3]-x[2]*a[0][2]-x[1]*a[0][1]）/a[0][0];
printf（"方程组的解为：
\n"）;
for（i=0;i<5;i++）
{
printf（"x[%d]=%12.6f\n",i,x[i]）;
}
}
Ø矩阵直接三角分解法
算法：
将方程组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
voidmain（void）
{
inti,j;
doublea[4][4]={{1,2,-12,8},{5,4,7,-2},{-3,7,9,5},{6,-12,-8,3}};
doublel[4][4],b[4]={27,4,11,49},y[4],x[4];
printf（"直接三角分解法求方程组的解：
\n"）;
for（i=0;i<4;i++）
{
l[i][0]=a[i][0];
l[0][i]=a[0][i];
}
l[1][1]=a[1][1]-l[1][0]*a[0][1];l[1][2]=a[1][2]-l[1][0]*a[0][2];l[1][3]=a[1][3]-l[1][0]*a[0][3];
l[2][1]=（a[2][1]-l[2][0]*l[0][1]）/l[1][1];l[2][2]=a[2][2]-l[2][0]*l[0][2]-l[2][1]*l[1][2];l[2][3]=a[2][3]-l[2][0]*l[0][3]-l[2][1]*l[1][3];
l[3][1]=（a[3][1]-l[3][0]*l[0][1]）/l[1][1];l[3][2]=（a[3][2]-l[3][0]*l[0][2]-l[3][1]*l[1][2]）/l[2][2];l[3][3]=a[3][3]-l[3][0]*l[0][3]-l[3][1]*l[1][3]-l[3][2]*l[2][3];
printf（"LU合并矩阵：
\n"）;
for（i=0;i<4;i++）
{
for（j=0;j<4;j++）
{
printf（"%12.6f",l[i][j]）;
}
printf（"\n"）;
}
y[0]=b[0];
y[1]=b[1]-y[0]*l[1][0];
y[2]=b[2]-y[0]*l[2][0]-y[1]*l[2][1];
y[3]=b[3]-y[0]*l[3][0]-y[1]*l[3][1]-y[2]*l[3][2];
printf（"Y矩阵：
\n"）;
for（i=0;i<4;i++）
printf（"y[%d]=%12.6f\n",i,y[i]）;
x[3]=y[3]/l[3][3];
x[2]=（y[2]-x[3]*l[2][3]）/l[2][2];
x[1]=（y[1]-x[3]*l[1][3]-x[2]*l[1][2]）/l[1][1];
x[0]=（y[0]-x[3]*l[0][3]-x[2]*l[0][2]-x[1]*l[0][1]）/l[0][0];
printf（"方程组的解：
\n"）;
for（i=0;i<4;i++）
printf（"x[%d]=%12.6f\n",i,x[i]）;
}
结果如下：
Ø迭代法
雅可比迭代法
算法：
设方程组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
#definee0.000001
voidmain（void）
{
floata,b,c,x[3];
inti;
printf（"Jacobi迭代法求方程组：
\n"）;
printf（"输入X1,X2,X3的初始值,以“，”间隔：
\n"）;
scanf（"%f,%f,%f",&a,&b,&c）;
for（i=0;;i++）
{
x[0]=（8-2*b-c）/5;
x[1]=（21-2*a+3*c）/8;
x[2]=（1-a+3*b）/（-6）;
if（fabs（x[0]-a）break;
else
a=x[0];b=x[1];c=x[2];
if（i>200）
{
printf（"发散\n"）;
break;
}
}
printf（"迭代%d次\n",i+3）;
printf（"方程组的解为：
\n"）;
for（i=0;i<3;i++）
{
printf（"x[%d]=%f\n",i,x[i]）;
}
}
结果如下：
Ø高斯-塞尔德迭代法
算法:
设方程组Ax=b的系数矩阵的对角线元素,
M为迭代次数容许的最大值,ε为容许误差
①取初始向量
令k=0。
②对i=1,2,…,n计算
③如果
，则输出
结束；否则执行④。
④如果
则不收敛，终止程序；否则
转②。
程序与实例
例5用高斯-塞尔德迭代法解方程组
结果为
X[0]=3.000000
X[1]=2.000000
X[2]=1.000000
程序如下：
#include
#include
#definee0.000001
voidmain（void）
{
inti;
floata[3][4]={{8,-3,2,20},{4,11,-1,33},{6,3,12,36}},x[3];
floatt,b,c;
printf（"高斯--赛德尔法解方程组\n"）;
printf（"输入X1，X2，X3的初始值，以逗号间隔：
\n"）;
scanf（"%f,%f,%f",&x[0],&x[1],&x[2]）;
for（i=0;;i++）
{
t=x[0];b=x[1];c=x[2];
x[0]=（a[0][3]-a[0][1]*x[1]-a[0][2]*x[2]）/a[0][0];
x[1]=（a[1][3]-a[1][0]*x[0]-a[1][2]*x[2]）/a[1][1];
x[2]=（a[2][3]-a[2][0]*x[0]-a[2][1]*x[1]）/a[2][2];
if（fabs（x[0]-t）break;
else
continue;
if（i>200）
{
printf（"发散\n"）;
break;
}
}
printf（"迭代%d次\n",i+1）;
for（i=0;i<3;i++）
printf（"x[%d]=%12.6f\n",i,x[i]）;
}
结果如下：
例6用雅可比迭代法解方程组
迭代4次得解
若用高斯-塞尔德迭代法则发散。
结果如下：
用高斯-塞尔德迭代法解方程组
迭代84次得解
，若用雅克比迭代法则发散。
结果如下：
_
升级会员 