lup5.m

分块五对角线性方程组的求解

%求解 AX=F
%利用分块五对角矩阵的LU分解
%本程序来自宁静志远网易博客gdgzzch.blog.163.com
%算法来自《求解块五对角方程组的新算法》 卢学飞 徐仲
%此算法非常精妙，只可惜审稿编辑数学素质太差，文章有多处错误。
% A=[a1,c2,e3,0,0,0;b2,a2,c3,e4,0;d3,b3,a3,c4,e5;0,d4,b4,a4,c5,e6]
% L=[s1,l1,i,0,0,0;0,s2,l2,I,0,0;0,0,s3,l3,I,0]
% U=[e1,0,0,0,0;t1,e2,0,0,0;u1,t2,e3,0,0]
function x=luP5(A,B,C,D,E,F,r,n)
S=cell(n,1); %每个cell里边是r*r矩阵
L=cell(n,1);
T=cell(n,1);
U=cell(n,1);
P=cell(n+2,1);
Q=cell(n+2,1);
R=cell(n+2,1);
G=cell(n+2,1);
H=cell(n+2,1);
K=cell(n+2,1);
I=eye(r); %r阶单位阵
O=zeros(r); %r阶零阵
S{1}=I;
E{1}=I;
L{1}=I;
T{1}=I;
S{2}=I;
E{2}=I;

U{1}=A{1}-S{1}*E{1}-L{1}*T{1};
L{2}=(B{2}-S{2}*T{1})/U{1};
T{2}=C{2}-L{1}*E{2};
U{2}=A{2}-S{2}*E{2}-L{2}*T{2};
P{1}=O;
P{2}=O;
Q{1}=-I;
Q{2}=O;
R{1}=O;
R{2}=-I;
T{n+1}=O;
E{n+1}=O;
E{n+2}=O;
for i=3:n
    S{i}=D{i}/U{i-2};
    L{i}=(B{i}-S{i}*T{i-1})/U{i-1};
    T{i}=C{i}-L{i-1}*E{i};
    U{i}=A{i}-S{i}*E{i}-L{i}*T{i};
end

for i=1:n
    P{i+2}=F{i}-S{i}*P{i}-L{i}*P{i+1};
    Q{i+2}=-S{i}*Q{i}-L{i}*Q{i+1}; 
    R{i+2}=-S{i}*R{i}-L{i}*R{i+1};
end

G{n+2}=zeros(r,r);
G{n+1}=zeros(r,r);
H{n+2}=zeros(r,r);
H{n+1}=zeros(r,r);
K{n+2}=zeros(r,r);
K{n+1}=zeros(r,r);
for i=n:-1:1
    Un=inv(U{i}); %矩阵的逆
    G{i}=Un*(P{i+2}-T{i+1}*G{i+1}-E{i+2}*G{i+2});
    H{i}=Un*(Q{i+2}-T{i+1}*H{i+1}-E{i+2}*H{i+2});
    K{i}=Un*(R{i+2}-T{i+1}*K{i+1}-E{i+2}*K{i+2});
end

M=inv(E{1}*H{1}+I);
y2=(T{1}*G{1}+E{2}*G{2})-(T{1}*H{1}+E{2}*H{2})*M*E{1}*G{1};
N=((T{1}*K{1}+E{2}*K{2}+I)-(T{1}*H{1}+E{2}*H{2})*M*E{1}*K{1});
y2=N\y2;
y1=M*(E{1}*(G{1}-K{1}*y2));

x=cell(n,1);
for i=1:n
    x{i}=G{i}-H{i}*y1-K{i}*y2;
end