Doolittle分解法matlab程序的实现

Doolittle分解方法是一种用于解决线性方程组的数值方法，它可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，从而可以方便地求解线性方程组。

在本文中，我们将介绍Doolittle分解方法的原理和实现过程，并用编程语言实现该方法来解方程组。

一、Doolittle分解方法原理1.1 Doolittle分解方法是一种LU分解的特例，它将一个矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。

其中，L 的主对角线元素全为1，U的主对角线以上的元素全为0。

这样的分解可以方便地求解线性方程组Ax=b，其中b是一个已知的列向量。

1.2 Doolittle分解方法的具体实现过程是通过高斯消元法来实现的。

将矩阵A分解为一个下三角矩阵L和一个上三角矩阵U，然后通过回代法求解线性方程组Ax=b。

具体来说，我们首先将矩阵A分解为L 和U，然后用L和U的乘积代替原来的矩阵A，将原来的线性方程组Ax=b变为LUx=b，然后通过两次回代法求解线性方程组Ly=b和Ux=y，最终得到线性方程组的解x。

1.3 Doolittle分解方法的优点是可以方便地求解多个方程组，因为一旦矩阵A被分解为L和U，就可以通过多次回代法来求解不同的线性方程组，而不需要重新分解矩阵A。

1.4 Doolittle分解方法的缺点是需要对原始的矩阵A进行分解，这需要一定的计算量，特别是对于比较大的矩阵来说。

Doolittle分解方法在实际应用中往往需要结合其他数值方法来提高求解线性方程组的效率。

二、Doolittle分解方法的实现过程2.1 我们需要定义一个函数来实现Doolittle分解。

该函数的输入是一个矩阵A，输出是矩阵A的下三角矩阵L和上三角矩阵U。

2.2 接下来，我们需要通过高斯消元法来实现Doolittle分解。

具体来说，我们首先对矩阵A进行行变换和列变换，使得矩阵A的主对角线元素非零，然后逐步消去矩阵A的非主对角线元素，得到下三角矩阵L和上三角矩阵U。

一种判断 LR 分解是否存在的优化算法沈家铭【摘要】使用定理直接判断方阵A是否存在LR分解的计算难度系数为O（n3），文中在此基础上提出了一个改进算法，将计算难度降为 O（n2）。

给出了设计思路及推导方法，并用Matlab实现了两种算法，通过例题验证了计算效率的提升。

%Since the calculating difficulty coefficient of LR decomposition is O (n3 ) , by using the theorem to judge whether phalanx A exists ,here weput forward an improved algorithm to reduce the calculating difficulty coefficient to O(n2 ) .Design and derivation method are given ,and two algorithms are implemented with matlab .The calculation efficiency improvement is verified with examples .【期刊名称】《长春工业大学学报（自然科学版）》【年(卷),期】2014(000)005【总页数】4页(P593-596)【关键词】LR分解;方阵;算法;Matlab【作者】沈家铭【作者单位】四川大学数学学院，四川成都 610065【正文语种】中文【中图分类】O151.21在有应用背景的数学问题中，很多求解问题最终都可归结为线性方程组的求解。

因此，解线性方程组在科学与工程计算中有着十分重要的作用。

对其系数矩阵进行三角分解，是一种解线性方程组的有效方法，LR分解是求三角分解的一个强有力的算法,受到了人们的极大关注。

在不进行行列置换的前提下，为了判断方阵A的LR 分解是否存在，是进行LR分解的前提[1-4]。

一、引言数值计算方法又称数值分析，是研究适合计算机求解的各种数学问题的近似方法及其理论。

它的内容包括函数逼近、数值微分与积分、非线性方程（组）的数值解、数值代数、常微分与偏微分方程数值解等。

这门课程起着承上启下的作用，承上是使线性代数、高等数学中的原理得以应用，方法得以实现，启下是为后续课程中数学问题的建模和求解提供方法，是高等理工科院校的重要基础课程。

如今，数值计算、理论研究及物理实验并列成为当今世界科学活动的三种主要方式。

为众多的科学与工程问题提供计算方法，提高计算的可靠性、有效性和精确性，是《数值计算方法》这门课程研究的主要内容。

在长期的教学实践中体会到在《数值计算方法》课程中做好理论内容的传授和学生实践能力的培养这两个环节非常重要。

如何合理的利用计算机软件进行有效地教与学是值得探讨的一个课题。

本文以具体教学为例，介绍了MATLAB软件在提高《数值计算方法》课堂教学质量中的具体使用。

