Análisis Numérico I 4° Semestre

Matemáticas

Análisis Numérico I

4° Semestre

Unidad 3. Sistemas de ecuaciones lineales

Clave

050920624/ 060920624

Universidad Abierta y a Distancia de México


UnADM | DCEIT | MT | MANU1 2

Índice

Unidad 3. Sistemas de ecuaciones lineales ............................................................................................................ 3

Presentación de la unidad .................................................................................................................................... 3

Competencia específica ........................................................................................................................................ 3

Logros de la unidad .............................................................................................................................................. 3

3.1. Sistemas de ecuaciones lineales ...................................................................................................................... 4

3.1.1. Matrices........................................................................................................................................................... 6

3.1.2. Matrices y sistemas de ecuaciones lineales ................................................................................................... 14

3.2. Normas de vectores y de matrices ............................................................................................................... 19

3.2.1. Solución de un sistema de ecuaciones lineales ............................................................................................. 28

3.2.2. Eliminación de Gauss .................................................................................................................................... 29

3.2.3. Método iterativo de Jacobi ............................................................................................................................ 37

3.3. Método iterativo de Gauss-Seidel ................................................................................................................ 42

3.4. Regresión lineal ............................................................................................................................................ 46

3.4.1. El modelo de regresión lineal ....................................................................................................................... 47

3.4.2. El modelo de regresión lineal ....................................................................................................................... 49

Cierre de la unidad ............................................................................................................................................. 51

Para saber más: .................................................................................................................................................. 51

Referencias Bibliográficas .................................................................................................................................. 52

Figura 1. Números ......................................................................................................................................... 3

Figura 2: Visualización de las distintas normas aplicadas a los vectores que cumplen 𝒙𝒑 = 𝟏................. 21

Figura 3. Ejemplo regresión lineal a una función del tipo y=mx+b............................................................ 47

Tabla 1. Operaciones elementales con sistemas de ecuaciones. ................................................................. 16

Tabla 2. Iteraciones hasta la convergencia ................................................................................................. 40

Tabla 3. Iteraciones ..................................................................................................................................... 42






Presentación de la unidad

Una de las aplicaciones básicas de cualquier computadora (o algoritmo) es la de resolver algún

sistema de ecuaciones a partir de un conjunto de puntos observados y sus resultados.

Los sistemas de ecuaciones lineales son la primera representación de problemas más complejos

que sólo resolver una ecuación, ya que representan restricciones en más de una dimensión y es a

través de un sistema de ecuaciones en las que la modelación matemática tiene una conexión

directa con el mundo que observamos diariamente. Es decir, es una de las actividades

matemáticas en las que tenemos mucha práctica desde los primeros niveles escolares en donde

representamos algún problema cotidiano en lenguaje matemático. De esta forma aprenderás

métodos para resolver sistemática y eficientemente estos sistemas con algún método formal.

Competencia específica

Utilizar matrices para representar y solucionar sistemas de ecuaciones lineales mediante su

representación matricial.

Logros de la unidad

Figura 1. Números

Aplicar métodos para resolver sistemas de

ecuaciones en su representación matricial.



3.1. Sistemas de ecuaciones lineales

La modelación matemática implica representar un fenómeno observado en un lenguaje

matemático, para expresarlo, en una ecuación o en un conjunto de ecuaciones que representen el

comportamiento general del fenómeno observado.

La solución de ecuaciones de una sola variable se hace a través de métodos específicos para

encontrar sus raíces. Esas ecuaciones representan toda una gama de problemas y modelos

asociados a ellos, así como los métodos diseñados para resolverlos. Este tipo de problemas es

insuficiente para modelar fenómenos más complejos.

Una parte de los fenómenos observados en disciplinas como la ingeniería y la economía se

modelan a través de sistemas de ecuaciones lineales (sistemas lineales). Son lineales porque

todas las variables están expresadas en términos lineales (es decir, el exponente de todas las

variables es 1 y no hay multiplicaciones entre los términos).

Cabe mencionar que herramientas de modelación no lineal de fenómenos usan la Teoría de

sistema de ecuaciones lineales para aproximar problemas de esta naturaleza no lineal mediante

un equivalente lineal.

Considera el siguiente problema:

Dos veces la edad de Juan más la edad de Pedro es igual a 10.

Este problema así planteado está subdeterminado ya que tenemos 2 variables y sólo una

ecuación

2𝑥 + 𝑦 = 10

Si agregamos el hecho de que la edad de Juan más tres veces la de Daniel da 40, ahora tenemos 3

incógnitas y sólo dos ecuaciones.



2𝑥 + 𝑦 = 10

𝑥 + 3𝑧 = 40

Si consideramos que la edad de Juan, más la de Pedro más dos veces la de Daniel es 30 entonces

ya tenemos tres ecuaciones con tres incógnitas lo que nos determina un sistema de ecuaciones de

la siguiente forma:

Definiendo la edad de Juan como 𝑥, la de Pedro como 𝑦 y la de Daniel como 𝑧 tenemos que:

2𝑥 + 1𝑦 = 10

𝑥 + 3𝑧 = 40

𝑥 + 𝑦 + 2𝑧 = 30

Al resolver este sistema de ecuaciones despejando y sustituyendo variables obtenemos que la

solución sea:

𝑥 = 4, 𝑦 = 2, 𝑧 = 12

El determinar las edades de forma lúdica como se hace en este juego revela varios aspectos que

caracterizan el modelado de un fenómeno a través de un sistema lineal de ecuaciones. Por

ejemplo, nos muestra las relaciones que guardan las variables entre sí, que es una relación lineal.

Otro aspecto de un sistema de ecuaciones como los que revisaremos es que contiene toda la

información necesaria para determinar su solución. (Un caso que estudiaremos más adelante son

los problemas de regresión que son mucho más comunes en la práctica y tienen la característica

de ser sistemas sobredeterminados, lo que significa que existen más datos que podrían

determinar un resultado pero, llegado el momento, examinaremos la técnica más popular para

determinar su solución)

Diseñar un método que no dependa de analizar e interpretar el sistema para determinar qué

ecuación podemos empezar a despejar y sustituir es complicado, para evitar esto usaremos una



representación matricial de nuestros sistemas lineales y emplearemos el álgebra definida para

describir métodos que nos ayuden a resolver sistemas lineales.

Por ejemplo, el sistema anterior lo podemos expresar como

𝐴𝑥 = 𝑏

Para nuestro problema particular esta ecuación toma los siguientes valores:

[2 1 01 0 31 1 2

] [𝑥𝑦𝑧] = [

104030

]

Donde el vector solución es el vector columna 𝑥𝑡= (4,2,12). Por notación nos referiremos a las

matrices con letras en tipografía común (aún cuando también sean vectores) y a los vectores con

una línea encima del símbolo o con letras en negritas.

3.1.1. Matrices

Las matrices son arreglos bidimensionales de números en filas y columnas. Las entradas de estos

arreglos usualmente toman valores en ℝ pero bien podrían hacerlo en el campo de los complejos.

Una matriz de 𝑛 filas por 𝑚 columnas (𝑛 × 𝑚) toma la siguiente forma

𝐴 = |

𝑎1,1 … 𝑎1,𝑚

⋮ ⋱ ⋮𝑎𝑛,1 … 𝑎𝑛,𝑚

|

Donde la entrada en la fila 𝑖 y la columna 𝑗 de la matriz A la denotaremos por 𝑎𝑖𝑗.

La definición de una matriz de dimensión (𝑛 × 𝑚) es:

𝐴 = [𝑎𝑖𝑗]𝑛×𝑚 𝑡. 𝑞. 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚 (1)



El conjunto de las matrices es un campo vectorial sobre ℝ ya que con las operaciones binarias de

{+,∙} (suma y multiplicación por escalar) podemos comprobar que cumple las siguientes

propiedades:

Sean 𝑋, 𝑌 dos matrices de 𝑛 × 𝑚.

1) 𝑋 + 𝑌 = 𝑌 + 𝑋(Conmutatividad)

2) (𝑋 + 𝑌) + 𝑍 = 𝑋 + (𝑌 + 𝑍) (Asociatividad)

3) ∃0 𝑡. 𝑞. 0 + 𝑋 = 𝑋 + 0 (Elemento neutro aditivo)

4) 𝑋 − 𝑋 = 𝑋 + (−𝑋) = 0 (Elemento opuesto a la suma)

5) (𝑎 + 𝑏)𝑋 = 𝑎𝑋 + 𝑏𝑋 (Distributividad de los escalares)

6) 𝑎(𝑋 + 𝑌) = 𝑎𝑋 + 𝑎𝑌 (Distributividad sobre los vectores)

7) 𝑎(𝑏𝑋) = (𝑎𝑏)𝑋 (Asociatividad de los escalares).

. En esta sección repasaremos algunas propiedades y operaciones para las matrices que nos

servirán para como herramientas para resolver sistemas lineales.

Sean 𝐴 = [𝑎𝑖𝑗]𝑛𝑥𝑚, 𝐵 = [𝑏𝑖𝑗]𝑛𝑥𝑚

dos matrices y 𝛼 ∈ ℝ.

Igualdad: 𝐴, 𝐵 son iguales si y solo si son iguales elemento a elemento, es decir

𝐴 = 𝐵 ⟺ 𝑎𝑖𝑗 = 𝑏𝑖𝑗 𝑝𝑎𝑟𝑎 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚

Suma: La suma de 𝐴, 𝐵 se realiza también elemento a elemento, es decir

𝐴 + 𝐵 = [𝑎𝑖𝑗 + 𝑏𝑖𝑗]𝑛×𝑚 𝑝𝑎𝑟𝑎 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚

Multiplicación por escalar: La matriz que resulta de multiplicar 𝛼𝐴 se hace multiplicando cada

entrada de 𝐴 por el escalar 𝛼.

𝛼𝐴 = 𝛼 ∙ [𝑎𝑖𝑗]𝑛×𝑚= [𝛼 ∙ 𝑎𝑖𝑗]𝑛×𝑚

𝑝𝑎𝑟𝑎 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚



Combinación lineal: Para hacer la combinación lineal 𝐴 + 𝛼𝐵 se puede deducir de las

definiciones anteriores:

𝐴 + 𝛼𝐵 = [𝑎𝑖𝑗 + 𝛼𝑏𝑖𝑗]𝑛×𝑚 𝑝𝑎𝑟𝑎 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚

Multiplicación de matrices: Sean 𝐴 = [𝑎𝑖𝑘]𝑛×𝑝 y 𝐵 = [𝑏𝑘𝑗]𝑝×𝑚 para 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚 y

1 ≤ 𝑘 ≤ 𝑝.

Entonces la multiplicación 𝐶 = 𝐴𝐵 se define por

𝐶 = [𝑎𝑖𝑘]𝑛×𝑝 ∙ [𝑏𝑘𝑗]𝑝×𝑚= [𝑐𝑖𝑗]𝑛×𝑚

