Taller Sistemas De Ecuaciones

Taller de Métodos Numéricos Sistema de Ecuaciones AL Ing. Matilde Montealegre Madero, MSc. Fecha de envio: Febrero 11, 2019 Para realizar en clase de Febrero 14 de 2019

SISTEMAS DE ECUACIONES ALGEBRAICAS LINEALES A través de un sistema de ecuaciones algebraicas lineales se determinan los valores de X1, X2, X3,…. Xn que simultáneamente satisfacen una serie de ecuaciones tal que: f1 ( X 1 , X 2 , X 3 ,... X n ) = 0 f 2 ( X 1 , X 2 , X 3 ,... X n ) = 0 … f n ( X 1 , X 2 , X 3 ,... X n ) = 0 ; La forma general para representar ecuaciones algebraicas lineales es mediante una matriz cuadrada de coeficientes, [A], un vector de constantes, {C} y un vector de soluciones {X}: a11 X 1 + a12 X 2 + ...a1n X n = C1 a21 X 1 + a22 X 2 + ...a2 n X n = C2 an1 X 1 + an 2 X 2 + ...ann X n = Cn Luego, se resuelve para [A]{X}={C}. La convergencia de los métodos para solucionar ecuaciones algebraicas lineales se garantiza si los elementos de la diagonal principal son todos diferentes de cero y la matriz [A] es diagonalmente dominante es decir: (3.1) aii >  aij , para j =1,2,.. n excluyendo j = i.

NOTA IMPORTANTE: En caso que el sistema no cumpla con la condición de ser diagonalmente dominante, intercambiar las ecuaciones de tal modo que cumplan la condición de convergencia, como se indica a continuación: Para el siguiente sistema de ecuaciones

−5 x1 4 x1

− x2

6 x1

+8 x2

+12 x3 = 60

ecuación 1

− x3 = −2

ecuación 2

= 45

ecuación 3

presenta una matriz de coeficientes [A] que no cumple la condición de convergencia necesaria para dar solución por los métodos aquí descritos para solucionar ecuaciones algebraicas lineales. Un elemento de la diagonal principal es igual a cero y para j =1,2,.. n excluyendo j = i se tiene aii <  aij ; por tanto la matriz [A] no es diagonalmente dominante.

−5 0 12 ሾ‫ܣ‬ሿ = ൭ 4 −1 −1൱ 6 8 0 Para aplicar los métodos, se intercambian las ecuaciones de tal modo que cumplan la condición de convergencia.

4 x1

− x2

6 x1

+8 x2

−5 x1

− x3 = −2

ecuación 2

= 45

ecuación 3

+12 x3 = 60

ecuación 1

Quedando la matriz de coeficientes [A]

4 −1 −1 ሾ‫ܣ‬ሿ = ൭ 6 8 0൱ −5 0 12

1. MÉTODO DE JACOBI De acuerdo con la forma general para un sistema de n ecuaciones, se puede resolver para X1 usando la primera ecuación, para X2 usando la segunda hasta Xn usando la n-esima ecuación. Luego, se tiene: C1 − a12 X 2 − a13 X 3 − ... − a1n X n (3.2.a) a11 C − a21 X 1 − a23 X 3 − ... − a2 n X n (3.2.b) X2 = 2 a22 ….. C − an1 X 1 − an 2 X 2 − ... − ann −1 X n−1 (3.2.n) Xn = n ann Generalizando se obtiene:   n  1  X i =  Ci −  aij X j  (3.3) aii  j =1  j <>i   El proceso de solución se comienza asumiendo valores para las X’s. Una primera aproximación es asumiendo que todas las X’s son iguales a cero. Estos valores son sustituidos en las ecuaciónes (3.2a) a (3.2n) para obtener un nuevo valor para X1=C1/a11, X2=C2/a22,….,Xn=Cn/ann. De esta primera iteración se obtiene un primer aproximado (aproximación nueva) para las X’s que es utilizado en la siguiente iteración para obtener nuevos resultados. Cada nueva iteración se basa en los resultados de la anterior: X1 =