二、MATLAB软件引入的必要性MATLAB是美国MathWorks公司自上世纪80年代中期推出的数学软件，其优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。

在欧美等高校，MATLAB已经成为线性代数、自动控制理论、数理统计、数值计算方法、动态系统仿真等高级课程的基本教学工具。

以前的《数值计算方法》课程常采用FORTRAN或者C语言进行教学和实验，要求学生既要对算法有充分的了解，又要熟练掌握这两种语言的语法和编程技巧，导致学生和教师将大量的时间和精力都花在烦琐的程序编写以及对各种结果的绘图上，学习效果往往令人不满意。

正如FORTRAN和C等高级语使人们摆脱了需要直接对计算机硬件资源进行操作一样，被称为第四代计算机语言的MATLAB，以其简洁的、更符合人们思维习惯的代码以及强大的绘图能力备受青睐。

《数值计算方法》课程内容多、课时少，如果运用传统教学方法，有些内容得不到细致地讲解，易使学生产生厌学情绪，收不到良好的教学效果。

计算方法上机作业(丁军)部门： xxx时间： xxx整理范文，仅供参考，可下载自行编辑实验一方程求根一、目的和意义通过matlab程序的编写，熟悉使用二分法、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。

巩固课堂和书本上所学的知识、加强实践能力、提高解决实际计算问题的水平、启发创新思想。

