Doolittle分解法matlab程序的实现
编程实现doolittle分解方法解方程组
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 分解是否存在的优化算法
一种判断 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]。
计算方法3-4
( 2) 找r2，使 a r2 2 max a i 2 ，
2 i n
对调2 r2 行 .
消元：用a 22把a i 2消为0 ( i 3,4, , n) : ai 2 第2行 a 第i行，则 22 ai 2 a ij a 2 j a a ij 22 （i 3,4， , n；j 2,3， , n 1 ）
3 3
Gauss 列主元消去法：
优点 ------ 计算结果更可靠； 缺点 ------ 挑主元花机时更多, 有变动，程序复杂。 次序
x1 ,, xn
用Matlab实现选列主元Gauss消去法解线性方程组 在Matlab程序编辑器中输入：
function x=nagauss2(a,b,flag) %a为系数矩阵；b为右 端列向量；flag若为0，则显示中间过程，否则不显示 if nargin<3,flag=0;end 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>k,t=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);
消 元： 用a11将ai 1 ( i 2,, n)化为零；
ai 1 把 a 第1行，加到第i 行。 11
《数值计算方法》教学中的MATLAB应用的方法研究
一、引言数值计算方法又称数值分析，是研究适合计算机求解的各种数学问题的近似方法及其理论。
它的内容包括函数逼近、数值微分与积分、非线性方程（组）的数值解、数值代数、常微分与偏微分方程数值解等。
这门课程起着承上启下的作用，承上是使线性代数、高等数学中的原理得以应用，方法得以实现，启下是为后续课程中数学问题的建模和求解提供方法，是高等理工科院校的重要基础课程。
如今，数值计算、理论研究及物理实验并列成为当今世界科学活动的三种主要方式。
为众多的科学与工程问题提供计算方法，提高计算的可靠性、有效性和精确性，是《数值计算方法》这门课程研究的主要内容。
在长期的教学实践中体会到在《数值计算方法》课程中做好理论内容的传授和学生实践能力的培养这两个环节非常重要。
如何合理的利用计算机软件进行有效地教与学是值得探讨的一个课题。
本文以具体教学为例，介绍了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分解
U(1,i)=A(1,i);endfor Nhomakorabea=2:n
L(r,1)=A(r,1)/U(1,1);
End
for i=2:n
for j=i:n
z=0;
for r=1:i-1
z=z+L(i,r)*U(r,j);
end
U(i,j)=A(i,j)-z;
end
if abs(U(i,i))<eps
flag='failure'
return;
end
for k=i+1:n
m=0;
for q=1:i-1
m=m+L(k,q)*U(q,i);
end
L(k,i)=(A(k,i)-m)/U(i,i);
end
end
L
U
end
2、在command window窗口中输入以下程序代码指令，回车一下即可。
Doolittle([3 2 2 4;7 3 2 7;1 3 6 9;2 5 9 2])
姓名学号专业班级
课程名称计算矩阵实验名称doolittle分解
实验日期
同组人员指导教师
得分
一、实验目的
1．了解三角分解与doolittle分解的基本概念及其性质；
2．掌握doolittle分解算法及matlab实现；
3．学习、掌握MATLAB软件有关的命令。
二、实验准备
熟悉doolittle分解，了解MATLAB矩阵的输出。
三、实验内容
利用doolittle分解下面矩阵
四、实验步骤
1、在实验时，双击matlab软件，首先建立M文件。具体操作：File→new→M-file→粘贴代码→保存。
matlab doolittle分解法代码
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程序
数值计算⽅法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)的解。
3.2.2 矩阵的doolittle分解
1
... ...
ln3 ... 1
因此再由Ux y 的解便得到 Ax b的解:
xn
yn unn
n
yr urj x j
xr
jr1
urr
u11
u12
u13
...
u22 u23 ...
u1,n
u2,n
U
...
...
u u n1,n1
n1,n
unn
r n 1, n 2,..., 2,1
解
a11
A
a21 a31 a41
a12 a22 a32 a42
a13 a23 a33 a43
a14 a24 a44
a15 a25 a35 a45
2 3 1 4
10 4 2 14
0 12
3 9
3 13 4 13
10 5 2 7
2 10 0 3 10
r
1
3 2 1
4 2
12 3
13 4
5
2
2
2
14
9
13
7
u1 j a1 j
li1
ai1 u11
y1 b1
r 1
yr br lrj y j j 1
r1
urj arj lrkukj k 1
2
l21 l31 l41
u22 l32 l42
u23 a33 a43
u24 a34 a44
y2
b3 b4
u11 u12 u13 u14
r
3
l21 l31 l41
u22 l32 l42
u23 u33 l43
u24 u34 a44
Doolittle分解
矩阵数值分析实验报告专业信息与计算科学班级学号姓名指导教师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。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
function [L,U,Y,X]=qw2014210705(A,b) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%切记！先将该文件重命名为doolittle.m, 然后存入到 bin文件中。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A为输入方阵
%然后在写入函数[L,U]=doolittle(A)即可获得B矩阵的Doolittle分解。
N=size(A);%获取A的行数和列数
if(N(1,1)~=N(1,2))%判断是否为方阵
disp('您输入的矩阵不是方阵，请重新输入');
return;
end
n=N(1,1);%获取A的维数
Y=zeros(n,1);
X=zeros(n,1);
A1=A;
YY=0;%判断是否顺序主子式不等于0，如果存在任何一个顺序主子式为0，则该值变为1. for i=1:n
if(det(A1(1:i,1:i))==0)
YY=1;
end
end
if(YY==1)
disp('您输入的矩阵不能进行Doolittle分解')
return;
end
L=eye(n);
U=zeros(n);
for i=1:n%求解L的第一列，U的第一行
U(1,i)=A1(1,i);
if(i>=2)
L(i,1)=A1(i,1)/U(1,1);
end
end
for k=2:n%以一行一列的方式，计算U和L阵中的其他元素
for j=k:n
s=0;
for r=1:k-1
s=s+L(k,r)*U(r,j);
end
U(k,j)=A1(k,j)-s;
end
for m=k+1:n
s=0;
for r=1:k-1
s=s+L(m,r)*U(r,k);
end
L(m,k)=(A1(m,k)-s)/U(k,k);
end
end
%求解过程
for i=1:n
s=0;
for j=1:n
s=s+Y(j,1)*L(i,j);
end
Y(i,1)=(b(i,1)-s)/L(i,i);
end
nn=n;
for i=1:n
s=0;
for j=1:n
s=s+X(j,1)*U(nn,j);
end
X(nn,1)=(Y(nn,1)-s)/U(nn,nn);
nn=nn-1;
end。