五对角矩阵LU分解求解Ax=B

=*

function [x]=FivediagLU(A,a,b,c,C,B)

%A,a,b,c,C依序为5从左到右5个对角线元素向量(其中不足n在对应元素补上0)，B为方程组右端向量

n=length(b);

%初始化：L,l,u,r,R对应A,a,b,c,C；

L=zeros(n,1);l=zeros(n,1);u=zeros(n,1);r=zeros(n,1);R=zeros(n,1);

%LU分解

u(1)=b(1);r(1)=c(1);%第一行

R(1:n-2)=C(1:n-2);

l(2)=a(2)/u(1);L(3)=A(3)/u(1);%第一列

u(2)=b(2)-l(2)*r(1);r(2)=c(2)-l(2)*R(1);%第二行

for k=3:n-1

    if u(k-1)==0

disp('h该矩阵为奇异矩阵！');

        return;

    end

    %计算L中的列

    l(k)=(a(k)-L(k)*r(k-2))/u(k-1);

    L(k+1)=A(k+1)/u(k-1);

    %计算U中的行

    u(k)=b(k)-L(k)*R(k-2)-l(k)*r(k-1);

    r(k)=c(k)-l(k)*R(k-1);  

end

l(n)=(a(n)-L(n)*r(n-2))/u(n-1);%最后一列

u(n)=b(n)-L(n)*R(n-2)-l(n)*r(n-1);%最后一行

%追赶法求解

%Ly=B求y

y(1)=B(1);

y(2)=B(2)-y(1)*l(2);

for k=3:n

    y(k)=B(k)-y(k-1)*l(k)-y(k-2)*L(k);

end

%Ux=y求x

if u(n)==0

disp('该矩阵为奇异矩阵！');

    return;

end

x(n)=y(n)/u(n);

x(n-1)=(y(n-1)-x(n)*r(n-1))/u(n-1);

for k=n-2:-1:1

    x(k)=(y(k)-x(k+1)*r(k)-x(k+2)*R(k))/u(k);

end

end