Método de Jacobi

O método de Jacobi é um algoritmo para resolver sistemas de equações lineares. Trata-se de uma versão simplificada do algoritmo de valores próprios de Jacobi. O método tem o nome do matemático Alemão Carl Gustav Jakob Jacobi.

O método iterativo de Jacobi é um método clássico que data do final do século XVIII. Técnicas iterativas são raramente utilizadas para solucionar sistemas lineares de pequenas dimensões, já que o tempo requerido para obter um mínimo de precisão ultrapassa o requerido pelas técnicas diretas como a eliminação gaussiana. Contudo, para sistemas grandes, com grande porcentagem de entradas de zero (sistemas esparsos), essas técnicas aparecem como alternativas mais eficientes. Sistemas esparsos de grande porte frequentemente surgem na análise de circuitos, na solução numérica de problemas de valor de limite e equações diferenciais parciais.^[1]

Descrição

Dada uma matriz quadrada de ${\textstyle n}$ equações lineares:

A\mathbf {x} =\mathbf {b}

em que:

A={\begin{bmatrix}a_{11}&a_{12}&\cdots &a_{1n}\\a_{21}&a_{22}&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}},\qquad \mathbf {x} ={\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}},\qquad \mathbf {b} ={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.

Então A pode ser decomposto num componente diagonal D e o resto R:

A=D R\qquad {\text{Onde}}\qquad D={\begin{bmatrix}a_{11}&0&\cdots &0\\0&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &a_{nn}\end{bmatrix}},\qquad R={\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\a_{21}&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &0\end{bmatrix}}

O sistema de equações lineares pode ser reescrito como:

D\mathbf {x} =\mathbf {b} -R\mathbf {x}

O método de Jacobi é um método iterativo que resolve o membro esquerdo da expressão em ordem a x ao usar o método resultante da iteração anterior no membro direito. Analiticamente, isto pode ser escrito como:

\mathbf {x} ^{(k 1)}=D^{-1}(\mathbf {b} -R\mathbf {x} ^{(k)})

Ou, equivalentemente:

x_{i}^{(k 1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j\neq i}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\ldots ,n.

Nota-se que a computação de $x_{i}^{k 1}$ é feita com base em todos os valores obtidos em iterações anteriores. Ao contrário do método de Gauss-Seidel, não se pode reescrever $x_{i}^{(k)}$ com $x_{i}^{(k 1)}$ , pois esse valor é necessário durante a continuação da computação. Esta é a diferença mais significativa entre o método de Jacobi e o método de Gauss-Seidel e é o motivo que o torna um algoritmo paralelo.

Algoritmo

MATLAB / FreeMat

X = input('Insira o vetor chute inicial ');

A = input('Insira a matriz dos coeficientes do sistema ');

b = input('Insira o vetor com os termos constantes ');

m = input('Insira o número máximo de iterações');

E = input('Insira a tolerância ');

n = input('Insira o número de equações ');

parar = 1;

ni = 0;

while (parar == 1 && ni < m)
    for i=1:n
        soma = 0;
        for j=1:n
            if(j~=i)
                soma = soma   A(i,j)*X(j)/A(i,i);
            end
        end
        x(i) = (b(i)/A(i,i)) - soma;
    end
    if (abs(norm(x) - norm(X))< E)
        parar = 0; 
    else
        X=x;
    end
    ni = ni   1;
end

disp('Resposta do sistema');

disp(X');

disp('Número de iterações utilizadas: ');

disp(ni);

Scilab

Na sequência, apresentamos um código Scilab que implementa uma versão simples do método de Jacobi. Os dados de entrada são: A matriz dos coeficientes, b vetor dos termos constantes, x0 vetor aproximação inicial, TOL tolerância, N número máximo de iterações. A saída é a calculada aproximação da solução de Ax = b.

function [x] = jacobi(A,b,x0,TOL,N)
   
    //preliminar
    nlin = size(A,1); //num. de linhas
    x = zeros(nlin,1); //cria x
   
    //iterações
    k = 1;
    while (k <= N)
       
        //iteração de Jacobi
        for i = 1:nlin
            x(i) = 0;
            for j = [1:i-1,i 1:nlin]
                x(i) = x(i) - A(i,j)*x0(j);
            end
            x(i) = (x(i)   b(i))/A(i,i);
        end
       
        //critério de parada
        if (norm(x-x0,'inf') < TOL) then
            return x;
        end
       
        //prepara nova iteração
        k = k 1;
        x0 = x;
    end
   
    //num. de iter. excedido
    disp('Método falhou.')
endfunction

JavaScript

Abaixo é apresentado a resolução utilizando a linguagem JavaScript.

function norm(Matrix) {
    //Requer implementação ou inclusão de biblitecas de terceiros
}
/*
* A = Matriz
* b = Resultado da matriz
*/
function gaussJacobi(A, b) {
	var X = new Array();
	var x = new Array();
	for (var k = 0; k < b.length; k  )
	{
		X[k] = 0; //Math.floor((Math.random() * 10000)   1);
	}

	var E = 0.00001;//Precisão
	var m = 1000;//Numero máximo de iterações
	var ni = 0;//Contador de iterações
	var continuar = true;

	while (continuar && ni < m) {
	    for (var i=0; i < b.length; i  ) {
	        soma = 0;
	        for (var j = 0; j < b.length; j  ) {
	            if (j != i) {
	                soma = soma   A[i][j]*X[j]/A[i][i];
	            }
	            x[i] = (b[i]/A[i][i]) - soma;
	        }
		}
	    if (Math.abs(norm(x) - norm(X)) < E) {
	        continuar = false;
		} else {
	        X=x.slice(0);
		}
	    ni = ni   1;
	}
	return X;
}

Java

A seguir apresenta-se o algoritmo na linguagem Java com o vetor padrão de aproximação inicial zero. Parâmetros: matriz → a matriz; b → o vetor dos resultados; e → a precisão; numIter → o número máximo de iterações. Esta função retorna um vetor onde cada posição corresponde a uma variável. Ex.: resultado[0] = x, resultado[1] = y, resultado[2] = z...

public static double[] jacobi(double[][] matriz, double[] b, double e, int numIter) {
	
	double[] x0 = new double[matriz.length];    //vetor dos resultados anteriores
	double[] x = x0.clone();                    //vetor dos resultados atuais
	double erro = 100, soma;
	int k = 0;      //A quantidade de iterações feitas

    //Zera o vetor x0
	for (int i = 0; i < x0.length; i  ) {
		x0[i] = 0;
	}

	while (erro >= e && k <= numIter) {
		for (int i = 0; i < matriz.length; i  ) {
			soma = 0;
			for (int j = 0; j < matriz[0].length; j  ) {
				if (i != j) {
					soma  = (matriz[i][j] * x0[j]) / matriz[i][i];
				}
				x[i] = (b[i] / matriz[i][i]) - soma;
			}
		}
		erro = calcErro(x, x0);
		x0 = x.clone();
		k  ;
	}
	
	return x;
	
}
	

//Faz o cálculo do erro
private static double calcErro(double[] a, double[] b) {
	double result[] = new double[a.length];
	for (int i = 0; i < a.length; i  ) {
		result[i] = Math.abs(a[i] - b[i]);
	}
	int cont = 0;            //um contador
	double maior = 0;       //o maior número contido no vetor
	while (cont < a.length) {
		maior = Math.max(maior, result[cont]);
		cont  ;
	}
	return maior;
}

Convergência

A iteração de Jacobi tem a forma:

$\mathbf {x} ^{(k 1)}=T\mathbf {x} ^{(k)} \mathbf {c}$

onde, ${\textstyle T=D^{-1}R}$ e $\mathbf {c} =D^{-1}b$ . Pode-se mostrar que, independentemente da aproximação inicial $\mathbf {x^{(0)}}$ , a sequência gerada $\{\mathbf {x} ^{(k)}\}$ converge para a solução de $A\mathbf {x} =\mathbf {b}$ se, e somente se, o raio espectral da matriz iteração $T$ for menor que $1$ .^[2] Ou seja, as iterações de Jacobi convergem para a solução do sistema se, e somente se, o raio espectral de ${\textstyle D^{-1}R}$ for menor que ${\textstyle 1}$ , i.e.:

$\rho (D^{-1}R)<1$ .

Consequentemente, o método de Jacobi é convergente sempre que a matriz ${\textstyle A}$ seja estrita ou irredutivelmente uma matriz estritamente diagonal dominante.

Recentemente, uma técnica de duplo ciclo foi introduzida para forçar a convergência para a solução correta do algoritmo de Jacobi mesmo quando as condições suficientes para convergência não são verificadas. A técnica de duplo ciclo funciona para matrizes positivas definidas ou dependentes de colunas.^[3]

Exemplos

1. Um sistema linear $Ax=b$ é dado por

A={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}},b={\begin{bmatrix}11\\13\\\end{bmatrix}}

e

x^{(0)}={\begin{bmatrix}1.1\\2.3\\\end{bmatrix}}.

Pretende-se usar a equação $x^{(k)}=Tx^{(k-1)} C$ . Agora é necessário encontrar $D^{-1}$ a matriz inversa dos valores da diagonal A através duma decomposição $A=LDU$

D^{-1}={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}},L={\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}}

