数值计算-线性方程组求解（2）-追赶法解三对角矩阵-MATLAB实现

问题描述1

  在用差分法求解二阶常微分方程的边值问题、热传导问题及三次样条插值函数的求解问题中，都会遇到下面形式的阶数较高的三对角方程组AX=f，即

⎡⎣⎢⎢⎢⎢⎢⎢⎢b1a2⋱c1b2an−1⋱c2bn−1⋱ancn−1bn⎤⎦⎥⎥⎥⎥⎥⎥⎥×⎡⎣⎢⎢⎢⎢⎢⎢⎢x1x2⋮xn−1xn⎤⎦⎥⎥⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢f1f2⋮fn−1fn⎤⎦⎥⎥⎥⎥⎥⎥⎥

其中系数矩阵A的元素满足条件(主对角占优条件)

⎧⎩⎨⎪⎪|b1|>|c1|>0,|bi|≥|ai|+|ci|andaici≠0(i=2,3,⋯,n−1),|bn|>|an|>0.

  此处，可以将A分解为如下形式

A=LU=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢α1γ2α2⋱⋱γn−1αn−1γnαn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥×⎡⎣⎢⎢⎢⎢⎢⎢⎢1β11β2⋱⋱1βn−11⎤⎦⎥⎥⎥⎥⎥⎥⎥

  由矩阵乘法原则可解出待定数αi,βi和γi

⎧⎩⎨⎪⎪⎪⎪⎪⎪α1=b1,γi=ai(i=2,3,⋯,n),βi−1=ci−1/αi−1(i=2,3,⋯,n),αi=bi−aiβi−1(i=2,3,⋯,n),

  由上式可知βi的递推公式

{β1=c1/b1,βi=ci/(bi−aiβi−1)(i=2,3,⋯,n−1)

并且只要算出了βi，其他待定数αi与γi均可通过已知数ai，bi与βi来表示.

算法归纳

1∘实现A=LU分解，即按照上述递推算式计算β1,β2,⋯βn；
2∘求解方程组LY=f，相应的递推算式是

{y1=f1/b1,yi=(fi−aiyi−1)/(bi−aiβi−1)(i=2,3,⋯,n)

3∘求解方程组UX=Y，相应的递推算式是

{xn=fn,xi=yi−βixi+1(i=n−1,n−2,⋯,2,1)

  由于计算β1→β2→⋯→βn及y1→y2→⋯→yn的过程称为追的过程，计算方程组的解xn→xn−1→⋯→x1的过程称为赶的过程，故上述求解方程组的方法称为追赶法

算法实现

% FILENAME: Chase.mfunction X = Chase( A, f )% here A must be a tridiagonal matrix    [N, ~] = size(A);    a = zeros(N-1);    b = zeros(N);    c = zeros(N-1);    % Step 1    for i = 1 : N        for j = 1 : N            if i == j                b(i) = A(i, i);            elseif i == j - 1                c(i) = A(i, j);            elseif i == j + 1                a(i) = A(i, j);            end        end    end    Beta = zeros(N-1, 1);    Beta(1) = c(1)/b(1);    for i = 2 : (N-1)        Beta(i) = c(i) / (b(i)-a(i) * Beta(i-1));    end    % Step 2    Y = zeros(N, 1);    Y(1) = f(1) / b(1);    for i = 2 : N        Y(i) = (f(i) -a(i)*Y(i-1)) / (b(i)-a(i)*Beta(i-1));    end    % Step 3    X = zeros(N, 1);    X(N) = Y(N);    for i = (N-1):-1:1        X(i) = Y(i) - Beta(i) * X(i + 1);    endend

% FILENAME: Chase_Run.mA = [2 1 0 0; 1/2 2 1/2 0; 0 1/2 2 1/2; 0 0 1 2];b = [-1/2 0 0 0]';fprintf('A = \n');disp(A);fprintf('b = \n');disp(b);X = Chase(A, b); % here A must be a tridiagonal matrixfprintf('X = \n');disp(X);

测试结果展示

>> Chase_RunA =     2.0000    1.0000         0         0    0.5000    2.0000    0.5000         0         0    0.5000    2.0000    0.5000         0         0    1.0000    2.0000b =    -0.5000         0         0         0X =    -0.2889    0.0778   -0.0222    0.0111>> A*X       % A*X==bans =   -0.5000    0.0000    0.0000         0

---

1. 问题描述完全参考于《数值计算方法（第三版）》 ↩