数值代数实验报告
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1.谈谈你对该算法的理解：（简单谈一下你是如何理解该算法的？）
先对84阶矩阵进行LU分解，通过Gauss消元法
对下三角形方程组利用前代法解出y，在对上三角方程组
用回代法解出x….
2.实验内容
function [ L,U ] = LUfac( A )
for k=1:n-1
A(k+1:n,k)=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
end
L=tril(A,0);
for i=1:n
L(i,i)=1;
end
U=triu(A,0);
End //进行LU分解
function [ b ] = TSL( L,b )
n=size(L,1);
for j=1:n-1
b(j)=b(j)/L(j,j);
b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);
end
b(n)=b(n)/L(n,n);
end //利用前代法解出y
function [ b ] = TSU( U,b )
n=size(U,1);
for j=n:-1:2
b(j)=b(j)/U(j,j);
b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j);
end
b(1)=b(1)/U(1,1);
end //利用回代法解出x
主函数程序
A=eye(84);
A=6*A;
for i=2:84
A(i,i-1)=8;
A(i-1,i)=1;
End //生成84阶的矩阵A
b=ones(84,1);
b=b*15;
b(1)=7;
b(84)=14;
[L,U]=LUfac(A);//调用函数LUfac对矩阵A进行分解
y=TSL(L,b);//调用函数TSL求解y
x=TSU(U,y); //调用函数TSU求解X
经过matlab…有
x’
ans =
1.0e+008 *
Columns 1 through 7
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 14
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 15 through 21
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 22 through 28
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 29 through 35
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 36 through 42
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 43 through 49
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 50 through 56
0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 Columns 57 through 63
-0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 Columns 64 through 70
0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 Columns 71 through 77
-0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 Columns 78 through 84
0.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684
function [L,U,P]=Lufac(A)
n=size(A,1);P=eye(n,n);
for i=1:n-1
[r,m]=max(abs(A(i:n,i)));
m=m+i-1;
A([i,m],:)=A([m,i],:);
P([i,m])=P([m,i]);
if A(i,i)~=0
A(i+1:n,i)=A(i+1:n,i)/ A(i,i);
A(i+1:n,i+1:n)= A(i+1:n,i+1:n)-A(i+1:n,i)*A(i,i+1:n); end
end
U=triu(A);
L=tril(A,-1)+eye(n,n);
end
A=eye(84);
A=6*A;
for i=2:84
A(i,i-1)=8;
A(i-1,i)=1;
end
b=ones(84,1);
b=b*15;
b(1)=7;
b(84)=14;
[L,U,P]=Lufac(A);
b=P*b;
y=TSL(L,b);
x=TSU(U,y);
结果如下：