= ∑ 𝑎𝑖𝑘𝑏𝑘𝑗

𝑝

𝑘=1

= 𝑎𝑖1𝑏1𝑗 + 𝑎𝑖2𝑏2𝑗 + ⋯+ 𝑎𝑖𝑝𝑏𝑝𝑗

Donde 1 ≤ 𝑖 ≤ 𝑛, 1 ≤ 𝑗 ≤ 𝑚.

De la expresión anterior observamos que una matriz multiplicada nueva queda definida por la

cantidad de filas del primer factor y por la cantidad de columnas en el segundo factor donde cada

entrada en la matriz nueva es el producto punto de la fila 𝑖 − é𝑠𝑖𝑚𝑎 de la primera matriz por la

columna 𝑗 − é𝑠𝑖𝑚𝑎 de la segunda matriz.

Ejemplo: Si 𝐴 = [3 0 21 2 00 1 1

] y 𝐵 = [2 10 11 0

] entonces 𝐶 = 𝐴𝐵 = [8 32 31 1

].

Para fijarnos con detalle la entrada 𝑐2,1 = ∑ 𝑎𝑖𝑘𝑏𝑘𝑗3𝑘=1 = 𝑎𝑖1𝑏1𝑗 + 𝑎𝑖2𝑏2𝑗 + 𝑎𝑖3𝑏3𝑗 = 𝑎2,1𝑏1,1 +

𝑎2,2𝑏2,1 + 𝑎2,3𝑏3,1 = 1 ∙ 2 + 2 ∙ 0 + 0 ∙ 1 = (1,2,0) ∙ (2,0,1) = 2.

El resto de las entradas se calcula de la misma forma.

Elementos destacados: En el conjunto de las matrices existen un par de elementos que pueden

ser el neutro aditivo o la matriz 𝟎 que se define como una matriz con 0 en todas sus entradas

0 = [0𝑖𝑗]𝑛×𝑚



Y la matriz identidad (𝑰𝒏), que a diferencia del caso anterior sólo tiene 1 en su diagonal y por lo

tanto es una matriz cuadrada (de dimensión 𝑛 × 𝑛).

𝐼𝑛 = [1𝑖𝑖]𝑛×𝑛

Esta matriz cumple con:

𝐼𝑛𝐴 = 𝐴 para todas las matrices𝐴𝑛×𝑝

𝐵𝐼𝑛 = 𝐵 para todas las matrices 𝐵𝑚×𝑛

Es importante señalar que la conmutatividad de la multiplicación es posible únicamente para el

espacio de matrices cuadradas.

Diagonal y matrices triangulares: Sea 𝐴 = [𝑎𝑖𝑗]𝑛×𝑛 una matriz cuadrada de orden n.

Denominamos como la diagonal de 𝐴 (o traza de A) a todos los elementos que cumplen que 𝑖 =

𝑗.

𝑑𝑖𝑎𝑔(𝐴) = {𝑎1,1, 𝑎2,2, 𝑎3,3, … , 𝑎𝑛,𝑛}

Si sucede que todas las entradas de 𝐴 = [𝑎𝑖𝑗 = 0]𝑛×𝑛

tal que 𝑖 > 𝑗 entonces 𝐴 se denomina

como matriz triangular superior

𝐴 = [

𝑎1,1 … 𝑎𝑛,𝑛

0 𝑎2,2 ⋮

0 0 𝑎3,3

]

Y en el caso de que las entradas de 𝐴 = [𝑎𝑖𝑗 = 0]𝑛×𝑛

tal que 𝑖 ≤ 𝑗 entonces 𝐴 queda

identificada como matriz triangular inferior.

𝐴 =

[ 𝑎1,1 = 0 0 0 0

𝑎2,1

⋮

𝑎2,2 = 0𝑎𝑖𝑗

0⋱

00

𝑎𝑛,𝑛 … 𝑎𝑛,𝑛−1 𝑎𝑛,𝑛 = 0]



Invertibilidad: Decimos que 𝐴 es invertible si existe (recuerda que no todas las matrices tienen

inversa) 𝐵 tal que

𝐴𝐵 = 𝐼 (2)

Decimos que 𝐵 es la inversa de 𝐴.

También denotamos a la inversa de 𝐴 como 𝐴−1. Sabemos que:

(𝐴−1)−1 = 𝐴

Más aún si 𝐴 y 𝐵 son invertibles cuadradas y del mismo orden, entonces su producto es

invertible y se cumple que:

(𝐴𝐵)−1 = 𝐵−1𝐴−1 (3)

Estos aspectos del espacio vectorial de las matrices se han demostrado en el curso de Álgebra

lineal, en este apartado únicamente las recordamos para utilizarlas en los métodos que

desarrollaremos.

Transpuesta: La transpuesta de una matriz 𝐴 = [𝑎𝑖𝑗]𝑛×𝑚 se denota por 𝐴𝑡 y significa

intercambiar líneas por columnas, entonces la transpuesta nos queda como:

𝐴𝑇 = [𝑎𝑗𝑖]𝑚×𝑛

Decimos que una matriz 𝐴 es singular si cumple alguna de las siguientes propiedades

equivalentes:

1. 𝐴 no tiene inversa, es decir, no existe matriz 𝑀 tal que 𝐴𝑀 = 𝑀𝐴 = 𝐼.

2. det(𝐴) = 0

3. rango(𝐴) < 𝑛 (el rango es el máximo número de filas o columnas linealmente

independientes que contiene la matriz)

4. 𝐴𝑧 = 0 p.a. vector z≠0



Decimos que una matriz es no singular en caso contrario.

Estas propiedades son fáciles de demostrar por lo que desarrollarás sus respectivas

demostraciones. Es necesario que tengas presente la definición de singularidad (y no

singularidad) ya que haremos referencia a ella continuamente a lo largo del texto.

En Octave representaremos matrices de la siguiente forma:

octave:1> A=[2 1 0; 1 0 3; 1 1 2]

A =

2 1 0

1 0 3

1 1 2

Es decir, usaremos corchetes para agrupar las entradas y el símbolo “;” (punto y coma) para

especificar donde termina un vector fila. La transpuesta de una matriz se obtiene con el operador

“’”. Por ejemplo

octave:2>B=A’

B =

2 1 1

1 0 1

0 3 2

La suma de matrices y multiplicación de escalares es intuitiva, se suman dos operandos siempre

y cuando cumplan con la definición, en cuanto a la multiplicación se aplica en el escalar que hay

que multiplicar por alguna matriz.

octave:3> C=A+B



C =

4 2 1

2 0 4

1 4 4

octave:4> D=4*A

D =

8 4 0

4 0 12

4 4 8

Para obtener la diagonal de una función usaremos la función diag que toma 1 o 2 argumentos. El

primero siempre es la matriz sobre la cual queremos obtener la diagonal, al no especificar ningún

otro parámetro entonces lo que obtenemos es 𝑡𝑟(𝐴).

octave:5>diag(D)

ans =

8

0

8

Pero si en el segundo parámetro especificamos algún entero ±𝑘 tal que 1 ≤ 𝑘 ≤ 𝑑𝑖𝑚 − 1 lo que

obtenemos son los elementos de 𝐷 tal que 𝑑𝑖𝑎𝑔(𝐷, 𝑘) = [𝑑𝑖,𝑖+𝑘] en el caso de que 𝑘 > 0, en el

caso de que 𝑘 < 0 lo que obtenemos es 𝑑𝑖𝑎𝑔(𝐷, 𝑘) = [𝑑𝑖+𝑘,𝑖].

octave:6>diag(D,1)

ans =

4

12

octave:10>diag(D,-2)

ans = 4



La multiplicación también tiene una aplicación directa con el operador “*” (asterisco) siempre y

cuando los operadores estén bien definidos. Por ejemplo.

octave:5> A=[1 1 3; 2 3 1]

A =

1 1 3

2 3 1

octave:6> B=[5 1; 2 0; 1 3]

B =

5 1

2 0

1 3

octave:8> A*B

ans =

10 10

17 5

Pero una diferencia es que Octave permite operaciones puntuales, es decir, tiene algunas

definiciones adaptadas para operar con matrices entrada a entrada siempre y cuando ambas

matrices tengan la misma dimensión. Por ejemplo, como 𝐶 y 𝐷 son ambas matrices de (3 × 3)

podemos realizar la suma puntual de ambas.

octave:9> C.+D

ans =

12 6 1



6 0 16

5 8 12

Como puedes ver es muy similar a la suma excepto porque antepusimos un punto al operador

suma. Lo mismo pasa con la multiplicación, exponenciación y división.

octave:10> C.*D

ans =

32 8 0

8 0 48

4 16 32

octave:11> C./D

ans =

0.50 0.50Inf

0.50 NaN 0.33

0.25 1.00 0.50

octave:12> C.^D

ans =

65536 16 1

16 1 16777216

1 256 65536

Es para el uso y operación de matrices para lo que se diseñaron Octave y MATLAB.

3.1.2. Matrices y sistemas de ecuaciones lineales

En este tema, utilizaremos la Teoría del álgebra de matrices para resolver nuestros sistemas de

ecuaciones. Retomando nuestro ejemplo anterior, el sistema de ecuaciones (sin el planteamiento)

es el siguiente:

(1) 2𝑥 + 1𝑦 = 10

(2) 𝑥 + 3𝑧 = 40

(3) 𝑥 + 𝑦 + 2𝑧 = 30

Cada una de estas ecuaciones plantea algún conjunto de restricciones observadas directamente

sobre el problema, para que quede más claro, se han enumeradolas ecuaciones para referirme a

ellas en esta subsección. Se observa que la ecuación más difícil para obtener alguna variable

despejada es la (3) por lo que es poco recomendable empezar a despejar cualquier variable a



partir de de esa ecuación. En este caso empezaremos por despejar 𝑧 de la ecuación (2) e 𝑦 de la

ecuación (1) para obtener las siguientes dos ecuaciones nuevas

(4) 𝑦 = 10 − 2𝑥

(5) 𝑧 =40 − 𝑥

3

Así 𝑦 y 𝑧 dependen exclusivamente de 𝑥. Ahora usaremos estas dos ecuaciones en (3) para

obtener

(6) 𝑥 + 10 − 2𝑥 +80 − 2𝑥

3= 30

Para eliminar el denominador 3 que aparece en (6) multiplicamos por 3 toda la ecuación ya que

esto la mantiene inalterada

(7) 3𝑥 + 30 − 6𝑥 + 80 − 2𝑥 = 90

Esta ecuación depende únicamente de 𝑥 por lo que el despeje es inmediato.

(8) 𝑥 =20

5= 4

Sustituyendo (8) en (4) y (5) obtenemos

(9) 𝑦 = 2

(10) 𝑧 = 12