b5E2RGbCAP二、计算公式(1>迭代法(2)牛顿法三、结构程序设计代码(1>迭代法k=1。

Tol=0.5*10^(-6>。

%误差限p0=1.0。

%初值设为1N=1000。

%最大迭代步数while k<=Np=((10-p0^3>/4>^(1/2>。

%迭代法迭代公式if abs(p-p0><Tol %当误差满足要求时，停止迭代break。

endk=k+1。

p0=p。

Enddisp(['p=',num2str(p,16>]>。

disp(['k=',num2str(k>]>。

%以16位有效数字输出p，输出整数kp1EanqFDPw(2>牛顿法k=1。

Tol=0.5*10^(-6>。

p0=1.0。

N=1000。

while k<=Np=p0-(p0^3+4*p0^2-10>/(3*p0^2+8*p0>。

%牛顿法迭代法公式if abs(p-p0><Tolbreakendk=k+1。

p0=p。

enddisp(['p=',num2str(p,16>]>。

disp(['k=',num2str(k>]>。

DXDiTa9E3d四、结果及讨论1、结果(1>迭代法(2)牛顿法2、讨论不动点迭代法是求解非线性方程的主要方法，而牛顿法应用广泛，是著名的非线性方程求根法，两种方法有如下不同：RTCrpUDGiT (1>牛顿法求解方程的单根时具有二阶收敛性，从结果可知，牛顿法比不动点迭代法的收敛速度快，能以较少的迭代达到理想的值。

matlab doolittle分解法代码在MATLAB中，Doolittle分解法（也称为LU分解）是一种常用的求解线性方程组的方法。

以下是一个使用MATLAB实现Doolittle分解法的示例代码：```MATLAB生成一个随机矩阵AA = rand(3, 3);初始化L和U矩阵为0L = eye(3);U = eye(3);采用Doolittle分解法分解矩阵Afor i = 1:3for j = 1:3if i == jU(i, j) = A(i, j);elseL(i, j) = A(i, j);endend寻找最大值并交换行max_idx = find(max(A(:, i)));L(:, i) = L(:, i) / A(max_idx, i);U(:, i) = U(:, i) / A(max_idx, i);end打印L和U矩阵disp('L矩阵：');disp(L);disp('U矩阵：');disp(U);验证分解是否正确A_original = A;L_original = eye(3);U_original = eye(3);[L_original, U_original] = doolittle(A_original);打印原始的L和U矩阵disp('原始的L矩阵：');disp(L_original);disp('原始的U矩阵：');disp(U_original);计算分解后的矩阵乘积与原矩阵是否相等result = L_original * U_original;disp('分解后的矩阵：');disp(result);检查分解后的矩阵是否与原矩阵相等isequal(A_original, result)```以上代码首先生成一个3x3的随机矩阵A，然后使用Doolittle分解法对其进行分解，并打印出分解后的L和U矩阵。

数值计算⽅法matlab程序function [x0,k]=bisect1(fun1,a,b,ep)if nargin<4ep=1e-5;endfa=feval(fun1,a);fb=feval(fun1,b);if fa*fb>0x0=[fa,fb];k=0;return;endk=1;while abs(b-a)/2>epx=(a+b)/2;fx=feval(fun1,x);if fx*fa<0b=x;fb=fx;elsea=x;fa=fx;k=k+1;endendx0=(a+b)/2;>> fun1=inline('x^3-x-1');>> [x0,k]=bisect1(fun1,1.3,1.4,1e-4)x0 =1.3247k =7>>N=500;endif nargin<3ep=1e-5;endx=x0;x0=x+2*ep;while abs(x-x0)>ep & kx0=x;x=feval(fun1,x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fun1=inline('(x+1)^(1/3)');>> [x0,k]=iterate1(fun1,1.5)x0 =1.3247k =7>> fun1=inline('x^3-1');>> [x0,k]=iterate1(fun1,1.5)x0 =Infk =9>>Steffesen加速迭代（简单迭代法的加速）function [x0,k]=steffesen1(fun1,x0,ep,N) if nargin<4N=500;endx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & kx0=x;y=feval(fun1,x0);z=feval(fun1,y);x=x0-(y-x0)^2/(z-2*y+x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fun1=inline('(x+1)^(1/3)');>> [x0,k]=steffesen1(fun1,1.5)x0 =1.3247k =3>> fun1=inline('x^3-1');>> [x0,k]=steffesen1(fun1,1.5)x0 =1.3247k =6Newton迭代function [x0,k]=Newton7(fname,dfname,x0,ep,N) if nargin<5N=500;endendx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & kx0=x;x=x0-feval(fname,x0)/feval(dfname,x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fname=inline('x-cos(x)');>> dfname=inline('1+sin(x)');>> [x0,k]=Newton7(fname,dfname,pi/4,1e-8) x0 =0.7391k =4⾮线性⽅程求根的Matlab函数调⽤举例：1.求多项式的根：求f(x)=x^3-x-1=0的根：>> roots([1 0 -1 -1])ans =1.3247-0.6624 + 0.5623i-0.6624 - 0.5623i2.求⼀般函数的根>> fun=inline('x*sin(x^2-x-1)','x')fun =Inline function:fun(x) = x*sin(x^2-x-1)>> fplot(fun,[-2 0.1]);grid on-1.5956>> x=fzero(fun,[-1 -0.1])x =-0.6180[x,f,h]=fsolve(fun,-1.6)x =-1.5956f =1.4909e-009h =1（h>0表⽰收敛，h<0表⽰发散，h=0表⽰已达到设定的计算函数值的最⼤次数）第三章：线性⽅程组的数值解法1. ⾼斯消元法function [A,x]=gauss3(A,b)%本算法⽤顺序⾼斯消元法求解线性⽅程组n=length(b);A=[A,b];for k=1:n-1A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1)); A((k+1):n,k)=zeros(n-k,1);A;endx=zeros(n,1);%上⾯为消元过程x(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1x(k)=(A(k,n+1)-A(k,(k+1):n)*x((k+1:n)))/A(k,k);end%上⾯为回代过程>> A=[2 3 4;3 5 2;4 3 30];>> b=[6,5,32]'b =>> [A,x]=gauss3(A,b)A =2.00003.00004.0000 6.00000 0.5000 -4.0000 -4.00000 0 -2.0000 -4.0000x =-1382列选主元的⾼斯消元法：function [A,x]=gauss5(A,b)%本算法⽤列选主元的⾼斯消元法求解线性⽅程组n=length(b);A=[A,b];for k=1:n-1%选主元[ap,p]=max(abs(A(k:n,k)));p=p+k-1;if p>kt=A(k,:);A(k,:)=A(p,:);A(p,:)=t;end%消元A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1)); A((k+1):n,k)=zeros(n-k,1);end%回代x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);>> A=[2 3 4;3 5 2;4 3 30]; b=[6,5,32]';>> [A,x]=gauss5(A,b)A =4.0000 3.0000 30.0000 32.00000 2.7500 -20.5000 -19.00000 0 0.1818 0.3636x =-1382三⾓分解法：Doolittle 分解function [L,U]=doolittle1(A)n=length(A);U=zeros(n);L=eye(n);U(1,:)=A(1,:);L(2:n,1)=A(2:n,1)/U(1,1);for k=2:nU(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);L(k+1:n,k)=A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,n)/U(k,k); End y=zeros(n,1);x=y;y(1)=b(1);for i=2:ny(i)=b(i)-L(i,1:i-1)*y(1:i-1);endx(n)=y(n)/U(n,n);for i=n-1:-1:1x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);end>> A=[1 2 3;2 5 2 ;3 1 5];b=[14 18 20]';>> [L,U,x]=doolittle1(A,b)3 -8 1U =1 2 30 1 -40 0 -36x =2.83331.33332.8333平⽅根法：function [L,x]=choesky3(A,b)n=length(A);L=zeros(n);L(:,1)=A(:,1)/sqrt(A(1,1));for k=2:nL(k,k)=A(k,k)-L(k,1:k-1)*L(k,1:k-1)';L(k,k)=sqrt(L(k,k));for i=k+1:nL(i,k)=(A(i,k)-L(i,1:k-1)*L(k,1:k-1)')/L(k,k); endendy=zeros(n,1);x=y;y(1)=b(1)/L(1,1);for i=2:ny(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i);endx(n)=y(n)/L(n,n);for i=n-1:-1:1x(i)=(y(i)-L(i+1:n,i)'*x(i+1:n))/L(i,i);end>> A=[4 -1 1;-1 4.25 2.75;1 2.75 3.5]-1.0000 4.2500 2.75001.00002.75003.5000>> b=[4 6 7.25]'b =4.00006.00007.2500[L,x]=choesky3(A,b)L =2.0000 0 0-0.5000 2.0000 00.5000 1.5000 1.0000x =111>>迭代法求⽅程组的解Jacobi迭代法：function [x,k]=jacobi2(a,b,x0,ep,N)%本算法⽤Jacobi迭代求解ax=b,⽤分量形式n=length(b); k=0;if nargin<5N=500;endif nargin<4ep=1e-5;endif nargin<3x0=zeros(n,1);y=zeros(n,1);while norm(x-x0,inf)>ep & kk=k+1;x0=x;for i=1:ny(i)=b(i);for j=1:nif j~=iy(i)=y(i)-a(i,j)*x0(j);endendif abs(a(i,i))<1e-10|k==Nwarning('a(i,i) is too small');returnendy(i)=y(i)/a(i,i);endx=y;enda=[4 3 0;3 4 -1; 0 -1 4];b=[24 30 -24]'；[x,k]=jacobi2(a,b)x =3.00004.0000-5.0000k =59Gauss-seidel迭代法：function [x,k]=gaussseide2(a,b,x0,ep,N)%本算法⽤Gauss-seidel迭代求解ax=b,⽤分量形式n=length(b); k=0;if nargin<5N=500;endendif nargin<3x0=zeros(n,1);y=zeros(n,1);endx=x0;x0=x+2*ep;while norm(x-x0,inf)>ep & kk=k+1;x0=x;y=x;for i=1:nz(i)=b(i);for j=1:nif j~=iz(i)=z(i)-a(i,j)*x(j);endendif abs(a(i,i))<1e-10|k==Nwarning('a(i,i) is too small');returnendz(i)=z(i)/a(i,i);x(i)=z(i);endend[x,k]=gaussseide2(a,b)x =3.00004.0000-5.0000k =25最速下降法function [x,k]=zuisuxiajiang(A,b,x0,ep,N)N=500;endif nargin<4ep=1e-8;endif nargin<3x0=ones(n,1);endx=x0;x0=x+2*ep;r=b-A*x;d=r;k=0;while norm(x-x0,inf)>ep & kk=k+1;x0=x;lamda=(d'*d)/(d'*A*d);x=x0+lamda*d;r=b-A*x;d=r;endif k==Nwarning('已达最⼤迭代次数')end共轭梯度算法function [x,k]=gongertidufa(A,b,x0,ep,N) %本算法⽤共轭梯度算法求解正定⽅程组Ax=b,,n=length(b);if nargin<5N=500;endif nargin<4ep=1e-8;x0=x+2*ep;r=b-A*x;d=r;k=0;while norm(x-x0,inf)>ep & kx0=x;lamda=(r'*r)/(d'*A*d);r1=r;x=x0+lamda*d;r=b-A*x;beta=(r'*r)/(r1'*r1);d=r+beta*d;endif k==Nwarning('已达最⼤迭代次数') end常微分⽅程数值解function [x,y]=Euler1(fun,xspan,y0,h)%本算法⽤欧拉格式计算微分⽅程y'=f(x,y)的解。

矩阵数值分析实验报告专业信息与计算科学班级学号姓名指导教师Doolittle 分解法一、实验目的在Gauss 消元法中,对于n 阶方程组,应用消去发经过n-1步消元之后,得到一个与Ax=b 等价的代数线性方程组)1()1(--=n n b x A ,而且)1(-n A 为一个上三角矩阵.所以我们想是否能把矩阵A 分解成一个下三角阵与一个上三角阵的乘积 A=LR,其中L 为下三角阵,R 为上三角阵.就变成了两个三角形方程组⎩⎨⎧==yRx b Ly , 的求解问题。

二、算法思想Setp1：利用for 循环求出r[k][j]=a[k][j]-∑-=1k 1p ]kp [r ]kp [l ，l[ik]=（a[ik]-∑-=1k 1p ]ip [r ]ip [l ）/r[k][k]。

Step2：y[i]=b[i]-∑-=1i 1k ]k [y ]ik [l ,得出x[i]=(y[i]-∑+=n1i k ]k [x ]ik [r )/r[i][i].三、程序代码#include <stdio.h>#include <stdlib.h>#define N 10 //矩阵大小范围float getmx(float a[N][N], float x[N], int i, int n) {float mx = 0;int r;for(r=i+1; r<n; r++){mx += a[i][r] * x[r];}return mx;}float getmy(float a[N][N], float y[N], int i, int n) {float my = 0;int r;for(r=0; r<n; r++){if(i != r) my += a[i][r] * y[r];}return my;}float getx(float a[N][N], float b[N], float x[N], int i, int n) {float result;if(i==n-1) //计算最后一个x的值result = (float)(b[i]/a[n-1][n-1]);else //计算其他x值（对于公式中的求和部分，需要调用getmx()函数） result = (float)((b[i]-getmx(a,x,i,n))/a[i][i]);return result;}float gety(float a[N][N], float b[N], float y[N], int i, int n) {float result;if(i==0) //计算第一个y的值result = float(b[i]/a[i][i]);else //计算其他y值（对于公式中的求和部分，需要调用getmy()函数） result = float((b[i]-getmy(a,y,i,n))/a[i][i]);return result;}void main(){float l[N][N]={0}; //定义L矩阵float u[N][N]={0}; //定义U矩阵float y[N]={0}; //定义数组Yfloat x[N]={0}; //定义数组Xfloat a[N][N]={{0},{0},{0}}; //定义系数矩阵float b[N]={0}; //定义右端项float sum=0;int i,j,k;int n=0;int flag=1;//用户手工输入矩阵while(flag){printf("请输入系数矩阵的大小：");scanf("%d", &n);if(n>N){printf("矩阵过大！\n");continue;}flag=0;}printf("请输入系数矩阵值：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("a[%d][%d]: ", i, j);scanf("%f", &a[i][j]);}}printf("请输入右端项数组：\n");for(i=0; i<n; i++){printf("b[%d]: ", i);scanf("%f", &b[i]);}/*显示原始矩阵*/printf("原始矩阵：\n");for(i=0; i<n; i++){for(j=0; j<n; j++)printf("%0.3f ",a[i][j]);printf("\n");}printf("\n\n");/*初始化矩阵l*/for(i=0; i<n; i++){for(j=0; j<n; j++){if(i==j) l[i][j] = 1;}}/*开始LU分解*//*第一步：对矩阵U的首行进行计算*/for(i=0; i<n; i++){u[0][i] = (float)(a[0][i]/l[0][0]); }/*第二步：逐步进行LU分解*/for(i=0; i<n; i++){/*对“L列”进行计算*/for(j=i+1; j<n; j++){for(k=0,sum=0; k<n; k++){if(k != i) sum += l[j][k]*u[k][i];}l[j][i] = (float)((a[j][i]-sum)/u[i][i]); }/*对“U行”进行计算*/for(j=i+1; j<n; j++){for(k=0,sum=0; k<n; k++){if(k != i+1) sum += l[i+1][k]*u[k][j]; }u[i+1][j] = (float)((a[i+1][j]-sum));}}/*输出矩阵l*/printf("矩阵L：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%0.3f ", l[i][j]);}printf("\n");}/*输出矩阵u*/printf("\n矩阵U：\n");for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%0.3f ", u[i][j]);}printf("\n");}/*回代方式计算数组Y*/for(i=0; i<n; i++){y[i] = gety(l,b,y,i,n);}/*显示数组Y*/printf("\n\n数组Y：\n");for(i=0; i<n; i++){printf("y%d = %0.3f\n", i+1,y[i]); }/*回代方式计算数组X*/for(i=n-1; i>=0; i--){x[i] = getx(u,y,x,i,n);}/*显示数组X*/printf("\n\n数组X：\n");for(i=0; i<n; i++){printf("x%d = %0.3f\n", i+1,x[i]); }}四、运行结果五、参考文献[1]刑志栋，矩阵数值分析，陕西：陕西科学技术出版社， 2005。