线性代数方程组求解

线性代数方程组求解 The Standardization Office was revised on the afternoon of December 13, 2020
线性代数方程组求解
一、实验要求
编程求解方程组：
方程组1：
方程组2：
方程组3：
要求：
用C/C++语言实现如下函数：
1.bool lu(double* a, int* pivot, int n);
实现矩阵的LU分解。

pivot为输出参数，pivot[0,n) 中存放主元的位置排列。

函数成功时返回false ，否则返回true 。

2. bool guass(double const* lu, int const* p, double* b, int n); 求线代数方程组的解
设矩阵Lunxn 为某个矩阵anxn 的LU 分解，在内存中按行优先次序存放。

p[0,n)为LU 分解的主元排列。

b 为方程组Ax=b 的右端向量。

此函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。

函数成功时返回false ，否则返回true 。

3. void qr(double* a, double* d, int n);矩阵的QR 分解
假设数组anxn 在内存中按行优先次序存放。

此函数使用HouseHolder 变换将其就地进行QR 分解。

d 为输出参数，d [0,n) 中存放QR 分解的上三角对角线元素。

4. bool hshld(double const*qr, double const*d, double*b, int n); 求线代数方程组的解
设矩阵qrnxn 为某个矩阵anxn 的QR 分解，在内存中按行优先次序存放。

d [0,n) 为QR 分解的上三角对角线元素。

b 为方程组Ax=b 的右端向量。

函数计算方程组Ax=b 的解，并将结果存放在数组b[0,n)中。

函数成功时返回false ，否则返回true 。

二、问题分析
求解线性方程组Ax=b ，其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。

因此，在求解线性方程组的过程中，把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要，然而矩阵A 的变换通常有两种分解方法：LU 分解法和QR 分解法。

1、LU 分解法：
将A 分解为一个下三角矩阵L 和一个上三角矩阵U ，即：A=LU,
其中 L=⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎣⎡10010012121 n n l l l ， U=⎥
⎥
⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡nn n n u u u u u u 000
00222112
11 2、QR 分解法：
将A 分解为一个正交矩阵Q 和一个上三角矩阵R,即：A=QR
三、实验原理
解Ax=b 的问题就等价于要求解两个三角形方程组： ⑴ Ly=b,求y; ⑵ Ux=y,求x.
设A 为非奇异矩阵,且有分解式A=LU , L 为单位下三角阵,U 为上三角阵。

L,U 的元素可以有n 步直接计算定出。

