线性方程组AXB的数值计算方法实验.docx
《线性方程组AXB的数值计算方法实验.docx》由会员分享，可在线阅读，更多相关《线性方程组AXB的数值计算方法实验.docx（40页珍藏版）》请在冰豆网上搜索。
线性方程组AXB的数值计算方法实验
线性方程组AX=B的数值计算方法实验
【摘要】在自然科学与工程技术中很多问题的解决常常归结为解线性代数方程组。
例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组的问题，用差分法或者有限元法解常微分方程，偏微分方程边值问题等都导致求解线性方程组。
线性代数方面的计算方法就是研究求解线性方程组的一些数值解法与研究计算矩阵的特征值及特征向量的数值方法。
关于线性方程组的数值解法一般有两类：
直接法和迭代法。
关键字高斯消元法、三角分解法、高斯-赛德尔迭代、稀疏矩阵
1、实验目的
1.掌握高斯消元法、三角分解法、高斯—赛德尔迭代发的编程技巧。
2.掌握线性方程组AX=B的数值计算方法。
3.掌握矩阵的基本编程技巧。
2、实验原理
1.高斯消元法
数学上，高斯消元法是线性代数规划中的一个算法，可用来为线性方程组求解。
高斯（Gauss）夏鸥按法其实是将一般的线性方程组变换为三角形（上三角）方程组求解问题（消元法），只是步骤规范，便于编写计算机程序。
一般高斯消元法包括两过程：
先把方程组化为同解的上三角形方程组，再按相反顺序求解上三角方程组。
前者称为消去或消元过程，后者称回代过程。
消去过程实际上是对增广矩阵作行初等变换。
对一般的n阶方程组，消去过程分n-1步：
第一步消去
下方元素。
第二步消去
下方元素，......，第n-1步消去
下方元素。
即第k步将第k行的适当倍数加于其后各行，或可说是从k+1~n行减去第k行的适当倍数，使它们第k列元素变为零，而其余列元素减去第k行对应列元素的倍数。
2.三角分解法
三角分解法是将原正方（square）矩阵分解成一个上三角形矩阵或是排列（permuted）的上三角形矩阵和一个下三角形矩阵，这样的分解法又称为LU分解法。
它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵，和求解联立方程组。
不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。
3.高斯—赛德尔迭代
高斯－赛德尔迭代（Gauss–Seidelmethod）是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。
研究雅可比迭代法，我们发现在逐个求
的分量时，当计算到
时,分量
，......，
都已经求得，而仍用旧分量
，......，
计算
。
由于新计算出的分量比旧分量准确些，因此设想一旦新分量
，......，
求出，马上就用新分量
，......，
代替雅可比迭代法中
，......，
来求
这就是高斯-赛德尔（Gauss-Seidel）迭代法。
把矩阵A分解成
（6）其中
分别为
的主对角元除外的下三角和上三角部分,于是,方程组
（1）便可以写成
即
其中
（7）
以
为迭代矩阵构成的迭代法（公式）
（8）
称为高斯—塞德尔迭代法（公式）,用变量表示的形式为
（9）
4.稀疏矩阵
矩阵中非零元素的个数远远小于矩阵元素的总数，并且非零元素的分布没有规律，则称该矩阵为稀疏矩阵（sparsematrix）；与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对角矩阵），则称该矩阵为特殊矩阵。
常见于进行大量数据计算。
3、实验内容
1.P1081
许多科学应用包含的矩阵带有很多的零。
在实际情况中很重要的三角形线性方程组有如下形式：
d1x1+c1x2=b1
a1x1+d2x2+c2x3=b2
a2x2+d3x3+c3x4=b3
····
····
····
aN-2xN-2+dN-1xN-1+cN-1xN=bN-1
aN-1xN-1+dNxN=bN
构造一个程序求解三角形线性方程组。
可假定不需要行变换，而且可用第k行消去第k+1行的xk。
2.P1201
求解线性方程组AX=B，其中
A=
B=
使用三角分解法求解X。
3.P1202
求解线性方程组AX=B，其中A=[aij]N×N，,aij=ij-1；而且B=[bij]N×1，b11=N，当i≥时，bij=（iN-1）/（i-1）。
对N=3，7，11的情况分别求解。
精确解为X=[11…11]’。
对得到的结果与精确解的差异进行解释。
4.P1203
通过重复求解N各线性方程组
ACJ=EJ，其中J=1，2，…，N
来得到A-1,
则
A[C1C2…CN]=[E1E2…EN]
而且
A-1=[C1C2…CN]
保证对LU分解只计算一次!
5.P1293
设有如下三角线性方程组，而且系数矩阵具有严格对角优势：
d1x1+c1x2=b1
a1x1+d2x2+c2x3=b2
a2x2+d3x3+c3x4=b3
····
····
····
aN-2xN-2+dN-1xN-1+cN-1xN=bN-1
aN-1xN-1+dNxN=bN
（i）根据方程组
（1），式
（2）和式（3），设计一个算法来求解上述方程组。
算法必须有效地利用系数矩阵的稀疏性。
a11x1+a12x2+···a1jxj+···+a1NxN=b1
a21x1+a22x2+···a2jxj+···+a2NxN=b2
·····
·····
·························方程组
（1）
aj1x1+aj2x2+···ajjxj+···+ajNxN=bj
·····
·····
·····
aN1x1+aN2x2+···aNjxj+···+aNNxN=bN
=
，j=1，2，···，N···式
（2）
=
，j=1，2，···，N··式（3）
（ii）根据（i）中设计的算法构建一个C++程序，并求解下列上三角线性方程组
（a）4m1+m2=3（b）4m1+m2=1
m1+4m2+m3=3m1+4m2+m3=2
m2+4m3+m4=3m2+4m3+m4=1
m3+4m4+m5=3m3+4m4+m5=2
········
········
········
m48+4m49+m50=3m48+4m49+m50=1
m49+m50=3m49+m50=2
6.P1294
利用高斯—赛德尔迭代法求解下列带状方程。
12x1-2x2+x3=5
-2x1+12x2-2x3+x4=5
x1-2x2+12x3-2x4+x5=5
x2-2x3+12x4-2x5+x6=5
···
x46-2x47+12x48-2x49+x50=5
x47-2x48+12x49-2x50=5
x48-2x49+12x50=5
4、实验结果及分析
1.P1081
实验描述：
本次实验使用系数矩阵的第k行消去第k+1行的xk，消除方法为第k行减去第k-1行乘上系数ak-1/bk-1，待消至第N行时，求解出xN，并依次会带求出各xN-1至x1，为了检验结果的正确性使用上面的方程组组
（1）及方程组
（2）进行验证。
方程组
（1）的结果为
，方程组
（2）的结果为
。
算法流程图：
Y
实验结果：
结果截图
（1）
图1
输入矩阵
=
，输出结果为X=
，与预期结果一致
（2）
图2
输入矩阵
=
，输出结果为X=
，与预期结果一致
实验结论：
通过对系数矩阵的增广矩阵进行高斯消元和回带容易得到线性方程组的解，同时，利用这种方法可以求得矩阵的逆。
2.P1201
实验描述：
本次实验的解法为使用LU矩阵求解X，该解法的内容为将系数矩阵A分解为一上三角矩阵及一下三角矩阵，且有A=LU，之后由LY=B，UX=Y分别求解出Y，X。
实验结果：
图3
输入矩阵A=
，B=
，输出结果为X=
实验结论：
将X代入AX后结果与矩阵B一致，运行结果正确无误，该程序正确，且有X=
3.P1202
实验描述：
（1）本次实验仍使用三角矩阵求解矩阵X的值，求解方法与实验二大致相同；
（2）矩阵A编写函数buildA完成，buildA的输入为矩阵A的阶数N，输出结果为矩阵A的地址，对于矩阵A有A=[aij]N×N，,aij=ij-1；
（3）矩阵B编写函数buildB完成，buildB的输入为矩阵B的阶数N，输出结果为矩阵B的地址，对于矩阵B有B=[bij]N×1，b11=N，当i≥时，bij=（iN-1）/（i-1）。
实验结果：
（1）N=3时
图4
（2）N=7时
图5
（3）N=11时
图6
实验结论：
由图4，图5，图6与结果对比可知，程序运行结果与预期结果相一致，程序正确无误。
4.P1203
实验描述：
（1）本次实验的目的为求解A的逆矩阵A-1，求解方法为利用A[C1C2…CN]=[E1E2…EN]，A-1=[C1C2…CN]求解，故可将之分解为对于ACk=Ek，k=1，2，3，···，N中对于Ck的求解，之后另A-1=[C1C2…CN]；
（2）由于A-1=[C1C2…CN]，A[C1C2…CN]=[E1E2…EN]，固有[E1E2…EN]=I，故Ek=[a1j]，j=1，2，···，N，其中a1k=1，other=0；
（3）对于ACk=Ek，k=1，2，3，···，N的求解方法使用LU三角分解法求解，用此方法求解出各个Ek对应的Ck，最后以此构成A-1；
（4）求出LU后，应判断矩阵对角线上是否存在为0的元素，若存在，则A不存在逆矩阵；若不存在，则可求解逆矩阵A-1；
（5）上述方法中的LU分解只需要进行一次；
（6）对于程序的正确性使用矩阵
及矩阵
进行验证，其中，矩阵
不存在逆矩阵，矩阵
的逆矩阵为
实验结果：
（1）输入为矩阵
图7
（2）输入矩阵为
图8
实验结论：
如果系数矩阵能分解为LU的形式，其中L为下三角矩阵，U为上三角矩阵，通过对系数矩阵的分解再求解可应用简单的迭代进行求解x。
5.P1293
实验描述：
（1）本实验的目的为使用高斯—赛德尔迭代求解带状方程组；
（2）由于实验中的带状方程有极强的稀疏性和相似性，故编程时应考虑该矩阵的以上特点以减少运算量及运行时占用的内存；
（3）为了减少程序的运行次数，故选择式（3）作为运行的程序；
（4）应无明确的结束标志，故选择delta=
<1e-5作为结束迭代的标志，此时再进行一轮迭代后输出结果。
实验结果：
a）4m1+m2=3
m1+4m2+m3=3
m2+4m3+m4=3
m3+4m4+m5=3
····
····
····
m48+4m49+m50=3
m49+m50=3
图9
（b）4m1+m2=1
m1+4m2+m3=2
m2+4m3+m4=1
m3+4m4+m5=2
····
····
····
m48+4m49+m50=1
m49+m50=2
图10
实验结论：
求解带状线性方程组的解可使用高斯-赛德尔迭代法。
6.P1294
实验描述：
（1）本实验的目的为使用高斯—赛德尔迭代求解带状方程组；
（2）由于实验中的带状方程有极强的稀疏性和相似性，故编程时应考虑该矩阵的以上特点以减少运算量及运行时占用的内存；
（3）为了减少程序的运行次数，故选择式（3）作为运行的程序；
（4）应无明确的结束标志，故选择delta=
<1e-5作为结束迭代的标志，此时再进行一轮迭代后输出结果。
算法流程图：
实验结果：
表1.通过高斯赛德尔迭代得到的x的解
1~10
0.46379552381655
0.5372846051999656
0.5090229246013306
0.4982216344361741
0.4989418602397619
0.4999853512481308
0.5000887238901357
0.500015318846052
0.4999947932669753
0.4999978569134675
11~20
0.5000001084251992
0.5000002015766873
0.500000022610945
0.4999999862385722
0.499999995873979
0.500000*********8
0.5000000004390388
0.5000000000239392
0.499999999966016
0.4999999999925573
21~30
0.5000000000017429
0.5000000000009169
0.5
0.4999999999999205
0.5
0.5
0.4999999999999212
0.5
0.5000000000009176
0.5000000000017445
31~40
0.499999999992555
0.499999999966017
0.500000000023939
0.5000000004390391
0.500000*********6
0.499999995873979
0.4999999862385723
0.5000000226109451
0.5000002015766873
0.500000108425199
41~50
0.4999978569134675
0.4999947932669753
0.500015318846052
0.5000887238901357
0.4999853512481308
0.4989418602397619
0.4982216344361741
0.5090229246013306
0.5372846051999657
0.4637955238165501
图11
实验结论：
X的值如上图所示。
附件（代码）：
1.P1081
#include
#include
usingnamespacestd;
//此函数用于计算矩阵的解
float*uptrbk（float**A,intN）
{
floatc;
float*x=newfloat[N-1];//生成一维数组x[N]
intn;
for（n=1;n{
c=A[n][n-1]/A[n-1][n-1];
A[n][n-1]=0;
A[n][n]=A[n][n]-c*A[n-1][n];
A[n][N]=A[n][N]-c*A[n-1][N];
}
x[N-1]=A[N-1][N]/A[N-1][N-1];
for（n=N-2;n>=0;n--）
{
A[n][N]=A[n][N]-x[n+1]*A[n][n+1];
x[n]=A[n][N]/A[n][n];
}
returnx;//带回计算结果数组x[N]的地址x
}
//实验的main函数
intmain（）
{
intN;
cout<<"请输入矩阵的阶数：
";//输入矩阵的阶数用于生成动态矩阵
cin>>N;
inti,k;
float*uptrbk（float**,int）;
float*x;
float**A=newfloat*[N-1];//生成动态增广矩阵
for（i=0;iA[i]=newfloat[N];
cout<<"请输入增广矩阵的值"<for（i=0;ifor（k=0;k<=N;k++）
cin>>A[i][k];
x=uptrbk（A,N）;//计算矩阵的解列向量X
cout<<"x的值为：
"<for（i=0;icout<<"x["<system（"pause"）;
return0;
}
2.P1201
#include
#include
usingnamespacestd;
//实验的main函数
intmain（）
{
intN,i,k;
float*x;
cout<<"请输入矩阵的阶数：
";//输入矩阵的阶数，用于生成动态矩阵
cin>>N;
N=N-1;
float*B=newfloat[N];
float*lufact（float**,float*,int）;
float**A=newfloat*[N];//生成动态系数矩阵A
for（i=0;i<=N;i++）
A[i]=newfloat[N];
cout<<"请输入矩阵A的值：
"<阵的值
for（i=0;i<=N;i++）
for（k=0;k<=N;k++）
cin>>A[i][k];
cout<<"请输入矩阵B的值：
"<for（i=0;i<=N;i++）
cin>>B[i];//输入矩阵B的值
x=lufact（A,B,N）;
cout<<"x的值为a：
"<for（i=0;i<=N;i++）//输出AX=B解的列向量X
cout<<"x["<system（"pause"）;
return0;
}
//使用LU法求解X的函数，A为系数矩阵，B为计算结果，N为矩阵阶数
float*lufact（float**A,float*B,intN）
{
inti,j,k;
floatc;
float**U;//生成二维数组U
float*x=newfloat[N];//生产保存结果的列矩阵X
float*y=newfloat[N];//生产用于保存中间值的列矩阵Y
U=A;
for（i=0;i<=N;i++）//将A转换为LU
{
for（j=i+1;j<=N;j++）
{
c=U[j][i]/U[i][i];
for（k=i;k<=N;k++）
U[j][k]=U[j][k]-c*U[i][k];
U[j][i]=c;
}
}
for（i=0;i<=N;i++）//计算中间矩阵Y的值
{
for（j=0;jB[i]=B[i]-y[j]*U[i][j];
y[i]=B[i];
}
for（i=N;i>=0;i--）//计算解X的值
{
for（j=N;j>i;j--）
y[i]=y[i]-x[j]*U[i][j];
x[i]=y[i]/U[i][i];
}
returnx;
}
3.P1202
#include
#include
#include
usingnamespacestd;
//实验的main函数
intmain（）
{
intN,i;
double*x;
double**A;
double*B;
double**buildA（int）;
double*buildB（int）;
double*lufact（double**,double*,int）;
cout<<"请输入矩阵的阶数：
";//输入矩阵的阶数，用于生成动态矩阵
cin>>N;
N=N-1;
A=buildA（N）;
B=buildB（N）;
x=lufact（A,B,N）;
cout<<"x的值为：
"<for（i=0;i<=N;i++）//输出AX=B解的列矩阵X
cout<<"x["<system（"pause"）;
return0;
}
//使用LU法求解X的函数，A为系数矩阵，B为计算结果，N为矩阵阶数
double*lufact（double**A,double*B,intN）
{
inti,j,k;
doublec;
double**U;//生成二维数组U，用于储存L&U矩阵
double*x=newdouble[N];//生成保存结果的列矩阵X
double*y=newdouble[N];//生成用于保存中间值的列矩阵Y
U=A;
for（i=0;i<=N;i++）//将A转换为LU矩阵
{
for（j=i+1;j<=N;j++）
{
c=U[j][i]/U[i][i];
for（k=i;k<=N;k++）
U[j][k]=U[j][k]-c*U[i][k];
U[j][i]=c;
}
}
for（i=0;i<=N;i++）//计算中间矩阵Y的值
{
for（j=0;jB[i]=B[i]-y[j]*U[i][j];
y[i]=B[i];
}
for（i=N;i>=0;i--）//计算解X的值
{
for（j=N;j>i;j--）
y[i]=y[i]-x[j]*U[i][j];
x[i]=y[i]/U[i][i];
}
returnx;
}
//用于产生阶数为N的矩阵A的函数
double**buildA（intN）
{
inti,j;
double**A=newdouble*[N];//生成二维动态数组A[N-1][N-1]
for（i=0;i<=N;i++）
A[i]=newdouble[N];
for（i=0;i<=N;i++）//产生矩阵A的元素并存储到A[N-1][N-1]中
for（j=0;j<=N;j++）
A[i][j]=pow（double（i+1）,j）;
returnA;
}
//用于产生1*N的列矩阵B
double*buildB（intN）
{
inti;
double*B=
升级会员 