Descarga la aplicación para disfrutar aún más
Vista previa del material en texto
NOTAS DE CLASE Tema: Sistemas de Ecuaciones Lineales Métodos Directos Métodos Iterativos Mal condicionamiento Tema: Sistemas de ecuaciones No Lineales Métodos de Punto Fijo y Newton-Raphson 14/12/2020 Dada una matriz A (m x n) y un vector B (m x 1), encontrar el vector X (n x 1) que satisface la ecuación A X = B (m ecuaciones con n incógnitas) Condiciones suficientes para la existencia de solución única son - m = n - determinante de A 0 EXISTENCIA y UNICIDAD DE LA SOLUCIÓN A es regular <==> det (A) 0 <==> existe A-1 Recordar: TRANSFORMACIONES ELEMENTALES Cualquiera de las siguientes operaciones aplicadas a un sistema de ecuaciones lineales produce un sistema equivalente - INTERCAMBIO: el orden de las ecuaciones puede alterarse - ESCALAMIENTO: multiplicar una ecuación por una constante no nula - SUSTITUCIÓN: una ecuación puede ser reemplazada por la suma de ella misma más un múltiplo de otra ecuación A nxn X nx1 = B nx1 ELIMINACIÓN GAUSSIANA sistema de ecuaciones MATRIZ AMPLIADA ELIMINACIÓN GAUSSIANA sistema de ecuaciones MATRIZ AMPLIADA La primera fila se usa para eliminar los elementos de la primera columna que están por debajo de la diagonal principal. El elemento a11 es el pivote. Los valores de mk1 = ak1/a11 son los multiplicadores. pivote m21= 2 m31= 4 m41= -3 Ahora, la segunda fila se usa para eliminar los elementos de la segunda columna que están por debajo de la diagonal principal. El elemento a22 es el pivote. Los valores de mk2 = ak2/a22 son los multiplicadores. pivote m32= 1.5 m42= -1.75 Finalmente, restando la cuarta fila de la tercera multiplicada por m43 = a43/a33 =-1.9 se obtiene el sistema triangular superior. pivote m43= -1.9 Utilizando el algoritmo de SUSTITUCIÓN HACIA ATRÁS obtenemos IDEA ALGORITMO PASO 1 ELIMINACIÓN GAUSSIANA Transformar el sistema A X = B en un sistema equivalente A* X = B* tal que A* sea una matriz triangular superior. A matriz invertible de n x n, B vector n x 1. Almacenar todos los coeficientes en la matriz ampliada (A|B) PASO 2 Eliminar la incógnita xq en todas las filas desde la q+1 en adelante. PASO 3 Sustitución hacia atrás. Código F90 BÁSICO del método de Eliminación de Gauss simple program sist_ec_lin use kinds implicit none integer(ikind4) :: n, i real(rkind8), allocatable :: A(:,:) real(rkind8), allocatable :: B(:), x(:) open(unit=10,file="input.in",status="old") read(10,*) n allocate(A(n,n),B(n),x(n)) do i=1,n ! Lectura de la matriz de coeficientes y t.i. read(10,*) A(i,:), B(i) end do close(10) call elim_gauss_v0(n,A,B,x) write(*,*)"x=", x deallocate(A,B,x) end program sist_ec_lin Sample input file: input.in subroutine elim_gauss_v0(n,A,B,x) use kinds ! Eliminación de Gauss (versión preliminar) implicit none integer(ikind4) :: n, i, j real(rkind8), dimension(n,n) :: A real(rkind8), dimension(n) :: B, x real(rkind8), dimension(n,n+1) :: Amp real(rkind8) :: mult, sumi ! ----------------------------------------------------------------------------------------------------- do i=1, n ! Armamos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1)=B(:) do i=1, n-1 ! Reducción de la matriz ampliada if ( Amp(i,i) == 0 ) then write(*,*)"ERROR: cabecera nula" ; stop endif do j=i+1, n mult=Amp(j,i)/Amp(i,i) ; Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do end do x(n)=Amp(n,n+1)/Amp(n,n) ! Sustitución hacia atrás do i=n-1,1,-1 sumi=0._rkind8 do j=i+1, n sumi=sumi+Amp(i,j)*x(j) end do x(i)=(Amp(i,n+1)-sumi)/Amp(i,i) end do ! ----------------------------------------------------------------------------------------------------- end subroutine elim_gauss_v0 2 1 1 2 1 -1 0 Mejoras simples, necesarias para los códigos anteriores: - verificar en el programa principal que A.x=B: evaluación de || A.x – B || - más mejoras ? - implementar en la subrutina la permutación de filas cuando se encuentra cabecera nula y detectar posible incompatibilidad o indeterminación del sistema de ecuaciones subroutine elim_gauss(n,A,B,x) ! Eliminación de Gauss (SIN pivoteo parcial) use kinds implicit none integer(ikind4) :: n, i, j real(rkind8), dimension(n,n) :: A real(rkind8), dimension(n) :: B, x real(rkind8), dimension(n,n+1) :: Amp real(rkind8) :: mult, sumi real(rkind8), dimension(n+1) :: aux ! --------------------------------------------------------------------------------------------- do i=1, n !Armarmos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1)=B(:) do i=1, n-1 ! Reducción de la matriz ampliada if ( Amp(i,i) == 0 ) then do j=i+1, n if ( Amp(j,i) /= 0 ) then aux=Amp(i,:) ; Amp(i,:)=Amp(j,:) ; Amp(j,:)=aux exit elseif ( j == n) then write(*,*)"ERROR: no se puede aplicar el método" stop endif end do endif do j=i+1, n mult=Amp(j,i)/Amp(i,i) ; Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do end do x(n)=Amp(n,n+1)/Amp(n,n) ! Sustitución hacia atrás do i=n-1,1,-1 sumi=0. do j=i+1, n sumi=sumi+Amp(i,j)*x(j) end do x(i)=(Amp(i,n+1)-sumi)/Amp(i,i) end do ! --------------------------------------------------------------------------------------------- end subroutine elim_gauss program sist_ec_lin use kinds implicit none integer(ikind4) :: n, i real(rkind8), allocatable :: A(:,:) real(rkind8), allocatable :: B(:), x(:), test(:) open(unit=10,file="input.in",status="old") read(10,*) n allocate(A(n,n),B(n),x(n),test(n)) do i=1,n ! Lectura de la matriz de coeficientes y t.i. read(10,*) A(i,:), B(i) end do close(10) call elim_gauss(n,A,B,x) write(*,*)"x=", x test=matmul(A,x)-B write(*,*) "Norma de Ax-B=",norm2(test) deallocate(A,B,x,test) end program sist_ec_lin Sample input file: input.in 2 0 1 1 1 1 2 ESTRATEGIAS DE PIVOTEO PARA REDUCIR LOS ERRORES DE REDONDEO 1.133 x1+5 .281 x2=6 .414 24 .14 x1−1 .210 x2=22 .03 Los valores son la solución del sistema x1=x2=1 .000 Resolvemos el sistema haciendo las operaciones con 4 dígitos significativos, se obtiene: x1=0 .9956 , x2=1 .001 24 .14 x1−1 .210 x2=22 .03 1.133 x1+5 .281 x2=6 .414 Intercambiamos las ecuaciones del sistema anterior Resolvemos nuevamente utilizando 4 dígitos significativos, se obtiene: x1=1.000 , x 2=1 .000 Pivoteo trivial y pivoteo parcialPivoteo trivial y pivoteo parcial ELIMINACIÓN GAUSSIANA CON PIVOTEO PARCIAL IDEA ALGORITMO PASO 1 Permanecen sin cambios PASO 2 Encontrar la fila k “más grande”, es decir, |akq| = max {|aqq|, |aq+1q|, …, |an-1q|, |anq|}. Intercambiar la fila q-ésima con la fila k-ésima Eliminar la incógnita xq en todas las filas desde la q+1 en adelante. PASO 3 Permanece sin cambios PIVOTEO PARCIAL VS. TOTAL PIVOTEO TOTAL: En la etapa r-esima del proceso de eliminación, elegimos el elemento a pq =max{ | a ij | / i,j ≥r } como nuevo pivote. Para ello intercambiamos las filas r y p, y las columnas r y q de forma que situamos el elemento a pq en laposicion (r,r). Obviamente hemos tomado i,j ≥r para no perturbar los ceros que ya teníamos. Posteriormente continuamos la eliminación con el nuevo pivote PIVOTEO PARCIAL: Se busca solamente en la r-esima columna; es decir, se toma a pr = max { | a ij | / i ≥ r } como nuevo pivote. Para ello intercambiamos las filas r y p, continuando posteriormente el proceso de eliminación. En la práctica, el método de Gauss con pivoteo total puede consumir mucho tiempo computacional, pues para hallar el máximo en cada paso hay que buscar entre (m - r +1).(n - r +1) elementos. En el otro caso, además del ahorro de tiempo, las incógnitas de nuestro sistema no cambian de orden en el sistema reducido. Por ello, en general, es suficiente utilizar un pivoteo parcial. subroutine elim_gauss_pivoteo_parcial(n,A,B,x)! Eliminación de Gauss CON pivoteo parcial use kinds implicit none integer(ikind4) :: n, i, j, row_max_abs_elem real(rkind8), dimension(n,n) :: A real(rkind8), dimension(n) :: B, x real(rkind8), dimension(n,n+1) :: Amp real(rkind8) :: mult, sumi, max_abs_elem real(rkind8), dimension(n+1) :: aux ! --------------------------------------------------------------------------------------------- do i=1, n ! Armarmos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1)=B(:) do i=1, n-1 ! Reducción de la matriz ampliada max_abs_elem=abs(Amp(i,i)) ; row_max_abs_elem=i do j=i+1,n if (abs(Amp(j,i)) > max_abs_elem) then max_abs_elem=abs(Amp(j,i)) ; row_max_abs_elem=j endif end do aux=Amp(i,:) ; Amp(i,:)=Amp(row_max_abs_elem,:) ; Amp(row_max_abs_elem,:)=aux if ( Amp(i,i) == 0 ) then write(*,*)"ERROR: no se puede aplicar el método" ; stop end if do j=i+1, n mult=Amp(j,i)/Amp(i,i) ; Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do end do x(n)=Amp(n,n+1)/Amp(n,n) ! Sustitución hacia atrás do i=n-1,1,-1 sumi=0._rkind8 do j=i+1, n sumi=sumi+Amp(i,j)*x(j) end do x(i)=(Amp(i,n+1)-sumi)/Amp(i,i) end do ! --------------------------------------------------------------------------------------------- end subroutine elim_gauss_pivoteo_parcial ¿ Qué cambios hay que introducir a los algoritmos para implementar eliminación de Gauss-Jordan ? subroutine elim_gauss_jordan_pivoteo_parcial(n,A,B,x) ! Eliminación de Gauss-Jordan CON pivoteo parcial implicit none integer(4) :: n, i, j, row_max_abs_elem real(8), dimension(n,n) :: A real(8), dimension(n) :: B, x real(8), dimension(n,n+1) :: Amp real(8) :: mult, sumi, max_abs_elem real(8), dimension(n+1) :: aux do i=1, n !Armarmos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1)=B(:) do i=1, n ! Reducción de la matriz ampliada max_abs_elem=abs(Amp(i,i)) ; row_max_abs_elem=i do j=i+1,n if (abs(Amp(j,i)) > max_abs_elem) then max_abs_elem=abs(Amp(j,i)) ; row_max_abs_elem=j end if end do aux=Amp(i,:) ; Amp(i,:)=Amp(row_max_abs_elem,:) Amp(row_max_abs_elem,:)=aux if ( Amp(i,i) == 0 ) then write(*,*)"ERROR: no se puede aplicar el método" ; stop endif do j=1, n if (j == i) cycle mult=Amp(j,i)/Amp(i,i) Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do Amp(i,:)=Amp(i,:)/Amp(i,i) end do x(:)=Amp(:,n+1) end subroutine elim_gauss_jordan_pivoteo_parcial program sist_ec_lin implicit none integer(4) :: n, i real(8), allocatable :: A(:,:) real(8), allocatable :: B(:), x(:), test(:) open(unit=10,file="input.in",status="old") read(10,*) n allocate(A(n,n),B(n),x(n),test(n)) do i=1,n ! Lectura de la matriz de coeficientes y t.i. read(10,*) A(i,:), B(i) end do close(10) call elim_gauss(n,A,B,x) ! Eliminación de Gauss write(*,*)"Resultado Gauss, x=", x test=matmul(A,x)-B write(*,*) "Norma de Ax-B=",norm2(test) call elim_gauss_jordan_pivoteo_parcial(n,A,B,x) ! Elim. Gauss-Jordan write(*,*)"Resultado Gauss-Jordan, x=", x test=matmul(A,x)-B write(*,*) "Norma de Ax-B=",norm2(test) deallocate(A,B,x,test) end program sist_ec_lin DEBE SER SIN PIVOTEO PARCIAL !!!!! subroutine elim_gauss_jordan_pivoteo_parcial_msist(n,A,m,B,x) ! Esta rutina resuelve simultaneamente, m sistemas de ecuaciones nxn de la forma: Ax=b_i con i=1,2,...,m; ! todos con la misma matriz de coeficientes A (nxn) y los m vectores b_i son las columnas de la matriz B (nxm). ! Las soluciones de los m sistemas de ecuaciones son las m columnas de la matriz x (nxm) que es calculada. ! ! El método utilizado es Eliminación de Gauss-Jordan CON pivoteo parcial implicit none integer(4) :: n, i, j, row_max_abs_elem, m real(8), dimension(n,n) :: A real(8), dimension(n,m) :: B, x real(8), dimension(n,n+m) :: Amp real(8) :: mult, sumi, max_abs_elem real(8), dimension(n+m) :: aux do i=1, n !Armarmos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1:n+m)=B do i=1, n ! Reducción de la matriz ampliada max_abs_elem=abs(Amp(i,i)) ; row_max_abs_elem=i do j=i+1,n if (abs(Amp(j,i)) > max_abs_elem) then max_abs_elem=abs(Amp(j,i)) ; row_max_abs_elem=j end if end do aux=Amp(i,:) ; Amp(i,:)=Amp(row_max_abs_elem,:) ; Amp(row_max_abs_elem,:)=aux if ( Amp(i,i) == 0 ) then write(*,*)"ERROR: no se puede aplicar el método" ; stop end if do j=1, n if (j == i) cycle mult=Amp(j,i)/Amp(i,i) ; Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do Amp(i,:)=Amp(i,:)/Amp(i,i) end do x=Amp(:,n+1:n+m) end subroutine elim_gauss_jordan_pivoteo_parcial_msist subroutine elim_gauss_jordan_msist(n,A,m,B,x) ! Esta rutina resuelve simultaneamente, m sistemas de ecuaciones nxn de la forma: Ax=b_i con i=1,2,...,m; ! todos con la misma matriz de coeficientes A (nxn) y los m vectores b_i son las columnas de la matriz B (nxm). ! Las soluciones de los m sistemas de ecuaciones son las m columnas de la matriz x (nxm) que es calculada. ! Notar que si m=n y B es la matriz identidad nxn, la matriz x calculada es la matriz inversa de A, PERO EN ESTE CASO ! el método de eliminación de Gauss-Jordan utilizado DEBE SER SIN pivoteo parcial implicit none integer(4) :: n, i, j, row_max_abs_elem, m real(8), dimension(n,n) :: A real(8), dimension(n,m) :: B, x real(8), dimension(n,n+m) :: Amp real(8) :: mult, sumi, max_abs_elem real(8), dimension(n+m) :: aux do i=1, n !Armarmos la matriz ampliada Amp(:,i)=A(:,i) end do Amp(:,n+1:n+m)=B do i=1, n ! Reducción de la matriz ampliada max_abs_elem=abs(Amp(i,i)) ; row_max_abs_elem=i do j=i+1,n if (abs(Amp(j,i)) > max_abs_elem) then max_abs_elem=abs(Amp(j,i)) ; row_max_abs_elem=j end if end do aux=Amp(i,:) ; Amp(i,:)=Amp(row_max_abs_elem,:) ; Amp(row_max_abs_elem,:)=aux if ( Amp(i,i) == 0 ) then write(*,*)"ERROR: no se puede aplicar el método" ; stop end if do j=1, n if (j == i) cycle mult=Amp(j,i)/Amp(i,i) ; Amp(j,:)=Amp(j,:)-mult*Amp(i,:) end do Amp(i,:)=Amp(i,:)/Amp(i,i) end do x=Amp(:,n+1:n+m) end subroutine elim_gauss_jordan_msist SI QUEREMOS USARLO PARA CALCULAR LA MATRIZ INVERSA: SIN PIVOTEO PARCIAL CONSIDERACIONES NUMÉRICAS ERROR DE TRUNCAMIENTO: Ninguno. En la ausencia de error de redondeo se obtiene la respuesta exacta en un número finito de pasos. ERROR DE REDONDEO: Puede ser serio, se recomienda el uso Pivoteo Parcial para minimizar el efecto. ESFUERZO COMPUTACIONAL: Proceso de eliminación requiere de (n3) operaciones. Sustitución hacia atrás (n2) operaciones. Total de operaciones (n3). SISTEMAS MAL CONDICIONADOS (pequeños cambios en los elementos de A o de B provocan grandes cambios en X): los métodos numéricos son proclives a tener más errores de redondeo (veremos más adelante). FIN CLASE 1 de SISTEMAS DE ECUACIONES Extienden los métodos iterativos introducidos en el capítulo de “raíces de ecuaciones no lineales” a espacios de mayor dimensión. En particular consideraremos extensiones del método de punto fijo aplicado a sistemas de ecuaciones lineales. METODOS ITERATIVOS x = g(x) x: vector g: función vectorial de x MÉTODO ITERATIVO DE JACOBI 4 x1− x2+ x3= 7 4 x1− 8 x 2+ x3= −21 −2 x1+ x2+ 5 x3= 15 sistema de ecuaciones El sistema anterior puede ser escrito en la forma siguiente: x1= 7+x2−x3 4 x2= 21+4 x1+x3 8 x3= 15+2 x1−x2 5 x1 ( k ) = 7+x2 (k−1 ) −x3 ( k−1) 4 x2 ( k ) = 21+4 x1 ( k−1 )+x3 (k−1) 8 x3 ( k ) = 15+2x1 (k−1 )−x2 ( k−1) 5 Lo cual sugiere el siguiente proceso iterativo x1 (2 ) = 7+3 .375−3 4 =1 .84375 x2 (2 ) = 21+4∗1.75+3 8 =3 .875 x3 (2 ) = 15+2∗1 .75−3 .375 5 =3 .025 Tomando el vector aproximación inicial x(0) = ( 1 , 2 , 2 ) x(1) = ( 1.75 , 3.375 , 3 ) x1 (1 ) = 7+2−2 4 =1.75 x2 (1 ) = 21+4∗1+2 8 =3 .375 x3 (1 ) = 15+2∗1−2 5=3 x(2) = ( 1.84375 , 3.875 , 3.025 ) MÉTODO ITERATIVO DE JACOBI PASO 1 PASO 2 PASO 3 Sigo? Si ||x(k+1) − x(k)|| < tol, detengo el proceso y retorno el valor de x(k+1) como respuesta Actualizo k ← k+1. Retorno al paso 1 x(k+1 )=g (x(k ) )Fórmula de recurrencia MÉTODO ITERATIVO DE JACOBI x j (k+1 ) = b j−a j 1 x1 (k ) −.. .−a jj−1 x j−1 (k ) −a jj+1 x j+1 (k ) −.. .−a jN xN (k ) a jj Norma PASO 1 PASO 2 PASO 3 MÉTODO ITERATIVO DE GAUSS-SEIDEL x j (k+1 ) = b j−a j 1 x1 (k+1 ) − .. .−a jj−1 x j−1 (k+1) −a jj+1 x j+ 1 (k ) −. . .−a jN xN (k ) a jj Sigo? Si ||x(k+1) − x(k)|| < tol, detengo el proceso y retorno el valor de x(k+1) como respuesta Actualizo k ← k+1. Retorno al paso 1 x(k+1 )=g (x(k ) )Fórmula de recurrencia Norma CONDICION SUFICIENTE DE CONVERGENCIA Definición: Matriz de diagonal estrictamente dominante |aii|>∑ j= 1 n |aij| ∀ i=1,… ,n Teorema: Si la matriz de coeficientes A del sistema de ecuaciones es diagonal estrictamente dominante, entonces el proceso iterativo convergerá a la solución exacta, independientemente del vector inicial x(0) What does it mean by ill conditioned and well-conditioned system of equations? - A system of equations is considered to be well conditioned if a small change in the coefficient matrix or a small change in the right hand side results in a small change in the solution vector. - A system of equations is considered to be ill conditioned if a small change in the coefficient matrix or a small change in the right hand side results in a large change in the solution vector. Sistemas de ecuaciones algebraicas bien/mal condicionados Some properties of the norm of a matrix Para poder “confiar” en los 4 primeros dígitos significativos de un número, su error relativo debe ser ~ 10-4 (*) (*) Por ejemplo, el mínimo número que se redondea a 5.043 (manteniendo 4 cifras significativas) es 5.0425 y el máximo es 5.0434999999999. Es decir, x=5.0430±0.0005. En una aritmética de punto flotante, se llama épsilon de la máquina (ε mach ) al menor valor de una determinada máquina que cumple lo siguiente: 1.0 + ε mach > 1.0 ε mach es el número decimal más pequeño que, sumado a 1, la computadora nos arroja un valor > 1 es decir, que no es redondeado a 1. Representa la exactitud relativa de la aritmética del computador. La existencia del épsilon de la máquina es una consecuencia de la precisión finita de la aritmética en coma flotante. Sistemas de ecuaciones algebraicas NO lineales MÉTODOS ITERATIVOS PARA SISTEMAS NO LINEALES f 1 ( x,y )=x 2−2 x− y+0 .5 f 2 ( x,y )=x 2+4 y2−4 Queremos hallar la forma de resolver el sistemas de ecuaciones: f 1 ( x,y )=0 f 2 ( x,y )=0 Consideremos las funciones MÉTODOS ITERATIVOS PARA SISTEMAS NO LINEALES f1(x,y) y f2(x,y) definen implícitamente curvas en el plano f 1 ( x,y )=x 2−2 x− y+0 .5 f 2 ( x,y )=x 2+4 y2−4 Una solución del sistema es un punto (p,q) en el que ambas curvas se cruzan f1(p,q)=0 y f2(p,q)=0 es una parábola es una elipse Definición Un Punto Fijo del sistema de 2 ecuaciones Definición El método de iteración de punto fijo es es un punto (p,q) tal que p=g1(p,q) y q=g2(p,q). x=g1 ( x,y ) y=g2 ( x,y ) pk+1=g1( pk ,qk ) qk+1=g2( pk ,qk ) k=0,1,. .. METODO ITERATIVO DE PUNTO FIJO despejar x de f1 despejar y de f2 x= x2− y+0 .5 2 y= −x2−4 y2+8 y+4 8 Tomando (p0,q0) generamos una sucesión de valores {(pk+1,qk+1)} pk+1=g1( pk ,qk )= pk 2−qk+0 .5 2 qk+1=g2( pk ,qk )= −pk 2−4 qk 2+8 qk+4 8 (s.m.m. -8y a f2) f 1 ( x,y )=x 2−2 x− y+0 .5=0 f 2 ( x,y )=x 2+4 y2−4=0 EJEMPLO Sea el siguiente sistema no lineal CASO I Si usamos como punto inicial (p0,q0) = (0,1), entonces CONVERGE ! CASO II Si usamos como punto inicial (p0,q0) = (2,0), entonces DIVERGE ! Cambiando las generatrices… (s.m.m. -11y a f2) (s.m.m. -2x a f1) CONVERGE ! pk+1=g1( pk ,qk )= −pk 2+4 pk+qk−0 .5 2 qk+1=g2( pk ,qk )= −pk 2−4 qk 2+11qk+4 11 TEORÍA Cómo determinar las generatrices de forma que la sucesión generada sea convergente? Sean las funciones g1 y g2 y sus derivadas parciales continuas en una región del plano llamada D que contiene un punto fijo (p,q) Si se cumple que Entonces la sucesión generada a partir de (p0,q0) ϵ D converge a (p,q) (| ∂ g i ∂ x |+| ∂ gi ∂ y |< 1 ) ∀( x,y )∈D pk+1=g1( pk ,qk ) qk+1=g2( pk ,qk ) k=0,1,. .. i=1,2 PUNTO FIJO ATRACTIVO (Caso bidimensional) (Caso tridimensional) | ∂ g1 ∂ x ( p,q ) |+| ∂ g1 ∂ y (p,q ) | < 1 | ∂ g2 ∂ x (p,q ) |+| ∂ g2 ∂ y ( p,q )| < 1 | ∂ g1 ∂ x ( p,q,r ) |+| ∂ g1 ∂ y (p,q,r ) |+| ∂ g1 ∂ z (p,q,r ) |< 1 | ∂ g2 ∂ x (p,q,r )|+| ∂ g2 ∂ y ( p,q,r ) |+| ∂ g2 ∂ z ( p,q,r ) | < 1 | ∂ g3 ∂ x (p,q,r )|+| ∂ g3 ∂ y ( p,q,r ) |+| ∂ g3 ∂ z ( p,q,r ) | < 1 Si alguna de las condiciones anteriores NO SE CUMPLE el proceso iterativo puede DIVERGER o CONVERGER. PUNTO FIJO REPULSIVO Si se cumple que: Entonces la iteración de punto fijo genera una sucesión divergente (Caso bidimensional) (Caso tridimensional) | ∂ g1 ∂ x ( p,q ) |+| ∂ g1 ∂ y (p,q ) | > 1 | ∂ g2 ∂ x (p,q ) |+| ∂ g2 ∂ y ( p,q )| > 1 | ∂ g1 ∂ x ( p,q,r ) |+| ∂ g1 ∂ y (p,q,r ) |+| ∂ g1 ∂ z (p,q,r ) |> 1 | ∂ g2 ∂ x (p,q,r )|+| ∂ g2 ∂ y ( p,q,r ) |+| ∂ g2 ∂ z ( p,q,r ) | > 1 | ∂ g3 ∂ x (p,q,r )|+| ∂ g3 ∂ y ( p,q,r ) |+| ∂ g3 ∂ z ( p,q,r ) | > 1 Método iterativo de Gauss-Seidel pk+1=g1( pk ,qk ) qk+1=g2( pk+1 ,qk ) k=0,1, . .. En forma análoga al método iterativo de Gauss-Seidel para sistemas de ecuaciones lineales podemos proponer una mejora al método iterativo de Punto Fijo para sistemas de ecuaciones no lineales. PASO 1 Evaluar la función Pk+1=g (Pk )=[ g1(Pk ) g2(Pk ) . . . gn(Pk ) ] PASO 2 Sigo? Si, ‖Pk+1−Pk‖max< tol ‖Pk+1−Pk‖max ‖Pk+1‖max+eps < tol detengo el proceso y retorno el valor de Pk+1 como respuesta ó PASO 3 k k + 1; vuelvo a PASO 1 PuntoFijo.mMETODO ITERATIVO DE PUNTO FIJO IDEA ALGORITMO Resolver el sistema no lineal F (x) = 0 mediante una sucesión {Pk} que converge a la solución, a partir de un punto inicial P0. Evaluar la sucesión Xk+1= G(Xk) EJEMPLO f 1 ( x )=x12−2 x1−x 2+0 .5=0 f 2 ( x )=x12+4 x 22−4=0 Dado el siguiente sistema de 2 ecuaciones no lineales resolver utilizando el método de Newton-Raphson F (X )=[ x 12 −2 x1−x2+0 .5 x 12 +4 x 22 −4 ]; J (X )=[ 2 x1−2 −1 2 x1 8 x2 ] F (2 .00,0 .25 )=[0 .250 .25 ]; J (2.00,0 .25 )=[ 2 .0 −1 .0 4 .0 2 .0 ] Tomando como estimación inicial (2.00, 0.25) [2 .0 −1 . 04 .0 2 .0 ] [ ΔPP1ΔPP2 ]=−[ 0. 25 0 .25 ] ΔPP=[ ΔPP1ΔPP2 ]=[ −0 .09375 0 .0625 ] Resolviendo el sistema lineal obtenemos donde Así el punto siguiente es P1=P0+ΔPP=[2 .000 .25 ]+[ −0 .09375 0 .0625 ]=[ 1.90625 0 .3125 ] de manera similar P2=[1 .9006910 .311213 ] P3=[1.9006770 .311219 ] MÉTODO DE NEWTON-RAPHSON PARA SISTEMAS DE ECUACIONES NO LINEALES IDEA ALGORITMO Resolver el sistema no lineal F (x) = 0 mediante una sucesión {x(k)} que esperamos converja a la solución, a partir de un punto inicial x(0). La sucesión {x(k)} es obtenida linealizando F(x). Evaluar la sucesión PASO 1 Evaluar la función vectorial PASO 2 Evaluar la matriz Jacobiana PASO 3 Cálculo el siguiente punto x(k+1) resolviendo el SISTEMA LINEAL (e.g. para k=0): PASO 4 detengo el proceso y retorno el valor de x(k+1) como respuesta De lo contrario, vuelvo al paso 1 hasta que converge o hasta alcanzar un cierto número máximo de iteraciones permitidas Sigo ? Si, || x(k+1)-x(k) || max < tol FIN NOTAS DE CLASE Tema: Sistemas de ecuaciones algebraicas Diapositiva 1 Diapositiva 2 Diapositiva 3 Diapositiva 4 Diapositiva 5 Diapositiva 6 Diapositiva 7 Diapositiva 8 Diapositiva 9 Diapositiva 10 Diapositiva 11 Diapositiva 12 Diapositiva 13 Diapositiva 14 Diapositiva15 Diapositiva 16 Diapositiva 17 Diapositiva 18 Diapositiva 19 Diapositiva 20 Diapositiva 21 Diapositiva 22 Diapositiva 45 Diapositiva 46 Diapositiva 47 Diapositiva 48 Diapositiva 49 Diapositiva 50 Diapositiva 51 Diapositiva 52 Diapositiva 53 Diapositiva 54 Diapositiva 55 Diapositiva 56 Diapositiva 57 Diapositiva 58 Diapositiva 59 Diapositiva 60 Diapositiva 61 Diapositiva 62 Diapositiva 63 Diapositiva 64 Diapositiva 65 Diapositiva 66 Diapositiva 67 Diapositiva 68 Diapositiva 69 Diapositiva 70 Diapositiva 71 Diapositiva 72 Diapositiva 73 Diapositiva 74 Diapositiva 75 Diapositiva 76 Diapositiva 77 Diapositiva 78 Diapositiva 79 Diapositiva 80 Diapositiva 81 Diapositiva 82 Diapositiva 83 Diapositiva 87 Diapositiva 88 Diapositiva 89 Diapositiva 90 Diapositiva 91 Diapositiva 92 Diapositiva 93 Diapositiva 94 Diapositiva 95 Diapositiva 96
Compartir