X inuevo = X ianterior

De donde:

nuevo

Xi

  n 1  anterior  = Ci −  aij X j  aii  j =1 j <> i  

(3.4)

Por cada iteración y para cada una de las X’s, se determina el error relativo aproximado utilizando la fórmula: X inuevo − X ianterior ε aprox (i ) = *100% (3.5) X inuevo El proceso se repite hasta que |εaprox | < εt.

2. MÉTODO DE GAUSS-SEIDEL Es el método de Jacobi modificado. Acelera la convergencia del método. Por cada iteración utiliza el valor de X, previamente calculado. Por tanto, de la fórmula generalizada (3.3) se obtienen los resultados de los nuevos valores de X’s usando   n 1  nuevo nuevo  (3.6) Xi = Ci −  aij X j  aii  j =1 j <>i   En la primera iteración se asume X2, X3, …Xn igual a cero para el cálculo de X1. Luego, para el cálculo de X2, se toma el valor calculado de X1 y los demás valores de X3, …Xn =0. Para el cálculo de X3, se usan los valores calculados de X1, X2 y los demás valores de X4, …Xn =0. Para el cálculo de Xn, todos los valores anteriores de X1, X2,….Xn-1 han sido calculados previamente, por tanto se utilizan estos valores. Terminada la primera iteración se calculan los errores para cada X, utilizando la fórmula dada por la Ec. 3.5 y se itera hasta que |εaprox | < εt .

3. MEJORAMIENTO DE LA CONVERGENCIA USANDO RELAJACIÓN Se basa en una pequeña modificación del Gauss-Seidel y está diseñado para mejorar la convergencia del método de Gauss-Seidel. Este método utiliza un factor de peso o de ponderación, λ, al cual se asigna un valor entre 0 y 2 para calcular: X inuevo = λ X inuevo + (1 − λ ) X ianterior

(3.7)

Sustituyendo la Ec. 3.7 en la fórmula de Gauss-Seidel, Ec. 3.6, se calcula para cada X usando:    n anterior para j > i , 1  nuevo   + (1 − λ ) X anterior nuevo para j i  El valor de λ, se determina de acuerdo con el problema físico especifico. Si λ = 1 , entonces 1 - λ = 0, por tanto el resultado es el método de Gauss-Seidel. Si 0 < λ < 1 , se tiene relajación por debajo, que es usada generalmente para hacer converger un sistema no convergente. Para valores de 1 < λ < 2, existe sobrerelajación, que es usada para acelerar la convergencia de un sistema convergente. Esto es típico de problemas que resultan de un proceso físico. Es muy importante seleccionar el valor adecuado de un λ para acelerar el proceso, pues un valor óptimo de dicho parámetro produce una convergencia extraordinariamente rápida. La estimación de este parámetro está en función del tamaño de la malla y de la heterogeneidad del sistema en estudio.

EJERCICIOS PROPUESTOS Pruebe con los sistemas siguientes:

I.

II.

III. Recuerde: Reorganizar cuando sea necesario para que la matriz sea Diagonalmente dominante. Prepare las subrutinas para solucionar sistemas de ecuaciones utilizando los métodos de Jacobi, Gauss-Seidel y Relajación. Cada método se debe ejecutar en la hoja activa que lleve su nombre usando la instrucción: Ejemplo: Sheets("jacobi").Activate para el caso del método Jacobi.

Utilice la instrucción de borrado para cada método dependiendo del tamaño del sistema de ecuaciones (n) haciendo llamado a la subrutina borrado() Sub borrado() Dim n As Integer, rango As Range, i As Integer 'n tamano de la matriz A n = Cells(1, 2).Value 'ciclo for para borrar 10 columnas a partir de la columna n+4 For i = 1 To 10

'selecciona la columna Range("a3" & ":c" & n + 2).Columns(n + 3 + i).Select 'asigna a la variable tipo range lo seleccionado Set rango = Selection Selection.ClearContents 'limpia la columna seleccionada Next i End Sub Organice la información en el libro de Excel tal como aparece en la siguiente imagen y Utilice los botones para el llamado de cada uno de los métodos.