e

U={\begin{bmatrix}0&-3\\0&0\\\end{bmatrix}}.

Para se encontrar T recorre-se à equação $T=D^{-1}(L U)$ .

T={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}\left\{{\begin{bmatrix}0&0\\-5&0\\\end{bmatrix}} {\begin{bmatrix}0&-3\\0&0\\\end{bmatrix}}\right\}={\begin{bmatrix}0&-1.5\\-0.714&0\\\end{bmatrix}}.

Agora, encontrado T, torna-se necessário obter C recorrendo à expressão $C=D^{-1}b$ .

C={\begin{bmatrix}1/2&0\\0&1/7\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}5.5\\1.857\\\end{bmatrix}}.

Assim, agora teremos $x^{(1)}=Tx^{(0)} C$ .

x^{(1)}={\begin{bmatrix}0&-1.5\\-0.714&0\\\end{bmatrix}}{\begin{bmatrix}1.1\\2.3\\\end{bmatrix}} {\begin{bmatrix}5.5\\1.857\\\end{bmatrix}}={\begin{bmatrix}2.05\\1.07\\\end{bmatrix}}.

E agora podemos encontrar $x^{(2)}$ .

x^{(2)}={\begin{bmatrix}0&-1.5\\-0.714&0\\\end{bmatrix}}{\begin{bmatrix}2.05\\1.07\\\end{bmatrix}} {\begin{bmatrix}5.5\\1.857\\\end{bmatrix}}={\begin{bmatrix}3.893\\0.3927\\\end{bmatrix}}.

