Thomas' algorithm in MATLAB with exhaustive explanation

%{

Thomas algorithm implementation on MATLAB with explaining comments on usage.

See for details and ask questions: http://joydevel.blogspot.com/2012/07/matlab.html


Тестовый расчет
---------------
Трехдиагональная матрица системы:
1 2 0 0|8
2 3 4 0|17
0 5 6 7|35
0 0 8 9|26
Результат вычисления
x1 = 2
x2 = 3
x3 = 1
x4 = 2

Проверка:
1*2 + 2*3 + 0*1 + 0*2 = 8
2*2 + 3*3 + 4*1 + 0*2 = 17
0*2 + 5*3 + 6*1 + 7*2 = 35
0*2 + 0*3 + 8*1 + 9*2 = 26

USAGE:
  [X, Correct, Stable] = thomas(A,B,C,D)
Заметим, что необходимо передавать -b в параметрах, т.к. исходная формулировка задачи выглядит так:
-b1  c1   0    0  | d1
a1  -b2   c2   0  | d2
0    a2  -b3   c3 | d3
0    0    a3  -b4 | d4
Аналогично для любой размерности задачи.
TEST:
  [[2 3 1 2], true, false] = thomas([2,5,8],[-1,-3,-6,-9],[2,4,7],[8,17,35,26])

%}


function [ X, Correct, Stable ] = thomas( A, B, C, D )
% THOMAS Решение трехдиагональной СЛАУ

N = numel(D); % размерность матрицы

% Начальные значения переменных
X = zeros(1, N);
Correct = true;
Stable = true;
Test = false;

% Прямой ход — преобразуем к наддиагональной матрице
P = zeros(1, N-1);
Q = zeros(1, N);

P(1) = C(1)/B(1);
Q(1) = -D(1)/B(1);
for i=2:1:N
    t = (B(i)-A(i-1)*P(i-1));
    if (t == 0)
        Correct = false;
        X = -1;
    end
    if i < N
        P(i) = C(i)/t;
        if (abs(P(i))>=1)
            Stable = false;
            X = -1;
        end
    end
    Q(i) = (A(i-1)*Q(i-1)-D(i))/t;
end

% Обратный ход
X(N) = Q(N);
for i=N-1:-1:1
    X(i) = P(i)*X(i+1) + Q(i);
end

end