SHARE
TWEET

Thomas' algorithm in MATLAB with exhaustive explanation

v_ignatyev Jan 15th, 2013 235 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %{
  2.  
  3. Thomas algorithm implementation on MATLAB with explaining comments on usage.
  4.  
  5. See for details and ask questions: http://joydevel.blogspot.com/2012/07/matlab.html
  6.  
  7.  
  8. Тестовый расчет
  9. ---------------
  10. Трехдиагональная матрица системы:
  11. 1 2 0 0|8
  12. 2 3 4 0|17
  13. 0 5 6 7|35
  14. 0 0 8 9|26
  15. Результат вычисления
  16. x1 = 2
  17. x2 = 3
  18. x3 = 1
  19. x4 = 2
  20.  
  21. Проверка:
  22. 1*2 + 2*3 + 0*1 + 0*2 = 8
  23. 2*2 + 3*3 + 4*1 + 0*2 = 17
  24. 0*2 + 5*3 + 6*1 + 7*2 = 35
  25. 0*2 + 0*3 + 8*1 + 9*2 = 26
  26.  
  27. USAGE:
  28.   [X, Correct, Stable] = thomas(A,B,C,D)
  29. Заметим, что необходимо передавать -b в параметрах, т.к. исходная формулировка задачи выглядит так:
  30. -b1  c1   0    0  | d1
  31. a1  -b2   c2   0  | d2
  32. 0    a2  -b3   c3 | d3
  33. 0    0    a3  -b4 | d4
  34. Аналогично для любой размерности задачи.
  35. TEST:
  36.   [[2 3 1 2], true, false] = thomas([2,5,8],[-1,-3,-6,-9],[2,4,7],[8,17,35,26])
  37.  
  38. %}
  39.  
  40.  
  41. function [ X, Correct, Stable ] = thomas( A, B, C, D )
  42. % THOMAS Решение трехдиагональной СЛАУ
  43.  
  44. N = numel(D); % размерность матрицы
  45.  
  46. % Начальные значения переменных
  47. X = zeros(1, N);
  48. Correct = true;
  49. Stable = true;
  50. Test = false;
  51.  
  52. % Прямой ход — преобразуем к наддиагональной матрице
  53. P = zeros(1, N-1);
  54. Q = zeros(1, N);
  55.  
  56. P(1) = C(1)/B(1);
  57. Q(1) = -D(1)/B(1);
  58. for i=2:1:N
  59.     t = (B(i)-A(i-1)*P(i-1));
  60.     if (t == 0)
  61.         Correct = false;
  62.         X = -1;
  63.     end
  64.     if i < N
  65.         P(i) = C(i)/t;
  66.         if (abs(P(i))>=1)
  67.             Stable = false;
  68.             X = -1;
  69.         end
  70.     end
  71.     Q(i) = (A(i-1)*Q(i-1)-D(i))/t;
  72. end
  73.  
  74. % Обратный ход
  75. X(N) = Q(N);
  76. for i=N-1:-1:1
  77.     X(i) = P(i)*X(i+1) + Q(i);
  78. end
  79.  
  80. end
RAW Paste Data
Top