Laboratorio de Computación Científica – Curso 2022/2023

Práctica 3. Sistemas de Ecuaciones.

Todos los scripts y funciones se copiarán en un Word de resultados (o en un Live Script de
Matlab), atendiendo a la Sección y número del ejercicio en el que nos encontremos.
Para aquellos scripts o funciones que sean llamados desde consola de comandos se copiará,
además, la llamada realizada, así como el resultado generado, tanto numérico como figuras.

1. Matlab incorpora funciones para conocer el condicionamiento de un sistema de ecuaciones

a través de las funciones intrínsecas cond, rcond o condest.

Haciendo uso de las tres funciones anteriores, escribe una función llamada condicion.m que
evalúe los valores de los tres números de condición, el determinante de la matriz y su inversa,
así como la solución del sistema (aplica el operador \). Esta función tendrá como variables de
entrada la matriz A del sistema de ecuaciones y el vector independiente b.

a) Aplícalo al siguiente sistema de ecuaciones:

0. 832 𝑥1 − 0. 448 𝑥2 = 1

0. 784 𝑥1 + 0. 421 𝑥2 = 0

b) Haz lo mismo con el siguiente sistema de ecuaciones:

A = [1.985, - 1.358; 0.953, -0.652]; b = [2.212; 1.062].

¿Qué pasa si se modifica b(2) por 1.061? ¿Por qué?

c) Haz lo mismo con el siguiente sistema de ecuaciones: A = [1 2; 2 4], b = [0; 1].

2. Crea una función que permita resolver un sistema de ecuaciones triangular inferior por el
método de sustituciones progresivas. Llama a la función sp.m y diséñala de tal modo que
reciba como parámetros de entrada la matriz de coeficientes A y un vector columna que
contenga los términos independientes b, y devuelva como salida un vector columna con las
soluciones del sistema x. Esto es A x = b, siendo A triangular inferior. Aplica la función creada
para obtener la solución del sistema:

𝑥1 = 1

0. 5𝑥1 – 2𝑥2 = – 3. 5

– 𝑥1 + 0. 5𝑥2 + 3𝑥3 = 9

3. Matlab suministra una función propia lu.m para el cálculo de la factorización LU, con pivoteo
de filas. La sintaxis es, [L,U,P]=lu(A); donde L es una matriz triangular inferior, U una matriz
triangular superior y P es la matriz de permutaciones que da cuenta del cambio de orden
aplicado a las ecuaciones del sistema, P*A=LU. Aplícalo para resolver el siguiente sistema de
ecuaciones

0. 5𝑥1 + 2𝑥2 + 2. 5𝑥3 + 3𝑥4 = 24

0. 7𝑥1 + 5. 2𝑥2 − 3𝑥3 + 𝑥4 = 6. 1

1
 Laboratorio de Computación Científica – Curso 2022/2023

0. 8𝑥1 − 6𝑥2 + 3. 4𝑥3 − 2𝑥4 =− 9

2. 1𝑥1 + 3. 2𝑥2 − 4. 5𝑥3 + 2. 3𝑥4 = 4. 2

4. Sean las siguientes sentencias que determinan las matrices estrictamente inferior L,
estrictamente superior U y diagonal D de una matriz A:

L=A-triu(A)

U=A-tril(A)

D=diag(diag(A))

En el método de Jacobi la expresión matricial para determinar el vector solución x viene dada
por:

x=inv (D)*(b-(L+U) *x0) siendo b el vector independiente y x0 el vector solución inicial.

Construye una función jacobi1.m que determine la primera solución x, partiendo de (0,0,0)
como aproximación inicial de la solución. Determina la primera solución del sistema de
ecuaciones:

3𝑥 + 𝑦 + 𝑧 = 4

2𝑥 + 5𝑦 + 1𝑧 = −1

− 𝑥 + 𝑦 + 3𝑧 = 4

5. Escribe una función jacobi.m, con estructura x=jacobi(A, b, x0, tol), de modo que resuelva
un sistema de ecuaciones usando el método de Jacobi.

a) Aplícalo al sistema del ejercicio anterior hasta conseguir un error inferior (usa el comando
Matlab norm) a la tolerancia tol y partiendo de un valor x0 como aproximación inicial a la
solución. Usa el comando while y la forma matricial del método de Jacobi. Toma una tolerancia
de 0.001. Las sucesivas soluciones se deben ver en pantalla en forma de tabla, tal que, en la
primera columna aparezca la iteración, en las tres siguientes columnas los elementos del vector
solución x correspondientes a esa iteración s y en la quinta columna el error entre dos
soluciones consecutivas x(s) y x(s-1). Recuerda usar el comando norm para esto último.

b) Si aplicas tu programa jacobi.m, al sistema de ecuaciones:

𝑥 + 2𝑦 + 3𝑧 = 4
4𝑥 + 5𝑦 + 6𝑧 =− 1

7𝑥 + 8𝑦 + 9𝑧 = 4.

partiendo de (0,0,0) como aproximación inicial de la solución verás que Jacobi no converge ya
que la matriz del sistema es singular. Para detectar estos casos, añade una condición a tu
función jacobi.m para que cuando el determinante de la matriz del sistema sea menor que
3*eps nos diga que no hay solución y se corte la ejecución del programa, por ejemplo, con la
orden return.

2
 Laboratorio de Computación Científica – Curso 2022/2023

c) Aplica ahora tu nuevo programa de Jacobi al siguiente sistema

2𝑥 − 2𝑦 = 4

2𝑥 + 3𝑦 − 𝑧 =− 1

− 5𝑥 + 2𝑧 = 4

Como verás, el método no converge aunque la matriz del sistema tenga determinante no nulo.