Agora que se está na posse de X matrizes é possível avaliar a convergência para aproximar as soluções.

2. Efetuando os cálculos com três casas decimais, determine a solução do sistema de equações lineares $AX=B$ por meio do método de Jacobi, onde :

A={\begin{bmatrix}8&1&-1\\1&-7&2\\2&1&9\\\end{bmatrix}}B={\begin{bmatrix}8\\-4\\12\\\end{bmatrix}}

Nesse caso tem-se:

$x_{1}^{(k 1)}=1-0.125x_{2}^{(k)} 0,125x_{3}^{(k)}$

$x_{2}^{(k 1)}=0.571 0.143x_{1}^{(k)} 0.286x_{3}^{(k)}$

$x_{3}^{(k 1)}=1.333-0.222x_{1}^{(k)}-0.111x_{2}^{(k)}$ ,

E os resultados de aplicação do método de Jacobi são mostrados na tabela:

k	1	2	3	4	5	6
$x_{1}^{k}$	1,000	1,095	0,994	0,992	1,001	1,000
$x_{2}^{k}$	0,571	1,095	1,027	0,990	1,000	1,000
$x_{3}^{k}$	1,333	1,047	0,968	0,998	1,001	1,000

Notemos que a partir da 5 interação o método se estabiliza, com um erro de aproximadamente 0,04%.

Critério de convergência (critério das linhas)

Seja o sistema linear $Ax=b$ e seja: α $_{k}$ = $\sum _{j=1,j\neq k}^{n}$ ${\frac {|a_{kj}|}{|a_{kk}|}}$ .

Se α= máx α $_{k}$ <1, sendo 1≤k≤n, então o método de Jacobi gera uma sequência { $x(k)$ } convergente para a solução do sistema dado, independente da escolha da aproximação inicial, $x(0)$ .

Exemplo: Analisando a matriz A do sistema linear:

A={\begin{bmatrix}10&2&1\\1&5&1\\2&3&10\\\end{bmatrix}}

,

Temos:

α $_{1}$ = ${\frac {2 1}{10}}$ = ${\frac {3}{10}}$ = 0.3 < 1;

α $_{2}$ = ${\frac {1 1}{5}}$ =0.4 < 1;

α $_{3}$ = ${\frac {2 3}{10}}$ = 0.5 < 1

E então máx α $_{k}$ = 0.5 < 1, sendo 1≤k≤3 donde, pelo critério das linhas, temos garantia de convergência para o método de Jacobi. Sempre que o critério de linhas não for satisfeito, devemos tentar uma permutação de linhas e/ou colunas de forma a obtermos uma disposição para a qual a matriz dos coeficientes satisfaça o critério das linhas.^[4]

Referências

↑ Cálculo numérico: características matemáticas. Sperandio,Décio. Mendes,João Teixeira.Moken e Silva,Luiz Henry. São Paulo,Pearson,2003.
↑ Burden, Richard L. (2008). Análise Numérica. [S.l.: s.n.] ISBN 8522106010
↑ Fixing the convergence of the Gaussian belief propagation algorithm. J. K. Johnson, D. Bickson and D. Dolev. In the International symposium on information theory (ISIT), July 2009. http://arxiv.org/abs/0901.4192
↑ Cálculo numérico: aspectos teóricos e computacionais. Ruggiero, Márcia A. Gomes e Lopes,Vera lúcia da Rocha. 2ª Ed. São Paulo, Editora: Pearson, 1996

Ver também

Ligações externas

Black, Noel; Moore, Shirley; and Weisstein, Eric W. «Jacobi method». MathWorld (em inglês)
«Jacobi Method from www.math-linux.com»
«Module for Jacobi and Gauss–Seidel Iteration»
«Numerical matrix inversion»

[1] Cálculo numérico: características matemáticas. Sperandio,Décio. Mendes,João Teixeira.Moken e Silva,Luiz Henry. São Paulo,Pearson,2003.

[2] Burden, Richard L. (2008). Análise Numérica. [S.l.: s.n.] ISBN 8522106010

[Bickson-3] Fixing the convergence of the Gaussian belief propagation algorithm. J. K. Johnson, D. Bickson and D. Dolev. In the International symposium on information theory (ISIT), July 2009. http://arxiv.org/abs/0901.4192

[4] Cálculo numérico: aspectos teóricos e computacionais. Ruggiero, Márcia A. Gomes e Lopes,Vera lúcia da Rocha. 2ª Ed. São Paulo, Editora: Pearson, 1996

[1]

[2]

[3]

[4]