Método da matriz tridiagonal

Cálculo Numérico - Versão Python

4.5 Método da matriz tridiagonal

O método da matriz tridiagonal ou algoritmo de Thomas⁶ ou ainda TDMA (do inglês tridiagonal matrix algorithm) é o caso particular da eliminação gaussiana aplicada a matrizes tridiagonais.

Uma matriz tridiagonal é uma matriz quadrada cujos únicos elementos não nulos estão na diagonal principal e nas diagonais imediatamente acima e abaixo da principal. Um sistema tridiagonal é um sistema de equações lineares cuja matriz associada é tridiagonal, conforme a seguir: $\begin{array}{rcl} [\begin{matrix} b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1} \\ d_{2} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}] . & (4.88) \end{array}$

Observamos que não é necessário armazenar todos os $n^{2}$ elementos da matriz em memória, sendo suficiente armazenar os vetores $a_{n}$ , $b_{n}$ e $c_{n}$ . Por conveniência, a partir daqui, definiremos os elementos inexistentes na matriz $a_{1}$ e $c_{n}$ como zero:

a_{1} = c_{n} = 0 .

(4.89)

O algoritmo para a solução do sistema tridiagonal (4.88) pelo algoritmo de Thomas é dada pelas seguintes expressões: $\begin{array}{rcl} c_{i}^{'} & = & \{\begin{matrix} \begin{matrix} \frac{c_{i}}{b_{i}}, & i = 1 \\ \frac{c_{i}}{b_{i} - a_{i} c_{i - 1}^{'}}, & i = 2, 3, \dots, n - 1 \end{matrix} \end{matrix} & (4.90) \\ e \\ d_{i}^{'} & = & \{\begin{matrix} \begin{matrix} \frac{d_{i}}{b_{i}}, & i = 1 \\ \frac{d_{i} - a_{i} d_{i - 1}^{'}}{b_{i} - a_{i} c_{i - 1}^{'}}, & i = 2, 3, \dots, n . \end{matrix} \end{matrix} & (4.91) \end{array}$

Finalmente a solução final é obtida por substituição reversa: $\begin{array}{rcl} x_{n} & = & d_{n}^{'} & (4.92) \\ x_{i} & = & d_{i}^{'} - c_{i}^{'} x_{i + 1}, i = n - 1, n - 2, \dots, 1 . & (4.93) \end{array}$

Teorema 4.5.1. A aplicação da eliminação gaussiana sem pivotamento ao sistema (4.88) produz o algoritmo dado em (4.90) e (4.92).

Demonstração. O primeiro passo consiste em dividir todos os elementos da primeira linha de (4.88) por

b_{1}

\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}], & (4.94) \end{array}

onde $c_{1}^{'} = \frac{c_{1}}{b_{1}}$ e $d_{1}^{'} = \frac{d_{1}}{b_{1}} .$

O segundo passo consiste em substituir a segunda linha por ela mesma subtraída da linha 1 multiplicada por $a_{2}$ ( $l_{2} \leftarrow l_{2} - a_{2} l_{1}$ ): $\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ 0 & b_{2} - a_{2} c_{1}^{'} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2} - a_{2} d_{1}^{'} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}] . & (4.95) \end{array}$

Em seguida, dividimos a segunda linha por $b_{2} - a_{2} c_{1}^{'}$ , a fim de normalizar a diagonal principal: $\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ 0 & 1 & c_{2}^{'} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2}^{'} \\ d_{3} \\ ⋮ \\ d_{n} \end{matrix}] . & (4.96) \end{array}$

onde $c_{2}^{'} = \frac{c_{2}}{b_{2} - a_{2} c_{1}^{'}}$ e $d_{2}^{'} = \frac{d_{2} - a_{2} d_{1}^{'}}{b_{2} - a_{2} c_{1}^{'}}$ .

O próximo passo consiste em substituir a terceira linha por ela mesma subtraída da linha 2 multiplicada por $a_{3}$ ( $l_{3} \leftarrow l_{3} - a_{3} l_{2}$ ): $\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ 0 & 1 & c_{2}^{'} \\ 0 & b_{3} - a_{3} c_{2}^{'} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2}^{'} \\ d_{3} - a_{3} d_{2}^{'} \\ ⋮ \\ d_{n} \end{matrix}] . & (4.97) \end{array}$

