求matlab LU 分解程序

matlab有多种LU分解程序下面算一种：

function [L,U]=myLU(A)

%实现对矩阵A的LU分解，L为下三角矩阵

A

[n,n]=size(A)

L=zeros(n,n)

U=zeros(n,n)

for i=1:n

L(i,i)=1

end

for k=1:n

for j=k:n

U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)')

end

for i=k+1:n

L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k)

end

end

用法，在控制台输入

A=[1 2 3 -4-3 -4 -12 132 10 0 -34 14 9 -13]

当G接近奇异时，有的奇异值较小，此时由于W中的元素相差过大，导致条件数极大，逆矩阵的计算误差较大，方程的解极不稳定。为了解决这个问题，维根斯(Wiggins)提出去掉较小的奇异值，用r×r矩阵We来代替W，从而有较稳定的广义逆[9]:

地球物理反演教程

其中:

地球物理反演教程

这样方程Gm=d就有广义逆解:

地球物理反演教程

下面以三个例子来说明如何使用奇异值分解程序[1]。

地球物理反演教程

下面程序段是维根斯法的具体实现步骤。用到两个子程序:奇异值分解子程序svdcmp和回代解线性方程组子程序svbksb。这两个子程序可以从Fortran power station 4.0中获得源程序。

首先调用svdcmp子程序进行奇异值分解，找到最大奇异值wmax然后设定最小奇异值和最大奇异值的比值界限ε(下面的程序令ε=10-6)，从而设定了最小奇异值wmin=wmax*ε，如果原来矩阵W中有奇异值小于wmin，则令它为零最后调用子程序svbksb用近似的奇异值矩阵回代解线性方程组。

! 例1，2，3

usemsimsl

! parameter(mp=3，np=2)! 例3

! parameter(mp=2，np=3)! 例2

parameter(mp=1，np=2)! 例1

integerm，n，np，mp

reala(mp，np)，X(np)，d(mp)，w(np)，v(np，np)，wmax，wmin

dataa/1，1/，d/2/ ! 例1

! dataa/1，0，1，0，0，1/，d/3，3/ ! 例2

! dataa/1，1，0，1，0，1/，d/3，1，1/! 例3

m=mp

n=np

! 对a进行奇异值分解

callsvdcmp(a，m，n，mp，np，w，v)

! SUBROUTINEsvdcmp(a，m，n，mp，np，w，v)

! INTEGERm，mp，n，np，NMAX

! REALa(mp，np)，v(np，np)，w(np)

! PARAMETER(NMAX=500)

!findmaximumsingularvalue

wmax=0.0

do13k=1，np

if(w(k).gt.wmax)wmax=w(k)

13 continue

! define"small"

wmin=wmax*(1.0e-6)

!zerothe"small"singularvalues

do14k=1，np

if(w(k).lt.wmin)w(k)=0.0

14 continue

write(*，*)'u'

doi=1，m

write(*，*)(a(i，j)，j=1，n)

enddo

! 回代解方程

callsvbksb(a，w，v，m，n，mp，np，d，X)

!这时a就是u

!SUBROUTINEsvbksb(u，w，v，m，n，mp，np，b，x)

!INTEGERm，mp，n，np，NMAX

!REALb(mp)，u(mp，np)，v(np，np)，w(np)，x(np)

!PARAMETER(NMAX=500)

write(*，*)'v'

doi=1，np

write(*，*)(v(i，j)，j=1，n)

enddo

write(*，*)'w'

doi=1，np

write(*，*)w(i)

enddo

write(*，*)'x'

write(*，*)(x(i)，i=1，np)

end

用以上方法解式(4.1)就可以直接获得模型的解。当然一般来说写不出那样的数据方程，而是用泰勒公式近似的，所以用以上方法解第3章的最小二乘法的线性方程组公式(3.14)，获得模型修改量，使反演的迭代过程可以稳定进行下去。

应用奇异值分解法并不要求系数矩阵为方阵，当G的维数M＜N时，方程组Gm=d欠定，奇异值分解法得到的是最小长度解。当G的维数M＞N时，方程组超定，得到最小方差解[1]。

[例1] m1+m2=2

这是一个欠定方程组。未知数个数大于方程个数，所以有无穷多个解，但是用奇异值分解法可以得到唯一解，这个解为最小长度解(解向量的长度最小)。

用上面程序解出m1=m2=0.999999≈1，从图4.1可知，符合[例1]方程的解有无穷多个，但是从原点到方程的直线的距离最小的点只有一个，这个点就是最小长度解。

图4.1 [例1]最小长度解示意图[1]

[例2]

这也是一个欠定方程组。用奇异值分解法也可以得到唯一最小长度解。

[例2]中第一个方程代表图4.2中FDCB所在的平面。第二个方程代表EDF所在的平面。这两个平面的相交线为FD所在的直线。FD直线上所有的点都是[例2]方程组的解。用奇异值分解法可以解出最小长度解(A点到原点的距离最小)。

图4.2 [例2]最小长度解示意图[1]

图4.3 [例3]最小方差解示意图

[例3]

地球物理反演教程

这是一个超定方程组，有2个未知数，3个方程，并且是一个矛盾方程组，因为无法同时满足3个方程，理论上这样的方程应该是无解的，但是奇异值分解法可以得到唯一解，即“最小方差解”。

图4.3所示黑点为奇异值分解法所得的“最小方差解”。它不满足[例3]中任何一个方程，但是它到每个方程所代表的直线的距离的平方和最小。

易读，模块化程序设计利于程序功能分块，加强了程序可阅读性，是良好的程序设计风格

利于排查错误，容易确定错误出现的位置

利于代码复用，相同功能的程序段写入函数，可以简单反复调用

灵活，利用调用函数的参数，实现程序段在不同数据基础上的运算，利用返回值也可以非常简单地获取焦点数据结果

---