matlab 追赶法

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 
--------------------- 
作者：zhangchao3322218 
来源：CSDN 
原文：https://blog.csdn.net/zhangchao3322218/article/details/7412681 
版权声明：本文为博主原创文章，转载请附上博文链接！

%追赶法求解三对角线性方程组，Ax=b，A用一维数组a，c，d存储。
function [L,U,x]=crout(a,c,d,b)%数组a存储三角矩阵A的主对角线元素，c、d存储主对角线上边下边带宽为1的元素
    n=length(a);
    n1=length(c);
    n2=length(d);
    %错误检查
    if n1~=n2%存储矩阵的数组维数错误
        error('MATLAB:Crout:不是三对角矩阵，参数数组中元素个数错误.');
    elseif n~=n1+1
        error('MATLAB:Crout:不是三对角矩阵，参数数组中元素个数错误.');
    end
   
    %初始化
    L=zeros(n);%生成n*n的全零矩阵
    U=zeros(n);
    p=1:n;
    q=1:n-1;
    x=1:n;
    y=1:n;
   
    %追赶法程序主体
    p(1)=a(1);
    for i=1:n-1
        q(i)=c(i)/p(i);
        p(i+1)=a(i+1)-d(i)*q(i);%d的下标改为1到n-1
    end
    %正解y
    y(1)=b(1)/p(1);%用x存储y
    for i=2:n
        y(i)=(b(i)-d(i-1)*y(i-1))/p(i);
    end
    %倒解x
    x(n)=y(n);
    for i=(n-1):-1:1
        x(i)=y(i)-q(i)*x(i+1);
    end
    %L,U矩阵
    for i=1:n
        L(i,i)=p(i);
        U(i,i)=1;
    end
    for i=1:n-1
        L(i+1,i)=d(i);
        U(i,i+1)=q(i);
end %end of function
--------------------- 
作者：love密密 
来源：CSDN 
原文：https://blog.csdn.net/jingmiaa/article/details/49660795 
版权声明：本文为博主原创文章，转载请附上博文链接！