LU分解法求线性方程组算法流程图
- 格式:doc
- 大小:28.61 KB
- 文档页数:1
报告内容:用LU分解法求解线性方程组一、报告目的学会用LU分解法解线性方程组,并且为了保证乘子的绝对值小于1,需要对主元数值很小的时候进行方程变换。
二、报告内容1 :测试基本线性方程:可见程序运行结果正确。
2 :测试主元有接近于0的线性方程:可见程序运行结果正确。
3:测试矩阵维度大于5情况:可见程序运行结果正确。
三、源程序function X=LU1(A,B)B=B';A=[A';B]',n=length(B');X=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);endendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=B(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*X(j);endX(i)=(y(i)-ly1)/U(i,i);end四、报告分析与心得MATLAB编程和C编程一样需要仔细的逻辑和十分的细心,自己还有很多需要学习的地方。
LU 和QR 法解线性方程组一、 问题描述求解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡20116126384102785124⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡4321x x x x ==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----3772, 要求:1、编写用三角(LU )分解法求解线性方程组;2、编写用正交三角(QR )分解法求解线性方程组。
二、问题分析求解线性方程组Ax=b ,其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵,从而简化方程组的求解。
因此,在求解线性方程组的过程中,把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要,然而矩阵A 的变换通常有两种分解方法:LU 分解法和QR 分解法。
1、LU 分解法:将A 分解为一个下三角矩阵L 和一个上三角矩阵U ,即:A=LU,其中 L=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡10010012121n n l l l , U=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn n n u u u u u u 0000022211211 2、QR 分解法:将A 分解为一个正交矩阵Q 和一个上三角矩阵R,即:A=QR三、实验原理1、LU 分解法解Ax=b 的问题就等价于要求解两个三角形方程组: ⑴ Ly=b,求y; ⑵ Ux=y,求x.设A 为非奇异矩阵,且有分解式A=LU , L 为单位下三角阵,U 为上三角阵。
L,U 的元素可以有n 步直接计算定出。
用直接三角分解法解Ax=b (要求A 的所有顺序主子式都不为零)的计算公式:① ),,2,1(n i a u li li ==,11/u a l il il = ,i=2,3,…,n. 计算U 的第r 行,L 的第r 列元素(i=2,3,…,n): ② ∑-=-=11r k ki rkri ri u la u , i=r,r+1,…,n;③ rr r k kr ikir ir u u la l /)(11∑-=-= , i=r+1,…,n,且r ≠n.求解Ly=b ,Ux=y 的计算公式;④:,3,2,,1111n i y l b y b y i k k ik i i =-==∑-=⑤.1,,2,1,/)(,/1--=-==∑+=n n i u x uy x u y x ii ni k k iki i nn n n2、QR 分解法四、实验步骤1、LU 分解法1>将矩阵A 保存进计算机中,再定义2个空矩阵L ,U 以便保存求出的三角矩阵的值。
Fortran LU分解法是一种用于解线性方程组的算法,它可以将一个方阵分解为一个下三角矩阵和一个上三角矩阵的乘积。
这种分解方法可以用于高斯消元法、最小二乘法和特征值计算等多种应用中。
在Fortran中,LU分解可以通过调用预先定义好的子程序实现。
下面是一个示例程序,演示如何使用Fortran LU分解法解一个线性方程组:program LUdecompositionimplicit noneinteger, parameter :: n = 3real, dimension(n,n) :: A, L, Ureal, dimension(n) :: b, xinteger :: i, j, k! Input matrix A and vector bA = reshape([2.0, 1.0, 1.0, 6.0, 2.0, 1.0, -2.0, -1.0, 3.0], shape(A))b = [1.0, 2.0, 3.0]! LU decomposition of Acall lu(n, A, L, U)! Solve Ax=b using forward and backward substitutiondo i = 1, nx(i) = b(i) / L(i,i)end dodo i = n, 1, -1do j = i+1, nx(i) = x(i) - U(i,j) * x(j)end dox(i) = x(i) / U(i,i)end do! Output solution xprint *, "Solution: ", xend program LUdecomposition在上面的程序中,我们首先定义了一个3x3的矩阵A和一个3x1的向量b。
然后,我们调用预先定义的子程序lu对矩阵A进行LU分解,得到下三角矩阵L和上三角矩阵U。
接下来,我们使用前向和后向代入法解方程组Ax=b,得到解向量x。
LU和QR分解法解线性方程组LU 和QR 法解线性方程组一、 问题描述求解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡20116126384102785124⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡4321x x x x ==⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----3772,要求:1、编写用三角(LU )分解法求解线性方程组;2、编写用正交三角(QR )分解法求解线性方程组。
二、问题分析求解线性方程组Ax=b ,其实质就是把它的系数矩阵A 通过各种变换成一个下三角或上三角矩阵,从而简化方程组的求解。
因此,在求解线性方程组的过程中,把系数矩阵A 变换成上三角或下三角矩阵显得尤为重要,然而矩阵A矩阵其余的元素,每次都要三次循环,求下一个元素需要上一个结果。
3>先由公式④ ,Ly=b:,3,2,,1111n i y l b y b y i k kiki i =-==∑-= 求出y ,因为L 为下三角矩阵,所以由第一行开始求y.4>再由公式⑤,Ux=y.1,,2,1,/)(,/1--=-==∑+=n n i u x u y x u y x ii ni k k ik i i nn n n求出x, 因为U 为上三角矩阵,所以由最后一行开始求x. 2、QR 分解法四、程序流程图1、LU 分解法开始输入系数矩阵A,常数项b及ndet←1K=1,n-1,1调选列主元子程序i=k+1,n,1aik←aik/akk(即求乘数mik)j=k+1,n,1aik←aij-aik *akjjbi←bi←aik*bkidet←akkdetkbn←bn/am(即求出xn)i=n-1,1,-1s←0j=i+1,n,1s←s+aij*bjjbi←(bi-s)/aij输出b及det结束idet←amndet回代过程消元过程由主程序转来d ←akk,p ←k i=k+1,n,1|aik|>|d|d ←aik,p ←ii d=0P=k j=k,n,1t ←apj,apj ←akj,akj ←t jt ←bp,bp ←bk,bk ←tdet ←det 返回主程序返回主程序输出失败标志,|A|=0det ←0结束选列主元≤=≠=2、QR 分解法开始输入数据A,b矩阵A的转置Householder变换转置矩阵相乘输出上三角阵R=H3*H2*H1*A输入正交阵Q=-(H3*H2*H1)T解上三角阵输出结果结束五、实验结果1、LU分解法2、QR分解法六、实验总结为了求解线性方程组,我们通常需要一定的解法。
三角矩阵的lu分解在线性代数中, LU分解(LU Decomposition)是矩阵分解的一种,可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积(有时是它们和一个置换矩阵的乘积)。
LU分解主要应用在数值分析中,用来解线性方程、求反矩阵或计算行列式。
什么是LU分解如果有一个矩阵A,将A表示成下三角矩阵L和上三角矩阵U的乘积,称为A的LU分解。
更进一步,我们希望下三角矩阵的对角元素都为1:一旦完成了LU分解,解线性方程组就会容易得多。
LU分解的步骤对于满秩矩阵A来说,通过左乘一个消元矩阵,可以得到一个上三角矩阵U。
可以看到,L实际上就是消元矩阵的逆。
容易知道二阶矩阵的逆:现在假设A是一个3×3矩阵,在不考虑行交换的情况下,通过消元得到上三角矩阵的过程是:LU 分解的前提并非所有矩阵都能进行LU分解,能够LU分解的矩阵需要满足以下三个条件:矩阵是方阵(LU分解主要是针对方阵);矩阵是可逆的,也就是该矩阵是满秩矩阵,每一行都是独立向量;消元过程中没有0主元出现,也就是消元过程中不能出现行交换的初等变换。
LU分解的意义LU分解的意义在于求解大型方程组。
一个方程组可以简化为Ax = b的形式,其中A是n阶方阵,x是未知数组成的向量,b是n×1矩阵,例如:以往求解的方式有两种,一是高斯消元法,二是对A求逆,使得x = A-1b。
第二种方式远比消元法复杂,先看一下消元法的计算量。
假设A是n阶满秩方阵,如果不写成增广矩阵,即不考虑 b,那么第一次消元达到的效果是:其中方块是A原来的元素,0是达到的效果,三角是经过消元运算后改变的元素。
以第二行为例,为了使第一个元素为0,需要让第一行乘以某个数(第一行n个元素,共进行了n次乘法运算),再将第一行和第二行相加或相减(第二行n个数与第一行的n个数相加,共进行了n次加法运算)。
如果把一组乘法和加法看成一次运算,那么第二行的消元共进行了n次运算;共有n-1行需要类似运算,所以第一次消元共进行了n(n - 1) ≈ n2次运算。
LU 分解法解线性方程组一. 实验目的1.熟练运用已学计算方法求解方程组2.加深对计算方法技巧,选择正确的计算方法来求解各种方程组3.培养使用电子计算机进行科学计算和解决问题的能力二.实验环境VC++6.0 语言:C++三.实验内容编写一个解线性方程组的L U -算法。
算法包括:1. 矩阵的输入输出2. 判断矩阵是否奇异3. 进行分解4. 输出L 、U5. 回带解出y 、x四.实验原理1.若一个线性方程组系数矩阵为n 阶方阵A 且各阶顺序主子式均不为零,则A 的L U -分解存在且唯一。
2.在满足1的条件下可推导得出以下公式 (1):11()j ij ij ik ij ijk l a l u u -==-∑11j j i j i j kk ik u a lu -==-∑(2): 11i i i ik kk y b l y -==-∑1()ni i ik k iik i x y u x u =+=-∑3.公式(1)用于求解矩阵L 、U ,公式(2)用于回带求解y 、x 。
从公式中可以看出:L 对角线上元素为1,U 第一行与A 第一行相同。
4.L U-分解的具体过程和计算顺序如下:(1)第一步分解:1111u a =(2)第二步分解:212111l a u =1212u a=22222112u a l u =-(3)第三步分解:313111l a =3232311222()l a l u =-1313u a = 23232113u a l u =-333331133223u a l u l u =--(n).第n 步分解:依次计算:1.1..n n n in nnl l u u -⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩1112121222121.1.......1n n n n nn u u u l uu l l u ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⨯⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦五.实验流程图1.程序总的流程图2.考虑主要目的这里只给出L U 分解的流程图:六.VC++源程序#include<iostream>using namespace std;double Fun(int n1,double a1[11][11]); //声明求矩阵行列式的递归函数函数double Fun(int n1,double a1[11][11]){int i_1,j_1,c; //c为数组b的行double b[11][11]; //用于存放余子式int p=0,q=0;int sum=0;if(n1==1)return a1[1][1];for(i_1=1;i_1<=n1;i_1++){for(c=1;c<=n1-1;c++){if(c<i_1)p=0;elsep=1;for(j_1=1;j_1<=n1-1;j_1++){b[c][j_1]=a1[c+p][j_1+1];}}if((i_1-1)%2==0)q=1;elseq=(-1);sum=sum+a1[i_1][1]*q*Fun(n1-1,b);}return sum;}//主程序,用于矩阵输入输出,判断,分解,回带void main(){int n,i,j,k;double temp=1.0;double a[11][11]={0};double L[11][11]={0};double b[11]={1};double x[11]={1};double m=0.0;double y[11]={1};double temp1=0.0;double temp2=0.0;//完成系数矩阵与右端项的输入与输出cout<<"!!!请输入待解方程组未知数的个数: n "<<endl; cin>>n;cout<<"!!!请输入待解方程组的系数"<<endl;for(i=1;i<=n;i++){for(int j=1;j<=n;j++){cin>>a[i][j];}}cout<<"!!!请输入待解方程组的系数矩阵如下:"<<endl;for(i=1;i<=n;i++){for(j=1;j<=n;j++){cout<<a[i][j]<<" ";}cout<<endl;}cout<<"!!!请输入待解方程组的右端项"<<endl;//初始化LUfor(int p=1;p<=n;p++){L[p][p]=1;}for(p=1;p<=n;p++){}cout<<endl;for(i=1;i<=n;i++){cin>>b[i];}cout<<"b="<<endl;for(i=1;i<=n;i++){cout<<b[i]<<endl;}for(i=2;i<=n;i++){temp=Fun(i,a);if(temp==0){cout<<endl<<endl;cout<<endl<<"output:temp ="<<temp<<endl;cout<<"矩阵奇异";break;}}cout<<endl;//LU分解if(temp!=0){for(i=1;i<=n;i++){for(j=1;j<=i-1;j++){if(a[j][j]==0)cout<<"请重新输入系数矩阵与右端项:"<<endl;elsem=0.0;for(k=1;k<=j-1;k++){m=m+L[i][k]*U[k][j];}L[i][j]=(a[i][j]-m)/U[j][j];}for(j=1;j<=i;j++){m=0.0;for(k=1;k<=j-1;k++){m=m+L[j][k]*U[k][i];}U[j][i] = a[j][i] - m;}}}if(temp!=0){cout<<"L="<<endl;for(i=1;i<=n;i++){for(j=1;j<=n;j++){cout<<L[i][j]<<" ";}cout<<endl<<endl;}cout<<endl<<endl<<"U="<<endl; for(i=1;i<=n;i++){for(j=1;j<=n;j++){cout<<U[i][j]<<" ";}cout<<endl<<endl;}//回带求解y[]cout<<"y="<<endl;for(int s=1;s<=n;s++){for(int t=1;t<=s-1;t++){temp1=temp1+L[s][t]*y[t];}y[s]=b[s]-temp1;cout<<" "<<y[s]<<" "<<endl;temp1=0.0;}cout<<endl;//回带求解x[]x[1]=y[1];for(s=n;s>=1;s--){for(int t=s+1;t<=n;t++){temp2=temp2+U[s][t]*x[t];}x[s]=(y[s]-temp2)/U[s][s];temp2=0.0;}cout<<"x="<<endl;for(i=1;i<=n;i++)cout<<" "<<x[i]<<endl;}}七.程序执行结果1.系数矩阵主子式有为0时,随便输入一个矩阵结果如下:2. 系数矩阵各阶主子式不为0时,输入课本70页习题12所示矩阵系数5791068109710895765A⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦与右端项2618229b⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,顺序执行程序,结果经matlab验证无误,结果如下:11 / 12验证:七.实验优缺点此算法可以一定程度上判断是否可解并解出其解,但由于未添加主元选取的步骤,导致对于有的方程组会判断失误并认为方程无解。
数值分析5LU分解法LU分解法是一种常用的数值分析方法,用于解线性方程组。
本文将详细介绍LU分解法的原理、算法步骤、优缺点以及应用领域,以期能够全面地掌握这一方法。
一、LU分解法原理LU分解法是将一个方程组的系数矩阵A分解为两个矩阵L和U的乘积的形式,其中L是下三角矩阵,U是上三角矩阵,通过分解可以简化方程组的求解过程。
LU分解法的基本思想是将原始方程组Ax=b分解为Ly=b和Ux=y两个方程组,其中L和U是通过A分解得到的矩阵。
二、算法步骤1.首先,将系数矩阵A分解为两个矩阵L和U。
L是下三角矩阵,主对角线元素均为1,而U是上三角矩阵。
2.然后,将原始方程组Ax=b转化为Ly=b,求解y的值。
3.最后,将解y代入Ux=y,求解x的值,即可得到方程组的解。
三、算法优缺点1.优点:LU分解法将原始方程组的系数矩阵分解为两个形式简单的矩阵,简化了方程组的求解过程。
对于重复使用系数矩阵A的情况,只需要进行一次LU分解,然后根据新的b值求解新方程组,提高了计算效率。
2.缺点:LU分解法需要进行矩阵分解计算,计算量较大,因此对于规模较大的方程组计算效率较低。
此外,当系数矩阵A存在奇异性或病态时,LU分解法可能会失败。
四、应用领域LU分解法在科学计算领域有着广泛的应用,特别是在求解线性方程组方面。
例如,在工程领域中,常需要通过数值方法求解复杂的结构力学问题,此时可以使用LU分解法求解由有限元方法离散得到的大规模线性方程组。
另外,LU分解法还可以用于解非线性方程组、求逆矩阵、计算矩阵的行列式等。
总结:LU分解法是一种常用的数值分析方法,用于求解线性方程组。
通过将系数矩阵A分解为两个矩阵L和U的乘积形式,可以简化方程组的求解过程。
LU分解法的优点是提高了方程组的求解效率,适用于重复使用系数矩阵A的情况。
然而,LU分解法也存在一定的缺点,如计算量较大、对奇异性和病态问题的处理较为困难。
LU分解法在科学计算领域有着广泛的应用,可以用于求解工程问题中的大规模线性方程组,解非线性方程组,求逆矩阵等。
实验六高斯列主元消去法和直接三角分解法(LU分解)一、实验名称:分别用高斯列主元消去法和直接三角分解法(LU分解)求方程组的解系数矩阵:10 7 8 7 常向量:107 5 6 5 88 6 10 9 67 5 9 10 7精确解为:(-60,102,-27,16)二、试验目的:分别用高斯列主元消去法和直接三角分解法(LU分解)求方程组的解,比较二者不同的特点。
三、算法描述:2、直接三角分解法(LU分解)四、源程序:1、高斯列主元消去法#include <stdio.h>void main(){float a[4][4]={{10,7,8,7},{7,5,6,5},{8,6,10,9},{7,5,9,10}}, y[4],c[4][4],x[4],d[4],m,b; int i,n,j,f;printf("请输入右端项:\n");for(i=0;i<=3;i++)scanf("%f",&y[i]);for(n=0;n<=2;n++){ m=a[n][n]; f=n;for(i=(n+1);i<=3;i++){if(m<a[i][n]){m=a[i][n]; f=i;}}if(f!=n){for(j=0;j<=3;j++){c[n][j]=a[n][j]; a[n][j]=a[f][j];a[f][j]=c[n][j];}d[n]=y[n]; y[n]=y[f]; y[f]=d[n];}for(i=(n+1);i<=3;i++){b=-a[i][n]/a[n][n];for(j=0;j<=3;j++)a[i][j]=a[n][j]*b+a[i][j];y[i]=y[n]*b+y[i];}}x[3]=y[3]/a[3][3];x[2]=(y[2]-a[2][3]*x[3])/a[2][2];x[1]=(y[1]-a[1][3]*x[3]-a[1][2]*x[2])/a[1][1];x[0]=(y[0]-a[0][3]*x[3]-a[0][2]*x[2]-a[0][1]*x[1])/a[0][0];printf("x1的值为%f\nx2的值为%f\nx3的值为%f\nx4的值为%f\n",x[0],x[1],x[2],x[3]);}2、直接三角分解法(LU分解)#include <stdio.h>void main (){float a[4][4]={{10,7,8,7},{7,5,6,5},{8,6,10,9},{7,5,9,10}},y[4], l[4][4],x[4],u[4][4],b[4]; int i,n,j;printf("请输入右端项:");for(i=0;i<=2;i++)scanf ("%f",&b[i]);for (i=0;i<=3;i++){ u[0][i]=a[0][i];l[i][0]=a[i][0]/u[0][0];}for (n=1;n<=2;n++){for(j=0;j<=(3-n);j++){u[n][n+j]=a[n][n+j];for (i=0;i<n;i++)u[n][n+j]=u[n][n+j]-l[n+j][i]*u[i][n]; }if (n=1){l[2][1]=(a[2][1]-l[2][0]*u[0][1])/u[1][1]; l[3][1]=(a[3][1]-l[3][0]*u[0][1])/u[1][1]; } else if (n=2){l[3][2]=(float)(a[3][2]-l[3][0]*u[0][2]-l[3][1]*u[1][2])/u[2][2]; } }for (n=0;n<=3;n++){y[n]=b[n];for (i=0;i<n;i++)y[n]=y[n]-l[n][i]*y[i]; }for (n=3;n>=0;n--){x[n]=y[n];for (i=3;i>n;i--)x[n]=(x[n]-u[n][i]*x[i]);x[n]=x[n]/u[n][n]; }for(i=0;i<=3;i++)printf("x%d的的值为:%f\n",i+1,x[i]); }五、输出结果1、高斯列主元消去法(1)请输入右端项:10 8 6 7X1的值为-59.999584X1的值为101.999306X1的值为-26.999817X1的值为15.999890(2)请输入右端项:5 6 7 8X1的值为-98.999306X1的值为163.998840X1的值为-40.999695X1的值为24.999817六、对算法的理解和改进改变右端项会对结果产生明显的影响,高斯列主元消去法仅考虑依次按列选取主元素,然后按行使之变到主元位置,再进行消去运算,消元结果冲掉A,计算解X冲掉常数项b,则在计算过程中由于右端项的不同解必然不同。
matlab求解线性⽅程组之LU分解线性代数中的⼀个核⼼思想就是矩阵分解,既将⼀个复杂的矩阵分解为更简单的矩阵的乘积。
常见的有如下分解:LU分解:A=LU,A是m×n矩阵,L是m×m下三⾓矩阵,U是m×n阶梯形矩阵QR分解:秩分解:A=CD , A是m×n矩阵,C是m×4矩阵,D是4×n矩阵。
奇异值分解:A=UDV T谱分解:在求解线性⽅程组中,⼀个核⼼的问题就是矩阵的LU分解,我们将⼀个矩阵A分解为两个更加简单的矩阵的复合LU,其中L是下三⾓矩阵,U是阶梯形矩阵。
下三⾓矩阵和上三⾓矩阵具有⾮常良好的性质:Lx=y 或者Ux=y 很容易求解。
问题1.对于任意的矩阵A,是否存在LU分解?定理:如果A⾏等价于阶梯形矩阵U,那么(E n E n-1......E1)A=U,其中的E i,i=1,2,.....,n是⾼斯消去矩阵,他们都是下三⾓矩阵,并且都可逆。
这个定理告诉我们三件事:1.并不是所有的矩阵都有LU分解的。
2.A=LU=(E n E n-1......E1)-1U=(E1-1E2-1.....E n-1)U。
3.这个定理还给出了求解矩阵A-1的⼀种⽅法。
数值算法1.Gauss消去⽤Gauss消去法将矩阵A⾏变换为U:⽤Gauss消去矩阵将A⾏变换为U:数值算法2.Gauss-jardon过程和Gauss-jardon基本⼀致,之不多在选择完最⼤元之后,将其化为1,这样就可以通过乘以⼀个倍数来消去其他⾏了。
选择主元当对某⼀列进⾏Gauss消去时,⼀般都是选择这⼀列中绝对值最⼤的⼀个元素作为主元,当然这会进⾏⾏交换。
其好处有⼀下⼏点:1.在Gauss会代的过程中,不会出现除数为0的情况。
2.减少误差传播,这主要是因为乘数⼩于等于1.(为何乘数⼩于等于1,如果选择这⼀列中绝对值最⼤的⼀个元素作为主元,我们假设这个元素是a,那么乘数等于-b/a,此时|b/a|<=1)。