高斯消去法lu分解法的优缺点

高斯消去法lu分解法的优点：高斯消元法的算法复杂度是O(n3)；这就是说，如果系数矩阵的是n × n，那么高斯消元法所需要的计算量大约与n3成比例，高斯消元法可用在任何域中。

高斯消去法lu分解法的缺点：高斯消元法对于一些矩阵来说是稳定的。对于普遍的矩阵来说，高斯消元法在应用上通常也是稳定的，不过亦有例外。高斯消去法是线性代数中的一个算法，可用来为线性方程组求解，求出矩阵的秩，以及求出可逆方阵的逆矩阵。

算法

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。

这正是所谓的杜尔里特算法：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方) 左右。

这里向你推荐一下克鲁特算法（其实就是对高斯列主元消元法进行优化，使之更适合于计算机编程），首先将矩阵A进行LU分解（将系数矩阵分解成一个上三角矩阵和一个下三角矩阵），分解的过程中用到了隐式的主元寻找法，同时利用克鲁特算法可以将两个n*n矩阵压缩到一个n*n矩阵中，大大节省了存储空间提高了计算速度。

方程可化为L*U*x=B，令U*x=y --->L*y=B

然后利用回代先求y，再利用y求x

因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算，所以它的计算速度应该能够达到高斯列主元消元法的三倍，但原理与其基本一致。

而且我在程序中使用了动态数组方便你今后进行扩展。

以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

如果你需要详细的理论讲解我可以将这两本书和源程序发给你，上面的论述相当详细足够你答辩用的了，我的邮箱hu_hu605@163.com

计算结果：

A矩阵:

2 2 5

3 4 7

1 3 3

B矩阵:

5

6

5

解矩阵:

x 1=-7

x 2=0.333333

x 3=3.66667

Press any key to continue

#include <cmath>

#include <iostream>

#include <iomanip>

#include <cstdlib>

#include <functional>

#include <vector>

#include <algorithm>

using namespace std

#define TINY 1.0e-20 //A small number.

#define N 3

void ludcmp(vector<vector<float>>&a, int n, vector<int>&indx, float &d)//对矩阵进行LU分解

void lubksb(vector<vector<float>>&a, int n, vector<int>&indx, vector<float>&b)//对矩阵进行前向和后向回代

void root(vector<vector<float>>&x,vector<float>&col)//解方程结果保存在y中

void iniv(vector<vector<float>>&x,vector<float>line,int n)//对二维动态数组进行初始化

void main()

{

int i,j,n=N//输入矩阵的维数

float A[N][N]={{2,2,5},{3,4,7},{1,3,3}}//左边A矩阵

float B[N]={5,6,5}//右边B矩阵

vector<vector<float>>x//建立动态二维数组存放A，保证你的程序进行扩展时只改A，B，N

vector<float>line

vector<float>y(n)//建立动态数组存放B

iniv(x,line,n)

y.clear()

for(i=0i<ni++)//将A赋给x，B赋给y

{

y.push_back(B[i])

for(j=0j<nj++)

{

x[i].push_back(A[i][j])

}

}

cout<<"A矩阵:"<<endl

for(i=0i<ni++)

{

for(j=0j<nj++)

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" "

}

cout<<endl

}

cout<<"B矩阵:"<<endl

for(i=0i<ni++)

{

cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl

}

root(x,y)//求根

cout<<"解矩阵:"<<endl

for(i=0i<ni++)

{

cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl

}

cout<<endl

}

void root(vector<vector<float>>&x,vector<float>&col)

{

int n=x.size(),i=0,j=0

vector<int>index(n)//用于记录寻找主元素过程中对矩阵的初等变换

index.clear()

float m=1.0//记录变换方式，此程序中无用

ludcmp(x,n,index,m)//进行LU分解

lubksb(x,n,index,col)//根据分解结果进行回带

}

//以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序

//如果你需要详细的理论讲解我可以将这两本书和源程序发给你，我的邮箱hu_hu605@163.com

void ludcmp(vector<vector<float>>&a, int n, vector<int>&indx, float &d)