A fim de normalizar o elemento da diagonal da terceira linha, dividimos toda a linha por $d_{3} - a_{3} d_{2}^{'}$ : $\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ 0 & 1 & c_{2}^{'} \\ 0 & 1 & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{n} & b_{n} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2}^{'} \\ d_{3}^{'} \\ ⋮ \\ d_{n} \end{matrix}] . & (4.98) \end{array}$

Este procedimento é realizado até que se atinja a última linha e temos o seguinte sistema: $\begin{array}{rcl} [\begin{matrix} 1 & c_{1}^{'} \\ 0 & 1 & c_{2}^{'} \\ 0 & 1 & ⋱ \\ ⋱ & ⋱ & c_{n - 1}^{'} \\ 0 & 1 \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ x_{n} \end{matrix}] = [\begin{matrix} d_{1}^{'} \\ d_{2}^{'} \\ d_{3}^{'} \\ ⋮ \\ d_{n}^{'} \end{matrix}] . & (4.99) \end{array}$

Neste estágio, podemos encontrar os $x_{n}$ através de substituição reversa, isto é: a última linha diz

x_{n} = d_{n}^{'} .

(4.100)

A penúltima linha diz

x_{n - 1} + c_{n - 1}^{'} x_{n} = d_{n - 1}^{'} \Rightarrow x_{n - 1} = d_{n - 1}^{'} - c_{n - 1}^{'} x_{n} .

(4.101)

Esse mesmo procedimento aplicado à linha $i = 1, \dots n - 1$ , nos dá

x_{i} = d_{i}^{'} - c_{i}^{'} x_{i + 1} .

(4.102)

$■$

Exemplo 4.5.1. Considere a resolução do seguinte sistema tridiagonal pelo algoritmo de Thomas: $\begin{array}{rcl} [\begin{matrix} 2 & 1 & 0 & 0 & 0 \\ 1 & 2 & 1 & 0 & 0 \\ 0 & 1 & 2 & 1 & 0 \\ 0 & 0 & 1 & 2 & 1 \\ 0 & 0 & 0 & 1 & 2 \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \end{matrix}] = [\begin{matrix} 4 \\ 4 \\ 0 \\ 0 \\ 2 \end{matrix}] . & (4.103) \end{array}$

Primeiramente identificamos os vetores $a$ , $b$ , $c$ e $d$ : $\begin{array}{rcl} a & = & (0, 1, 1, 1, 1) & (4.104) \\ b & = & (2, 2, 2, 2, 2) & (4.105) \\ c & = & (1, 1, 1, 1, 0) & (4.106) \\ d & = & (4, 4, 0, 0, 2) & (4.107) \end{array}$

Agora, calculamos os vetores $c^{'}$ e $d^{'}$ : $\begin{array}{rcl} c_{1}^{'} & = & \frac{c_{1}}{b_{1}} = \frac{1}{2} & (4.108) \\ c_{2}^{'} & = & \frac{c_{2}}{b_{2} - a_{2} c_{1}^{'}} = \frac{1}{2 - 1 \cdot \frac{1}{2}} = \frac{2}{3} & (4.109) \\ c_{3}^{'} & = & \frac{c_{3}}{b_{3} - a_{3} c_{2}^{'}} = \frac{1}{2 - 1 \cdot \frac{2}{3}} = \frac{3}{4} & (4.110) \\ c_{4}^{'} & = & \frac{c_{4}}{b_{4} - a_{4} c_{3}^{'}} = \frac{1}{2 - 1 \cdot \frac{3}{4}} = \frac{4}{5} & (4.111) \\ d_{1}^{'} & = & \frac{d_{1}}{b_{1}} = \frac{4}{2} = 2 & (4.112) \\ d_{2}^{'} & = & \frac{d_{2} - a_{2} d_{1}^{'}}{b_{2} - a_{2} c_{1}^{'}} = \frac{4 - 1 \cdot 2}{2 - 1 \cdot \frac{1}{2}} = \frac{4}{3} & (4.113) \\ d_{3}^{'} & = & \frac{d_{3} - a_{3} d_{2}^{'}}{b_{3} - a_{3} c_{2}^{'}} = \frac{0 - 1 \cdot \frac{4}{3}}{2 - 1 \cdot \frac{2}{3}} = - 1 & (4.114) \\ d_{4}^{'} & = & \frac{d_{4} - a_{4} d_{3}^{'}}{b_{4} - a_{4} c_{3}^{'}} = \frac{0 - 1 \cdot (- 1)}{2 - 1 \cdot \frac{3}{4}} = \frac{4}{5} & (4.115) \\ d_{5}^{'} & = & \frac{d_{5} - a_{5} d_{4}^{'}}{b_{5} - a_{5} c_{4}^{'}} = \frac{2 - 1 \cdot \frac{4}{5}}{2 - 1 \cdot \frac{4}{5}} = 1 & (4.116) \end{array}$

Finalmente, calculamos o vetor $x$ : $\begin{array}{rcl} x_{5} & = & d_{5}^{'} = 1 & (4.117) \\ x_{4} & = & d_{4}^{'} - c_{4}^{'} \cdot x_{5} = \frac{4}{5} - \frac{4}{5} \cdot 1 = 0 & (4.118) \\ x_{3} & = & d_{3}^{'} - c_{3}^{'} \cdot x_{4} = - 1 - \frac{3}{4} \cdot 0 = - 1 & (4.119) \\ x_{2} & = & d_{2}^{'} - c_{2}^{'} \cdot x_{3} = \frac{4}{3} - \frac{2}{3} \cdot (- 1) = 2 & (4.120) \\ x_{1} & = & d_{1}^{'} - c_{1}^{'} \cdot x_{2} = 2 - \frac{1}{2} \cdot 2 = 1 & (4.121) \end{array}$

E assim, obtemos o vetor $x = [1, 2, - 1, 0, 1]$ .

Código Python: Método da matriz tridiagonal

import numpy as np

def TDMA(a,b,c,d):
    #preliminares
    a = a.astype(’double’)
    b = b.astype(’double’)
    c = c.astype(’double’)
    d = d.astype(’double’)

    #recupera ordem do sistema
    n=np.shape(a)[0]

    #inicializa vetores auxiliares
    cl=np.zeros(n)
    dl=np.zeros(n)
    x=np.zeros(n)

    #calcula cl e dl
    cl[0]=c[0]/b[0]
    for i in np.arange(1,n-1,1):
       cl[i]=c[i]/(b[i]-a[i]*cl[i-1])

    dl[0]=d[0]/b[0]
    for i in np.arange(1,n,1):
       dl[i]=(d[i]-a[i]*dl[i-1])/(b[i]-a[i]*cl[i-1])

    #faz substituicao reversa para obter a solucao x
    x[n-1]=dl[n-1]
    for i in np.arange(n-2,-1,-1):
       x[i]=dl[i]-cl[i]*x[i+1]

    return x

Nesse código, usou-se cl e dl para denotar $c^{'}$ e $d^{'}$ . Observe que se for desnecessário preservar os valores originais dos vetores $c$ e $d$ , eles podem, com economia de memória e simplicidade de código, ser sobrescritos pelos vetores $c^{'}$ e $d^{'}$ , respectivamente. Eis uma nova implementação:

import numpy as np

def TDMA(a,b,c,d):
    #preliminares
    a = a.astype(’double’)
    b = b.astype(’double’)
    c = c.astype(’double’)
    d = d.astype(’double’)

    #recupera ordem do sistema
    n=np.shape(a)[0]

    #inicializa vetor x
    x=np.zeros(n)

    #calcula cl e dl sobrescrevendo-os em c e d
    c[0]=c[0]/b[0]
    for i in np.arange(1,n-1,1):
       c[i]=c[i]/(b[i]-a[i]*c[i-1])

    d[0]=d[0]/b[0]
    for i in np.arange(1,n,1):
       d[i]=(d[i]-a[i]*d[i-1])/(b[i]-a[i]*c[i-1])

    #faz substituicao reversa para obter a solucao x
    x[n-1]=d[n-1]
    for i in np.arange(n-2,-1,-1):
       x[i]=d[i]-c[i]*x[i+1]

    return x

A solução do sistema do Exemplo 4.5.1 pode ser obtida através dos seguintes comandos:

>>>a=np.array([1,1,1,1,1])
>>>b=np.array([2,2,2,2,2])
>>>c=np.array([1,1,1,1,1])
>>>d=np.array([4,4,0,0,2])
>>>TDMA(a,b,c,d)

E 4.5.1. Considere o problema linear tridiagonal dado por $\begin{array}{rcl} [\begin{matrix} 5 & 4 & 0 & 0 & 0 & 0 \\ 1 & 3 & 1 & 0 & 0 & 0 \\ 0 & 2 & 4 & 1 & 0 & 0 \\ 0 & 0 & 1 & 2 & 1 & 0 \\ 0 & 0 & 0 & 2 & 3 & 2 \\ 0 & 0 & 0 & 0 & 1 & 2 \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ x_{6} \end{matrix}] = [\begin{matrix} 13 \\ 10 \\ 20 \\ 16 \\ 35 \\ 17 \end{matrix}] . & (4.122) \end{array}$

Identifique os vetores $a$ , $b$ , $c$ e $d$ relativos ao algoritmo da matriz tridiagonal. Depois resolva o sistema usando o computador.

Resposta.

\begin{array}{rcl} a & = & (0, 1, 2, 1, 2, 1) & (4.123) \\ b & = & (5, 3, 4, 2, 3, 2) & (4.124) \\ c & = & (4, 1, 1, 1, 2, 0) & (4.125) \\ d & = & (13, 10, 20, 16, 35, 17) & (4.126) \\ x & = & (1, 2, 3, 4, 5, 6) & (4.127) \end{array}

E 4.5.2. Considere o seguinte sistema de equações lineares: $\begin{array}{rcl} x_{1} - x_{2} & = & 0 \\ - x_{j - 1} + 5 x_{j} - x_{j + 1} & = & cos (j ∕ 10), 2 \leq j \leq 10 \\ x_{11} & = & x_{10} ∕ 2 & (4.128) \end{array}$

Identifique os vetores $a$ , $b$ , $c$ e $d$ relativos ao algoritmo da matriz tridiagonal no sistema linear dado. Depois resolva o sistema usando o computador. Veja também Exercício 4.7.4

Resposta.

\begin{array}{rcl} a & = & (0, - 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1 ∕ 2) & (4.129) \\ b & = & (1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1) & (4.130) \\ c & = & (- 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1, - 1, 0) & (4.131) \\ d & = & (0, cos (2 ∕ 10), cos (3 ∕ 10), cos (4 ∕ 10), cos (5 ∕ 10), & (4.132) \\ cos (6 ∕ 10), cos (7 ∕ 10), cos (8 ∕ 10), cos (9 ∕ 10), cos (1), 0) & (4.133) \\ x & = & (0, 324295, 0, 324295, 0, 317115, 0, 305943, 0, 291539, & (4.134) \\ 0, 274169, 0, 253971, 0, 230846, 0, 20355, 0, 165301, 0, 082650) & (4.135) \end{array}

Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 19/8/2020 às 17:36:36.

Recursos

Projeto

IME - UFRGS
- Página do IME
- Página da UFRGS

UFRGS - IME - Recursos Educacionais Abertos de Matemática. Contato: reamat@ufrgs.br.

Informe erros ou edite você mesmo!