matlab求解三对角方程组

在求解偏微分方程时，最经典的就是热传导方程的求解。
求解方法普遍采用隐格式。差分之后很显然会得到下图所示的三对角矩阵：
求解三对角矩阵有很多方法，下面只列举两种方法
1、直接求解
行如上面的矩阵，编写代码如下：

n = n0-1;
a2 = sparse(2:n,1:n-1,ones(1,n-1)*(-r),n,n);
a1 = sparse(1:n,1:n,(1+2*r)*ones(1,n),n,n);
a = a1+a2+a2';
for i = 1:n-1
	b(i,1) = u(i,j)
end
b(1) = u(1,j)+rg1;
b(n-1) = u(n-1,j)+rg2;
x = a\b;

2、追赶法

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

为了验证两种方法的效果，编写算例：

%% 稀疏矩阵的算例
n = 10;
a2 =  sparse(2:n,1:n-1,ones(1,n-1),n,n)
a1 = sparse(1:n,1:n,4*ones(1,n),n,n);
a = a1+a2+a2';
b = ones(n,1);
x1 = a\b;

a4 = 4*ones(1,n);
a3 =  ones(1,n-1);
[L,U,x2]=crout(a4,a3,a3',b);

x = x1'-x2

计算结果发现：
x =

1.0e-16 *

     0         0         0         0    0.2776   -0.2776         0         0         0         0

结果几乎一样，所以在实际求解中完全可以采用直接求解法。
PS:我的零食小屋开业啦！有想吃零食的欢迎找我！
有辣辣的，干果还有甜品和水果哦！vx:li1360347cn