{

int i,imax,j,k

float big=0,dum=0,sum=0,temp=0

vector<float>vv(n)

vv.clear()

d=1.0

for (i=0i<ni++)

{

big=0.0

for (j=0j<nj++)

if ((temp=fabs(a[i][j])) >big)

big=temp

vv[i]=1.0/big

}

for (j=0j<nj++)

{

for (i=0i<ji++)

{

sum=a[i][j]

for (k=0k<ik++)

sum -= a[i][k]*a[k][j]

a[i][j]=sum

}

big=0.0

for (i=ji<ni++)

{

sum=a[i][j]

for (k=0k<jk++)

sum -= a[i][k]*a[k][j]

a[i][j]=sum

if ( (dum=vv[i]*fabs(sum)) >= big)

{

big=dum

imax=i

}

}

if (j != imax)

{

for (k=0k<nk++)

{

dum=a[imax][k]

a[imax][k]=a[j][k]

a[j][k]=dum

}

d = -(d)

vv[imax]=vv[j]

}

indx[j]=imax

if (a[j][j] == 0.0)

a[j][j]=TINY

if (j != n)

{

dum=1.0/(a[j][j])

for (i=j+1i<ni++)

a[i][j] *= dum

}

}

}

void lubksb(vector<vector<float>>&a, int n, vector<int>&indx, vector<float>&b)

{

int i,ii=0,ip,j

float sum

for(i=0i<ni++)//按LU分解时寻找主元所进行的初等变换进行反边变换。

{

ip=indx[i]

sum=b[ip]

b[ip]=b[i]

b[i]=sum

}

sum=0

for (i=1i<ni++)

{

sum=0

for(j=0j<ij++)

{

sum+=a[i][j]*b[j]

}

b[i]=b[i]-sum

}

b[n-1]=b[n-1]/a[n-1][n-1]

for (i=n-2i>=0i--)

{

sum=0

for(j=i+1j<nj++)

{

sum+=a[i][j]*b[j]

}

b[i]=(b[i]-sum)/a[i][i]

}

}

void iniv(vector<vector<float>>&x,vector<float>line,int n)

{

int i,j

for(i=0i<ni++)

{

x.push_back(line)

for(j=0j<nj++)

{

x[i].clear()

}

}

}

我们以方程组2x1 + 6x2 - x3 = -12

5x1 - x2 +2x3 = 29

-3x1 - 4x2 + x3 = 5

为例 来说明楼主自己把方程组化为矩阵形式。以下为源代码 。

#include <stdio.h>

#include <stdlib.h>

#include <malloc.h>

#include <math.h>

int GS(int,double**,double *,double)

double **TwoArrayAlloc(int,int)

void TwoArrayFree(double **)

int main(void)

{

int i,n

double ep,**a,*b

n = 3

ep = 1e-4

a = TwoArrayAlloc(n,n)

b = (double *)calloc(n,sizeof(double))

if(b == NULL)

{

printf("memory get error\n")

exit(1)

}

a[0][0]= 2a[0][1]= 6a[0][2]= -1

a[1][0]= 5a[1][1]=-1a[1][2]= 2

a[2][0]=-3a[2][1]=-4a[2][2]= 1

b[0] = -12b[1] = 29 b[2] = 5

if(!GS(n,a,b,ep))

{

printf("can't solve it with GS elimination\n")

exit(0)

}

printf("The solution of equations is as follows：\n")

for(i=0i<3i++)

{

printf("x%d = %.2f\n",i,b[i])

}

TwoArrayFree(a)

free(b)

return 0

}

int GS(n,a,b,ep)

int n

double **a

double *b

double ep

{

int i,j,k,l

double t

for(k=1k<=nk++)

{

for(l=kl<=nl++)

if(fabs(a[l-1][k-1])>ep)

break

else if(l==n)

return(0)

if(l!=k)

{

for(j=kj<=nj++)

{

t = a[k-1][j-1]

a[k-1][j-1] =a[l-1][j-1]

a[l-1][j-1] =t

}

t=b[k-1]

b[k-1]=b[l-1]

b[l-1]=t

}

t=1/a[k-1][k-1]

for(j=k+1j<=nj++)

a[k-1][j-1]=t*a[k-1][j-1]

b[k-1]*=t

for(i=k+1i<=ni++)

{

for(j=k+1j<=nj++)

a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1]

b[i-1]-=a[i-1][k-1]*b[k-1]

}

}

for(i=n-1i>=1i--)

for(j=i+1j<=nj++)

b[i-1]-=a[i-1][j-1]*b[j-1]

return(1)

}

double **TwoArrayAlloc(int r,int c)

{

double *x,**y

int n

x=(double *)calloc(r*c,sizeof(double))

y=(double **)calloc(r,sizeof(double*))

for(n=0n<=r-1++n)

{

y[n]=&x[c*n]

}

return y

}

void TwoArrayFree(double **x)

{

free(x[0])

free(x)

}

---