El proceso realizado se basa en nuestra experiencia y capacidad informática que tenemos. De esa

forma evaluamos cuál de las primeras ecuaciones es la más adecuada para empezar nuestros

despejes, (desarrollar este proceso en una computadora es todavía imposible, aunque cada vez

más cercano), aunque es probable que exista algún algoritmo discernible para obtener

automáticamente la solución.

Para resolver sistemas lineales como el que acabamos de presentar se sugiere tomar en cuenta:

a) Las variables consideradas serán todas aquellas que jueguen algún papel en



cualquiera de las ecuaciones del sistema. Si alguna variable no está mencionada en

alguna ecuación se considerará que su coeficiente es 0.

b) Las ecuaciones pueden intercambiarse libremente. Es claro que no importa en qué

orden leamos las ecuaciones éstas representan el mismo problema.

c) Cualquier ecuación puede ser multiplicada por algún 𝛼 ∈ ℝ, 𝛼 ≠ 0. Esto nos

reemplaza la ecuación en una ecuación equivalente pero no la anula (anular una

ecuación significaría tener un sistema distinto que describa otra cosa).

d) Podemos reemplazar la ecuación por una suma entre ella y cualquier multiplicación

por un escalar de cualquier otra ecuación en el sistema. Esta operación tampoco

altera el resultado ya que sólo opera con las restricciones impuestas por el sistema

en cada una de sus dimensiones.

Tabla 1. Operaciones elementales con sistemas de ecuaciones.

La forma en que resolvimos nuestro ejemplo fue desarrollando estas operaciones. Ahora

consideremos otro sistema:

𝑒𝑐1: 𝑤 + 𝑥 + 3𝑧 = 4

𝑒𝑐2: 2𝑤 + 𝑥 − 𝑦 + 𝑧 = 1

𝑒𝑐3: 3𝑤 − 𝑥 − 𝑦 + 2𝑧 = −3

𝑒𝑐4: − 𝑤 + 2𝑥 + 3𝑦 − 𝑧 = 4

Para resolver este sistema de ecuaciones vamos a aplicar las reglas de la Tabla 1:

• Primero sumaremos a la 𝑒𝑐 2 el producto escalar −2 ∙ 𝑒𝑐1. Esto significa que a la 𝑒𝑐1 la

multiplicaremos por el escalar -2:

−2 ∙ 𝑒𝑐1 = −2𝑤 − 2𝑥 − 6𝑧 = −8

• Al sumarla término a término con 𝑒𝑐2 nos queda:

𝑒𝑐2∗: 0𝑤 − 𝑥 − 𝑦 − 5𝑧 = −7

• Para eliminar el término 𝑤 de 𝑒𝑐3 podemos sumarle la ecuación uno multiplicada por el

escalar 3, es decir sumarle −3 ∙ 𝑒𝑐1

−3 ∙ 𝑒𝑐1 = −3𝑤 − 3𝑥 − 9𝑧 = −12

• Al sumarla con 𝑒𝑐3 obtenemos:



𝑒𝑐3∗: 0𝑤 − 4𝑥 − 𝑦 − 7𝑧 = −15

• Y finalmente si a la ecuación 4 le sumamos la ecuación 1 obtenemos:

𝑒𝑐4∗: 0𝑤 + 3𝑥 + 3𝑦 + 2𝑧 = 8

De esta forma hemos eliminado todas las variables 𝑤 de todas las ecuaciones después de la

ecuación 1.

Aplicaremos un procedimiento similar para la variable 𝑥 en todas las ecuaciones después de la

ecuación 2.

Los asteriscos usados en las ecuaciones anteriores tenían la función de distinguirlas de la

ecuación original indicando que son equivalentes, pero no idénticas,

Para no sobrecargar la notación cada ecuación se llamará como la ecuación lo indique (sin el

asterisco).

Nuestro sistema equivalente ahora tiene la siguiente forma:

𝑒𝑐1: 𝑤 + 𝑥 + 3𝑧 = 4

𝑒𝑐2 : − 𝑥 − 𝑦 − 5𝑧 = −7

𝑒𝑐3 : − 4𝑥 − 𝑦 − 7𝑧 = −15

𝑒𝑐4: 3𝑥 + 3𝑦 + 2𝑧 = 8

Para eliminar la variable 𝑥 de la ecuación 3 podemos sumarle −4 ∙ 𝑒𝑐2 con lo que nos queda:

−4 ∙ 𝑒𝑐2: 4𝑥 + 4𝑦 + 20𝑧 = 28

Al sumarla con 𝑒𝑐3 obtenemos:

𝑒𝑐3∗: 3𝑦 + 13𝑧 = 13

Para eliminar la variable 𝑥 de la última ecuación podemos multiplicar 𝑒𝑐2 por 3 y sumarla a 𝑒𝑐4.

3 ∙ 𝑒𝑐2 : − 3𝑥 − 3𝑦 − 15𝑧 = −21

Entonces nuestra nueva ecuación 4 nos queda como:

𝑒𝑐4∗ : − 13𝑧 = −13

Nuestro nuevo sistema de ecuaciones equivalentes queda como:

𝑒𝑐1: 𝑤 + 𝑥 + 3𝑧 = 4

𝑒𝑐2 : − 𝑥 − 𝑦 − 5𝑧 = −7

𝑒𝑐3: 3𝑦 + 13𝑧 = 13

𝑒𝑐4 : − 13𝑧 = 29



Como puedes observarm lo que hicimos fue sumarle a la ecuación 𝑗 − é𝑠𝑖𝑚𝑎 (𝑗 > 𝑖) un múltiplo

de la ecuación 𝑖 − é𝑠𝑖𝑚𝑎 de tal forma que los coeficientes de las variables 𝑖 − é𝑠𝑖𝑚𝑎𝑠 en cada

ecuación debajo de la ecuación 𝑖 − é𝑠𝑖𝑚𝑎 se anulan.

Es decir, usando las operaciones de la Tabla 1, que nos garantizan que nuestro sistema

transformado va a ser equivalente al original obtenemos que seva a obtener coeficiente cero en

las variables debajo de la diagonal.

𝑒𝑐1: 𝑤 + 𝑥 + 0𝑦 + 3𝑧 = 4

𝑒𝑐2: − 𝑥 − 𝑦 − 5𝑧 = −7

𝑒𝑐3: 3𝑦 + 13𝑧 = 13

𝑒𝑐4: − 13𝑧 = −13

Concluimos que fue conveniente realizar estas operaciones porque ahora nuestro sistema tiene

una forma triangular (triangular superior para ser específicos) con la característica de que cada

ecuación de abajo hacia arriba tiene una variable indeterminada más que la anterior y 𝑒𝑐4 sólo

tiene una variable que puede ser despejada inmediatamente y de forma concatenada podemos ir

despejando cada variable desde 𝑒𝑐4 hasta la 𝑒𝑐1.

La solución es {𝑤 = −1, 𝑥 = 2, 𝑦 = 0, 𝑧 = 1}.

Si separamos los coeficientes de cada ecuación y ocupamos la operación de multiplicación de

matrices podemos reconstruir el sistema original.

[

1 1 0 3 2 1 −1 1 3−1

−1 2

−1 3

2−1

] [

𝑤𝑥𝑦𝑧

] = [

4 1−3 4

]

Podemos ocupar el conjunto de operaciones que definimos para nuestro sistema de ecuaciones y

aplicarlo al mismo sistema, pero en su representación matricial y únicamente operar con los

puros coeficientes.

El sistema equivalente simplificado es



[

1 1 0 3 0 −1 −1 −5 0 0

0 0

3 0

13−13

] [

𝑤𝑥𝑦𝑧

] = [

4 −7 13−13

]

3.2. Normas de vectores y de matrices

Nuestra definición de error dependía de la magnitud de los datos observados: sin embargo,

necesitamos seguir cuantificando el error y la sensibilidad numérica de nuestros problemas por lo

que es importante establecer la forma de medir el tamaño de las matrices para poder aplicar

nuestros conceptos de condicionamiento, tomando en cuenta que las matrices son elementos de

un espacio vectorial.

Normas de Vectores

Una Norma de Vector 𝒗 en un espacio vectorial 𝑉/𝐾 es una función que va del campo 𝑉a los

reales no negativos.

𝜑: 𝑉 ⟶ ℝ+

Debe cumplir las siguientes condiciones

a) 𝜑(𝒗) = 0 ⟺ 𝒗 = 𝟎 𝑐𝑜𝑛 𝟎 ∈ 𝑉

b) Para 𝜆 ∈ 𝐾. 𝜑(𝜆𝒗) = |𝜆| 𝜑(𝒗)

c) 𝜑(𝒗𝟏 + 𝒗𝟐) ≤ 𝜑(𝒗𝟏) + 𝜑(𝒗𝟐) (Desigualdad del triángulo)

La norma más común es la Euclidiana (‖ ∙ ‖) que se usa para determinar el tamaño de vectores

𝒙 ∈ 𝑅𝑛 definida por

‖𝒙‖2 = (∑|𝑥𝑖|2

𝑖=1

)

12

= (𝑥12 + ⋯+ 𝑥𝑛

2)12

Pero ese es un caso particular para la potencia 2.

Esta definición se puede generalizar con la expresión conocida como Norma de Hölder:



‖𝒙‖𝑝 = (∑|𝑥𝑖|𝑝

𝑛

𝑖=1

)

1𝑝

= (𝑥1𝑝 + ⋯+ 𝑥𝑛

𝑝)1𝑝 𝑝 ≥ 1 (4)

Donde cuando 𝑝 = 2 obtenemos la Norma Euclidiana,

Hay casos significativos por analizar:

• 𝑝 = 1: También conocida en las Ciencias de la Computación como la Norma

Manhattan. Esto se debe a la referencia de que en la ciudad antes mencionada, la

distancia entre dos puntos se hace contando simplemente la cantidad de cuadras que la

separan.

• 𝑝 = 2: Es el caso de la Norma Euclidiana.

• 𝑝 = ∞: También conocida como Norma Infinita y corresponde al caso en el que

sacamos raíz infinita a la suma de todas las entradas elevadas al infinito. Es decir,

corresponde al caso en que:

lim𝑝→∞

‖𝒙‖𝑝 = max1≤𝑖≤𝑛

|𝑥𝑖|

Ésta es una norma práctica ya que significa sacar el valor absoluto a todas las entradas para

después escoger el valor más grande.

Identificaremos a cada norma con un subíndice que indique su orden. En general para cualquier

vector 𝑥 ∈ ℝ𝑛de estas tres normas podemos decir que:

‖𝑥‖1 ≥ ‖𝑥‖2 ≥ ‖𝑥‖∞

En la figura 1 puedes observar cómo se ven cada una de estas tres normas aplicadas a los

diferentes vectores del círculo unitario:

