Tridiagonal matrix algorithm

In <a href="/facts/Numerical_linear_algebra/jnS6Tteb">numerical linear algebra</a>, the tridiagonal matrix algorithm, also known as the Thomas algorithm (named after <a href="/facts/Llewellyn_Thomas/cAzzHNCg">Llewellyn Thomas</a>), is a simplified form of <a href="/facts/Gaussian_elimination/23Rvf4ht">Gaussian elimination</a> that can be used to solve <a href="/facts/Tridiagonal_matrix/1SK5A0rP">tridiagonal systems of equations</a>. A tridiagonal system for n unknowns may be written as

a
          
            i
          
        
        
          x
          
            i
            −
            1
          
        
        +
        
          b
          
            i
          
        
        
          x
          
            i
          
        
        +
        
          c
          
            i
          
        
        
          x
          
            i
            +
            1
          
        
        =
        
          d
          
            i
          
        
        ,
      
    
    {\displaystyle a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},}

where 
 
 
 
 
 a
 
 1
 
 
 =
 0
 
 
 {\displaystyle a_{1}=0}
 
 and 
 
 
 
 
 c
 
 n
 
 
 =
 0
 
 
 {\displaystyle c_{n}=0}
 
.

[
            
              
                
                  
                    b
                    
                      1
                    
                  
                
                
                  
                    c
                    
                      1
                    
                  
                
                
                
                
                  0
                
              
              
                
                  
                    a
                    
                      2
                    
                  
                
                
                  
                    b
                    
                      2
                    
                  
                
                
                  
                    c
                    
                      2
                    
                  
                
                
                
              
              
                
                
                  
                    a
                    
                      3
                    
                  
                
                
                  
                    b
                    
                      3
                    
                  
                
                
                  ⋱
                
                
              
              
                
                
                
                  ⋱
                
                
                  ⋱
                
                
                  
                    c
                    
                      n
                      −
                      1
                    
                  
                
              
              
                
                  0
                
                
                
                
                  
                    a
                    
                      n
                    
                  
                
                
                  
                    b
                    
                      n
                    
                  
                
              
            
            ]
          
        
        
          
            [
            
              
                
                  
                    x
                    
                      1
                    
                  
                
              
              
                
                  
                    x
                    
                      2
                    
                  
                
              
              
                
                  
                    x
                    
                      3
                    
                  
                
              
              
                
                  ⋮
                
              
              
                
                  
                    x
                    
                      n
                    
                  
                
              
            
            ]
          
        
        =
        
          
            [
            
              
                
                  
                    d
                    
                      1
                    
                  
                
              
              
                
                  
                    d
                    
                      2
                    
                  
                
              
              
                
                  
                    d
                    
                      3
                    
                  
                
              
              
                
                  ⋮
                
              
              
                
                  
                    d
                    
                      n
                    
                  
                
              
            
            ]
          
        
        .
      
    
    {\displaystyle {\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}}.}

For such systems, the solution can be obtained in 
 
 
 
 O
 (
 n
 )
 
 
 {\displaystyle O(n)}
 
 operations instead of 
 
 
 
 O
 (
 
 n
 
 3
 
 
 )
 
 
 {\displaystyle O(n^{3})}
 
 required by <a href="/facts/Gaussian_elimination/23Rvf4ht">Gaussian elimination</a>. A first sweep eliminates the 
 
 
 
 
 a
 
 i
 
 
 
 
 {\displaystyle a_{i}}
 
's, and then an (abbreviated) backward substitution produces the solution. Examples of such matrices commonly arise from the discretization of 1D <a href="/facts/Poisson_equation/Pr8zkCXB">Poisson equation</a> and natural cubic <a href="/facts/Spline_interpolation/AowdpQwm">spline interpolation</a>.
Thomas' algorithm is not <a href="/facts/Numerical_stability/jRDdo8zV">stable</a> in general, but is so in several special cases, such as when the matrix is <a href="/facts/Diagonally_dominant/qbRQV7CV">diagonally dominant</a> (either by rows or columns) or <a href="/facts/Symmetric_positive_definite/Hbr6AuPS">symmetric positive definite</a>; for a more precise characterization of stability of Thomas' algorithm, see Higham Theorem 9.12. If stability is required in the general case, <a href="/facts/Gaussian_elimination/23Rvf4ht">Gaussian elimination</a> with <a href="/facts/Partial_pivoting/Lo9drS0e">partial pivoting</a> (GEPP) is recommended instead.

Tridiagonal matrix algorithm open-in-new

Tridiagonal matrix algorithm