lu分解法c语言
① 对矩阵x进行QR分解和LU分解,QR分解和LU分解是什么意思呢
你好!
LU分解是矩阵的三角分解,产生一个上三角矩阵和一个下三角矩阵。
QR分解是矩阵的正交分解。
如果对你有帮助,望采纳。
② 矩阵的LU分解该怎么具体做。亲,。举个例子吧
LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。
(A,E)~r2-2r1
1 -2 3 1 0 0
0 8 -4 -2 1 0
0 1 1 0 0 1
r3-1/8r2
1 -2 3 1 0 0
0 8 -4 -2 1 0
0 0 3/2 1/4 -1/8 1
该矩阵记做(U,P)
求出矩阵P=
1 0 0
-2 1 0
1/4 -1/8 1
的逆P-1=
1 0 0
2 1 0
0 1/8 1
因为PA=U,所以A=P-1U=LU,L=P-1
③ 怎样用LU分解法解线性方程组
Ax=B,改写成Ly=B,Ux=y的方程组。就相当于将A=LU分解成了两个矩阵。称为矩阵A的三角分解,或LU分解。如果L为单位下三角阵,则叫Doolittle分解,若U为单位上三角阵,则叫Crout分解。只要A的各顺序主子式不为零,则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积。
•设Ax=b,A=LU,则Ax=LUx=b
于是令Ux=y,则Ly=b
这样原来方程能化为两个简单方程组
在线性代数中, LU分解(LU Decomposition)是矩阵分解的一种,可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。LU分解主要应用在数值分析中,用来解线性方程、求反矩阵或计算行列式。
(3)lu分解法c语言扩展阅读:
相关算法:
LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵,其变换矩阵就是一个单位下三角矩阵。
这正是所谓的杜尔里特算法:从下至上地对矩阵A做初等行变换,将对角线左下方的元素变成零,然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵,这一系列单位下三角矩阵的乘积的逆就是L矩阵,它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方) 左右。
④ c语言能解决方程吗
我学过《数值分析》,李庆扬老师主编,施普辛格出版社。老师上课时虽然没有讲用C语言实现,但是给了流程图。我学过的用C语言解方程主要是牛顿迭代法,主要解非线性方程;解线性方程组的主要是LU分解法和列主元高斯消去法,这些东西都是我们的前辈早就研究出来的,但是现在可以用计算机方便地实现。
⑤ 计算方法LU分解法求解线性方程组 在线等
其中L为下三角矩阵,U为上三角矩阵。
例如,4×4矩阵A的情况,(1)式如下:
(2)
可以用如(1)式分解来求解线性方程组
(3)
首先求解向量y使得
(4)
然后再来求解
(5)
此拆分方法的优点在于求解一个三角形方程组相当容易,这样,(4)式可用向前替代过程求解,如下:
(6)
(5)式可用回代过程求解,这与(2)式~(3)式一样,
(7)
(6)式和(7)式共需执行N2次内层循环(对每个右端项b),每个内层循环包括一次乘法和一次加法。如果有N个右端项,它们是单位列向量(在求矩阵逆时就是这种情况),考虑这些零元素可把(6)式的总执行次数从N3/2减少到N3/6,而(7)式的执行次数不变,仍为N3/2。
注意:一点对A进行了LU分解,就可以一次求解所有要解的右端项。
算法实现:
首先,写出(1)式或(2)式的第i,j分量。它总是一个和式,开始部分形式如下:
和式中的项数依赖于i和j中较小的数。事实上有三种形式:
(8,9,10)
显然,(8)~(10)式共有N2个方程,而要求N2+N个未知的α和β(因对角线的未知元素有两套),既然未知数的个数比方程个数多,就人为指定N各位指数,然后再来求解其他的未知数。事实上,总是令
(11)
有一个算法称为Crout算法,它仅按某种次序排列方程,就能容易的求出(8)式~(11)式的N2+N各方程中的所有α和β。步骤如下:
设
,即(11)式
对每个j=0,1,2,...,N-1进行以下两步:
第一步,对每个i=0,1,...,j用(8)式、(9)式和(11)式来解βij,即
(12)
第二步,对每个i=j+1,j+2,...,N-1用(10)式来求解αij,即
(13)
在求解下一个j之前要保证进行了以上两步。
如果按上述过程进行几次迭代后,就会发现(12)式和(13)式右端的α和β在需要时已经得到,还会发现,每一个aij仅被使用一次就不再使用了。这意味着分解是“同址”进行的。简言之Crout算法得到的矩阵是混合矩阵,对本例排列如下:
注:不是把矩阵A分解成LU形式,而是将其按行置换的方式分解。
Crout算法的精妙之处:
l (12)式,在i=j(最后一次应用)时,与(13)式(除后者还要做一次除法外)是完全一样的,这两种情况要求和的上线都是k=j-1(=i-1)。这意味着,不必费心去考虑对角线元素βjj是否会正落在对角线上,也不必考虑该列中,它下面的某个元素(未做除法的)αij,i=j+1,j+2,...,N-1是否会提升成为对角线元素β。
l 它首先找到每行的最大元素,而后(在找最大主元时)乘以一个比例系数,这就实现了“隐式主元法”。
运行示例:
Origin coefficient matrix:
| 0.0 2.0 0.0 1.0 |
| 2.0 2.0 3.0 2.0 |
| 4.0 -3.0 0.0 1.0 |
| 6.0 1.0 -6.0 -5.0 |
-----------------------------------------------
LU mixed matrix:
| 6.0 1.0 -6.0 -5.0 |
| 0.0 2.0 0.0 1.0 |
| 0.3333333333333333 0.8333333333333334 5.0 2.833333333333333 |
| 0.6666666666666666 -1.8333333333333333 0.8 3.8999999999999995 |
-----------------------------------------------
Origin left-hand vector b:
| 0.0 |
| -2.0 |
| -7.0 |
| 6.0 |
-----------------------------------------------
Final solution vector:
| -0.5000000000000003 |
| 1.0000000000000002 |
| 0.33333333333333337 |
| -2.0000000000000004 |
-----------------------------------------------
示例程序:
package com.nc4nr.chapter02.lu;
public class LU ...{
// 4 * 4 coefficient matrix a
double[][] a = ...{
...{0.0, 2.0, 0.0, 1.0},
...{2.0, 2.0, 3.0, 2.0},
...{4.0, -3.0, 0.0, 1.0},
...{6.0, 1.0, -6.0, -5.0}
};
// 4 * 1 coefficient matrix b
double[] b = ...{
0.0,
-2.0,
-7.0,
6.0
};
int anrow = 4;
int[] indx = new int[anrow];
int parity = 1;
private void lucmp() ...{
final double tiny = 1.0e-20;
int imax = 0, n = anrow;
double big, m, sum, temp;
double[] vv = new double[n];
System.out.println("Origin coefficient matrix:");
output(a,4);
for (int i = 0; i < n; i++) ...{
big = 0.0;
for (int j = 0; j < n; j++) ...{
if ((temp = Math.abs(a[i][j])) > big) big = temp;
}
if (big == 0.0) System.out.println("lu: singular matrix in lucmp.");
vv[i] = 1.0/big;
}
for (int j = 0; j < n; j++) ...{
for (int i = 0; i < j; i++) ...{
sum = a[i][j];
for (int k = 0; k < i; k++) sum -= a[i][k]*a[k][j];
a[i][j] = sum;
}
big = 0.0;
for (int i = j; i < n; i++) ...{
sum = a[i][j];
for (int k = 0; k < j; k++) sum -= a[i][k]*a[k][j];
a[i][j] = sum;
if ((m = vv[i]*Math.abs(sum)) >= big) ...{
big = m;
imax = i;
}
}
if (j != imax) ...{
for(int k = 0; k < n; k++) ...{
m = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = m;
}
parity = -parity;
m = vv[imax];
vv[imax] = vv[j];
vv[j] = m;
}
indx[j] = imax;
if (a[j][j] == 0.0) a[j][j] = tiny;
if (j != n - 1) ...{
m = 1.0/a[j][j];
for (int i = j+1; i < n; i++) a[i][j] *= m;
}
}
System.out.println("LU mixed matrix:");
output(a,4);
}
private void lubksb() ...{
double sum;
int n = anrow, ii = 0;
System.out.println("Origin left-hand vector b:");
output(b,4);
for (int i = 0; i < n; i++) ...{
int ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii != 0)
for (int j = ii - 1; j < i; j++) sum -= a[i][j]*b[j];
else if (sum != 0.0)
ii = i + 1;
b[i] = sum;
}
for (int i = n-1; i >= 0; i--) ...{
sum = b[i];
for(int j = i + 1; j < n; j++) sum -= a[i][j]*b[j];
b[i] = sum / a[i][i];
}
System.out.println("Final solution vector:");
output(b,4);
}
private void output(double a[][], int anrow) ...{
for (int i = 0; i < anrow; i++) ...{
System.out.println(" | " + a[i][0] + " " +
a[i][1] + " " +
a[i][2] + " " +
a[i][3] + " | ");
}
System.out.println("-----------------------------------------------");
}
private void output(double[] b, int bnrow) ...{
for (int i = 0; i < bnrow; i++) ...{
System.out.println(" | " + b[i] + " | ");
}
System.out.println("-----------------------------------------------");
}
public LU() ...{
lucmp(); // 分解
lubksb(); // 回代
}
public static void main(String[] args) ...{
new LU();
}
}
⑥ 急球c语言!矩阵直接分解法(lu分解法)
以四维为例。系数在M中,常数项在N中
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void main()
{
int i,j,m,n;
double TM=0,TMm=0,TN=0,TNn=0;
int NN = 4;
double L[4][4]= {{1,0,0,0}, {0,1,0,0}, {0,0,1,0},{0,0,0,1}, };
// M coef 4x4, N const 4x1
double M[4][4]={4,-1,0,2,-1,4,-1,0,0,-1,4,-1,2,0,-1,4};
double N[4][1]={-1,-7,9,0};
for(i=1;i<NN;i++) {
for(j=0;j<i;j++) {
for(m=0;m<j;m++) {
TM=M[i][m] * L[j][m] + TMm;
TMm=TM;
TM=0;
};
M[i][j]=M[i][j]-TMm;
TMm=0;
L[i][j]=M[i][j]/M[j][j];
};
for(n=0;n<i;n++) {
TN=M[i][n]*L[i][n]+TNn;
TNn=TN;
TN=0;
}
M[i][i]=M[i][i]-TNn;
TNn=0;
};
// solve
for(i=0;i<NN;i++) {
for(m=0;m<i;m++) {
TM=L[i][m]*N[m][0]+TMm;
TMm=TM;
TM=0;
};
N[i][0]=N[i][0]-TMm;
TMm=0;
};
for(i=0;i<NN;i++) {
N[i][0]=N[i][0]/M[i][i];
};
for(i=NN-1;i>=0;i--) {
for(m=i+1;m<NN;m++) {
TM=L[m][i]*N[m][0]+TMm;
TMm=TM;
TM=0;
};
N[i][0]=N[i][0]-TMm;
TMm=0;
};
for(i=0;i<NN;i++) printf("%0.2f\n",N[i][0]);
}
⑦ 求C语言课程设计:用高斯列主元消元法解线性方程组
这里向你推荐一下克鲁特算法(其实就是对高斯列主元消元法进行优化,使之更适合于计算机编程),首先将矩阵A进行LU分解(将系数矩阵分解成一个上三角矩阵和一个下三角矩阵),分解的过程中用到了隐式的主元寻找法,同时利用克鲁特算法可以将两个n*n矩阵压缩到一个n*n矩阵中,大大节省了存储空间提高了计算速度。
方程可化为L*U*x=B,令U*x=y --->L*y=B
然后利用回代先求y,再利用y求x
因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算,所以它的计算速度应该能够达到高斯列主元消元法的三倍,但原理与其基本一致。
而且我在程序中使用了动态数组方便你今后进行扩展。
以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写,LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序
如果你需要详细的理论讲解我可以将这两本书和源程序发给你,上面的论述相当详细足够你答辩用的了,我的邮箱[email protected]
计算结果:
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=0;i<n;i++)//将A赋给x,B赋给y
{
y.push_back(B[i]);
for(j=0;j<n;j++)
{
x[i].push_back(A[i][j]);
}
}
cout<<"A矩阵:"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";
}
cout<<endl;
}
cout<<"B矩阵:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;
}
root(x,y);//求根
cout<<"解矩阵:"<<endl;
for(i=0;i<n;i++)
{
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语言数值计算法方法大全》第二章的程序
//如果你需要详细的理论讲解我可以将这两本书和源程序发给你,我的邮箱[email protected]
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)
{
int i,imax,j,k;
float big=0,m=0,sum=0,temp=0;
vector<float> vv(n);
vv.clear();
d=1.0;
for (i=0;i<n;i++)
{
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big)
big=temp;
vv[i]=1.0/big;
}
for (j=0;j<n;j++)
{
for (i=0;i<j;i++)
{
sum=a[i][j];
for (k=0;k<i;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++)
{
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (m=vv[i]*fabs(sum)) >= big)
{
big=m;
imax=i;
}
}
if (j != imax)
{
for (k=0;k<n;k++)
{
m=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=m;
}
d = -(d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0)
a[j][j]=TINY;
if (j != n)
{
m=1.0/(a[j][j]);
for (i=j+1;i<n;i++)
a[i][j] *= m;
}
}
}
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)
{
int i,ii=0,ip,j;
float sum;
for(i=0;i<n;i++)//按LU分解时寻找主元所进行的初等变换进行反边变换。
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
b[i]=sum;
}
sum=0;
for (i=1;i<n;i++)
{
sum=0;
for(j=0;j<i;j++)
{
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-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<n;j++)
{
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=0;i<n;i++)
{
x.push_back(line);
for(j=0;j<n;j++)
{
x[i].clear();
}
}
}