d) En los métodos iterativos es posible previamente determinar la convergencia del método

estudiando las propiedades de la matriz de convergencia, que relaciona dos soluciones
consecutivas. Ello puede hacerse incorporando al programa jacobi.m las siguientes sentencias:
>> J=-inv(D)*(L+U); %J es la matriz de convergencia del método de Jacobi, que relaciona dos soluciones
consecutivas

>>autovalores=eig(J); % eig calcula los autovalores de una matriz (lo verás en álgebra)

>>abs_auto=abs(eig(J));

>>radio_espectral=max(abs(eig(J)))

>>if radio_espectral>1

>> disp('El metodo no converge. No da la solución del sistema')

>> return

>>end

Escribe tu programa final jacobi.m, las ejecuciones para los tres sistemas de ecuaciones y los
resultados.

6. Escribir una función x=jacobi_amor(A,b,x0,w,tol) que resuelva un sistema de ecuaciones por

el método iterativo Jacobi amortiguado, hasta conseguir un error de convergencia inferior a la
tolerancia tol, partiendo de un valor x0 como aproximación inicial a la solución y w como peso
(factor de amortiguamiento).

Nota: Con añadir en tu programa jacobi.m la sentencia de amortiguamiento tras la sentencia

que define el método de Jacobi es suficiente. Además, debe tenerse en cuenta que la matriz
de convergencia cambia y ya no es la misma que la del método de Jacobi. ¿Cuál será ahora la
matriz de convergencia?

a) Aplica tu función al siguiente sistema de ecuaciones con una tolerancia de 0.001, un peso de
0.9 y vector solución inicial [0;0;0]. Te tiene que salir un radio espectral de 0.9463, por tanto el
método converge.

2𝑥 − 2𝑦 = 4

2𝑥 + 3𝑦 − 𝑧 =− 1

− 5𝑥 + 2𝑧 = 4

3
 Laboratorio de Computación Científica – Curso 2022/2023

b) Si la solución inicial fuera x0=[3; 1; 10], ¿cuántas iteraciones se realizan hasta llegar a la
solución deseada? ¿Cuál es tu conclusión?

c) Si el peso valiera 1.1, ¿cuál será el valor del radio espectral para el sistema de ecuaciones
anterior?

7. Resolución de sistemas de ecuaciones por Gauss-Seidel.

a) Escribe una función de la forma xs1 = gseidel1(A,b,xs) que realice una iteración del método
de Gauss-Seidel y aplícalo al siguiente sistema de ecuaciones con xs = [0; 0; 0]:

2𝑥 − 2𝑦 = 4

2𝑥 + 3𝑦 − 𝑧 =− 1

5𝑥 + 2𝑧 = 4.

b) Escribe una función gseidel.m que resuelva de manera iterativa el sistema de ecuaciones
anterior por el método de Gauss-Seidel y que llame a la función escrita en el apartado anterior.
Considera un valor de tolerancia para la solución de 0.001. No olvides incluir las
consideraciones sobre el valor del determinante de A y de la matriz de convergencia.

4
 Laboratorio de Computación Científica – Curso 2022/2023

EJERCICIO APLICADO. Cálculo de tensiones en una estructura.

Determinar las fuerzas estáticas que aparecen en la siguiente estructura. La figura muestra la
numeración de los nodos, fuerzas (𝐹𝑖𝑗 =− 𝐹𝑗𝑖) y reacciones en los puntos de apoyo (𝑅𝑖). Se

somete a unas cargas de 500 y 250 N en los nodos 4 y 5, según indican las flechas que
representan dichas cargas.

A partir de esta imagen se pueden plantear fácilmente las ecuaciones de equilibrio en cada
nodo (izquierda), lo que da un sistema de 10 ecuaciones y 10 incógnitas (𝐹𝑖𝑗, 𝑅𝑖𝑥, 𝑅𝑖𝑦), los
términos independientes corresponden a las cargas externas en los nodos 4 y 5. Este sistema se
puede modelar en forma matricial (derecha), donde A es la matriz de coeficientes:
fila=ecuación y columna=fuerza (variable a calcular).

En otras palabras, este problema se puede representar de forma matricial A·x=b, donde:
● A es una matriz que viene determinada por la estructura del puente.
● b es un vector columna que representar las cargas a las que está sometido el puente.
● x es un vector columna que incluye las tensiones y reacciones que se producen en el
puente [R1x; R1y; R3y; F12; F23; F14; F24; F25; F35; F45].

Ejercicio:

a) Comprobar que la matriz A está bien construida utilizando situaciones de carga que
generen simetría y sean fáciles de interpretar: carga vertical en nodo 2 o carga horizontal
en nodo 3.

5
 Laboratorio de Computación Científica – Curso 2022/2023

Utilizar el operador que proporciona Matlab para resolver sistemas e interpretar los
resultados de estas situaciones. Finalmente resolver el problema para el vector de cargas
propuesto en la figura.
b) Resolver el sistema por descomposición LU, utilizando las funciones lu, sp y sr.

Para resolver el problema utilizando métodos iterativos es necesario reordenar el sistema

para que sea diagonal dominante, para ello se puede utilizar la matriz P del apartado b). Si
no se reordena, los métodos no convergen.

c) Representar la evolución del radio espectral del método de Jacobi Amortiguado frente al
peso.

d) Resolver el sistema por Jacobi Amortiguado, con peso 0.3 y tolerancia 1, partiendo desde
la solución cero.

e) ¿Cuándo interesará utilizar Jacobi Amortiguado frente a LU? Partiendo de la solución

anterior, vuelve a resolver el sistema para una pequeña variación de la carga y compara los
tiempos.