A. COMO COLOCAR BOTONES? Renombre la hoja 1 de Excel como “llamadometodos” y coloque los botones para ejecutar cada método siguiendo los pasos 1, 2 y 3 como se muestra en la imagen.

B. Haga el llamado desde cada botón al método que corresponda como sigue: (DOBLE CLICK sobre el botón en modo de diseño). Private Sub cmdgaussseidel_Click() Call GAUSSSEIDEL End Sub Private Sub cmdjacobi_Click() Call JACOBI End Sub Private Sub Relajacion_Click() Call RELAJACION End Sub

C.  Por cada uno de los métodos determine: 1. Numero de iteraciones 2. Resultados por cada iteración. Imprima resultados como se indicó en clase.

Para el método de relajación ejecute el programa para diferentes valores de λ . Recuerde que 0 < λ < 2 . 1. Grafique λ vs Numero de iteraciones 2. Seleccione el mejor λ . Analice resultados.

4. ALGORITMO DE THOMAS Este se usa continuamente para resolver sistemas de ecuaciones lineales donde la matriz resultante es una matriz tridiagonal. El algoritmo de Thomas es básicamente una descomposición LU. Es un método directo, es decir, no requiere iteraciones. Cuando se modela el flujo de calor o de algunos fluidos en una dimensión donde se tiene un sistema de n ecuaciones con n incógnitas a saber: X1, X2,... Xn, el sistema resultante se puede expresar en la forma matricial dada por: Tu = d (3.9)  b1 c1  u1   d1        a2 b2 c2  u2   d 2    u3   d3  a3 b3 c3   =  . .   .   .   an −1 bn −1 cn −1  un −1   d n −1       an bn    un   d n  Aplicando descomposición LU: T = T LTu (3.10) Substituyendo la Ec. (3.10) en la Ec. (3.9); TL y= d siendo: y=Tuu Paso 1: Factorización Después de aplicar la eliminación Gaussiana, las matrices TL y Tu se definen como: 0 0 0  1 β1 0 0 0 0   α1 0 0 0 1 β   0 0 0  0 0 0 2   a2 α 2 0   0 1 β3 0 0  a3 α 3 0 0 0 Tu =  TL =    . . . .  . . . .      0 1 β n −1  an −1 α n −1 0      0 1  an α n    Los coeficientes α y ß, se obtienen de la siguiente forma: c β1 = 1 α 1 = b1

α1

Para i = 2, 3, 4, ..., n

βi =

c i −1

α i −1

α i = bi - ai β i

Paso 2: Sustitución Progresiva y1 = d 1 / α 1

y i = (d i - a i yi −1 ) /α i

Para i = 2, 3, 4, ..., n

Paso 3: Solución regresiva un = yn Puesto que un se conoce, se puede ahora hallar un-1, un-2, un-3,…., u2, u1

u i = y i - β i +1u i+1

Para i = n-1, n-2,..., 1.

Algoritmo Thomas() matrizA() ‘Matriz de coeficientes a() ‘vector de términos por debajo de la diagonal b() ‘vector de términos de la diagonal principal c() ‘vector de terminos por encima de la diagonal d() ‘vector de términos independientes u() ‘vector resultados INICIO ‘entrada de datos a la matriz y asignación de valores para a(),b(),c() y d() DatosThomas () Calculathomas() ‘realiza los caclculos para obtener el vector resultados u() ImprimirThomas() ‘muestra resultados FIN Option Explicit Sub Thomas() Dim matrizA() As Double, a() As Double, b() As Double, c() As Double Dim d() As Double, n As Integer, u() As Double n = Cells(6, 3).Value 'redimensiona las variables acorde con el tamaño del sistema ReDim Preserve matrizA(n, n) ReDim Preserve a(n) ReDim Preserve b(n) ReDim Preserve c(n) ReDim Preserve d(n) ReDim Preserve u(n) Call DatosThomas(matrizA, a, b, c, d, n, u) Call CalculaThomas(a, b, c, u, d, n) Call ImprimirThomas(matrizA, a, b, c, d, n, u) End Sub Sub DatosThomas(matrizA, a, b, c, d, n, u) 'Entrada de datos e inicialización de variables Dim i As Integer, j As Integer For i = 1 To n For j = 1 To n matrizA(i, j) = Cells(8 + i, j).Value Next j d(i) = Cells(8 + i, 5).Value u(i) = 0 a(i) = 0 b(i) = 0 c(i) = 0 Next i 'El vector a esta por debajo de la diagonal luego se cumple j=i-1 'El vector b es la diagonal luego, i es igual a j 'El vector c esta por encima de la diagonal, luego se cumple j=i+1 'De la matrizA, se coloca la informacion en los vectores a, b, c For i = 1 To n