用直接三角分解法解Ax=b （要求A 的所有顺序主子式都不为零）的计算公式： ① ),,2,1(n i a u li li ==,11/u a l il il = ,i=2,3,…，n. 计算U 的第r 行，L 的第r 列元素(i=2,3,…,n)： ② ∑-=-=1
1r k ki rk ri ri u l a u ， i=r,r+1,…,n;
③ rr r k kr ik ir ir u u l a l /)(1
1
∑-=-= , i=r+1,…,n,且r ≠n.
求解Ly=b ，Ux=y 的计算公式；
④
:
,3,2,,
1
1
11n i y l b y b y i k k ik i i =-==∑-=
⑤
.
1,,2,1,/)(,/1
--=-
==∑+=n n i u x u y x u y x ii n
i k k ik i i nn n n
四、实验步骤
1>将矩阵A 保存进计算机中，再定义2个空矩阵L ，U 以便保存求出的三角矩阵的值。

利用公式①，②，③将矩阵A 分解为LU ，L 为单位下三角阵,U 为上三角阵。

2>可知计算方法有三层循环。

先通过公式①计算出U 矩阵的第一行元素li u 和L 矩阵的第一列元素il l 。

再根据公式②和③，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。

3>先由公式④ ，Ly=b
:
,3,2,,1
1
11n i y l b y b y i k k ik i i =-==∑-=
求出y ，因为L 为下三角矩阵，所以由第一行开始求y. 4>再由公式⑤，Ux=y
.
1,,2,1,/)(,/1
--=-
==∑+=n n i u x u y x u y x ii n
i k k ik i i nn n n
求出x, 因为U 为上三角矩阵，所以由最后一行开始求x.
五、程序流程图
1、LU 分解法
回代过
程
2、QR分解法
六、实验结果
1、LU分解法
方程组1 ：
方程组2：
方程组3：
2、QR分解法
方程组1：
方程组2：
方程组3：
七、实验总结
为了求解线性方程组，我们通常需要一定的解法。

其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。

在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。

1、三角分解法
三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。

它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。

不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。

2、 QR分解法
QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。

在编写这两个程序过程中，起初遇到不少麻烦！虽然课上老师反复重复着：“算法不难的，It's so easy!”但是当自己实际操作时，感觉并不是那么容易。

毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。

每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。

回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。

附源代码:
#include <iostream>
#include <>
#include <>
#include <>
#include <>
using namespace std;
bool lu(double* a, int* pivot, int n);6f\t",B[i]); }
cout<<"\nGuass解法的误差为:\n";
for (int i=0; i<n; i++)
{
if(i==3)
cout<<endl;
printf("%.16f\t",B[i]-1);
expct=expct + B[i]-1;
}
printf("\n误差期望值为%.16f\t\n",expct/6); cout<<endl;
}
else
{
cout<<"矩阵QR分解: \n";
qr(A,D,n);6f\t",B[i]);
}
cout<<"\nHouseholder解法的误差为:\n";
for (int i=0; i<n; i++)
{
if(i==3)
cout<<endl;
printf("%.16f\t",B[i]-1);
expct=expct + B[i]-1;
}
printf("\n误差期望值为%.16f\t\n",expct/6); cout<<endl;
}
return 0 ;
}
bool lu(double* a, int* pivot, int n)//矩阵LU分解{
int i,j,k;
double max,temp;
max = 0;
for (i=0; i<n-1; i++)//依次对第i列进行处理
{
// 选出i列的主元,记录主元位置
max = fabs(a[n*i + i]);
pivot[i]=i;
for(j=i+1; j<n; j++)//j表示第j行
{
if( fabs(a[n*j + i])>max)
{
max= fabs(a[n*j + i]) ;
pivot[i]=j;
}
}
// 对第i列进行行变换，使得主元在对角线上
if(pivot[i]!=i)
{
for(j=i; j<n; j++)//ij与pivot[i]j换只用对上三角进行处理
{
temp=a[n*i + j];
a[n*i + j]=a[n*pivot[i]+ j];
a[n*pivot[i]+ j]=temp;
}
}
for(j=i+1; j<n; j++)//Pi 部分下三角L
a[n*j + i]=a[n*j+i]/a[n*i+i];
for(j=i+1; j<n; j++)//计算上三角U
for(k=i+1; k<n; k++)
a[n*j + k]=a[n*j+k]- a[n*j+i]*a[n*i+k];
}
//计算下三角 L
for(i=0; i<n-2; i++)//i行k列
for(k=i+1; k<n-1;k++)
{
temp=a[n*pivot[k] + i];
a[n*pivot[k] + i]=a[k*n + i];
a[k*n + i]=temp;
}
return false ;
}
bool guass(double const* lu, int const* p, double* b, int n)//求线性代数方程组的解
{
double temp;
//按qivot对b行变换，与LU匹配
for(i=0; i<n-1; i++) //貌似错误在这里哦
{
temp = b[p[i]];
b[p[i]] = b[i];
b[i]=temp;
}
//Ly=b,将y的内容放入b
for(i=0; i<n; i++)
for(j=0; j<i; j++)
b[i]=b[i]-lu[n*i+j]*b[j];
//Uy=x,将x的内容放入b
for(i=n-1; i>=0; i--)
{
for(j=n-1; j>i; j--)
b[i]=b[i]-lu[n*i+j]*b[j];
b[i]=b[i]/lu[n*i+i];
}
return false;
}
void qr(double* a, double* d, int n) //矩阵的QR分解{
int i,j,l,k;
double tem,m;
double *temp;
temp = (double *)malloc(sizeof(double)*n);
for (i=0; i<n-1; i++)//依次对第i列进行处理
{
//获得tem值
m = 0 ;
for(j=i; j<n; j++)//j表示第j行
m = m +a[n*j + i]*a[n*j + i ];
if(a[n*i + i ]>0)
m = -sqrt(m);
else
m = sqrt(m);
//获得temp放入矩阵，并存主元d
tem = 0 ;
d[i] = m ;
a[n*i +i] = a[n*i +i] - m;
for(j=i; j<=n-1; j++)
tem=tem + a[n*j +i]*a[n*j +i];
tem= sqrt(tem);
for(j=i; j<=n-1; j++)
a[n*j +i] = a[n*j +i]/tem ;
// 调整矩阵
for(k=i+1;k<n;k++)
{
for(j=i; j<n; j++)
{
tem = 0 ;
for(l=i; l<n; l++)
tem =tem + a[n*j + i]*a[n*l + i]*a[n*l + k];
temp[j] = a[j*n+k] - 2*tem;
}
for(j=i; j<n; j++)
a[j*n+k] = temp[j];
}
}
d[n-1] = a[(n-1)*n+n-1];
}
bool householder(double const*qr, double const*d, double*b, int n)//求线性代数方程组的解
{
int i,j,l;
double rem;
double *temp;
temp = (double *)malloc(sizeof(double)*n);
for(i=0; i<n-1; i++)
{
for(j=i; j<n; j++)
{
rem = 0;
for(l=i;l<n; l++)
rem = rem + qr[l*n+i]*qr[j*n+i]*b[l];
temp[j] = b[j] - 2*rem;
}
for(j=i; j<n; j++)
b[j] = temp[j];
}
for(j=n-1; j>-1; j--)
{
for(l=n-1; l!=j;--l)
b[j] =b[j] - b[l]*qr[j*n+l];
b[j] = b[j] /d[j];
}
return false; }。