LU和QR法解线性方程组
一、 问题描述
求解方程组==
,
要求：
1、编写用三角（LU）分解法求解线性方程组；
2、编写用正交三角（QR）分解法求解线性方程组。
2、问题分析
求解线性方程组Ax=b，其实质就是把它的系数矩阵A通过各种变换成一个下三角或上三角矩阵，从而简化方程组的求解。因此，在求解线性方程组的过程中，把系数矩阵A变换成上三角或下三角矩阵显得尤为重要，然而矩阵A的变换通常有两种分解方法：LU分解法和QR分解法。
1、LU分解法：
将A分解为一个下三角矩阵L和一个上三角矩阵U，即：A=LU,
其中 L=， U=
2、QR分解法：
将A分解为一个正交矩阵Q和一个上三角矩阵R,即：A=QR
三、实验原理
1、LU分解法
解Ax=b 的问题就等价于要求解两个三角形方程组： ⑴ Ly=b,求y; ⑵ Ux=y,求x.
设A为非奇异矩阵,且有分解式A=LU, L为单位下三角阵,U为上三角阵。
L,U的元素可以有n步直接计算定出。用直接三角分解法解Ax=b（要求A的所有顺序主子式都不为零）的计算公式：
① ,
,i=2,3,…，n.计算U的第r行，L的第r列元素(i=2,3,…,n)： ②
， i=r,r+1,…,n;
③ , i=r+1,…,n,且r≠n.求解Ly=b，Ux=y的计算公式；
④
⑤
2、QR分解法
4、实验步骤
1、LU分解法
1>将矩阵A保存进计算机中，再定义2个空矩阵L，U以便保存求出的三角矩阵的值。
利用公式①，②，③将矩阵A分解为LU，L为单位下三角阵,U为上三角阵。
2>可知计算方法有三层循环。
先通过公式①计算出U矩阵的第一行元素 和L矩阵的第一列元素
。
再根据公式②和③，和上次的出的值，求出矩阵其余的元素，每次都要三次循环，求下一个元素需要上一个结果。
3>先由公式④ ，Ly=b
求出y，因为L为下三角矩阵，所以由第一行开始求y.
4>再由公式⑤，Ux=y
求出x, 因为U为上三角矩阵，所以由最后一行开始求x.
2、QR分解法
5、程序流程图
1、LU分解法
word/media/image16_1.png
word/media/image17_1.png
2、QR分解法
word/media/image18_1.png
6、实验结果
1、LU分解法
2、QR分解法
7、实验总结
为了求解线性方程组，我们通常需要一定的解法。其中一种解法就是通过矩阵的三角分解来实现的，属于求解线性方程组的直接法。在不考虑舍入误差下，直接法可以用有限的运算得到精确解，因此主要适用于求解中小型稠密的线性方程组。
1、三角分解法
三角分解法是将A矩阵分解成一个上三角形矩阵U和一个下三角形矩阵L，这样的分解法又称为LU分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程，求反矩阵和求解联立方程组。不过要注意这种分解法所得到的上下三角形矩阵并非唯一，还可找到数个不同 的一对上下三角形矩阵，此两三角形矩阵相乘也会得到原矩阵。
2、 QR分解法
QR分解法是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为QR分解法。
在编写这两个程序过程中，起初遇到不少麻烦！虽然课上老师反复重复着：“算法不难的，It's so easy!”但是当自己实际操作时，感觉并不是那么容易。毕竟是要把实际的数学问题转化为计算机能够识别的编程算法，所以在编写程序之前我们仔细认真的把所求解的问题逐一进行详细的分析，最终转化为程序段。每当遇到问题时，大家或许有些郁闷，但最终还是静下心来反复仔细的琢磨，一一排除了错误，最终完成了本次实验。
回头一想原来编个程序其实也没有想象的那么复杂，只要思路清晰，逐步分析，就可以慢慢搞定了。
通过这次实验，让我们认知到团队的作用力度是不容忽视的，以后不管干任何时都要注重团队合作，遇到不懂得不明白的大家一起讨论，越讨论越清晰，愈接近最优答案。这样不管干什么都能事半功倍。庆幸自己有这么个团队，也明白大家一起劳动的果实最珍贵。
附 源代码:
LR分解法：
#include
void input_A(); //输入矩阵A
void input_b(); //输入矩阵b
void output_x(float x[4]); //输出方程组的根
void main()
{
int i,k,r;
float A[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}},b[4]={-2,-7,-7,-3},x[4],u[4][4],l[4][4],y[4]; //给定的方程组
// input_A();
// input_b();
//显示矩阵A、b
printf("矩阵A[4][4]:\n");
for(i=0;i<4;i++)
{
for(k=0;k<4;k++)
printf("%-10f",A[i][k]);
printf("\n");
}
for(i=0;i<4;i++)
u[0][i]=A[0][i];
for(i=0;i<4;i++)
{
l[i][0]=A[i][0]/u[0][0];
l[i][i]=1;
}
for(r=1;r<4;r++) //计算矩阵L和U
{
for(i=r;i<4;i++)
{
float t=0;
for(k=0;k
t=t+l[r][k]*u[k][i];
u[r][i]=A[r][i]-t;
}
for(i=r;i<4;i++)
{
float t1=0;
for(k=0;k
t1+=l[i][k]*u[k][r];
l[i][r]=(A[i][r]-t1)/u[r][r];
}
}
y[0]=b[0];
for(i=1;i<4;i++)
{
float t=0;
for(k=0;k
t+=l[i][k]*y[k];
y[i]=b[i]-t;
}
for(i=3;i>=0;i--)
{
float t=0;
for(k=i+1;k<4;k++)
t+=u[i][k]*x[k];
x[i]=(y[i]-t)/u[i][i];
}
printf("*****************************\n");
output_x(x);
}
//输入矩阵A
void input_A()
{
float A[4][4];
int i,j;
printf("input matrix A[4][4]:\n");
for(i=0;i<4;i++)
for(j=0;j<4;j++)
scanf("%f",&A[i][j]);
}
//输入矩阵b
void input_b()
{
float b[4];
int i;
printf("input matrix b[4]:\n");
for(i=0;i<4;i++)
scanf("%f",&b[i]);
}
//输出方程的根x
void output_x(float x[4])
{
int i;
printf("方程组的根为:\n");
for(i=0;i<4;i++)
printf("x[%d]=%-13f",i+1,x[i]);
printf("\n");
}
QR分解法：
#include
#include
#define N 4
void matrix_time(double A[][N],double B[][N],double C[][N],int n); //两个矩阵相乘,结果存在矩阵C[][N]
void matrix_vec(double A[][N],double B[N],double C[N],int n); //矩阵和向量相乘,结果存在向量C[N]
double vec_value(double A[],int n); //求向量的模
void vec_time(double a[],double H[][N],int n); //两个向量相乘得一个矩阵；
void householder(double *a,double H[][N],int n, int m); //求解Householder矩阵函数
void matrix_turn(double A[][N],int n); //求矩阵装置
void solution(double A[][N],double b[],int n); //求解上三角方程的解
void QR(double A[][N],double b[],int n);
void main()
{
double A[4][4]={{4,2,1,5},{8,7,2,10},{4,8,3,6},{12,6,11,20}};
double b[4]={-2,-7,-7,-3};
int i,j;
int n=4;
printf("待求解的方程组系数矩阵A为：\n");
for(i=0;i
{
for(j=0;j
printf("%-10f",A[i][j]);
printf("\n");
}
QR(A,b,n);
}
void matrix_time(double A[][N],double B[][N],double C[][N],int n) //两个矩阵相乘，结果存在矩阵C[][N]
{
int i,j,k;
for(i=0;i
for(j=0;j
{
C[i][j]=0;
for(k=0;k
C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
}
void matrix_vec(double A[][N],double B[N],double C[N],int n) //矩阵和向量相乘,结果存在向量C[N]
{
int i,j;
for(i=0;i
{
C[i]=0;
for(j=0;j
C[i]=C[i]+A[i][j]*B[j];
}
}
double vec_value(double A[],int n) //求向量的模
{
double Sum=0;
int i;
for(i=0;i
Sum=Sum+A[i]*A[i];
Sum=sqrt(Sum);
return Sum;
}
void vec_time(double a[],double H[][N],int n) //两个向量相乘得一个矩阵；
{
int i,j;
for(i=0;i
for(j=0;j
H[i][j]=a[i]*a[j];
}
void householder(double *a,double H[][N],int n, int m) //计算Householder矩阵
{
int i,j;
double first; //存放首元素
double vec_v; //存放向量的模
double a1[N]={0};
for(i=0;i
for(j=0;j
{
if(i==j)
H[i][j]=1;
else
H[i][j]=0;
}
first=a[m]; //取矩阵A转置的第m行（取矩阵A的第m列数）
vec_v=vec_value(&a[m],n-m); //计算m列元素所构成向量的模
a[m]=a[m]-vec_v; //w的分子部分
vec_v=vec_value(&a[m],n-m); //w的分母部分
vec_v=vec_v*vec_v;
for(i=m;i
for(j=m;j
H[i][j]+=-a[i]*a[j]*2/vec_v;
}
void matrix_turn(double A[][N],int n) //求矩阵的转置
{
double a[N][N]={0};
int i,j;
for(i=0;i
for(j=0;j
a[i][j]=A[i][j];
for(i=0;i
for(j=0;j
A[i][j]=a[j][i];
}
void solution(double A[][N],double b[],int n) //求解上三角方程的解
{
int i,j;
double x[N]={0};
double sum=0;
for(i=n-1;i>=0;i--)
{
for(j=n-1;j>i;j--)
sum+=A[i][j]*x[j];
x[i]=(b[i]-sum)/A[i][i];
sum=0;
}
printf("矩阵的QR分解求解结果为：\n");
for(i=0;i
printf("X[%d]=%-10f\n",i+1,x[i]);
}
void QR(double A[][N],double b[],int n)
{
int i,j,k,t;
double H1[N][N]={0},H2[N][N]={0},H3[N][N]={0}; //H1,H2存放相邻两次的Householder矩阵，H3为中间变量矩阵
double temp[N][N]={0};
double Q[N][N]={0},R[N][N]={0},Q1[N][N]={0};
double b1[N]={0};
for(i=0;i
for(j=0;j
temp[i][j]=A[i][j];
for(i=0;i
{
H2[i][i]=1; //单位阵
}
matrix_turn(temp,n); //矩阵A的转置 （方便后续取A的某一列元素）
for(i=0;i
{
householder(temp[i],H1,n,i); //得到Householder矩阵H1
matrix_time(H1,A,Q,n); //矩阵H1和A相乘放在Q中
for(k=0;k
for(t=0;t
A[k][t]=Q[k][t];
for(k=0;k
for(t=0;t
temp[k][t]=A[k][t];
matrix_turn(temp,n); //对矩阵转置
matrix_time(H1,H2,H3,n);
for(k=0;k
for(t=0;t
H2[k][t]=H3[k][t];//H2是一次变换与A相乘后的矩阵（H1*A）
}
for(k=0;k
for(t=0;t
R[k][t]=A[k][t];
printf("*********************************\n");
printf("系数矩阵A进行QR分解后的R矩阵为：\n");
for(i=0;i
{
for(j=0;j
printf("%-10f",R[i][j]);
printf("\n");
}
for(k=0;k
for(t=0;t
Q[k][t]=H2[k][t]; //此时Q为最终正交矩阵Q的转置
matrix_vec(Q,b,b1,n);//Q矩阵为正交阵，故Q的逆就等于Q装置，求解Q的逆
matrix_turn(Q,n); //转置，求解Q
printf("*********************************\n");
printf("系数矩阵A进行QR分解后的Q矩阵为：\n");
for(i=0;i
{
for(j=0;j
printf("%-10f",Q[i][j]);
printf("\n");
}
solution(R,b1,n); //求解上三角方程
}
本文来源：https://www.2haoxitong.net/k/doc/5e2c5f8c680203d8ce2f24d5.html
文档为doc格式