在方程Ax=d中,只需要输入A和d即可,如果不是三对角阵就会报错。
下面是代码
function [x] = chase(A,d)
%求解线性方程组 Ax=d, 其中 A是三对角阵,如果输入的不是三对角矩阵就会报错
%A是系数矩阵
%d是方程组的右端向量
if ~isequal(tril(A,-1)-diag(diag(tril(A,-1),-1),-1),zeros(size(A)))
error('输入矩阵不是一个三角矩阵');
end
if ~isequal(triu(A,1)-diag(diag(triu(A,1),1),1),zeros(size(A)))
error('输入矩阵不是一个三角矩阵');
end
a=[0;diag(tril(A,-1),-1)];%下对角线
b=diag(A);%中对角线
c=[diag(triu(A,1),1);0];%上对角线
l=zeros(size(a,1),1);%求L
u=zeros(size(b,1),1);%求U
n=size(b,1);%矩阵的维度
x=zeros(n,1);
y=zeros(n,1);
u(1)=b(1);
for i=2:n
l(i)=a(i)/u(i-1);
u(i)=b(i)-l(i)*c(i-1);
end
y(1)=d(1);
for i=2:n
y(i)=d(i)-l(i)*y(i-1);
end
x(n)=y(n)/u(n);
for i=n-1:-1:1
x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
end