For j = 1 To n If i = j Then b(i) = matrizA(i, j) 'vector diagonal End If If j = i + 1 Then c(i) = matrizA(i, j) 'vector por encima de la diagonal End If If j = i - 1 Then a(j + 1) = matrizA(i, j) 'vector por debajo de la diagonal a partir de la posición a(2) End If Next j Next i End Sub Sub CalculaThomas(a, b, c, u, d, n) Dim alpha() As Double, beta() As Double, y() As Double Dim i As Integer ReDim Preserve alpha(n) ReDim Preserve beta(n) ReDim Preserve y(n) 'Factorización y sustitución regresiva alpha(1) = b(1) beta(1) = c(1) / alpha(1) y(1) = d(1) / alpha(1) For i = 2 To n beta(i) = c(i - 1) / alpha(i - 1) alpha(i) = b(i) - (a(i) * beta(i)) y(i) = (d(i) - (a(i) * y(i - 1))) / alpha(i) Next i 'solución Regresiva u(n) = y(n) For i = n - 1 To 1 Step -1 u(i) = y(i) - (beta(i + 1) * u(i + 1)) Next i End Sub Sub ImprimirThomas(matrizA, a, b, c, d, n, u) Dim i As Integer, j As Integer For i = 1 To n For j = 1 To n Cells(21 + i, j).Value = matrizA(i, j) Next j Cells(21 + i, 5).Value = d(i) Next i For i = 1 To n Cells(14, i).Value = a(i) Cells(15, i).Value = b(i) Cells(16, i).Value = c(i) Cells(17, i).Value = d(i) Next i For i = 1 To n Cells(6 + i, 8).Value = u(i)

Next i End Sub Resolver usando el siguiente Ejemplo: Aplicación de Leyes de Kirchhoff para obtener el sistema de ecuaciones Considere la malla que consiste de nueve resistencias y una batería, como se muestra en la figura.

Obtener el sistema de ecuaciones lineal y resolver por el algoritmo de Thomas. Para las corrientes i1, i2 , i3 , dadas las resistencias y la batería r1 = 10, r2 = 10, r3 = 10, r4 = 10, r5 = 20, r6 = 20, r7 = 10, r8 = 10, r9 = 10, e1 = 36 Al aplicar las leyes de Kirchoff se tiene:

(r1 + r4 + r7 )i1 (−r4 )i1

− r4 i2

= e1

(r4 + r2 + r5 + r8 )i2 − r5i2

(r5 + r3 + r6 + r9 )i3 = 0

Reemplazando:

Dado que la matrizA resulta una matriz tridiagonal se resuelve por el algoritmo de Thomas cuya solución es i1 = 1.3 , i2 = 0.3 , i3 = 0.1

−10i2

= 36

−10i1 50i2

−20i3 = 0

30i1

−20i2

−r5i3 = 0

+60i3 = 0

Los vectores usados en el cálculo son: vector a vector b vector c vector d

0 30 -10 36

-10 50 -20 0

-20 60 0 0

Ejemplo 2: Compruebe el algoritmo de Thomas con el siguiente problema:  1 2 0  u1   3        3 4 5  u1  =  2        0 6 7  u1   1 

u1 = 0.545  Solucion = u2 = 1.227  u3 = -0.909