%直接三角分解法（1）

function [x,y,L,U]=nalu(a,b) n=length(a);

x=zeros(n,1);y=zeros(n,1); U=zeros(n,n);L=eye(n,n); U(1,:)=a(1,:);

L(2:n,1)=a(2:n,1)/U(1,1); for k=2:n

U(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,k))/U(k,k); end

for i=2:n

y(1,1)=b(1,1);

y(i,1)=b(i,1)-L(i,1:i-1)*y(1:i-1,1); end

y(:,1);

for i=n-1:-1:1

x(n,1)=y(n,1)/U(n,n);

x(i,1)=(y(i,1)-U(i,i+1:n)*x(i+1:n,1))/U(i,i); end

x(:,1);

clear all;

clc;

A=[1,2,3;2,5,2;3,1,5]; b=[14;18;20];

[x,y,L,U]=nalu(A,b);

function [x,y,L,U]=sanjiao(a,b) n=length(a); x=zeros(n,1); y=zeros(n,1); L=eye(n,n); U=zeros(n,n);

%L,U·??a

U(1,:)=a(1,:);

L(2:n,1)=a(2:n,1)/U(1,1); for j=2:n

U(j,j:n)=a(j,j:n)-L(j,1:j-1)*U(1:j-1,j:n);

L(j+1:n,j)=(a(j+1:n,j)-L(j+1:n,1:j-1)*U(1:j-1,j))/U(j,j); end

%?ó?a￡?áíUx=y,Ly=b y(1,1)=b(1,1); for i=2:n

y(i,1)=b(i,1)-L(i,1:i-1)*y(1:i-1,1); end

x(n,1)=y(n,1)/U(n,n); for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1:n)*x(i+1:n,1))/U(i,i); end

clear all;

a=[1,2,3;2,5,2;3,1,5]; b=[14;18;20];

[x,y,L,U]=sanjiao(a,b); x y