Tridiagonal matris algoritması

Vikipedi, özgür ansiklopedi
Atla: kullan, ara

Tridiagonal matris algoritması' (b.b.d. TDMA), Thomas algoritması olarak da bilinmektedir (Llewellyn Thomas ardından isimlendirilmiştir), sayısal lineer cebirde tridiagonal denklem sistemlerini çözmek için kullanılan basitleştirilmiş bir Gauss eleme yöntemidir.

n bilinmeyenli bir tridiagonal sistem şöyle yazılabilir:


a_i x_{i - 1}  + b_i x_i  + c_i x_{i + 1}  = d_i , \,\!

 a_1  = 0\, ve  c_n = 0\, 'dır.

Bu algoritma köşegen olarak baskın (b.b.d. "diagonally dominant") sistemlere uygulanabilir. Bu sistem matris formunda aşağıdaki gibi belirtilebilir:


\begin{bmatrix}
   {b_1} & {c_1} & {   } & {   } & { 0 } \\
   {a_2} & {b_2} & {c_2} & {   } & {   } \\
   {   } & {a_3} & {b_3} & \ddots & {   } \\
   {   } & {   } & \ddots & \ddots & {c_{n-1}}\\
   { 0 } & {   } & {   } & {a_n} & {b_n}\\
\end{bmatrix}
\begin{bmatrix}
   {x_1 }  \\
   {x_2 }  \\
   {x_3 }  \\
   \vdots   \\
   {x_n }  \\
\end{bmatrix}
=
\begin{bmatrix}
   {d_1 }  \\
   {d_2 }  \\
   {d_3 }  \\
   \vdots   \\
   {d_n }  \\
\end{bmatrix}
.


Bu tür sistemler için çözüm, O(n) mertebesinde aşama sayısı ile O(n^3) mertebesinde aşama sayısı gerektiren Gauss eleme yöntemi kullanılmadan elde edilebilir. İlk aşamada bütün a_i'ler elenir ve ardından geriye doğru yerine koyma yöntemi ile sonuca ulaşılır. Bu tip matrislerle çoğunlukla bir boyutlu Poisson denkleminin ayrıklaştırılmasında (ör. bir boyutlu difüzyon problemi) ve kübik eğri interpolasyonunda karşılaşılır.

Yöntem[değiştir | kaynağı değiştir]

İlk aşama matris katsayılarının aşağıdaki gibi düzenlenmesinden oluşur. Düzenlenmiş yeni katsayılar ayraç ile işaretlenmiştir:

c'_i =
\begin{cases}
\begin{array}{lcl}
  \cfrac{c_i}{b_i}                  & ; & i = 1 \\
  \cfrac{c_i}{b_i - c'_{i - 1} a_i} & ; & i = 2, 3, \dots, n-1 \\
\end{array}
\end{cases}
\,

ve:

d'_i =
\begin{cases}
\begin{array}{lcl}
  \cfrac{d_i}{b_i}                  & ; & i = 1 \\
  \cfrac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i} & ; & i = 2, 3, \dots, n. \\
\end{array}
\end{cases}
\,

Bu aşama ileri süpürme (b.b.d. "forward sweep") olarak adlandırılır. Çözüm ise geriye doğru yerine koyma yöntemi ile elde edilir:

x_n = d'_n\,
x_i = d'_i - c'_i x_{i + 1} \qquad ; \ i = n - 1, n - 2, \ldots, 1.

C içinde uygulaması[değiştir | kaynağı değiştir]

void solveMatrix (int n, double *a, double *b, double *c, double *v, double *x)
{
	/**
	 * n - number of equations
	 * a - sub-diagonal (means it is the diagonal below the main diagonal) -- indexed from 1..n-1
	 * b - the main diagonal
	 * c - sup-diagonal (means it is the diagonal above the main diagonal) -- indexed from 0..n-2
	 * v - right part
	 * x - the answer
	 */
	for (int i = 1; i < n; i++)
	{
		double m = a[i]/b[i-1];
                b[i] = b[i] - m * c[i - 1];
		v[i] = v[i] - m*v[i-1];
	}
 
	x[n-1] = v[n-1]/b[n-1];
 
	for (int i = n - 2; i >= 0; --i)
		x[i] = (v[i] - c[i] * x[i+1]) / b[i];
}

MATLAB içinde uygulaması[değiştir | kaynağı değiştir]

function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
n = length(b); % n is the number of rows
 
% Modify the first-row coefficients
c(1) = c(1) / b(1);    % Division by zero risk.
d(1) = d(1) / b(1);    % Division by zero would imply a singular matrix.
 
for i = 2:n-1
    temp = b(i) - a(i) * c(i-1);
    c(i) = c(i) / temp;
    d(i) = (d(i) - a(i) * d(i-1))/temp;
end
 
d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));
 
% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1
    x(i) = d(i) - c(i) * x(i + 1);
end
end

Fortran 90 içinde uygulaması[değiştir | kaynağı değiştir]

      subroutine solve_tridiag(a,b,c,v,x,n)
      implicit none
!	 a - sub-diagonal (means it is the diagonal below the main diagonal)
!	 b - the main diagonal
!	 c - sup-diagonal (means it is the diagonal above the main diagonal)
!	 v - right part
!	 x - the answer
!	 n - number of equations
 
        integer,intent(in) :: n
        real(8),dimension(n),intent(in) :: a,b,c,v
        real(8),dimension(n),intent(out) :: x
        real(8),dimension(n) :: bp,vp
        real(8) :: m
        integer i
 
! Make copies of the b and v variables so that they are unaltered by this sub
        bp(1) = b(1)
        vp(1) = v(1)
 
        !The first pass (setting coefficients):
firstpass: do i = 2,n
         m = a(i)/bp(i-1)
         bp(i) = b(i) - m*c(i-1)
         vp(i) = v(i) - m*vp(i-1)
        end do firstpass
 
         x(n) = vp(n)/bp(n)
        !The second pass (back-substition)
backsub:do i = n-1, 1, -1
          x(i) = (vp(i) - c(i)*x(i+1))/bp(i)
        end do backsub
 
    end subroutine solve_tridiag

Kaynaklar[değiştir | kaynağı değiştir]

Dış bağlantılar[değiştir | kaynağı değiştir]