En el caso de 𝑝 = 2 la figura es el círculo, en el caso de 𝑝 = 1 obtenemos cuatro líneas que unen

las intersecciones de círculo con los distintos ejes, esto se obtiene con los diferentes valores que



puede tomar la ecuación de la recta |𝑥| + |𝑦| = 𝑘 (para los cuadrantes la ecuación de las rectas

son 𝑦 = 1 − 𝑥 para el cuadrante I, 𝑦 = 1 + 𝑥 para el cuadrante II,−𝑦 = −1 − 𝑥 para el

cuadrante III y – 𝑦 = −1 + 𝑥 para el cuadrante IV) y en el caso de 𝑝 = ∞ se conforma un

cuadrado tomando la entrada más grande para cada vector.

Las normas para distintos valores de 𝑝 son esencialmente equivalentes pero su uso puede

depender de la aplicación específica, por ejemplo, las normas para 𝑝 = 1 ó𝑝 = ∞ se pueden usar

para realizar análisis de sensibilidad por su facilidad de aplicación, mientras que la Norma

Euclidiana (bajo ciertas transformaciones) es útil dentro de los diferentes métodos que requieran

calcular el tamaño de un vector. Cabe mencionar que esta norma también tiene la propiedad de

ser diferenciable para todo vector 𝒙.

Normas de Matrices

El caso de los vectores es intuitivo en el sentido que podemos asociar (al menos para ℝ3) a cada

vector una flecha como su representación, cuya magnitud queda determinada por qué tan

grande es dicha flecha.

En el caso de las matrices bidimensionales una representación tan directa no es tan fácil y

requerimos definir un concepto de norma.

Figura 2: Visualización de las distintas normas aplicadas a los vectores que cumplen ‖𝒙‖𝒑 = 𝟏



Una definición de norma que estaremos usando es la de norma subordinada a la norma de un

vector que se define como sigue:

‖𝐴𝑛×𝑛‖ = max𝐱≠0

‖𝐴𝒙‖

‖𝒙‖ (5)

Donde 𝒙 ∈ ℝ𝑛, decimos que esta norma está subordinada a la norma de vector‖𝒙‖. Es

importante aclarar que el término 𝐴𝒙 es la multiplicación de la matriz 𝐴 por el vector 𝒙, por lo

tanto, lo que tenemos es un vector en ℝ𝑛 entonces la cantidad ‖𝐴‖ es el valor más largo

determinado por la norma de ‖𝐴𝒙‖ normalizada por todos los vectores 𝒙 distintos de cero.

Los ejemplos más comunes son la Norma de matrices 1 y la Norma infinito, que se verán más

adelante.

Intuitivamente la norma de una matriz, considerando que éstas son transformaciones lineales

sobre vectores, mide qué tanto deforma esta matriz a cualquier vector en términos de la norma en

la que está el vector.

Tomando en cuenta que ya definimos normas para vectores podemos aplicarlas en las matrices.

Por ejemplo la norma de matriz subordinada a la distancia Manhattan:

‖𝐴‖1 = max𝑗

∑|𝑎𝑖𝑗|

𝑖=1

Es equivalente a la suma máxima sobre las columnas de 𝐴. O bien la norma de matrices

subordinada a la norma infinita:

‖𝐴‖∞ = max𝑖

∑|𝑎𝑖𝑗|

𝑗=1

Es equivalente a la suma máxima sobre las filas de A. El caso de la norma de matrices

subordinada a la norma euclidiana es mucho menos intuitivo y corresponde a la raíz del eigen



valor (también conocidos como valores propios o autovalores y son los valores 𝜆 definidos como

𝑑𝑒𝑡(𝐴 − 𝜆𝐼) = 0) más grande del producto 𝐴𝑡𝐴.

A esta definición se le conoce como forma espectral:

‖𝐴‖2 = (𝜌(𝐴𝑇𝐴))12 = 𝜎𝑀𝐴𝑋(𝐴)

Donde el radio espectral es la cantidad:

𝜌(𝐴) = max{|𝜆|: 𝑑𝑒𝑡(𝐴 − 𝜆𝐼) = 0}

El valor 𝜎𝑀𝐴𝑋(𝐴) es el valor singular más grande de 𝐴.

Otra norma que tiene similitudes con la Norma Euclidiana es la Norma de Frobenius y se

define como sigue:

‖𝐴‖𝐹2 = (∑∑|𝑎𝑖𝑗|

2𝑚

𝑗=1

𝑚

𝑖=1

)

12

= 𝑡𝑟(𝐴𝑇𝐴)12

D onde 𝑡𝑟(𝐴) es la traza de la matriz A.

Estas normas deben satisfacer las siguientes propiedades.

Suponiendo 𝐴 y 𝐵 matrices:

a) ‖𝐴‖ > 0 si 𝐴 ≠ 0.

b) ‖𝛼𝐴‖=|𝛼| ∙ ‖𝐴‖ donde 𝛼 ∈ ℝ

c) ‖𝐴 + 𝐵‖ ≤ ‖𝐴‖ + ‖𝐵‖

d) ‖𝐴𝐵‖ ≤ ‖𝐴‖ ∙ ‖𝐵‖

e) ‖𝐴𝒙‖ ≤ ‖𝐴‖‖𝒙‖ para cualquier vector 𝒙

Ejemplos.



Determinar la norma de la matriz asociada a la norma-1 (Norma Manhattan) y a la Norma

Infinita de las siguientes matrices:

𝐴 = [1 4 12 6 00 7 1

] , 𝐵 = [1 0.2 0.03

0.2 0.2 0.030.03 0.03 0.03

]

Para la primera matriz la norma-1 ‖𝐴‖1 = max{1 + 2,4 + 6 + 7,1 + 1} = max{3,17,2} = 17 y

la Norma Infinita es ‖𝐴‖∞ = max{1 + 4 + 1,2 + 6,7 + 1} = max {6,8,8} = 8.

La Norma de Frobenius es:‖𝐴‖𝐹 = √(1 + 4) + (16 + 26 + 49) + (1 + 1) =

√5 + 101 + 2 = √108 = 10.392

Para la segunda matriz la norma-1:

‖𝐴‖1 = max{1 + 0.2 + 0.03, 0.2 + 0.2 + 0.03, 0.03 + 0.03 + 0.03}

= max{1.23, 0.43, 0.09} = 1.23

y la norma infinita es ‖𝐴‖∞ = max{1 + 0.2 + 0.03, 0.2 + 0.2 + 0.03,0.03 + 0.03 + 0.03}

= max{1.23, 0.43, 0.09} = 1.23.

La Norma de Frobenius es:

‖𝐴‖𝐹 = √(12 + 0.22 + 0.032) + (0.22 + 0.22 + 0.032) + 3(0.032)

= √1 + 3 ∗ 0.22 + 5 ∗ 0.032 + √1 + 3(0.04) + 5(0.009909)

Condición de una Matriz

Estas definiciones tienen por objetivo analizar qué tanto deforma la transformación lineal al

vector original y qué tanto afectan pequeñas perturbaciones al resultado del problema.

Necesitamos resolver el siguiente sistema:



𝐴𝒙 = 𝑏

Como 𝑏 está sujeto a pequeñas perturbaciones 𝛿𝑏 entonces en vez de resolver para𝒙 queremos

resolver para �̂�.

Los sistemas que queremos resolver son:

𝐴𝒙 = 𝑏 y 𝐴�̂� = 𝑏 + 𝛿𝑏

es decir𝐴(�̂� − 𝒙) = 𝛿𝑏

y multiplicando por la inversa de 𝐴 (esto se puede porque la matriz no es singular) observamos

que el cambio en la entrada, es decir, el error en los datos está dado por:

(�̂� − 𝒙) = 𝐴−1𝛿𝑏

El tamaño de dicho error, para alguna norma, es:

‖�̂� − 𝒙‖ = ‖𝐴−1𝛿𝑏‖ ≤ ‖𝐴−1‖ ∙ ‖𝛿𝑏‖

La desigualdad la podemos establecer por la propiedad d). Esta medida nos indica el error

absoluto en los datos de entrada, pero también el error relativo así que para esto vamos a tomar la

siguiente relación:

‖𝑏‖ = ‖𝐴𝒙‖ ≤ ‖𝐴‖ ∙ ‖𝒙‖

Donde podemos observar que:

1

‖𝒙‖≤

‖𝐴‖

‖𝑏‖

Multiplicando ambos lados de esta desigualdad tenemos que:

‖𝒙 − 𝒙‖

‖𝒙‖≤ ‖𝐴‖ ∙ ‖𝐴−1‖

‖𝛿𝑏‖

‖𝑏‖

Es decir, observamos la relación que guarda el error relativo en la entrada con el miembro

izquierdo de la desigualdad respecto del error relativo en la salida en el miembro derecho de la



desigualdad relacionados por el número ‖𝐴‖ ∙ ‖𝐴−1‖ al que denominaremos como el

condicionamiento o número de condición de A:

𝜅(𝐴) = ‖𝐴‖ ∙ ‖𝐴−1‖ (6)

Para alguna norma dada.

Esta definición nos permitirá ver qué tan sensible a perturbaciones son los problemas en los que

esté involucrada la matriz 𝐴.

Si el valor de 𝜅 ≫ 1 entonces diremos que 𝐴 está mal condicionada.

Ejemplo.

Para los casos anteriores 𝜅(𝐶) es, en el caso de la matriz 𝐴 y las normas 1, 2,∞ y F

respectivamente 𝜅(𝐴) = {31.167, 14.391, 15.333, 15.945}.

Para la matriz 𝐵 los condicionamientos respectivos son 𝜅(𝐵) =

{413.21, 354.47, 413.21, 358.30}.

Recuerda que el condicionamiento de una matriz es un solo número, pero en este caso, estoy

usando la notación para presentar los números de condiciones especificados de forma eficiente.

Al ver los número anteriores, comparando el condicionamiento equivalente entre cada matriz

para las distintas normas, podemos observar que en general la segunda matriz está peor

condicionada que la primera, esto quiere decir que si un problema involucra dicha matriz es más

sensible a pequeñas perturbaciones en sus datos.

En Octave las matrices se definen elemento a elemento separando cada fila por un “;” (punto y

coma)



octave-3.2.4.exe:1> A=[1 4 1; 2 6 0; 0 7 1]

A =

1 4 1

2 6 0

0 7 1

octave-3.2.4.exe:2> B=[ 1 0.2 0.03; 0.2 0.2 0.03; 0.03 0.03 0.03]

B =

1.000000 0.200000 0.030000

0.200000 0.200000 0.030000

0.030000 0.030000 0.030000

octave-3.2.4.exe:3>

Las normas las puedes obtener con la función norm y un parámetro para especificar la norma

específica que deseas:

octave-3.2.4.exe:3> norm(A,1)

ans = 17

octave-3.2.4.exe:4> norm(A,Inf)

