求解三对角方程组的追赶法(Matlab程序)

来源:互联网 发布:邮件 rss源码 编辑:程序博客网 时间:2024/05/16 17:30

clear all;clc;
fprintf('输入n:(10,20,30)\n');
n=input('');
n
a=zeros(1,n);b=zeros(1,n);c=zeros(1,n-1);
A=zeros(n,n);
for i=1:n
    for j=1:n
        if j==i
            A(i,j)=2;
        elseif (j-i==-1)||(j-i==1)
            A(i,j)=1;
        else
            A(i,j)=0;
        end
    end
end
fprintf('系数矩阵:');A
fprintf('主对角线元素:');b=diag(A)'
fprintf('第1条对角线元素:');c=diag(A,1)'
fprintf('第-1条对角线元素:');a(2:n)=c
d=zeros(1,n);d(1)=1;d(n)=((-1)^(n+1));
fprintf('给定的系数矩阵(转置后):');d
y=zeros(n,1);x=y;
u=zeros(1,n);l=u;%u矩阵U的主对角线元素,l矩阵L的第-1条对角线元素
%---------追过程--------
u(1)=b(1);y(1)=d(1);
for i=2:n
    l(i)=a(i)/u(i-1);
    u(i)=b(i)-l(i)*c(i-1);
    y(i)=d(i)-l(i)*y(i-1);
end
fprintf('中间解向量:');y
fprintf('单位下三角阵L的第-1条对角线元素:');l
fprintf('上三角阵U的主对角线元素:');u  
U=zeros(n);
L=eye(n);
for i=1:n-1
    L(i+1,i)=l(i+1);
end
fprintf('单位下三角阵:');L    %单位下三角阵

for i=1:n-1
    U(i,i)=u(i);
    U(i,i+1)=c(i);
end
U(n,n)=u(n);
fprintf('上三角阵:');U    %上三角阵
%---------赶过程--------
x(n)=y(n)/u(n);
for i=n-1:-1:1
    x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
fprintf('解向量:');x