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分解的Fortran子程序 希望可以有所帮助
!求解au=b,u
!n表示为方程维数
subroutine lu(a,b,n,u)
implicit real(8) (a-h,o-z)
real(8)::a(n,n),b(n),u(n),a_bak(n,n),b1(n),aL(n,n),aU(n,n),y(n)
!exchange rows
do i=1,n
tmpMax=0.d0
do ic=i,n
if(tmpMax<dabs(a(ic,i))) then
tmpMax=dabs(a(ic,i))
i_rec=ic
endif
enddo
if(i_rec.ne.i) then
do jc=i,n
tmp=a(i,jc)
a(i,jc)=a(i_rec,jc)
a(i_rec,jc)=tmp
enddo
tmp=b(i)
b(i)=b(i_rec)
b(i_rec)=tmp
endif
!decomposition
do j=i,n
tmp=0.d0
do k=1,i-1
tmp=tmp+aL(i,k)*aU(k,j)
enddo
aU(i,j)=a(i,j)-tmp
tmp=0.d0
do k=1,i-1
tmp=tmp+aL(j,k)*aU(k,i)
enddo
aL(j,i)=(a(j,i)-tmp)/aU(i,i)
enddo
enddo
!find answer
do i=1,n
tmp=0.d0
do j=1,i-1
tmp=tmp+aL(i,j)*y(j)
enddo
y(i)=b(i)-tmp
enddo
do i=n,1,-1
tmp=0.d0
do j=i+1,n
tmp=tmp+aU(i,j)*u(j)
enddo
u(i)=(y(i)-tmp)/aU(i,i)
enddo
end
其中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, dum, sum, temp
double[] vv = new double[n]
System.out.println("Origin coefficient matrix:")
output(a,4)
for (int i = 0i <ni++) ...{
big = 0.0
for (int j = 0j <nj++) ...{
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 = 0j <nj++) ...{
for (int i = 0i <ji++) ...{
sum = a[i][j]
for (int k = 0k <ik++) sum -= a[i][k]*a[k][j]
a[i][j] = sum
}
big = 0.0
for (int i = ji <ni++) ...{
sum = a[i][j]
for (int k = 0k <jk++) sum -= a[i][k]*a[k][j]
a[i][j] = sum
if ((dum = vv[i]*Math.abs(sum)) >= big) ...{
big = dum
imax = i
}
}
if (j != imax) ...{
for(int k = 0k <nk++) ...{
dum = a[imax][k]
a[imax][k] = a[j][k]
a[j][k] = dum
}
parity = -parity
dum = vv[imax]
vv[imax] = vv[j]
vv[j] = dum
}
indx[j] = imax
if (a[j][j] == 0.0) a[j][j] = tiny
if (j != n - 1) ...{
dum = 1.0/a[j][j]
for (int i = j+1i <ni++) a[i][j] *= dum
}
}
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 = 0i <ni++) ...{
int ip = indx[i]
sum = b[ip]
b[ip] = b[i]
if (ii != 0)
for (int j = ii - 1j <ij++) sum -= a[i][j]*b[j]
else if (sum != 0.0)
ii = i + 1
b[i] = sum
}
for (int i = n-1i >= 0i--) ...{
sum = b[i]
for(int j = i + 1j <nj++) 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 = 0i <anrowi++) ...{
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 = 0i <bnrowi++) ...{
System.out.println(" | " + b[i] + " | ")
}
System.out.println("-----------------------------------------------")
}
public LU() ...{
lucmp()// 分解
lubksb()// 回代
}
public static void main(String[] args) ...{
new LU()
}
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）