ans = 8

octave-3.2.4.exe:5> norm(A,2)

ans = 10.236

octave-3.2.4.exe:6> norm(A,'fro')

ans = 10.392

octave-3.2.4.exe:7>

Y para determinar el condicionamiento de una matriz puedes aplicar la definición recién

presentada o usar la función cond.



octave-3.2.4.exe:7>cond(A)

ans = 14.391

octave-3.2.4.exe:8> norm(pinv(A),2)*norm(A,2)

ans = 14.391

octave-3.2.4.exe:9>cond(B)

ans = 42.214

octave-3.2.4.exe:10> norm(pinv(B),2)*norm(B,2)

ans = 42.214

octave-3.2.4.exe:11>

La norma usada por omisión para calcular el condicionamiento es la Norma Euclidiana, la

función pinv se usa para calcular la pseudo inversa de una matriz y sólo resta concluir que la

matriz B está peor condicionada que la A, esto significa que un problema que haga uso de esa

matriz de coeficientes será más sensible a pequeños cambios en los datos de entrada.

3.2.1. Solución de un sistema de ecuaciones lineales

Una vez establecida una relación entre los sistemas de ecuaciones y el álgebra de matrices y

examinadas definiciones para operar con ellas, podemos decir que el problema equivalente a

resolver es:

𝐴𝑚,𝑛𝒙𝒏,𝟏 = 𝑏𝑛,1 (7)

Para algún vector columna 𝒙. donde A es la matriz de coeficientes del sistema de ecuaciones, a

𝒙 lo denominamos vector solución y al vector columna 𝑏 vector de soluciones u observaciones.

No podemos despejar𝒙 como en las ecuaciones de una sola variable, pues éste, es un sistema de

ecuaciones y por lo tanto el despeje va a requerir manipulación matricial.



Llo que podemos hacer es multiplicar ambos lados de la ecuación por una matriz no singular 𝑀

tal que:

𝑀𝐴 = 𝐼

𝑀𝐴𝒙 = 𝑀𝑏

Entonces nuestro sistema se transforma en:

𝒙 = 𝑀𝑏 (8)

Nos estamos refiriendo a los sistemas equivalentes con matrices que son transformaciones

lineales sobre los vectores columna 𝒙. Es decir, podemos conceptualizar a la expresión (7) como

la deformación lineal de la matriz 𝐴 sobre 𝒙 restringido por los valores del vector columna 𝑏.

Los métodos que veremos en esta sección son para sistemas que no están sobredeterminados ni

subdeterminado, es decir, donde hay tantas variables como valores observados y por lo tanto

tenemos matrices cuadradas. En los siguientes algoritmos veremos la forma para construir a 𝑀 y

poder satisfacer la expresión (8).

3.2.2. Eliminación de Gauss

Este es el método usado para resolver un sistema de ecuaciones empleando operaciones

elementales.

Para poder entender con más facilidad el Método de Gauss primero hay que definir algunas

formas matriciales.

La forma más sencilla de un sistema de ecuaciones es donde la ecuación tiene una sola variable

asignada a un valor observado:

𝑎11𝑥1 = 𝑣1

Por otro lado, si alguna ecuación tiene dos variables, una de las cuales es la variable anterior

hasta llegar a tener alguna ecuación con 𝑛 variables tenemos que:

𝑎21𝑥1 + 𝑎12𝑥2 = 𝑣2

⋮



𝑎𝑛1𝑥1 + ⋯+ 𝑎𝑛𝑛𝑥𝑛 = 𝑣𝑛

La facilidad de trabajar con estos sistemas es que la ecuación con una variable ya está

determinada y como la ecuación con dos variables contiene a esta variable entonces la podemos

ir sustituyendo hasta resolver el sistema completo.

Esto lo puedes revisar en la sección 3.1.1.

Al transformar el sistema en su forma matricial obtenemos una matriz 𝐴 con muchos ceros. Los

renglones de estas matrices pueden reacomodarse con las operaciones básicas sobre matrices sin

alterar el sistema de ecuaciones de tal forma que el renglón con más coeficientes sea el primero,

el segundo será el que tiene una variable menos y así hasta que el último renglón representará la

ecuación que sólo tiene una variable asociada a un valor observado.

La matriz resultante tendrá una gran cantidad de ceros en una zona de la matriz. Este tipo de

matrices se denominan triangulares superiores o inferiores, dependiendo de dónde se

encuentren todas las entradas de la matriz cuyo valor es 0.

A continuación describimos en qué consisten estas matrices, así como su importancia en la

resolución de un sistema de ecuaciones lineales.

Sistemas triangulares

Si una matriz de coeficientes 𝐿 tiene la siguiente forma:

𝐿 = [

𝑙11 0𝑙21 𝑙22

… 0… 0

⋮ ⋮𝑙𝑛1 𝑙𝑛2

⋱ ⋮… 𝑙𝑛𝑛

]



Es decir, si todas las entradas donde 𝑖 < 𝑗 tenemos que 𝑙𝑖𝑗 = 0 entonces a 𝐿 se le denomina

matriz triangular inferior. Si los elementos de la diagonal 𝑙𝑖𝑖 = 0 entonces decimos que es

estrictamente diagonal inferior.

Suponiendo que tenemos un sistema de ecuaciones en su representación matricial de la siguiente

forma:

𝐿𝒙 = 𝑏

Donde el vector columna 𝑏 = (𝑏1, … , 𝑏𝑛)𝑇, entonces el vector solución 𝒙 está determinado por

𝑥1 =𝑏1

𝑙11

𝑥𝑖 = (𝑏𝑖 − ∑ 𝑙𝑖𝑗𝑥𝑗

𝑛

𝑗=𝑖+1

) 𝑙𝑖𝑖⁄ , 𝑖 = 2,… , 𝑛

(9)

A este método se le denomina comúnmente Método de sustitución hacia adelante y puedes

observar en la segunda expresión de (9) que obtenemos el siguiente valor de 𝒙 ,que resulta

indeterminado una vez calculado los anteriores.

Análogamente si una matriz de coeficientes 𝑈 tiene la siguiente forma:

𝑈 = [

𝑢11 𝑢12

0 𝑢22

… 𝑢1𝑛

… 𝑢2𝑛

⋮ 00 0

⋱ ⋮0 𝑢𝑛𝑛

]

Es decir, si todas las entradas donde 𝑖 > 𝑗 tenemos que 𝑢𝑖𝑗 = 0 entonces a 𝑈 se le denomina

matriz triangular superior. Adicionalmente si los elementos de la diagonal 𝑢𝑖𝑖 = 0 entonces

decimos que es estrictamente diagonal superior.



Supón que tenemos un sistema de ecuaciones en su representación matricial de la siguiente

forma:

𝑈𝒙 = 𝑏

Donde el vector columna 𝑏 = (𝑏1, … , 𝑏𝑛)𝑇, entonces el vector solución 𝒙 está determinado por:

𝑥𝑛 =𝑏𝑛

𝑢𝑛𝑛

𝑥𝑖 = (𝑏𝑖 − ∑𝑢𝑖𝑗𝑥𝑗

𝑛

𝑗=1

) 𝑢𝑖𝑖⁄ , 𝑖 = 𝑛 − 1, … ,1

(10)

A este método se le denomina Método de sustitución hacia atrás.

Eliminación Gaussiana

Tener un sistema de ecuaciones como en las representaciones de expresiones (9) y (10) es útil,

pues el Método de Gauss aprovecha esta característica de los sistemas triangulares.

En el caso de tener un sistema lineal como el que sigue:

𝑒𝑐1: 𝑎11𝑥1 + ⋯+ 𝑎1𝑗𝑥𝑗 + ⋯+ 𝑎1𝑛𝑥𝑛 = 𝑏1

𝑒𝑐2: 𝑎21𝑥1 + ⋯+ 𝑎2𝑗𝑥𝑗 + ⋯+ 𝑎2𝑛𝑥𝑛 = 𝑏2

⋮

𝑒𝑐𝑁: 𝑎𝑛1𝑥1 + ⋯+ 𝑎𝑛𝑗𝑥𝑗 + ⋯+ 𝑎𝑛𝑛𝑥𝑛 = 𝑏𝑛

(11)

Cuya representación matricial es:

𝐴𝒙 = 𝑏

Donde 𝐴 no es singular.

Para obtener 𝒙 el Método de Gauss se divide en dos fases:



a) Fase de Conversión: A través de operaciones elementales convertimos a nuestro sistema

original en uno equivalente de la forma:

𝑈𝒙 = 𝑐

Esto lo haremos operando sobre las distintas ecuaciones de (11) de tal forma que la

ecuación 𝑖 − 𝑒𝑠𝑖𝑚𝑎 sea sustituida por alguna otra fila 𝑏𝑗 − é𝑠𝑖𝑚𝑎 denominada pivote

multiplicado por un número real 𝜆.

𝑒𝑐𝑖 ← 𝑒𝑐𝑖 + 𝜆𝑒𝑐𝑗

Esta operación es la parte medular del Método de Gauss ya que el valor de 𝜆 debe ser

calculado de tal forma que el primer coeficiente de 𝑒𝑐𝑖 se anule, por esto la ecuación 𝑒𝑐𝑗

cumple con que 𝑗 > 𝑖 y su primer coeficiente distinto de cero es un elemento de la

diagonal.

Supongamos que en algún paso de nuestra fase de sustitución tenemos una matriz 𝐴∗ que

no está completamente transformada en 𝑈, entonces debe tener la siguiente forma:

𝐴∗ = [

𝑢11 𝑢12

0 𝑎𝑖𝑖

… 𝑢1𝑛

… 𝑎𝑖𝑛

⋮ 𝑎𝑗𝑖

⋮ ⋮

⋱ 𝑎𝑗𝑛

⋮ ⋮

]

Lo que nosotros necesitamos es que 𝑎𝑘𝑖 = 0 a partir de algún múltiplo de 𝑎𝑖𝑖

0 = 𝑎𝑗𝑖 + 𝜆 𝑎𝑖𝑖

Entonces el factor por el que debemos multiplicar nuestra fila pivote debe cumplir:

𝜆𝑖 = −𝑎𝑗𝑖

𝑎𝑖𝑖 (12)

Recuerda que todas las operaciones (permutaciones y multiplicaciones) que realices

sobre la matriz 𝐴 las debes hacer sobre la matriz 𝑏 ya que de otra forma terminarás con



un sistema distinto y no equivalente que es lo que buscamos, es por eso que el sistema

final debe ser 𝑈𝒙 = 𝑐

b) Fase de Resolución. Una vez que tenemos el sistema original transformado a uno del tipo

𝑈𝒙 = 𝑏

Entonces resolvemos para 𝒙 usando sustitución hacia atrás.

Ejemplo:

Algún fenómeno está representado por el siguiente sistema:

