LU分解
矩阵分解
练习：
求矩阵??
???????
???-=1,1,2,42,1,4,26,12,3,32,4,4,2A 的LU 分解。
代码：
#include
#include
int lluu(n,a,l,u)
int n;
double a[],l[],u[];
{
int i,j,k,w,v ,ll,m,p;
for (k=0; k<=n-2; k++)
{
ll=k*n+k;
if (fabs(a[ll])+1.0==1.0)
{
printf("fail\n");
return(0);
}
for (i=k+1; i<=n-1; i++)
{
w=i*n+k;
a[w]=a[w]/a[ll];
}
for (i=k+1; i<=n-1; i++)
{
w=i*n+k;
for (j=k+1; j<=n-1; j++)
{
v=i*n+j;
a[v]=a[v]-a[w]*a[k*n+j];
}
}
printf("\tk=%d\n",k);
printf("\t%f\t%f\t%f\t%f\n",a[0],a[1],a[2],a[3]);
printf("\t%f\t%f\t%f\t%f\n",a[4],a[5],a[6],a[7]);
printf("\t%f\t%f\t%f\t%f\n",a[8],a[9],a[10],a[11]);
printf("\t%f\t%f\t%f\t%f\n",a[12],a[13],a[14],a[15]);
}
for (i=0; i<=n-1; i++)
{
for (j=0; j { w=i*n+j; l[w]=a[w]; u[w]=0.0; } w=i*n+i; l[w]=1.0; u[w]=a[w]; for (j=i+1; j<=n-1; j++) { w=i*n+j; l[w]=0.0; u[w]=a[w]; } } return(1); } //包含在文件LLUU.C中 #include main() { int i,j; double l[4][4],u[4][4]; double a[4][4]={ {2.0,4.0,4.0,2.0}, {3.0,3.0,12.0,6.0}, {2.0,4.0,-1.0,2.0}, {4.0,2.0,1.0,1.0}}; i=lluu(4,&a[0][0],&l[0][0],&u[0][0]); if (i!=0) { printf("MAT L IS:\n"); for (i=0; i<=3; i++) { for (j=0; j<=3; j++) printf("%13.7f ",l[i][j]); printf("\n"); } printf("\n"); printf("MAT U IS:\n"); for (i=0; i<=3; i++) { for (j=0; j<=3; j++) printf("%13.7f ",u[i][j]); printf("\n"); } printf("\n"); } } 运行结果： 研究生课程设计 (论文) 电力系统稳态分析Stability Analysis of Power Systems 教学单位自动化学院 姓名薛媛媛 学号 111101130 年级2011 级 专业电力系统及其自动化 指导教师杨伟 职称副教授 2011年12月6日 摘要：快速分解法是一种定雅克比法，形成系数矩阵'B 、''B 时忽略了支路电阻、对地导纳和理想变压器非标准变化，以及θP -迭代过程中节点电压的不同取值，从而避免每次迭代重新形成雅克比矩阵及其因子表，计算效率大幅提高。本文采用快速分解法中的XB 型算法，并基于MATLAB 软件仿真分析具体实例，发现具有较好的收敛性，特别适合在线计算。 关键词：电力系统 潮流计算 快速分解法 雅可比矩阵 0 引言 用牛顿-拉夫逊法计算潮流时，每次迭代都需要重新形成雅克比矩阵并分解因子表。为避免每次迭代重新形成雅克比矩阵及其因子表，20世纪70年代初Stott 提出了快速分解法。 快速分解法采用了一些假设： a 、电力系统有功功率主要受电压相角影响，无功功率主要受电压幅值影响。 b 、高压网线路的r< 潮流计算的快速分解法 摘要：本文采用快速分解法进行潮流计算，分析其基本理论，并使用MATLAB软件进行编程设计。最后运用实例进行验证。结果表明快速分解法具有较好的迭代速度。 关键词：潮流计算快速分解法 MATLAB编程，实例验证 1引言 潮流计算是电力系统分析最基本、最重要的计算，是电力系统运行、规划以及安全性、可靠性分析和优化的基础，也是各种电磁暂态和机电暂态分析的基础和出发点。潮流计算要求具有可靠的收敛性，占用内存少，计算速度快，调整和修改容易，使用灵活方便。各种算法的改进以及新算法的提出，很多都是为了使潮流计算能更好地满足计算要求。本文应用快速分解法进行潮流计算，并给出算例分析。 2潮流计算的快速分解法 研究表明，用牛顿-拉夫逊法计算潮流时，每次迭代都要重新形成雅可比矩阵，然后重新对它进行因子表分解并求解修正方程。为避免每次迭代重新形成雅可比矩阵及其因子表，人们研究用定雅可比矩阵取代随迭代过程不断变化的雅可比矩阵，这种方法叫定雅可比法。此外，人们还结合电力系统的物理特点，发展了各种版本的解耦潮流算法，20世纪70年代初提出的快速分解法是这一阶段的主要研究成果。 关于快速分解潮流算法，有三项里程碑意义的研究成果。其一是Stott在1974年发现的XB型算法；其二是Van Amerongen在1989年发现的BX型算法；其三是Monticelli等人在1990年所作的关于快速分解潮流算法收敛机理的理论阐述。这些研究工作不仅是电力系统计算方面的典范，也揭示了这样一个事实：工程上有效的方法一定有其深刻的理论来支持。 2.1 快速分解法的修正方程及迭代格式 将极坐标型定雅可比法的修正公式重写如下： §3 LU 分解法 ——Gauss 消去法的变形 知识预备： 1矩阵的初等行变换、初等矩阵及其逆、乘积 2矩阵的乘法 3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆 一、Gauss 消去法的矩阵解释 Gauss 消去法实质上是将矩阵A 分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。 第一轮消元：相当于对A （1） 左乘矩阵L 1，即 )2()1(1A A L = 其中 ) 1(11 ) 1(1 1)2()2(2 )2(2)2(22 ) 1(1) 1(12)1(11) 2(131211,,1001 011a a l a a a a a a a A l l l L i i nn n n n n =?????? ??????? ?=????? ?? ?????????---= 第二轮消元：对应于 )3()2(2A A L = 一般地 1,,2,1) 1()(-==+n k A A L k k k (1) 其中 n k k i a a l l l L k kk k ik ik nk k k k ,,2,1,,10111)() (1 ++==??? ???? ?? ? ?? ??? ???? ?--=+整个消元过程为 U A A L L L L n n n 记)(1221=-- ????? ? ? ?? ?? ? ? ?=nn n n u u u u u u 22211211………（2） 从而 U L U L L L L U L L L L A n n n n ?===---------1 112121111221)( 其中L 是单位下三角矩阵，即 ,1,,1,,3,2,,11 11)() (21323121??? ? ??-===????? ?? ?????????=n j n i l l l l l l L j jj j ij a a ij n n …（3） 【注】消元过程等价于A 分解成LU 的过程 回代过程是解上三角方程组的过程。 二、矩阵的三角分解 1、若将A 分解成L ?U ，即A=L ?U ，其中L 为单位下三角矩阵,U 为非奇异上三角矩阵,则称之为对A 的Doolittle 分解。 当A 的顺序主子式都不为零时，消元运算可进行，从而A 存在唯一的Doolittle 分解。 证明：若有两种分解，A=L 1U 1，A=L 2U 2，则必有L 1=L 2，U 1=U 2。 因为L 1U 1=L 2U 2，而且L 1，L 2都是单位下三角矩阵，U 1,U 2都是可逆上三角矩阵,所以有 112112--=U U L L 例：给定一4阶矩阵42158 7210 48366 8 4 9A ????? ?=?????? ，通过LU 分解求逆矩阵1A -。 解：算法过程为：1 1 1 1 ()A L U U L u l ----=?=?=?， 第一步：求LU 矩阵 设 0000 0102031011111213 202122222330 31 32 33330000000000 L U U U U L L U U U L U L L L U U L L L L U → ???? ?????????=↓ ????????????? ， 通过（4）~（7）式可逐步进行矩阵L 和U 中元素的计算，如下所示： 00112233000001010202030310203010203000 00 00 111110011212100213131003() 1,(U ) 4,2,1,5,(L )8462,1, 1.5, 4 4 4(U ) 7223,2210,10L L L L L U a U a U a U a a a a L L L U U U U a L U U a L U U a L U ============= === === ===-=-?==-=-?==-=计算的对角的第一行的第一列的第二行250, -?= 21212001113131300111 22222002211223232003211332323002311222 (L )11()(812)23115()(8 1.52)3 3 U 311202,615201,L 115()(4 1.512 3 L a L U U L a L U U U a L U L U U a L U L U L a L U L U U =-=?-?== -= ?-?= =--=-?-?==--=-?-?== --= ?-?- 的第二列（的第三行） （的第三列）33333003311332230) 1.25, U 59 1.550 1.2510.25; 3 U a L U L U L U ?==---=-?- ?-?=（的第四行） 经迭代计算，最后得到L 和U 矩阵为： 第二步：求L 和U 矩阵的逆u ，l 1 1 ,;u U l L --== （1）求U 矩阵的逆 1 1 0001020300010203111213 111213 2223222333334215000300000000210 0.25U U U U u u u u U U U u u u u U U u u U u --?? ????????????????===↑?? ???????????? ?? ?? 由式（9）可得矩阵U 的逆的各元素计算如下： 矩阵的LU 分解 一、题目 求一个4阶矩阵的LU 分解。（????? ???????=112 92 74 48 78 66 54 36 46 40 34 24 18 16 14 12 A ） 二、方法 Doolittle （杜里特尔）分解法 三、程序 jiangLU.M 的程序如下： function[L,U,flag]=jiangLU(A) [n,n]=size(A); L=eye(n);U=zeros(n);flag='OK'; for k=1:n for j=k:n z=0; for q=1:k-1 z=z+L(k,q)*U(q,j); end U(k,j)=A(k,j)-z; end if abs(U(k,k)) 48 74 92 112 >> [L,U,flag]=jiangLU(A) L = 1 0 0 0 2 1 0 0 3 2 1 0 4 3 2 1 U = 12 14 16 18 0 6 8 10 0 0 2 4 0 0 0 2 flag = OK 五、拓展 1、矩阵分解成LU形式是有条件的，首先矩阵必须是非奇异的矩阵，其次矩阵的全部顺序主子式非零的时候才能完全保证矩阵可分解成LU且分解唯一。不然，会有如下情况出现：>> A=[1 2 3;2 4 1;4 6 7] A = 1 2 3 2 4 1 4 6 7 >> [L,U,flag]=jiangLU(A) L = 1 0 0 2 1 0 4 0 1 U = 1 2 3 0 0 -5 0 0 0 flag = failure 2、Doolittle分解法是不选主元的三角分解法，若U（k,k）绝对值很小时，按此分解法可能引起舍入误差的积累。可以通过交换矩阵A的行实现矩阵PA的LU分解（前提是A为非奇异的矩阵）。 实验二：插值法 专业班级：信息与计算科学112班 姓名：马晶 学号：2011014856 （1）实验目的： a.熟悉MA TLAB 编程，学会用MATLAB 实现用插值法求多项式；掌握用画图命令显示所得的多项式的图形，根据所得图形分析问题 b. 学习插值方法及程序设计算法，体会多项式插值法和三次样条插值法的不同，比较两种插值哪个更精确？ （2）实验题目： 1. 已知函数在下列各点的值为 x i 0.2 0.4 0.6 0.8 1.0 F(x i ) 0.98 0.92 0.81 0.64 0.38 试用4次牛顿插值多项式P 4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值.用图给出{( x i,y i ), x i =0.2+0.08i,i=0,1,11,10}, P 4(x)及S(x). 2.在区间[-1,1]上分别取n=10﹑20用两组等距节点对龙格函数2 1()125f x x =+作多项式插值及三 次样条插值,分别画出插值函数及()f x 的图形. 3．下列数据点的插值 x 0 1 4 9 16 2 5 3 6 49 64 y 0 1 2 3 4 5 6 7 8 可以得到平方根函数的近似,在区间[0,64]上作图. (1)用这9个点作8次多项式插值 () 0010012010101()()[,]()[,,]()()[,,,]()n n n P x f x f x x x x f x x x x x x x f x x x x x x x -=+-+--+--L L L (2)用三次样条(第一边界条件)程序求().S x 从得到结果看在[0,64]上,哪个插值更精确;在区间[0,1]上,两种插值哪个更精确? （3）试验理论与理论基础： 1.多项式插值: 设在区间[a,b]上给定n+1个点 01n a x x x b ≤<<<≤L 上的函数值()(0,1,,),i i y f x i n ==L 求次数不超过n 的多项式,使(),0,1,,.i i y P x i n ==L 然后得到关于系数01,,,n a a a L 的n+1元线性方程组 010*********,, n n n n n n n n n a a x a x y a a x a x y a a x a x y ?+++=?+++=???? +++=?L L M L 该方程的系数矩阵为范德蒙德矩阵,通过带入数据可以直接得到唯一解,,,a a a L ,故可求得插值多项式 数值分析实验报告 （1） 学院：信息学院 班级：计算机0903班 姓名：王明强 学号：20092954 课题一 A.问题提出 给定下列几个不同类型的线性方程组，请用适当的方法求解 线性方程组1、设线性方程组 ?? ??????????????????? ??? ??? ?????--------------------------1368 2438 100 41202 9 13 726422123417911101610352431205362177586832337616244911315120130123122400 1 0563568000 0121 324?? ??? ???? ???? ???????????????????1098765432 1x x x x x x x x x x =???? ??????? ????? ????????????????-2119381346323125 x *= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 )T 2、设对称正定阵系数阵线方程组 ??????????????????????????----------------------1924336002141103520411144334 3104221812334161206538114140231212200420424?????????????? ????????????87654321x x x x x x x x = ??? ??? ??? ?????? ??? ?? ? ?????---4515229232060 x * = ( 1, -1, 0, 2, 1, -1, 0, 2 )T 3、三对角形线性方程组 (1)列主元素消去法求解线性方程： #include潮流计算的快速分解法程序
潮流计算的快速分解法
数值分析5 LU分解法
LU矩阵分解实例
矩阵的LU分解
LU分解和列主元消去法求解
数值分析高斯顺序消去法、列主元消去法LU分解法
列主元消去法和LU分解法(C语言)