• Author / Uploaded
• Anonymous UPRQTEVnHl

Citation preview

UNIVERSIDAD NACIONAL DE SAN AGUSTIN FACULTAD DE INGENIERIA CIVIL ESCUELA PROFESIONAL DE INGENIERIA SANITARIA

CURSO: MÉTODOS NUMÉRICOS Y PROGRAMACION TEMA: MÉTODO DE JACOBI ALUMNOS:  JAIME CÁCERES CAHUA  JAMES SEGALES SAICO

DOCENTE: DR. BIDER CALAPUJA

AREQUIPA-2016

Método de Jacobi

Marco Teórico En análisis numérico el método de Jacobi es un método iterativo, usado para resolver sistemas de ecuaciones lineales del tipo Ax  b

 a11 a  21 A   a31    am1

a12

a13

a22

a23

a32

a33

am 2

am3

 x1  x   2 X   x3      x   n

a1n  ... a2 n  ... a3n    ... amn  ...

b1  b   2 B  b3      b   m

El algoritmo toma su nombre del matemático alemán Carl Gustav Jakob. El método consiste en usar fórmulas como iteración de punto fijo. La base del método consiste en construir una sucesión convergente definida iterativamente. El límite de esta sucesión es precisamente la solución del sistema. A efectos prácticos si el algoritmo se detiene después de un número finito de pasos se llega a una aproximación al valor de x de la solución del sistema. En la iteración de Jacobi, se escoge una matriz Q que es diagonal y cuyos elementos diagonales son los mismos que los de la matriz A. La matriz Q toma la forma:

 a11 0  Q0    0

0

0

a22

0

0

a33

0

0

0  ... 0  ... 0    ... ann  ...

y la ecuación general se puede escribir como

Qx ( k )  (Q  A) x ( k 1)  b

Si denominamos R a la matriz A-Q:

0 a  21 R   a31    an1

a12

a13

0

a23

a32

0

an 2

an3

a1n  ... a2 n  ... a3n    ... 0 

...

la ecuación se puede reescribir como:

Qx ( k )   Rx ( k 1)  b El producto de la matriz Q por el vector columna x(k) será un vector columna. De modo análogo, el producto de la matriz R por el vector columna x(k-1)será también un vector columna. La expresión anterior, que es una ecuación vectorial, se puede expresar por n ecuaciones escalares (una para cada componente del vector). De este modo, podemos escribir, para un elemento i cualquiera y teniendo en cuenta que se trata de un producto matriz-vector: n

n

j 1

j 1

 qij x j (k )   rij x j (k 1)  bi Si tenemos en cuenta que en la matriz Q todos los elementos fuera de la diagonal son cero, en el primer miembro el único término no nulo del sumatorio es el que contiene el elemento diagonal qii, que es precisamente aii. Más aún, los elementos de la diagonal de Rson cero, por lo que podemos eliminar el término i=j en el sumatorio del segundo miembro. De acuerdo con lo dicho, la expresión anterior se puede reescribir como:

aii xi

(k )



n



j 1, j i

aijx ( k 1) b j

i

de donde despejando xi(k) obtenemos:

xi  (b i  (k )

n

a

j 1, j i

ijx j ( k 1)

) / aii

que es la expresión que nos proporciona las nuevas componentes del vector x(k) en función de vector anterior x(k-1) en la iteración de Jacobi. A continuacion se presenta un algoritmo para el método de Jacobi.

input n, (aij ), (bi ), ( xi ), M for

k  1, 2,...., Mdo for

i  1, 2,...., ndo ui  (bi   j  1, j  aij xj ) / aii n

end for

i  1, 2,...., ndo

xi  ui end

end output k , ( x j )

El método de Jacobi se basa en escribir el sistema de ecuaciones en la forma:

 x1  (b1  a21 x2  a31 x3... an1 xn ) / a11   x  (b  a x  a x a x ) / a  2 12 1 32 3... n 2 n 22   2 .    .  .     xn  (bn  a1n x1  a2 n x2... ) / ann 

Partimos de una aproximación inicial para las soluciones al sistema de ecuaciones y sustituimos estos valores en la ecuación. De esta forma, se genera una nueva aproximación a la solución del sistema, que en determinadas condiciones, es mejor que la aproximación inicial. Esta nueva aproximación se puede sustituir de nuevo en la parte derecha de la ecuación y así sucesivamente hasta obtener la convergencia.

Programa para resolver sistemas lineales utilizando el método de Jacobi function jacovi = jacovi(A,b,val_inicial,max_error,max_iter) try iter = 0; xk = val_inicial; D = diag(diag(A)); L = (tril(A) - diag(diag(A))).*(-1); U = (triu(A) - diag(diag(A))).*(-1); while(1) final = (inv(D)*b) + ((inv(D)*(L+U))*xk); iter=iter+1; fprintf('Iteracion: %d \n',iter) disp('Valores X k+1') disp(final) error = abs((final-xk)./xk) if ( (max_iter~=-1 & max_error==-1)| (max_iter=max_iter) break; else xk = final; end end if(max_iter