𝑒𝑐1: 𝑥1 + 𝑥2 + 3𝑥4 = 4

𝑒𝑐2: 2𝑥1 + 𝑥2 − 𝑥3 + 𝑥4 = 1

𝑒𝑐3: 3𝑥1 − 𝑥2 − 𝑥3 + 2𝑥4 = −3

𝑒𝑐4 : − 𝑥1 + 2𝑥2 + 3𝑥3 − 𝑥4 = 4

Este sistema en su forma matricial queda como:

[

1 1 0 32 1 −1 13 −1 −1 2

−1 2 3 −1

] [

𝑤𝑥𝑦𝑧

] = [

41

−34

]

Usamos a 𝑒𝑐1 como pivote y usando la expresión (12) para obtener los factores 𝜆𝑖 ∈

{−2

1, −

3

1,1

1} , 𝑖 = 2,3,4 para ser usados en las siguientes sustituciones 𝑒𝑐2 ← 𝑒𝑐2 −

2𝑒𝑐1,𝑒𝑐3 ← 𝑒𝑐3 − 3𝑒𝑐1,𝑒𝑐4 ← 𝑒𝑐4 + 𝑒𝑐1 con lo que nuestro sistema se transforma en:

[

1 1 0 30 −1 −1 −50 −4 −1 −70 3 3 2

] [

𝑤𝑥𝑦𝑧

] = [

4−7

−158

]



Volvemos a usar (12) para encontrar los nuevos factores 𝜆𝑖 ∈ {−4,3} para usarlos en nuestro

nuevo pivote 𝑒𝑐2 de la siguiente forma 𝑒𝑐3 ← 𝑒𝑐3 − 4𝑒𝑐2 y 𝑒𝑐4 ← 𝑒𝑐4 + 3𝑒𝑐2 entonces nuestro

sistema se transforma en:

[

1 1 0 30 −1 −1 −50 0 3 130 0 0 −13

] [

𝑤𝑥𝑦𝑧

] = [

4−713

−13

]

Con esto obtenemos un sistema triangular superior, ahora sólo resta resolver para

{𝑤, 𝑥, 𝑦, 𝑧} utilizando sustitución hacia atrás usando (9):

𝑧 =−13

−13= 1,

𝑦 = 13 − (13 ∙ 1) = 0,

𝑥 = (−7 − (−1 ∙ 0 − 5 ∙ 1)) −1⁄ = 2,

𝑤 = (4 − (1 ∙ 2 + 0 ∙ 0 + 3 ∙ 1)) 1⁄ = −1.

Ejemplo

Resolver el sistema para 𝒙 = (𝑥1, 𝑥2, 𝑥3) planteado en el siguiente sistema lineal:

𝑒𝑐1: 4𝑥1 − 2𝑥2 + 𝑥3 = 11

𝑒𝑐2 : − 2𝑥1 + 4𝑥2 − 2𝑥3 = −16

𝑒𝑐3: 𝑥1 − 2𝑥2 + 4𝑥3 = 17

En forma matricial este sistema queda como:

[4 −2 1

−2 4 −21 −2 4

] [𝑥1𝑥2

𝑥3

] = [11

−1617

]

Tomando 𝑒𝑐1 como nuestra ecuación pivote para operar en las ecuaciones 2 y 3 así como los

factores 𝜆𝑖 ∈ {1

2,1

4} respectivamente para cada operación transformamos nuestro sistema en:



[4 −2 10 3 −1.50 −1.5 3.75

] [𝑥1𝑥2

𝑥3

] = [11

−10.514.25

]

Ahora tomamos a 𝑒𝑐2 como nuestro pivote y usamos 𝜆 = −−1.5

3=

1

2 para hacer la operación

𝑒𝑐3 ← 𝑒𝑐3 + 𝜆𝑒𝑐2 con lo que transformamos nuestro sistema a:

[4 −2 10 3 −1.50 0 3

] [𝑥1𝑥2

𝑥3

] = [11

−10.59]

Así terminamos nuestra fase de eliminación para pasar a la fase de sustitución:

𝑥3 = 3,

𝑥2 = (−10.5 + 1.5 ∙ 3) 3⁄ = −2,

𝑥1 = (11 − (−2 ∙ −2 + 1 ∙ 3)) 4⁄ = 1

Este método es el más conocido y fácil de implementar computacionalmente en Octave. Para

resolver un sistema de ecuaciones puedes usar el operador “\” una vez que tengas definida la

matriz 𝐴 y 𝑏 para obtener el valor de 𝒙.

Este operador es transparente al usuario y usa el método para optimizar la resolución de sistemas

lineales, por ejemplo, descomponiendo 𝐴 = 𝐿𝑈 para alguna 𝐿 y 𝑈 particulares.

octave-3.2.4.exe:1> A=[1 1 0 3; 2 1 -1 1; 3 -1 -1 2; -1 2 3 -1]

A =

1 1 0 3

2 1 -1 1

3 -1 -1 2

-1 2 3 -1



octave-3.2.4.exe:2> b=[4;1;-3;4]

b =

4

1

-3

4

octave-3.2.4.exe:3> x=A\b

x =

-1.0000e+000

2.0000e+000

-8.6348e-017

1.0000e+000

octave-3.2.4.exe:4>

Observa que el valor 𝑥3 que en nuestro ejemplo es 0 para Octave

octave-3.2.4.exe:5> x(3)

ans = -8.6348e-017

octave-3.2.4.exe:6>

Cuando obtenemos algún valor 𝑣 × 10−17 lo consideraremos 0.

3.2.3. Método iterativo de Jacobi

El Método de Gauss no es el más adecuado para sistemas grandes, puede ser muy costoso, pero

podemos establecer una estrategia distinta que se base en una solución inicial aproximada y



sobre ella seguir aplicando el método una y otra vez hasta que bajo algún parámetro o umbral de

paro nos detengamos. A este tipo de métodos se les llama Métodos iterativos.

Un Método Iterativo𝑓 es un procedimiento en el que tomando la solución inicial aproximada

𝑥(0) la nueva solución 𝑚 + 1 se calcula a partir de la solución anterior calculada en la iteración,

esto es:

𝑥(𝑚+1) = 𝑓(𝑥(𝑚))

Hasta que se alcance el criterio de convergencia predefinido. las convenciones para detener el

algoritmo pueden ser variadas, pero, entre ellas, podemos observar que las más populares son:

a) La cantidad 𝑚 > 𝑀 de iteraciones máxima se alcanza

b) El error absoluto es ‖𝑥(𝑚+1)– 𝑥(𝑚)‖ < 휀 se rebasa y es menor que algún 휀 predefinido

tal que 0 < 휀 ≪ 1

c) El error relativo ‖𝑥(𝑚+1)– 𝑥(𝑚)‖ ‖𝑥(𝑚+1)‖⁄ < 휀 respecto de la solución anterior se

rebasa. Este será el criterio de convergencia que usaremos preferentemente.

El Método de Jacobi (también conocido como Método de los desplazamientos simultáneos).

Su implementación es extremadamente sencilla ya que los nuevos vectores solución se calculan a

partir de despejar cada variable 𝑥𝑖 a partir de las ecuaciones originales en el sistema, esto es:

𝑥𝑖(𝑚+1)

= (𝑏𝑖 − ∑𝑎𝑖𝑗𝑥𝑗(𝑚)

𝑗≠𝑖

) 𝑎𝑖𝑖⁄ , 𝑖 = 1,… , 𝑛 (13)

Lo que significa despejar cada variable de cada ecuación 𝑖 − é𝑠𝑖𝑚𝑎 para obtener la nueva

solución 𝒙(𝒎+𝟏).



De forma expresa esto significa que para aplicar el Método de Jacobi primero tenemos que

tomar nuestro sistema

𝑎11𝑥1 + ⋯+ 𝑎1𝑗𝑥𝑗 + ⋯+ 𝑎1𝑛𝑥𝑛 = 𝑏1

𝑎21𝑥1 + ⋯+ 𝑎2𝑗𝑥𝑗 + ⋯+ 𝑎2𝑛𝑥𝑛 = 𝑏2

⋮

𝑎𝑛1𝑥1 + ⋯+ 𝑎𝑛𝑗𝑥𝑗 + ⋯+ 𝑎𝑛𝑛𝑥𝑛 = 𝑏𝑛

Y después despejar 𝑥_𝑖 de cada una de las ecuaciones:

𝑥1 =1

𝑎11∙ (𝑏1 − 𝑎12𝑥2 − ⋯− 𝑎1𝑗𝑥𝑗 − ⋯− 𝑎1𝑛𝑥𝑛)

𝑥2 =1

𝑎22∙ (𝑏2 − 𝑎21𝑥1 − ⋯− 𝑎2𝑗𝑥𝑗 − ⋯− 𝑎2𝑛𝑥𝑛)

⋮

𝑥𝑛 =1

𝑎𝑛𝑛∙ (𝑏𝑛 − 𝑎𝑛1𝑥1 − ⋯− 𝑎𝑛𝑗𝑥𝑗 − ⋯− 𝑎𝑛𝑛−1𝑥𝑛−1)

En caso de que algún elemento de la diagonal sea 0 sólo tienes que intercambiar las filas para

evitar esto.

Ejemplo:

Considera el sistema

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

−3𝑥1 + 9𝑥2 + 𝑥3 = 2

2𝑥1 − 𝑥2 − 7𝑥3 = 3

Al despejar cada variable obtenemos:

𝑥1 =1

5∙ (−1 + 2𝑥2 − 3𝑥3)

𝑥2 =1

9∙ (2 + 3𝑥1 − 𝑥3)

𝑥3 = −1

7∙ (3 − 2𝑥1 + 𝑥2)



Y si tomamos como 𝒙(𝟎) = (0,0,0) nuestra primera iteración queda:

𝑥1(1)

=1

5∙ (−1 + 2 ∙ 0 − 3 ∙ 0) = −0.200

𝑥2(1)

=1

9∙ (2 + 3 ∙ 0 − 0) ≈ 0.222

𝑥3(1)

= −1

7∙ (3 − 2 ∙ 0 + 0) ≈ −0.429

A continuación puedes ver las iteraciones hasta alcanzar la convergencia en la siguiente tabla:

𝒙𝟏 𝒙𝟐 𝒙𝟐 𝒆𝒓𝒓𝒐𝒓𝒓𝒆𝒍 𝒊𝒈𝒖𝒂𝒍

0.000 0.000 0.000

1 -0.200 0.222 -0.429 0.522547783 no

2 0.146 0.203 -0.517 0.555519273 no

3 0.192 0.328 -0.416 0.543117008 no

4 0.181 0.332 -0.421 0.660847684 no

5 0.185 0.329 -0.424 0.661631653 no

6 0.186 0.331 -0.423 0.6604241 no

7 0.186 0.331 -0.423 0.662452308 si

8 0.186 0.331 -0.423 0.662492338 si

9 0.186 0.331 -0.423 0.662428382 si

10 0.186 0.331 -0.423 0.66245973 si

11 0.186 0.331 -0.423 0.662461787 si

12 0.186 0.331 -0.423 0.662460076 si

13 0.186 0.331 -0.423 0.662460527 si

14 0.186 0.331 -0.423 0.662460598 si

15 0.186 0.331 -0.423 0.662460561 si

16 0.186 0.331 -0.423 0.662460567 si

17 0.186 0.331 -0.423 0.662460568 si

18 0.186 0.331 -0.423 0.662460568 si

Tabla 2. Iteraciones hasta la convergencia



Lo que indican las columnas de esta tabla es iteración, primera, segunda y tercera entrada del

vector calculado, error relativo respecto de la solución anterior y una función que calcula si todas

las entradas fueron iguales a la solución de la iteración anterior con un umbral de 0.001.

La solución por la que optaremos es 𝒙(𝟏𝟖)y es que a partir de 𝑚 = 19, el error relativo se

estabilizó. Este será el criterio por el que optaremos para escoger la mejor solución, en caso de

que el error no se estabilice nunca (algo que bien podría suceder) entonces escogeremos la

iteración para la que el vector tuvo todas sus entradas iguales.

Nuestro resultado expresamente es 𝒙 = (0.186,331,−0.423).

Ejemplo:

Considera el sistema:

10𝑥1 + 𝑥2 + 𝑥3 = 12

𝑥1 + 10𝑥2 + 𝑥3 = 12

𝑥1 + 𝑥2 + 10𝑥3 = 12

Aplicando (13) para cada ecuación obtenemos:

𝑥1 =1

10∙ (12 − 𝑥2 − 𝑥3)

𝑥2 =1

10∙ (12 − 𝑥1 − 𝑥3)

𝑥3 =1

10∙ (12 − 𝑥1 − 𝑥2)

Nuestra primera aproximación también será 𝒙(𝟎) = (0,0,0) (Toma en cuenta que no

necesariamente tiene que ser así, pero en los ejemplos puede aplicarse esta aproximación inicial).

Aplicando simultáneamente la aproximación inicial a las expresiones anteriores obtenemos:

𝑥1(1)

=1

10∙ (12 − 0 − 0) = 1.2



𝑥2(1)

=1

10∙ (12 − 0 − 0) = 1.2

𝑥3(1)

=1

10∙ (12 − 0 − 0) = 1.2

Las diferentes iteraciones las puedes ver en la siguiente tabla que tiene las mismas

consideraciones que la tabla anterior

m x_1 x_2 x_2 error_rel igual

0.000 0.000 0.000

1 -0.200 0.156 -0.508 0.567624019 no

2 0.167 0.334 -0.429 0.6172092 no

3 0.191 0.333 -0.422 0.668266027 no

4 0.186 0.331 -0.423 0.664702009 si

5 0.186 0.331 -0.423 0.662406964 si

6 0.186 0.331 -0.423 0.66243197 si

7 0.186 0.331 -0.423 0.662461068 si

8 0.186 0.331 -0.423 0.662460934 si

9 0.186 0.331 -0.423 0.662460564 si

10 0.186 0.331 -0.423 0.662460563 si

11 0.186 0.331 -0.423 0.662460568 si

Tabla 3. Iteraciones

Nuestra solución es el vector 𝒙(𝟏𝟏) = (0.186,0.331,−423).

3.3. Método iterativo de Gauss-Seidel

Este método en particular se basa en la aplicación de un concepto denominado inversa

aproximada.

Una matriz 𝐶 es inversa aproximada de 𝐴 si:

‖𝐼 − 𝐶𝐴‖ < 1 (14)



Para alguna norma. Al aplicar un método iterativo del tipo

𝑔(𝒙) = 𝒙 − 𝐶𝑓(𝒙)

Donde 𝑓(𝒙) = 𝐴𝒙 − 𝑏 obtenemos que:

𝑔(𝒙) = 𝒙 − 𝐶(𝐴𝒙 − 𝑏) = 𝐶𝑏 + (𝐼 − 𝐶𝐴)𝒙 (15)

Al aplicar (15) a dos vectores obtenemos que:

𝑔(𝒙) − 𝑔(𝒚) = (𝐼 − 𝐶𝐴)(𝒙 − 𝒚)

Y si tomamos su norma vemos que se cumple que:

‖𝑔(𝒙) − 𝑔(𝒚)‖ ≤ ‖𝐼 − 𝐶𝐴‖ ∙ ‖𝒙 − 𝒚‖

Al considerar (14) podemos concluir que 𝑔 es un mapeo contrayente, lo que garantiza la

convergencia del Método de Jacobi.

El Método de Gauss-Seidel hace una modificación a este método ya que 𝐶 usualmente se toma

como 𝐷−1donde 𝐷 = 𝑑𝑖𝑎𝑔(𝐴), pero si tomamos a 𝐶 = 𝐿𝐴 + 𝐷 (donde 𝐿𝐴 es la matriz triangular

inferior de 𝐴) entonces podríamos esperar que el método converja más rápido, esto es porque

esperamos que

‖𝐼 − (𝐿 + 𝐷)−1𝐴‖ ≤ ‖𝐼 − 𝐷−1𝐴‖ < 1

Lo anterior no es totalmente cierto pero sí para algunos casos específicos de matrices como las

tridiagonales (aquellas con elementos alrededor de la diagonal), las dominantes por diagonal

(aquellas en que para cada fila el elemento de la diagonal es mayor que la suma del resto de los

elementos de esa misma fila) o en general para matrices ralas que es la aplicación principal de

estos métodos o bien cuando 𝐴 tiene elementos positivos en la diagonal y negativos o cero fuera

de ella.

En este método la forma general de iteración es como sigue:



𝒙(𝒎+𝟏) = 𝒙(𝒎) + (𝐿 + 𝐷)−1(𝑏 − 𝐴𝒙(𝒎))

(𝐿 + 𝐷)𝒙(𝒎+𝟏) = (𝐿 + 𝐷)𝒙(𝒎) + (𝑏 − 𝐴𝒙(𝒎))

(𝐿 + 𝐷)𝒙(𝒎+𝟏) = (𝐿 + 𝐷 − 𝐴)𝒙(𝒎) + 𝑏

𝐷𝒙(𝒎+𝟏) = −𝐿𝑥(𝑚+1) − 𝑈𝑥(𝑚) + 𝑏 ; 𝑈 = 𝐿 + 𝐷 − 𝐴

Es decir, la definición de la siguiente solución viene dada por elementos ya calculados en la

presente iteración (−𝑈𝑥(𝑚)). La fórmula explícita queda como:

𝑥𝑖(𝑚+1)

= (−∑𝑎𝑖𝑗𝑥𝑗(𝑚+1)

−

𝑗<𝑖

∑𝑎𝑖𝑗𝑥𝑗(𝑚)

+ 𝑏𝑖

𝑗>𝑖

) 𝑎𝑖𝑖⁄ (16)

Al Método de Gauss-Seidel también se le conoce como Método de los Desplazamientos

Sucesivos (a diferencia de Jacobi en el que los desplazamientos son simultáneos), esto quiere

decir, y a partir de la fórmula (16), que la aplicación del resultado recién computado se puede

aplicar inmediatamente en la siguiente entrada del vector 𝒙.

Ejemplo:

Calcularemos el resultado de

𝐴𝒙 = 𝑏

Para los ejemplos anteriores calculados con el Método de Jacobi. Entonces el sistema y su

respectivo despeje quedan igual, así como la condición inicial 𝒙(𝟎), lo que cambia es la primera

iteración la cual queda de la siguiente forma:

Calculamos la primera entrada del vector resultado

𝑥1(1)

=1

5∙ (−1 + 2 ∙ 0 − 3 ∙ 0) = −0.200

En cuanto tenemos 𝑥1 lo usamos para calcular 𝑥2



𝑥2(1)

=1

9∙ (2 − 3 ∙ 0.200 − 0) ≈ 0.156

Y hacemos lo mismo con 𝑥3 aprovechando los valores recién calculados de 𝑥1 y 𝑥2

𝑥3(1)

= −1

7∙ (3 + 2 ∙ 0.200 + 0.156) ≈ −0.508

El resto de las iteraciones las puedes ver en la siguiente tabla con las mismas consideraciones

que los ejemplos anteriores

𝒎 𝒙𝟏 𝒙𝟐 𝒙𝟐 𝒆𝒓𝒓𝒐𝒓𝒓𝒆𝒍 𝒊𝒈𝒖𝒂𝒍

0 0.000 0.000 0.000

1 -0.200 0.156 -0.508 0.567624019 no

2 0.167 0.334 -0.429 0.6172092 no

3 0.191 0.333 -0.422 0.668266027 no

4 0.186 0.331 -0.423 0.664702009 si

5 0.186 0.331 -0.423 0.662406964 si

6 0.186 0.331 -0.423 0.66243197 si

7 0.186 0.331 -0.423 0.662461068 si

8 0.186 0.331 -0.423 0.662460934 si

9 0.186 0.331 -0.423 0.662460564 si

10 0.186 0.331 -0.423 0.662460563 si

11 0.186 0.331 -0.423 0.662460568 si


Observa que tardó 8 iteraciones menos para llegar al resultado(0.186,0.331,−0.423).

Para el segundo ejemplo revisa las iteraciones en la siguiente tabla:

𝒎 𝒙_𝟏 𝒙_𝟐 𝒙_𝟐 𝒆𝒓𝒓𝒐𝒓_𝒓𝒆𝒍 𝒊𝒈𝒖𝒂𝒍

0.0000 0.0000 0.0000 0.0000



1 1.2000 1.0800 0.9720 1.8845 no

2 0.9948 1.0033 1.0002 2.0936 no

3 0.9996 1.0000 1.0000 2.0033 no

4 1.0000 1.0000 1.0000 2.0000 si

5 1.0000 1.0000 1.0000 2.0000 si


El resultado expreso al que llegamos en la 5ª iteración es (1,1,1).

3.4. Regresión lineal

Otro problema con el que nos encontramos al tratar con datos experimentales es el de hacer pasar

una curva que describa el comportamiento de la forma más cercana posible a los datos

observados.

Una de las formas de construir esta función es con los Métodos de Interpolación en los cuales

la función construida pasa exactamente por los puntos observados. Pero muchas veces esta

aproximación no es realista y esto se debe básicamente a la existencia de ruido y procesos que no

controlamos en la generación de datos

Muchas veces, hay que conformarnos con la curva que mejor describa la tendencia de los datos

observados o medidos, a este proceso de construcción se le conoce como Regresión y decimos

que es lineal cuando la curva que describa la tendencia general de los datos sea una recta.

En este tipo de problemas existen más datos observados que variables, es decir, tenemos sistemas

sobredeterminados, esta característica la usamos para aprovechar la información estadística con

una gran cantidad de datos.

En la figura 1 puedes ver la recta que mejor describe los datos observados marcados con (+) en la

gráfica de la figura 2.



Los datos son

𝑋 = {1.7, 1.6, 2.8, 5.6, 1.3, 2.2, 1.3, 1.1, 3.2, 1.5, 5.2, 4.6, 5.8, 3.0}

𝑌 = {3.7, 3.9, 6.7, 9.5, 3.4, 5.6, 3.7, 2.7, 5.5, 2.9, 10.7, 7.6, 11.8, 4.1}

Y la recta obedece a la ecuación:

𝑦 = 1.6699𝑥 + 0.96447

Figura 3. Ejemplo regresión lineal a una función del tipo y=mx+b

En esta unidad analizaremos los Métodos de regresión lineal por mínimos cuadrados usando

ecuaciones normales y el de regresión lineal múltiple.

3.4.1. El modelo de regresión lineal

Considerando que tenemos 𝑚 puntos observados de la forma (𝑡𝑖, 𝑦𝑖) el método de regresión

lineal intenta ajustar alguna combinación lineal de un vector con𝑛 parámetros 𝝓 = (𝜙1, … , 𝜙𝑛).

La combinación lineal que buscamos es de la forma:



𝑓(𝑡, 𝑥) = 𝑥1𝜙1(𝑡) + ⋯+ 𝑥𝑛𝜙𝑛(𝑡) (17)

Pero de tal forma que sea la mejor posible con el criterio a seguir de minimizar la distancia de

cada punto a la recta, es decir:

min∑(𝑦𝑖 − 𝑓(𝑡𝑖, 𝑥))2

𝑖=1

(18)

Donde las funciones 𝜙𝑖serán determinadas a priori, usualmente basándonos en la experiencia

propia a partir de la distribución de los datos para que el ajuste sea el mejor posible.

A este método se le denomina Ajuste Por Mínimos Cuadrados, es importante, encontrar un

vector de 𝑛 entradas usando 𝑚 puntos que serán evaluados en las funciones 𝜙𝑖(𝑡). Es decir,

vamos a trabajar con un sistema sobredeterminado de 𝑚 × 𝑛 entradas.

Es importante notar que en este caso no buscamos encontrar el sistema 𝐴𝒙 = 𝑏 sino más bien:

𝐴𝒙 ≈ 𝑏 (19)

Debido a a que estamos considerando un error intrínseco en la determinación de los datos

observados.

A partir de esta expresión tenemos que considerar el vector de residuos que cumple con:

𝐴𝑥 − 𝑏 = 𝑟

El tamaño de este vector 𝑟 por definición es:

‖𝑟‖22 = 𝑟𝑇𝑟

‖𝑟‖22 = (𝑏 − 𝐴𝒙)𝑇(𝑏 − 𝐴𝒙) = 𝑏𝑇𝑏 − 2𝒙𝑇𝐴𝑇𝑏 + 𝒙𝑻𝐴𝑇𝐴𝒙

El vector encontrado es el que queremos minimizar, para tal efecto necesitamos derivar esta

expresión respecto de cada dirección de 𝒙 e igualarla a 0.

𝜕

𝜕𝑥𝑖(𝑏2 − 2𝑥𝑖𝑎𝑗𝑖𝑏𝑖 + (𝑎∗)𝑖𝑗

2 𝑥𝑖2) = 0 ; 𝑎𝑗𝑖𝑎𝑖𝑗 = 𝑎∗



Lo que se traduce en:

2(𝑎∗)𝑖𝑗2 𝑥𝑖 − 2(𝑎𝑗𝑖𝑏𝑖) = 0

En forma matricial tenemos que esta expresión es:

𝐴𝑡𝐴𝒙 = 𝐴𝑡𝑏 (20)

A la expresión (20) la conocemos como el conjunto de ecuaciones normales, el cual nos

permitirá determinar el vector solución 𝒙 que buscamos.

3.4.2. El modelo de regresión lineal

El modelo de regresión lineal múltiple funciona bajo el mismo principio que el de regresión

lineal simple. Es decir, si tenemos n descriptores o variables independientes asociadas a una

variable dependiente suponemos que existe una combinación lineal que describe el

comportamiento de esta variable dependiente.

Cabe mencionar que los datos deben cumplir las siguientes características

a) Linealidad: Todas las variables funcionan bajo un modelo lineal 𝑦 = 𝑚𝑥 + 𝑏.

b) Homocedasticidad: Significa que la varianza de los errores que consideremos es la

misma (o aproximadamente) para todos.

c) Independencia: los errores aleatorios son independientes entre sí

d) Normalidad: Los errores se colocan a través de la distribución normal con media 0 y

alguna varianza fija 𝜎02

Nuestro problema es que consideramos que los datos observados (𝑦𝑖) respondan teóricamente a

un modelo lineal como:

𝑦𝑖 = 𝑏0 ∙ 1 + 𝑏1 ∙ 𝑥1𝑗 + ⋯+ 𝑏𝑛 ∙ 𝑥𝑛𝑗 + 𝜖𝑖



Después el ajuste de nuestros datos será �̂�𝑖 y lo que queremos minimizar para cada uno de estos

valores es:

𝜖𝑖 = (𝑦𝑖– �̂�𝑖)2 (21)

De forma matricial obtenemos:

�̂� = [

1 𝑥11 … 𝑥1𝑛

1 𝑥21 … 𝑥2𝑛

⋮ ⋮ ⋮ ⋮1 𝑥𝑚1 … 𝑥𝑚𝑛

] [

𝑏0

𝑏1

⋮𝑏𝑛−1

]

Considerando el error 𝜖 en la expresión (21) implica de forma matricial que:

𝜖 = [

𝑦1 − �̂�1

𝑦2 − �̂�2

⋮𝑦𝑛 − �̂�𝑛

] = 𝑌 − �̂�

Pero sabemos que

𝜖 = 𝑌 − 𝑋𝑏

r Considerando la definición de varianza residual sabemos que 𝜎2 = 𝝐2 = ‖𝝐‖𝑇 ∙ ‖𝝐‖. En

nuestro caso esto define una función que depende de b y que llamaremos 𝑟(𝑏).

𝑟(𝑏) = ‖𝝐‖𝑡 ∙ ‖𝝐‖ (22)

La expresión anterior es la que queremos minimizar para encontrar el vector 𝑏 que nos

representará la forma para las pendientes de nuestras rectas.

𝑟(𝒃) = (𝑌 − 𝑋𝒃)𝑇(𝑦 − 𝑋𝒃)

De la cual necesitamos:

𝜕

𝜕𝑏𝑟(𝒃) =

𝜕(𝑌 − 𝑋𝒃)𝑇(𝑌 − 𝑋𝒃)

𝜕𝑏= 0

𝜕

𝜕𝑏𝑟(𝒃) = −𝑋𝑇𝑌 − 𝑋𝑇𝑌 + 2𝑋𝑇𝑋𝑏 = 0 ⇒ 𝑋𝑇𝑋𝒃 = 𝑋𝑇𝑌

Al aplicar la transpuesta de 𝑋𝑇𝑋 en ambos lados del último término obtenemos que:

(𝑋𝑇𝑋)−1𝒃 = (𝑋𝑇𝑋)−1𝑋𝑇𝑌

Por lo tanto nuestro vector 𝒃 queda determinado por la expresión:

𝒃 = (𝑋𝑇𝑋)−1𝑋𝑇𝑌



Cierre de la unidad

En esta unidad hemos revisado los métodos numéricos más utilizados para resolver sistemas de

ecuaciones lineales. Revisamos adicionalmente el modelo de regresión lineal. Ambos temas

tienen una gran variedad de aplicaciones en las distintas áreas del conocimiento, tanto en la

ingeniería como en las ciencias sociales.

Para saber más:

Te recomendamos los siguientes links, para revisar información importante sobre integración

numérica, auxiliado con el software.

Weisstein, Eric W. (2013) "Matrix Norm."FromMathWorld--A Wolfram Web

Resource.http://mathworld.wolfram.com/MatrixNorm.html

Brookes, M., (2011) "The Matrix Reference Manual". [on line]

http://www.ee.imperial.ac.uk/hp/staff/dmb/matrix/intro.html

Black, Noel; Moore, Shirley; and Weisstein, Eric W. (2013) "Jacobi Method."From Math World-

-A Wolfram Web Resource.

http://mathworld.wolfram.com/JacobiMethod.html

Peter (2012) "Jacobi method" From CFD on line. [on line]

http://www.cfd-online.com/Wiki/Jacobi_method

Resource.http://mathworld.wolfram.com/JacobiMethod.html

Peter (2012) "Jacobi method" From CFD on line. [on line]


http://mathworld.wolfram.com/MatrixNorm.html

http://www.ee.imperial.ac.uk/hp/staff/dmb/matrix/intro.html







GraphPad Software (2012) [on line] http://graphpad.com/guides/prism/6/curve-fitting/

Referencias Bibliográficas

Burden, R. (2011) Análisis numérico (7ª edición) México: Cengage Learning.

Mathews, J., Fink, K. (2000). Métodos Numéricos con MATLAB. (3ª edición) Madrid, España.

Prentince Hall.

Álvarez, L, Martínez, A. (s.f.). Tema 3: resolución de sistemas lineales y no lineales. Recuperado

de:

http://www.dma.uvigo.es/~lino/Tema3.pdf

Ajuste. (s.f.) Recuperado de:

http://exa.unne.edu.ar/matematica/metodos/5-3-material-teorico/tema_6_ajuste_2011.pdf

Angarita, A. (2013). Apuntes de análisis numérico. Recuperado de:

http://sb46f5727470feb20.jimcontent.com/download/version/1449503459/module/8900648268/n

ame/Apuntes%20de%20Analisis%20Numerico.pdf

Jiménez López, V. Pallarés Ruiz, A. (s.f.) Métodos numéricos. Recuperado de:

http://www.um.es/docencia/vjimenez/ficheros/textos/metodosnumericos.pdf

http://www.dma.uvigo.es/~lino/Tema3.pdf

http://exa.unne.edu.ar/matematica/metodos/5-3-material-teorico/tema_6_ajuste_2011.pdf

http://sb46f5727470feb20.jimcontent.com/download/version/1449503459/module/8900648268/name/Apuntes%20de%20Analisis%20Numerico.pdf

http://sb46f5727470feb20.jimcontent.com/download/version/1449503459/module/8900648268/name/Apuntes%20de%20Analisis%20Numerico.pdf

http://www.um.es/docencia/vjimenez/ficheros/textos/metodosnumericos.pdf

Análisis Numérico I 4° Semestre

Documents

Transcript of Análisis Numérico I 4° Semestre