Metodos Numericos

Métodos Numéricos
para Ingenierı́a
usando Python
W ILLIAM L A C RUZ
Universidad Central de Venezuela
Facultad de Ingenierı́a
Escuela de Ingenierı́a Eléctrica
Departamento de Electrónica, Computación y Control
c 2019 W. La Cruz
Contenido
Contenido ii
Antes de leer y usar los programas de este libro iii
1 Herramientas Matemáticas Básicas 1

1.1 Conceptos Básicos del Álgebra Lineal . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Espacio Vectorial, Producto Escalar y Normas Vectoriales . . . . . . . 1
1.1.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.1.3 Normas Matriciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.4 Alcance y Rango de una Matriz . . . . . . . . . . . . . . . . . . . . . 12
1.1.5 Matrices Especiales . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1.6 Autovalores y Autovectores . . . . . . . . . . . . . . . . . . . . . . . 15
1.2 Conceptos Básicos del Cálculo . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.1 Funciones de una Variable Real . . . . . . . . . . . . . . . . . . . . . 16
1.2.2 Funciones de Varias Variables . . . . . . . . . . . . . . . . . . . . . . 21
1.3 Calculando con Números . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.3.1 Representación Numérica . . . . . . . . . . . . . . . . . . . . . . . . 24
1.3.2 Representación de Punto-Flotante . . . . . . . . . . . . . . . . . . . . 25
1.3.3 Eficiencia, Estabilidad y Condicionamiento . . . . . . . . . . . . . . . 31
1.3.4 Velocidad de Convergencia . . . . . . . . . . . . . . . . . . . . . . . 38
1.4 Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Bibliografı́a 43
ii
Antes de leer y usar los programas
de este libro
Todos los programas en este libro se escribieron en el lenguaje Python versión 3.6
(https://www.python.org/), para una computadora personal. Los programas del libro
conforman el módulo calnum, que está disponible en el sitio Web:
http://neutron.ing.ucv.ve/cruzw/metnum.html
A lo largo del libro se asume que todos los objetos del módulo calnum están cargados
previamente, a través de la instrucción:
>>> from calnum import ∗
con la cual se cargan todos los programas del módulo.
En el libro también se utilizan las siguientes librerı́as cientı́ficas y gráficas de Python:
• Numpy: http://www.numpy.org/
• Scipy: https://www.scipy.org/
• Mapplotlib: https://matplotlib.org/
Los objetos de la librerı́a Numpy se utilizan con el prefijo np (siguiendo la sintaxis de una
clase), ya que se asume que Numpy se ha cargado previamente con la instrucción:
>>> import numpy as np
iii
1 Herramientas Matemáticas
Básicas
En este libro se usan sistemáticamente conceptos matemáticos básicos que el lector

deberı́a saber, es posible que no los recuerde de forma inmediata. Emplearemos este
capı́tulo para recordar tales conceptos e introducir nuevos conceptos que pertenecen al
Análisis Numérico. Comenzamos condensando nociones que son tı́picas de los cursos de
Álgebra Lineal y Cálculo, reformulándolas de manera adecuada para su uso en el Cálculo
Cientı́fico. Solo comentaremos resultados fundamentales, por lo cual hemos evitado las
demostraciones de los mismos, ya que esto escapa del objetivo del libro. Finalizamos el
capı́tulo calculando con números, esto es, exploramos la aritmética en computadora. En
las siguientes referencias se pueden encontrar las demostraciones de cada uno de los teo-
remas que enunciaremos en este capı́tulo [1, 20, 12, 23, 26, 10, 4].
1.1 Conceptos Básicos del Álgebra Lineal

En esta sección haremos una revisión de los conceptos básicos del algebra lineal tales
como: vectores, espacio vectorial, producto escalar, norma, matrices, inversa y determi-
nante de una matriz, norma de una matriz, autovalores y autovectores.
1.1.1 Espacio Vectorial, Producto Escalar y Normas Vectoriales

El concepto de espacio vectorial es uno de los más importantes en las Matemáticas.
Definición 1.1. Un conjunto no vacı́o V de elementos x, y, z, . . . se llama espacio vecto-

rial, cuando satisface las siguientes condiciones:
I. Para cualesquiera dos elementos x, y ∈ V está definido unı́covamente un tercer
elemento w ∈ V , llamado suma de ellos y denotado por x + y tal que
1) x + y = y + x (conmutatividad),
2) x + (y + z) = (x + y) + z (asociatividad),
3) existe en V un elemento 0 tal que x + 0 = x para todo x ∈ V ,
4) para todo x ∈ V existe un elemento −x tal que x + (−x) = 0.
II. Para cualquier número α y cualquier elemento x ∈ V está definido el elemento
αx ∈ V de manera que:
1) α(βx) = (αβ)x, para cualquier número β,
2) 1 · x = x.
III. Las operaciones de adición y multiplicación están relacionadas entre sı́ mediante
las leyes distributivas:
1) (α + β)x = αx + βx, para cualesquiera números α y β,
2) α(x + y) = αx + αy, para cualquier número α.
1
2 1.1. CONCEPTOS BÁSICOS DEL ÁLGEBRA LINEAL
A los elementos de un espacio vectorial se les acostumbra llamar vectores. Veamos dos
ejemplos de espacios vectoriales, dejando como ejercicio para el lector la comprobación,
en cada uno de ellos, de los axiomas enunciados en la Definición 1.1.
Ejemplo 1.1.
1. El conjunto de los números reales, denotado por R, con las operaciones habituales
de adición y multiplicación representa un espacio vectorial.
2. El conjunto formado por las n-uplas de números reales, esto es, x = (x1 , x2 , . . . , xn ),
donde xi ∈ R para todo i = 1, 2, . . . , n, es un espacio vectorial. Éste se denota por
Rn . Los elementos de Rn también se pueden representar como un vector columna,
 
x1
 x2 
 
x =  . .
 .. 
xn
A lo largo del libro denotaremos un vector de Rn como una n-upla o un vector columna.
Definición 1.2. Decimos que un conjunto no vacı́o W de un espacio vectorial V es un

subespacio de V si, y sólo si W es un espacio vectorial. En particular, el conjunto W de
las combinaciones lineales de un sistema de p vectores, {v1 , . . . , vp }, es un subespacio
de V , llamado subespacio generado por el sistema de vectores y es denotado por
W = span(v1 , . . . , vp )
= {v = αi v1 + · · · + αp vp : αi ∈ R}.
Definición 1.3. Un sistema de vectores {v1 , . . . , vm } de un espacio vectorial V es lla-

mado linealmente independiente si la relación
αi v1 + · · · + αm vm = 0
con α1 , . . . , αm ∈ R implica que α1 = · · · = αm = 0. De lo contrario, el sistema se dice

linealmente dependiente, esto es, existen escalares α1 , . . . , αm no todos nulos tales que
αi v1 + · · · + αm vm = 0.
Definición 1.4. Si en un espacio vectorial V se pueden encontrar n vectores lineal-

mente independientes y cualesquiera n + 1 vectores de este espacio son linealmente
dependientes, entonces se dice que el espacio V tiene dimensión n. Se llama base de un
espacio vectorial V de n dimensiones, a todo conjunto de n vectores linealmente inde-
pendientes. Se dice que cualquier espacio vectorial que tenga una base con un número
finito de elementos es de dimensión finita. Otros espacios vectoriales se conocen como
de dimensión infinita.
CAPÍTULO 1. HERRAMIENTAS MATEMÁTICAS BÁSICAS 3
Es fácil comprobar que Rn es de dimensión n, además, una base de Rn , denominada

base canónica, es el conjunto de vectores e1 , e2 , . . . , en definidos por
ei = (0, 0, . . . , 0, 1 , 0, . . . , 0),
↑
i-ésima componente
es decir, todas las componentes de ei son nulas salvo la i-ésima que es igual a 1.
Pasemos ahora a definir producto escalar, con el cual podremos introducir el concepto
de vectores ortogonales o perpendiculares.
Definición 1.5. El producto escalar de dos vectores x e y del espacio vectorial Rn , es un

número real hx, yi que verifica las siguientes condiciones:
1) hx, yi = hy, xi,
2) hx, y + zi = hx, yi + hx, zi, para todo z ∈ V ,
3) hλx, yi = λhx, yi, para todo escalar λ,
4) hx, xi ≥ 0 y hx, xi = 0 sólo si x = 0.
En adelante, consideraremos el espacio vectorial Rn . El producto escalar en Rn que

utilizaremos a lo largo del libro se define como Se deja como
ejercicio verifi-
n
X car que el pro-
T
hx, yi ≡ x y = xi y i , ducto escalar
i=1 hx, yi ≡ xT y
en Rn satis-
para x = (x1 , x2 , . . . , xn ) e y = (y1 , y2 , . . . , yn ), donde xT se lee traspuesta de x, cuyo face los axio-
mas de la Defi-
concepto se introducirá en la Sección 1.1.2.
nición 1.5
Otro concepto importante en álgebra lineal es la norma o longitud de un vector. Segui-
damente mostramos la definición formal de norma.
Definición 1.6. Una norma vectorial, denotada por el sı́mbolo k · k, es una función real
definida en el espacio vectorial Rn , que satisface las siguientes propiedades:
1) kxk ≥ 0, para x ∈ Rn y kxk = 0 si y sólo si x = 0.
2) kαxk = |α| kxk, para x ∈ Rn y α ∈ R.
3) kx + yk ≤ kxk + kyk, para x, y ∈ Rn .
En palabras, las condiciones de la Definición 1.6 requieren que: 1) la norma de un

vector no nulo es positiva, 2) la norma de un vector escalado es igual a la norma del vector
escalada por la misma cantidad, y 3) la norma de la suma vectorial no excede la suma de
las normas de sus partes.
Las normas vectoriales más comúnmente usadas son los casos especiales de las normas
de Hölder: !1/p
n
X
p
kxkp = |xi | , 1 ≤ p ≤ ∞.
i=1
Los casos particulares de p = 1, p = 2 y p = ∞, son las normas más importantes en la

práctica y las mismas se definen a continuación. La bola cerrada {x ∈ Rn : kxk ≤ 1}
correspondiente de cada norma se muestra a la derecha para el caso n = 2.
n
X
kxk1 = |xi | = |x1 | + |x2 | + · · · + |xn |,
i=1
n
!1/2
X 1/2
kxk2 = x2i = x21 + x22 + · · · + x2n ,
i=1
kxk∞ = máx |xi |.

i=1,...,n
La norma-2 también se denomina norma Euclı́dea. Las normas que más emplearemos
en este libro son la norma-2 y su norma matricial inducida, que más adelante definiremos.
El concepto de ortogonalidad es importante en muchas aplicaciones, particularmente
en álgebra lineal.
Definición 1.7. Dos vectores x, y ∈ Rn , se dicen ortogonales si

n
X
xT y = xi yi = 0.
i=1
Un vector x para el cual kxk2 = 1 se dice normalizado. El conjunto de vectores

{x1 , x2 , . . . , xn } se dice ortonormal, si ellos son mutuamente ortogonales y normali-
zados, esto es, (
0, si i 6= j,
xTi xj =
1, si i = j.
Ejemplo 1.2. La base canónica de Rn es un conjunto ortonormal. En efecto, se tiene que
eTi ej = (0, . . . , 0, 1 , 0, . . . , 0)T (0, . . . , 0, 1 , 0, . . . , 0) = 0,

↑ ↑
i-ésima componente j-ésima componente
siempre que i 6= j; además, es fácil ver que eTi ei = 1.
El ángulo entre dos vectores de Rn se define de la siguiente manera.
Definición 1.8. El águlo θ entre dos vectores x e y de Rn está dado por

xT y
cos(θ) = .
kxk2 kyk2
De esta forma, dos vectores x e y son ortogonales si el ángulo entre ellos es θ = 90◦ .
Una propiedad importante de las normas de Hölder es la desigualdad de Hölder,
xT y ≤ kxkp kykq ,
donde
1 1
+ = 1.
p q
Un caso especial de la desigualdad de Hölder es la desigualdad de Cauchy-Schwarz,
xT y ≤ kxk2 kyk2 . (1.1)
En muchos métodos numéricos es necesario construir un sistema de vectores ortogo-

nales. El proceso de Gram-Schmidt se utiliza para generar un conjunto de vectores orto-
gonales (u ortonormales si se normaliza) a partir de un conjunto de vectores linealmente
independientes.
Teorema 1.1 (Proceso de Gram-Schmidt). Sea {v1 , . . . , vm } un sistema de vectores li-

nealmente independientes de Rn . Para i = 1, . . . , m se definen los vectores ui de manera
recursiva,
i−1
X
vi − (viT uj )uj
j=1
ui = .
i−1
X
vi − (vi uj )uj
j=1
2
Entonces, {u1 , . . . , um } es un sistema de vectores ortonormales.
1.1.2 Matrices
Una matriz es un concepto muy importante en álgebra lineal, que tiene muchas aplica-
ciones en ingenierı́a y en las ciencias aplicadas.
Definición 1.9. Una matriz real A = (aij ) de orden m × n es un arreglo rectangular

de números reales aij , i = 1, . . . , m, j = 1, . . . , n. El conjunto de todas las matrices
reales de orden m × n se denota por Rm×n . Las matrices de orden m × 1 se llaman
vectores columna o simplemente vectores; en cambio, las matrices de orden 1 × n se
llaman vectores fila.
Definición 1.10. La traspuesta de una matriz A ∈ Rm×n es la matriz AT ∈ Rn×m tal

que (AT )ij = aji .
En el siguiente ejemplo mostramos algunas matrices y sus traspuestas.

Ejemplo 1.3. Observe las matrices A ∈ R3×3 , B ∈ R1×4 y C ∈ R3×1 , con sus respectivas
traspuestas.
  
3 −1.1 0 4
A= 0 2 0.4 , B = 1 −3.4 1
2 9 , C = −0.001 ,
−1 1 2.1   2
  1
3 0 −1 −3.4
AT = −1.1 2 1 , B T =  
 1/2  , C T = 4 −0.001 2 .
0 0.4 2.1
9
Aprovechemos este ejemplo para comenzar a trabajar con Python. Las matrices A, B y
C pueden ser creadas y manipuladas en Python. Para ello, se utiliza la librerı́a Numpy.
Seguidamente se muestra la manera de hacerlo. Primero, en el Shell (escritorio de trabajo)
de Python, cargamos la librerı́a Numpy.
>>> import numpy as np
Con estas instrucciones se carga la librerı́a Numpy con el nombre np. Todas las funciones
disponibles en Numpy se utilizan con la siguiente sintaxis: np.funcion (siguiendo las reglas
de una clase). A lo largo del escrito utlizaremos el prefijo np para denotar la librerı́a Numpy,
la cual, por supuesto, debe cargarse previamente con el nombre np. Ahora bien, para crear
una matriz se pueden utilizar los comandos array o matrix de Numpy. Por ejemplo, las
matrices A, B y C dadas arriba se crean de la siguiente manera:
>>> A = np . a r r a y ( [ [ 3 , − 1 . 1 , 0 ] , [ 0 , 2 , 0 . 4 ] , [ − 1 , 1 , 2 . 1 ] ] )
>>> B = np . a r r a y ( [ 1 , − 3 . 4 , 0 . 5 , 9 ] )
>>> C = np . a r r a y ( [ [ 4 ] , [ − 0 . 0 0 1 ] , [ 2 ] ] )
>>> p r i n t (A)
[ [ 3 . −1.1 0 . ]
[ 0. 2. 0.4]
[ −1. 1. 2.1]]
>>> p r i n t (B)
[ 1 . −3.4 0 . 5 9 . ]
>>> p r i n t (C)
[ [ 4.00000000e+00]
[ −1.00000000e−03]
[ 2.00000000e+00]]
Para hallar la traspuesta de una matriz se utiliza uno de los comandos .transpose() o .T,
como se indica a continuación:
>>> A . T
a r r a y ( [ [ 3 . , 0 . , −1. ] ,
[ −1.1 , 2 . , 1 . ] ,
[ 0. , 0.4 , 2.1]])
>>> C . t r a n s p o s e ( )
a r r a y ( [ [ 4.00000000e+00, −1.00000000e −03, 2.00000000e +00]])
Definamos las principales operaciones con matrices.

Definición 1.11. Sean A = (aij ), B = (bij ) y C = (cij ) matrices.
• Suma: La suma de las matrices A ∈ Rm×n y B ∈ Rm×n , es una matriz C = A + B

tal que
cij = aij + bij , i = 1, . . . , m, j = 1, . . . , n.
• Multiplicación por un escalar: La multiplicación de un número real α por una

matriz A ∈ Rm×n , es una matriz C = αA tal que
cij = αaij , i = 1, . . . , m, j = 1, . . . , n.
• Multiplicación por otra matriz: La multiplicación de dos matrices A ∈ Rm×p y

B ∈ Rp×n , es una matriz C = AB tal que C ∈ Rm×n y
p
X
cij = aik bkj , i = 1, . . . , m, j = 1, . . . , n.
k=1
En general, el producto de matrices no es conmutativo.
El conjunto de las matrices reales Rm×n con las operaciones algebraicas dadas en la
definición anterior, conforma un espacio vectorial.
Ejemplo 1.4.
         
−1 3 2 5 1 8 −1 3 −4 12
 5 4 + −6 −2 = −1 2 , 4 5 4 =  20 16 ,
0 −2 8 6 8 4 0 −2 0 −8
    
3 −1 0 2 12
 0   
2 4 −6 = 20 . 
−1 1 2 8 8
En Python las operaciones: A + B, A − B, aC y AB, se realizan respectivamente como:

A+B, A-B, a*C y A.dot(B). Las operaciones anteriores se realizan de la siguiente manera:
>>> np . a r r a y ([[ −1 ,3] ,[5 ,4] ,[0 , −2]])+ np . a r r a y ( [ [ 2 , 5 ] , [ − 6 , − 2 ] , [ 8 , 6 ] ] )

array ([[ 1 , 8] ,
[−1, 2 ] ,
[ 8 , 4]])
>>> 4∗np . a r r a y ( [ [ − 1 , 3 ] , [ 5 , 4 ] , [ 0 , − 2 ] ] )
a r r a y ([[ −4 , 12] ,
[20 , 16] ,
[ 0 , −8]])
>>> np . a r r a y ( [ [ 3 , − 1 , 0 ] , [ 0 , 2 , 4 ] , [ − 1 , 1 , 2 ] ] ) . dot ( np . a r r a y ( [ [ 2 ] , [ − 6 ] , [ 8 ] ] ) )
array ([[12] ,
[20] ,
[ 8]])
Definición 1.12. Una matriz A ∈ Rm×n es cuadrada si tiene el mismo número de co-
lumnas y filas, es decir, si m = n. Cuando una matriz cuadrada A tiene la propiedad
que AT = A, se dice que A es simétrica.
Una matriz cuadrada importante es la matriz identidad de orden n dada por:

 
1 0 0 ··· 0
0 1 0 · · · 0
 
 
I = 0 0 1 · · · 0 ,
 .. .. .. . . .. 
. . . . .
0 0 0 ··· 1
en otras palabras, I = (δij ) es una matriz de orden n tal que

(
1, i = j,
δij =
0, i 6= j.
La matriz identidad satisface la igualdad IA = AI = A, para toda matriz cuadrada A de

orden n. Con el comando eye(n) de la librerı́a Numpy se crea la matriz identidad de orden
n. A continuación creamos la matriz identidad de orden 3.
>>> np . eye ( 3 )
array ([[ 1. , 0. , 0.] ,
[ 0. , 1. , 0.] ,
[ 0. , 0. , 1.]])
Definición 1.13. Una matriz cuadrada A de orden n se dice invertible o no singular, si

existe una matriz cuadrada B de orden n tal que BA = AB = I. B es llamada matriz
inversa de A y se denota por A−1 . Una matriz cuando no es invertible se denomina
singular.
Teorema 1.2. Si A y B son matrices no singulares de orden n, entonces:
i) (A−1 )−1 = A.
ii) (AB)−1 = B −1 A−1 .
Ahora definiremos un número asociado con la matriz cuadrada A = (aij ) de orden n,

llamado determinante de la matriz y denotado por det(A). El determinante de orden n
se define en términos de los determinantes de orden n − 1, junto con la condición que el
determinante de una matriz 1 × 1 de un sólo elemento, A = (a), está dado por
det(A) = a.
Esto significa que el determinante de una matriz de orden 2 × 2 está definido en términos
de determinantes de orden 1 × 1. El determinante de una matriz de orden 3 × 3 está
definido en términos de determinantes de orden 2 × 2, y ası́ sucesivamente. Para precisar,
introducimos las siguientes definiciones.
Definición 1.14. El determinante de una matrix de orden (n − 1) obtenida eliminando

la i-ésima fila y la j-ésima columna de una matrix A = (aij ) de orden n, es llamado el
menor de aij y es denotado por Mij . El número
Aij = (−1)i+j Mij ,
es llamado el cofactor de aij .
Definición 1.15. El determinante de una matriz cuadrada A de orden n está definido

por la expresión
Xn
det(A) = a1j A1j ,
j=1
conjuntamente con la definición del det(A) para n = 1. En otras palabras, el determi-

nante de A es la suma de los productos de los elementos en la primera fila de la matriz
por sus respectivos cofactores.
Teorema 1.3. Si A y B son matrices cuadradas de orden n, entonces:

i) det (AB) = det (BA) = det (A) det (B).
ii) det (AT ) = det (A).
iii) det (αA) = αn det (A).
iv) det (I) = 1.
v) A es no singular si, y sólo si det(A) 6= 0.
Definición 1.16. Sean A = (aij ) una matriz cuadrada y Aij el cofactor de aij . La
traspuesta de la matriz de cofactores se conoce como la adjunta de la matriz A, denotada
por adj A,
adj A = (Aji ).
Teorema 1.4. Si A es una matriz no singular, entonces la inversa de A está dada por
1
A−1 = adj A.
det A
Definición 1.17. La traza de una matriz cuadrada A = (aij ) de orden n es la suma de

n
P
los elementos de su diagonal, esto es, tr(A) = aii .
i=1
En el siguiente ejemplo ilustramos los conceptos vistos hasta ahora. Aprovechemos en

este ejemplo para introducir los comandos de Python que permiten calcular la inversa, el
determinante y la traza de una matriz.
Ejemplo 1.5. Consideremos la matriz

 
1 0 0
A = 0 2 1
0 −1 3
Se tiene que la inversa, la traza y el determinante de A son:

 
1 0 0
A−1 = 0 3/7 −1/7 , tr(A) = 6, y det(A) = 1 · det(A11 ) = 2 · det(3) − 1 · det(−1) = 7.
0 1/7 2/7
En Python, el determinante y la inversa de una matriz se calculan, respectivamente, con

los comandos det e inv del módulo linalg de la librerı́a Numpy, y La traza se calcula con el
comando trace de la librerı́a Numpy. Con la matriz A anterior se tiene:
>>> from numpy . l i n a l g import inv , d e t
>>> from numpy import t r a c e
>>> A = np . a r r a y ( [ [ 1 , 0 , 0 ] , [ 0 , 2 , 1 ] , [ 0 , − 1 , 3 ] ] )
>>> d e t (A)
7.0000000000000009
>>> t r a c e (A)
6
>>> i n v (A)
array ([[ 1. , 0. , 0. ],
[ 0. , 0.42857143 , −0.14285714] ,
[ 0. , 0.14285714 , 0.28571429]])
Observe que las instrucciones

from numpy . l i n a l g import inv , d e t
from numpy import t r a c e
cargan solamente las funciones det, inv y trace, sin cargar las otras funciones de la librerı́a.
En algunas aplicaciones, es necesario calcular la inversa de una matriz B la cual di-

fiere de una matriz singular A solamente en una perturbación de rango 1. La fórmula de
Sherman-Morrison permite calcular la inversa de B.
Teorema 1.5 (Fórmula de Sherman-Morrison). Si u y v son dos vectores de Rn y A es

una matriz no singular de orden n, entonces
1
(A − uvT )−1 = A−1 + α(A−1 uvT A−1 ), con α = , si vT A−1 u 6= 1.
(1 − vT A−1 u)
Es importante subrayar que son necesarias n! operaciones (multiplicaciones y sumas)

para calcular el determinante de una matriz de orden n empleando la definición recursiva.
Por ello el cálculo del determinante es prohibitivo cuando la dimensión de la matriz es
grande (n ≫ 1000); por ejemplo, si n = 100, entonces son necesarias aproximadamente
9.3326 × 10157 operaciones para calcular el determinante. Como el cálculo de la inversa de
una matriz puede involucrar el cálculo de determinantes, entonces obtener la inversa de
matrices grandes es computacionalmente costoso. De esta forma, si la inversa o el determi-

nante no están disponibles, mejor serı́a no calcularlos directamente usando su definición,
sino buscar métodos alternativos, más adelante hablaremos de ello.
Otra clase de matrices cuadradas importantes son las matrices definidas positivas.
Definición 1.18. Una matriz cuadrada A de orden n es definida positiva, si xT Ax > 0

para todo vector x ∈ Rn no nulo.
Ejemplo 1.6. La matriz A del Ejemplo 1.5 es definida positiva. En efecto, para todo
0 6= x ∈ R3 se tiene:
  
1 0 0 x1
T
x Ax = x1 x2 x3 0  2 1   x2  = x21 + 2x22 + 3x23 > 0.
0 −1 3 x2
Generalmente, verificar que una matriz es definida positiva a partir de la definición

no es una tarea sencilla, ya que la condición xT Ax > 0 se debe satisfacer para un vector
arbitrario x 6= 0. Más adelante estudiaremos algunos procedimientos para determinar
cuando una matriz es definida positiva.
1.1.3 Normas Matriciales

Una matriz de orden m × n se puede ver como un vector en el espacio nm-dimensional:
a partir de la matriz A = (aij ) se construye, por ejemplo, el vector
(a11 , a21 , . . . , am1 , . . . , a1n , a2n , . . . , amn ) ∈ Rmn .
Por tanto, cualquier norma vectorial en Rmn puede medir el “tamaño” de tal matriz. Sin
embargo, al tratar con el espacio vectorial de las matrices existen normas especiales que
son más útiles que las normas vectoriales definidas previamente. Éstas son las normas
inducidas.
Definición 1.19. Para una matriz A ∈ Rm×n , se define la norma de A como
kAxkp
kAkpq = máx .
x∈R n , x6=0 kxkq
La norma matricial k · kpq es inducida por las normas vectoriales k · kp y k · kq .
Las normas matriciales satisfacen las propiedades usuales de las normas, esto es,
1. kAk = 0, si y sólo si A = 0, para A ∈ Rm×n .
2. kαAk = |α| kAk, para A ∈ Rm×n y α ∈ R.
3. kA + Bk ≤ kAk + kBk, para A, B ∈ Rm×n .

El caso p = q es de interés muy particular y la norma asociada k · kpq es simplemente

denotada por k · kp y se llama norma p. Los casos más importantes son de nuevo los
asociados con p, q = 1, 2, ∞. Para p = 1 y p = ∞, se tienen las siguientes definiciones
alternativas:
m
X X n
kAk1 = máx |aij |, kAk∞ = máx |aij |.
j=1,...,n i=1,...,n
i=1 j=1
La norma matricial más importante que no es inducida por una norma vectorial es la
norma de Frobenius, definida por
 1/2
Xm X n
kAkF = tr(AT A) =  a2ij  , para A ∈ Rm×n .
i=1 j=1
Ejemplo 1.7. Para la matriz  

1 4 6
A = 3 −1 5
0 8 7
se tiene
kAk1 = máx {1 + 3 + 0, 4 + | − 1| + 8, 6 + 5 + 7} = 18,
kAk∞ = máx {1 + 4 + 6, 3 + | − 1| + 5, 0 + 8 + 7} = 15,
1713 √
kAk2 = , kAkF = 201.
130
El comando norm del módulo linalg de la librerı́a Numpy, se utiliza para calcular las normas
p y la norma de Frobenius. Para la matriz A,
>>> A = np . a r r a y ( [ [ 1 , 4 , 6 ] , [ 3 , − 1 , 5 ] , [ 0 , 8 , 7 ] ] )
se tiene que las normas calculadas anteriormente están dadas respectivamente por:
>>> np . l i n a l g . norm (A , 1 )
18.0
>>> np . l i n a l g . norm (A , np . i n f )
15.0
>>> np . l i n a l g . norm (A , 2 )
13.176932423207695
>>> np . l i n a l g . norm (A , ’ f r o ’ )
14.177446878757825
1.1.4 Alcance y Rango de una Matriz

Un subespacio importante asociado con una matriz A ∈ Rm×n es su alcance o espacio
de llegada, definido por
alc(A) = {y ∈ Rm : y = Ax, x ∈ Rn }.
El rango de una matriz A, denotado por rang(A), es la dimensión del alcance de A, es

decir, el número de columnas linealmente independientes de A.
Ejemplo 1.8. El rango de la matriz

 
1 4 6
A1 = 3 −1 5
0 8 7
es rang(A1 ) = 3. Ahora, el rango de la matriz

 
1 0 2
A2 = 0 1 0
1 1 2
es rang(A2 ) = 2. Con el comando matrix_rank del módulo linalg de la librerı́a Numpy, se

calcula el rango de una matriz en Python:
>>> A1 = np . a r r a y ( [ [ 1 , 4 , 6 ] , [ 3 , − 1 , 5 ] , [ 0 , 8 , 7 ] ] )
>>> A2 = np . a r r a y ( [ [ 1 , 0 , 2 ] , [ 0 , 1 , 0 ] , [ 1 , 1 , 2 ] ] )
>>> np . l i n a l g . m a t r i x r a n k (A1)
3
>>> np . l i n a l g . m a t r i x r a n k (A2)
2
1.1.5 Matrices Especiales

Algunas matrices tienen estructuras particulares que a menudo son convenientes para
propósitos computacionales. La siguiente lista, aunque incompleta, da una idea de estas
matrices especiales que desempeñan un papel importante en el análisis numérico y en
aplicaciones del cálculo cientı́fico.
Matriz Diagonal
Las matrices diagonales son matrices cuadradas de orden n de la forma
 
d11
 d22 
 
D= .. ,
 . 
dnn
donde dij = 0 para i 6= j. Utilizamos la notación D = diag(d11 , d22 , . . . , dnn ) para indicar
que D es una matriz diagonal.
El determinante de una matriz diagonal es el producto de los elementos de su diagonal.
Una matriz diagonal D tiene inversa si todos los elementos de su diagonal son distintos de
cero y  −1 
d11
 d−1 
 22 
D −1 =  .. .
 . 
d−1
nn
Matrices Triangulares
Una matriz cuadrada A de orden n se dice triangular superior si aij = 0 para i > j,
mientras que ésta se dice triangular inferior si aij = 0 para i < j. Seguidamente se muestra
la forma de las matrices triangular inferior y triangular superior, respectivamente.
   
l11 0 · · · 0 u11 u12 · · · u1n
 l21 l22 · · · 0   0 u22 · · · u1n 
   
L= . .. . . ..  , U =  .. .. .. ..  .
 .. . . .   . . . . 
ln1 ln2 · · · lnn 0 0 · · · unn
Recordemos algunas propiedades algebraicas de las matrices triangulares que son fáciles
de verificar.
• El determinante de una matriz triangular superior (inferior) es el producto de los

elementos de su diagonal;
• la inversa de una matriz triangular inferior (respectivamente, superior) también es

una matriz triangular inferior (respectivamente, superior);
• el producto de dos matrices triangulares inferiores (respectivamente, superiores) es

una matriz triangular inferior (respectivamente, superior).
Matrices Bidiagonales
Una matriz cuadrada B de orden n se dice bidiagonal superior si bij = 0 para j 6= i o j 6=

i+1, mientras ésta se dice bidiagonal inferior si bij = 0 para j 6= i o j 6= i−1. Seguidamente
se muestra la forma de las matrices bidiagonales inferior y superior, respectivamente.
   
a11 b11 b12
a21 a22   .. 
   b 22 . 
A= . .  , B= 


 .. ..   ..
. bn−1,n 
an,n−1 ann bnn
Matrices Tridiagonales
Una matriz cuadrada T de orden n se dice tridiagonal si tij = 0 para i, j tal que |j − i| >
1. Seguidamente se muestra la forma de una matriz tridiagonal.
 
t11 t12
 .. 
t21 t22 . 
T = 


 .. ..
. . tn−1,n 
tn,n−1 tnn
Matriz de Hessenberg superior
Una matrix A = (aij ) se dice matriz de Hessenberg superior si ella tiene ceros debajo
de la primera subdiagonal principal, es decir, aij = 0 para i > j + 1. A continuación
mostramos ejemplos de matrices de Hessenberg.

   
1 8 0 1 −2 1 −1 3
7 3 8 −6 0 7 3 2
   

A = 0 2 −4 5 1 , 
B = 0 2 5
 
0 0 7 3 4 0 0 6
0 0 0 −8 5 5×5 0 0 0 5×3
Matriz Permutación
Una matriz cuadrada se dice matriz permutación si sus columnas son una permutación
de las columnas de la matriz identidad. Seguidamente se muestran ejemplos.
     
1 0 0 0 0 0 0 1 0 0 0 1
0 0 1 0 0 0 1 0 0 1 0 0
P1 = 
0 1 0 0
 , P2 = 
0 1 0 0
 , P3 = 
0 0 1 0
 .
0 0 0 1 4×4 1 0 0 0 4×4 1 0 0 0 4×4
1.1.6 Autovalores y Autovectores

Sea A una matriz cuadrada de orden n. El número λ ∈ C se llama autovalor de A, si
existe un vector no nulo x ∈ Rn tal que Ax = λx. El vector x es el autovector asociado con
el autovalor λ y el conjunto de los autovalores de A se denomina espectro de A, denotado
por σ(A). Diremos que x e y son respectivamente un autovector derecho y un autovector
izquierdo, asociados al autovalor λ, si se cumple
Ax = λx, yT A = λyT .
El número λ es la solución de la ecuación caracterı́stica
pA (λ) = det(A − λI) = 0,
donde pA (λ) es el polinomio caracterı́stico. Como éste último es un polinomio de grado n

con respecto a λ, existen n autovalores de A no necesariamente distintos. Las siguientes
identidades se cumplen:
n
Y n
X
det(A) = λi , tr(A) = λi ,
i=1 i=1
y como det(AT − λI) = det((A − λI)T ) = det(A − λI), se concluye que σ(AT ) = σ(A).
De todo lo anterior se obtiene que una matriz es Qnsingular si y sólo si ésta tiene al menos
un autovalor nulo, dado que pA (0) = det(A) = i=1 λi = 0. Ahora, como A es una matriz
real, pA (λ) es un polinomio con coeficientes reales; por lo tanto, los autovalores complejos
de A deberán aparecer en pares complejos conjugados.
El módulo máximo de los autovalores de A se denomina radio espectral de A y es de-
notado por: ρ(A) = máxλ∈σ(A) |λ|. Seguidamente damos el Teorema de Cayley-Hamilton.
Teorema 1.6 (Cayley-Hamilton). Si p(λ) es el polinomio caracterı́stico de A, entonces

p(A) = 0.
16 1.2. CONCEPTOS BÁSICOS DEL CÁLCULO
1.2 Conceptos Básicos del Cálculo

Seguidamente realizamos una revisión de los conceptos básicos del Cálculo tales como
función, lı́mite, continuidad, derivada, integral, Teorema del Valor Medio, desarrollo de
Taylor, entre otros. A lo largo del libro utilizamos la notación habitual de la Teorı́a de
Conjuntos, además, empleamos los siguientes conjuntos: N: conjunto de los números
naturales, Z: conjunto de los números enteros, Q: conjunto de los números racionales,
R: conjunto de los números reales y C: conjunto de los números complejos.
1.2.1 Funciones de una Variable Real

Para introducir el concepto de función es necesario considerar el plano cartesiano o
plano xy, que consiste en formar el producto cartesiano de R × R = R2 y es representado
geométricamente a través de dos rectas perpendiculares de referencia (llamadas ejes coor-
denados), uno horizontal (llamado eje de las x o de las abscisas), y el otro vertical (el eje de
las y o de las ordenadas); su punto de intersección, se indica por 0 y se denomina origen.
El par de números (x, y) con x, y ∈ R, se denomina par ordenado.
Definición 1.20. Una función f es un conjunto de pares ordenados (x, y) ninguno de

los cuales tiene el mismo primer elemento. Escribiremos y = f (x) para indicar que y es
imagen x bajo la función f ; el sı́mbolo f (x) se lee “f de x”, “f en x” o “f evaluada en
x”. El conjunto de x donde f está bien definida (que no se produzca división por cero)
se denomina dominio y el conjunto de imágenes de f se denomina rango.
En el Cálculo elemental se tiene interés en considerar funciones en las que el dominio y

el rango son conjuntos de números reales. Estas funciones se llaman funciones de variable
real o más brevemente funciones reales y se pueden representar geométricamente mediante
una gráfica en el plano xy. Ahora bien, es conveniente distinguir entre funciones, es decir,
es importante saber cuando dos funciones son iguales. Diremos que dos funciones f y g
son iguales, si f y g tienen el mismo dominio y f (x) = g(x) para todo x del dominio de f .
Seguidamente damos uno ejemplos de funciones de uso común en el cálculo.
Ejemplo 1.9. 1. Funciones lineales. Una función f definida para x ∈ R mediante la
fórmula f (x) = ax + b, donde a, b ∈ R.
2. Funciones potenciales. Una función f definida para x ∈ R mediante la fórmula f (x) =
xn , donde n es un entero positivo.
3. Polinomios. Un polinomio p es una función definida para x ∈ R por una ecuación de la
forma
n
X
n
p(x) = a0 + a1 x + · · · + an x = ak xk ,
k=1
donde los números a0 , a1 , . . . , an son los coeficientes del polinomio y el entero no nega-
tivo n es su grado (si an 6= 0).
En la Figura 1.1 se muestran una parte de las gráficas de las funciones f (x) = 2,
g(x) = 79 x + 1 y p(x) = 12 x4 − x2 .
y y y
p(x) = 21 x4 − x2
2 f (x) = 2
g(x) = 79 x + 1
x x x
(a) f (x) = 2, función cons- (b) g(x) = 97 x + 1, función li- (c) p(x) = 12 x4 − x2 , polino-
tante neal mio de grado 4
Figura 1.1. Tipos de funciones
Para funciones f y g se pueden construir la funciones: suma f + g, producto f · g

(o f g) y cociente f /g. El dominio de cada una de esta nuevas funciones viene dado por
la intersección de los dominios de las funciones involucradas, además, la función f /g,
generalmente, no está definida en los números x tales que g(x) = 0.
Definición 1.21. Una función f se dice que es creciente en un conjunto S ⊂ R, si f (x) ≥

f (y) para x, y ∈ S con x < y. Si se verifica la desigualdad estricta f (x) < f (y) para todo
x < y en S se dice que la función es estrictamente creciente. Análogamente, una función
se dice decreciente en S si f (x) ≤ f (y) para x, y ∈ S con x < y. Si f (x) > f (y) para todo
x < y en S se dice que la función es estrictamente decreciente. Una función se denomina
monótona en S si es creciente en S o decreciente en S. Estrictamente monótona significa
que f , o es estrictamente creciente en S o es estrictamente decreciente en S.
Describamos ahora el concepto de continuidad, uno de los más importantes en el

Cálculo. Para ello es necesario introducir el concepto de lı́mite de una función.
Definición 1.22. Sean f una función definida de un intervalo (a, b) a R, y x0 ∈ (a, b).
Se dice que w ∈ R es un lı́mite de f cuando x tiende a x0 , si para todo ε > 0, existe
δ = δ(ε) tal que |f (x) − w| < ε siempre que 0 < |x − x0 | < δ. Para indicar que w es un
lı́mite de f cuando x tiende a x0 , se escribe
lı́m f (x) = w.
x→x0
Definición 1.23. Se dice que una función f es continua en x0 ∈ R si:

i) f está definida en x0 , y
ii) lı́m f (x) = f (x0 ).
x→x0
A continuación damos los teoremas fundamentales sobre lı́mites y algunas propiedades

de las funciones continuas.
Teorema 1.7. Sean f y g dos funciones reales tales que lı́m f (x) = A y lı́m g(x) = B.
x→x0 x→x0
Entonces:
i) lı́m [f (x) ± g(x)] = A ± B,
x→x0
ii) lı́m [f (x) · g(x)] = A · B,
x→x0
iii) lı́m [f (x)/g(x)] = A/B, si B 6= 0.
x→x0
Teorema 1.8. Sean f y g dos funciones reales continuas en el punto x0 . La suma f + g, la

diferencia f − g, y el producto f ·g son también continuas en x0 . Si g(x0 ) 6= 0, entonces f /g
también es continua en x0 . Además, si f es continua en x0 y g es continua en y0 = f (x0 ),
entonces la función compuesta h(x) = g(f (x)) es continua en x0 .
Teorema 1.9 (Teorema de Bolzano). Sea f continua en cada punto del intervalo cerrado
[a, b] y supongamos que f (a)f (b) < 0. Entonces, existe por lo menos un punto c ∈ (a, b) tal
que f (c) = 0.
Teorema 1.10 (Teorema del valor intermedio). Sea f continua en cada punto de un
intervalo [a, b]. Si x1 < x2 son dos puntos cualesquiera de [a, b] tales que f (x1 ) 6= f (x2 ),
entonces la función f toma todos los valores comprendidos entre f (x1 ) y f (x2 ) por lo menos
una vez en el intervalo (x1 , x2 ).
A continuación describimos los conceptos básicos del cálculo diferencial.
Definición 1.24. Sea f una función definida por lo menos en un intervalo abierto (a, b).
La derivada de f en x ∈ (a, b), denotada por f ′ (x), se define como
f (x + h) − f (x)
f ′ (x) = lı́m ,
h→0 h
con tal que el lı́mite exista.
Como el Teorema 1.7, que nos enseña a calcular el lı́mite de la suma, diferencia, pro-
ducto y cociente de dos funciones, el siguiente teorema nos proporciona un conjunto de
reglas para el cálculo de derivadas.
Teorema 1.11. Sean f y g dos funciones definidas en un conjunto común de números

reales S. Supongamos que f y g tienen derivadas en cada punto de S. Entonces, la suma
f + g, la diferencia f − g, el producto f ġ y el cociente f /g, tienen derivadas que se calculan
por las fórmulas:
i) (f ± g)′ = f ′ ± g ′ ,
ii) (f · g)′ = f ′ · g + g ′ · f ,
′
f f ′ · g − g′ · f
iii) = , en puntos x tales que g(x) 6= 0.
g g2
Además, si f = u ◦ v tal que las derivadas v ′ (x) y u′ (y) existen con y = v(x), entonces la
derivada f ′ (x) también existe y está dada por la fórmula
f ′ (x) = u′ (v(x)) · v ′ (x),
denominada regla de la cadena.
La derivada tiene muchos usos en el Cálculo, especialmente se puede utilizar para la

localización de los máximos y mı́nimos de las funciones.
Definición 1.25. Se dice que una función real f tiene un máximo absoluto en un con-
junto S ⊂ R, si existe por lo menos un punto c ∈ S tal que f (x) ≤ f (c), para todo
x ∈ S. Ahora, se dice que f tiene un máximo relativo en un punto c ∈ S, si existe un
cierto intervalo abierto I ⊂ R que contiene a c tal que f (x) ≤ f (c), para todo x ∈ I ∩ S.
Los conceptos de mı́nimo absoluto y mı́nimo relativo se definen del mismo modo con la
desigualdad invertida.
Teorema 1.12. Sea f una función definida en un intervalo abierto I, y supongamos que
f tiene un máximo relativo o un mı́nimo relativo en un punto c ∈ I. Si la derivada f ′ (c)
existe, entonces f ′ (c) = 0.
A continuación presentamos otros teoremas importantes del cálculo diferencial.
Teorema 1.13 (Teorema del Valor Medio para Derivadas). Si f es una función continua
en todo punto del intervalo [a, b] que tiene derivada en cada punto del intervalo (a, b),
entonces existe por lo menos un punto c ∈ (a, b) para el que
f (a) − f (b) = f (c)(b − a).
Teorema 1.14. Supongamos f es continua en un intervalo [a, b] y que existe la derivada

f ′ en todo punto del intervalo (a, b), excepto en un punto c ∈ (a, b).
i) Si f ′ (x) > 0 para todo x < c y f ′ (x) < 0 para todo x > c, entonces f tiene un máximo
relativo en c.
ii) Si f ′ (x) < 0 para todo x < c y f ′ (x) > 0 para todo x > c, entonces f tiene un mı́nimo
relativo en c.
Teorema 1.15. Sea c un punto crı́tico de f en un intervalo (a, b); esto es, supongamos que
a < c < b y f ′ (c) = 0. Supongamos también que la segunda derivada f ′′ existe en (a, b).
Entonces:
i) Si f ′′ (x) < 0 para todo x ∈ (a, b), entonces f tiene un máximo relativo en c.
ii) Si f ′′ (x) > 0 para todo x ∈ (a, b), entonces f tiene un mı́nimo relativo en c.
Seguidamente describimos algunos conceptos básicos del cálculo integral.
Definición 1.26. Una partición P de un intervalo [a, b] es un conjunto finito de puntos

x0 , x1 , . . . , xn tales que
a = x0 ≤ x1 ≤ · · · ≤ xn = b.
Escribimos
∆xi = xi − xi−1 , para i = 1, . . . , n.
Definición 1.27. Sea f una función acotada en el intervalo [a, b], esto es, existen
números m y M tales m ≤ f (x) ≤ M para x ∈ [a, b]. Para cada partición P de [a, b]
hacemos
n
X n
X
Mi = sup f (x), mi = inf f (x), U (P, f ) = Mi ∆xi , L(P, f ) = mi ∆xi .
xi−1 ≤x≤xi xi−1 ≤x≤xi
i i
Las integrales superior e inferior de Riemann se definen respectivamente por

Z b Z b
f (x) dx = inf U (P, f ) y f (x) dx = sup L(P, f ).
a a
Si las integrales superior e inferior son iguales, se dice que f es integrable en [a, b] y el
valor común de estas integrales se representa con el sı́mbolo
Z b
f (x) dx,
a
que se denomina integral de Riemann de f o simplemente integral de f en [a, b].
Los siguientes teoremas relacionan la diferenciación con la integración.
Teorema 1.16. Sea f integrable en el [a, b]. Para x ∈ [a, b] hagamos

Z x
F (x) = f (t) dt.
a
Entonces, F es continua en [a, b]; además, si f es continua en el algún punto z ∈ [a, b], la
derivada F ′ (z) existe y
F ′ (z) = f (z).
La función F del teorema anterior usualmente se denomina primitiva de f . Seguida-

mente damos el Teorema Fundamental del Cálculo.
Teorema 1.17 (Teorema Fundamental del Cálculo). Si f es integrable en el [a, b] y si

existe una función F derivable en [a, b] tal que F ′ = f , entonces
Z b
f (x) dx = F (b) − F (a).
a
El Teorema de Taylor juega un papel importante en la definición de métodos para la

aproximación de funciones, la resolución de ecuaciones diferenciales, entre otros.
Teorema 1.18 (Teorema de Taylor). Sean f una función que tiene derivada continua de
orden n + 1 en el intervalo [a, b] y x0 ∈ [a, b]. Entonces, para todo x ∈ [a, b] existe β = β(x)
tal que
Xn
f (k) (x0 )
f (x) = (x − x0 )k + En (x),
k!
k=0
f (n+1) (β)
donde En (x) = (n+1)! (x − x0 )n+1 .
1.2.2 Funciones de Varias Variables
Los conceptos de campo escalar, campo vectorial y transformación lineal, son muy
importantes para el estudio de métodos para la aproximación de funciones, sistemas de
ecuaciones lineales y no lineales y optimización.
Definición 1.28. Sea D un subconjunto de Rn . Una función f de D en R se llama

un campo escalar o una función real de n variables. La función f asigna a cada vector
x ∈ D ⊂ Rn un valor real f (x).
Definición 1.29. Una función f : D ⊂ Rn → Rm se llama función vectorial de varias

variables. Si n = m > 1, la función se denomina campo vectorial. La función f se puede
expresar como
f (x) = (f1 (x), . . . , fm (x)) para x ∈ D,
donde fi : D ⊂ Rn → R, para i = 1, . . . , m, se llaman funciones componentes de f .
Definición 1.30. Se dice que una función vectorial T : Rn → Rm es una transformación

lineal si
T (x1 + x2 ) = T (x1 ) + T (x2 ), T (αx) = αT (x),
para todo x, x1 , x2 ∈ Rn y todo escalar α.
El concepto de continuidad de una función de varias variables se define como una

generalización de la continuidad de una función de una variable.
Definición 1.31. Sean f : D ⊂ Rn → R y x0 ∈ D. Se dice que f es continua en x0 si

para ε > 0 existe δ = δ(ε) > 0 tal que |f (x) − f (x0 )| < ε, siempre que 0 < kx − x0 k < δ
con x ∈ D.
Definición 1.32. Una función vectorial f : D ⊂ Rn → Rm es continua en x0 ∈ D, si

cada una de sus funciones componentes es continua en x0 .
Definición 1.33. Sean f : D ⊂ Rn → Rm y x0 ∈ D. Se dice que f es diferenciable en x0

si existe una transformación lineal T : Rn → Rm tal que
kf (x0 + h) − f (x0 ) − T (h)k

lı́m = 0.
h→0 khk
La generalización del concepto de derivaba de una función vectorial dada arriba, no es

muy operativa. Seguidamente damos algunos conceptos que permiten calcular la derivada
de funciones vectoriales.
Definición 1.34. Sean f : D ⊂ Rn → R y x ∈ D. La derivada parcial del campo escalar

f con respecto a la j-ésima variable xj de x se define por
∂f f (x + ej ) − f (x)
(x) = lı́m ,
∂xj t→0 t
en caso de que el lı́mite exista.
Las derivadas parciales de una función vectorial f : Rn → Rm , se definen de manera

análoga. La derivada parcial de f con respecto a la j-ésima variable xj de x existe si y sólo
si existen las derivadas parciales de sus funciones componentes f1 , . . . , fm y, en este caso,

∂f ∂f1 ∂fm
(x) = (x), . . . , (x) .
∂xj ∂xj ∂xj
Esto último nos permite definir la matiz jacobiana o jacobiano de un campo vectorial.
Definición 1.35. Sean f : D ⊂ Rn → Rm y x ∈ D. Sean f1 , . . . , fm las funciones

componentes de f . El jacobiano del campo vectorial f en x está dado por
 
∂f1 ∂f1 ∂f1
(x) (x) · · · (x)
 ∂x1 ∂x2 ∂xn 
 
 
 ∂f2 ∂f2 ∂f2 
 (x) (x) · · · (x) 
 
f ′ (x) =  ∂x1 ∂x2 ∂xn 
 .. .. .. 
 .. 
 . . . . 
 
 ∂f ∂f ∂f 
m m m
(x) (x) · · · (x)
∂x1 ∂x2 ∂xn
en caso de que las derivadas parciales existan.
En el caso de un campo escalar el jacobiano tiene una sola fila y se le conoce con el
nombre de gradiente.
Definición 1.36. Sean f : D ⊂ Rn → R un campo escalar y x ∈ D. El gradiente de f en

x se define por
∂f ∂f ∂f
∇f (x) = (x), (x), . . . , (x)
∂x1 ∂x2 ∂xn
Como el gradiente es un campo vectorial, entonces su jacobiano es una matriz cuadrada

denominada Hessiana.
Definición 1.37. Sean f : D ⊂ Rn → R un campo escalar y x ∈ D. La Hessiana de f en

x se define por
 2 
∂ f1 ∂ 2 f1 ∂ 2 f1
 ∂x2 (x) ∂x1 ∂x2
(x) · · ·
∂xn ∂x1
(x)
 1 
 
 ∂f ∂ 2 f ∂ 2 f 
 2
(x)
2
(x) · · ·
2
(x) 
 2 
2 
∇ f (x) =  ∂x 1 ∂x 2 ∂x 2 ∂x n ∂x 2 

 .. .. .. .. 
 . . . . 
 
 
 ∂ 2 fm ∂ 2 fm ∂ 2 fm 
(x) (x) · · · 2
(x)
∂x1 ∂xn ∂x2 ∂xn ∂xn
A continuación damos el Teorema de Taylor para funciones de dos variables.
Teorema 1.19 (Teorema de Taylor para funciones de dos variables). Sean f(t, y) y sus
derivadas parciales hasta orden (n + 1) continuas en el rectángulo abierto
R = {(t, y) ∈ R2 : a < t < b, c < y < d},
y sea (t0 , y0 ) ∈ R. Entonces,

∂f ∂f
f (t, y) = f (t0 , y0 ) + (t − t0 ) (t0 , y0 ) + (y − y0 ) (t0 , y0 ) + · · ·
∂t ∂y
" n #
1 X n n−i i ∂nf
+ (t − t0 ) (y − y0 ) n−i i (t0 , y0 ) + En (t, y),
n! i ∂t ∂y
i=0
donde En (t, y) es el término de error que envuelve las derivadas parciales de f (t, y) de
orden (n + 1).
24 1.3. CALCULANDO CON NÚMEROS
1.3 Calculando con Números

Todo algoritmo numérico, no importa cuán complicado sea, siempre se reduce a rea-
lizar una secuencia de operaciones aritméticas simples. Prácticamente todos los cálculos
numéricos se realizan en algún tipo de computadora digital; por ello es necesario describir
los detalles de la aritmética en computadora. Primero, revisemos cómo se representan los
números en una computadora.
1.3.1 Representación Numérica

La representación o sistema decimal que usamos a diario se nos hace muy familiar y
trabajamos con éste fácilmente. La representación decimal es un sistema posicional que
usa la base 10. Esto es, el valor de un número depende de sus dı́gitos y de la posición
de los mismos con respecto a un punto fijo, el punto decimal. El valor de cada número
se multiplica por una potencia de 10, creciendo las potencias positivas a la izquierda del
punto decimal y decreciendo las potencias negativas a la derecha. Por ejemplo, el número
512.31 es representado como un número decimal por:
(512.31)10 = 5 × 102 + 1 × 101 + 2 × 100 + 3 × 10−1 + 1 × 10−2
El sistema decimal no es eficiente para ser usado internamente por una computadora.
En lugar de ello se utiliza, generalmente, el sistema de base 2 o binario. Un sistema
binario es semejante al sistema decimal, excepto que la base es 2, los dı́gitos binarios
(bits) permitidos son 0 y 1, y cada bit es multiplicado por una potencia de 2, dependiendo
de su posición relativa al punto binario. La conversión de binario a decimal y viceversa es
fácil:
(1001.001)2 = 1 × 23 + 0 × 22 + 0 × 21 + 0 × 2−1 + 0 × 2−2 + 1 × 2−3

1
= 8 + 1 + = (9.125)10
8
Una dificultad con esta representación, a menudo llamada de punto-fijo, es su inefi-
ciencia para representar números muy grandes o muy pequeños; por ejemplo, el número
0.00000015 tiene nueve dı́gitos, pero solamente dos de ellos ofrecen información, el resto
son todos ceros que son necesarios para escalar el número. Tal inconveniente puede
remediarse con la notación cientı́fica, cuyo objetivo principal es representar un número
en dos partes. La primera es un número decimal con representación de punto-fijo y,
la segunda parte, es un factor de escala que permite representar un rango de números.
Especı́ficamente, un número real x se representa en notación cientı́fica como
x = ± 0.d1 d2 · · · dn dn+1 · · · × β e
donde di ∈ {0, 1, . . . , β − 1}, para i = 1, 2, . . . , e es un entero positivo y β es la base

(comúnmente se usa β = 2 o β = 10). Los dı́gitos di empezando con el primer dı́gito
distinto de cero se denominan dı́gitos significativos. Si d1 es distinto de cero, se dice x está
normalizado. Por ejemplo, el número 0.00000015 puede ser representado como 0.15 × 10−6
o 0.015 × 10−5 (existen infinitas representaciones, se deja al lector encontrar otra). La pri-
mera representación es la notación cientı́fica normalizada de 0.00000015; en cambio la
segunda no lo es, además, se aprecia que 0.00000015 tiene dos dı́gitos significativos. A par-
tir de ahora supondremos, a menos que se indique lo contrario, que los números estarán
representados en notación cientı́fica normalizada con lo cual garantizamos la unicidad en

la representación. Una versión binaria de esta notación es usada comúnmente en compu-
tadoras y se describirá en detalle en la siguiente sección.
1.3.2 Representación de Punto-Flotante

La representación de punto-flotante de un número tiene dos partes: una mantisa m,
que representa los dı́gitos significativos y un exponente e, que da un factor de escala. Am-
bos m y e están codificados en binario con un número fijo de bits, y el factor de escala es
alguna potencia de 2. Dentro de estas especificaciones existen muchas posibilidades de es-
cogencia. En 1985 el Instituto de Ingenierı́a Eléctrica y Electrónica (conocido por sus siglas
IEEE, en inglés Institute of Electrical and Electronics Engineers) produjo un conjunto de
especificaciones y requerimientos para la aritmética de punto-flotante que son conocidos
como estándar IEEE. Actualmente, el esquema más usado para punto-flotante es el formato
doble IEEE o de doble precisión, que se muestra en la Figura 1.2.
s
|{z}|
e
{z }|
m
{z }
1 bit 11 bits 52 bits
Figura 1.2. Formato doble IEEE
El formato doble IEEE utiliza 64 bits (u 8 bytes); un bit es usado para el signo del
número s (1 cuando el número es negativo y 0 cuando es positivo), y 52 bits son dedi-
cados a la mantisa m. Los restantes 11 bits representan el exponente e. Un número x
representado por este arreglo se escribe como:
x = (−1)s × 1.(m)10 × 2E
El campo del exponente e es una representación sesgada de su valor actual E (entero en
base decimal); además, e y E están relacionados por
(e)10 = E + 1023.
La mantisa m es una fracción con el punto binario inmediatamente a su izquierda, por
ejemplo,
m = 1101000000000000000000000000000000000000000000000000
tiene un valor decimal
(m)10 = 2−1 + 2−2 + 2−4 = 0.8125
Ahora, e y s son considerados enteros. Pasemos a ver un ejemplo.
Ejemplo 1.10. Hallemos la representación del número −11.6 en el formato doble IEEE. Se
comienza expresando a −11.6 punto-flotante
−11.6 = (−1)1 × (1.4375) × 23
De allı́ se obtiene que s = 1. Ahora, para encontrar la representación binaria de la
parte fraccionaria, .4375, se procede de la siguiente manera. Debemos hallar los bits
a−1 , a−2 , a−3 , . . . tales que
.4375 = a−1 2−1 + a−2 2−2 + a−3 2−3 + · · ·
Para hallar a−1 , multiplicamos por 2 la ecuación anterior en ambos lados. Ası́ obtenemos:
0.8750 = a−1 + a−2 2−1 + a−2 2−2 + · · ·
luego, a−1 = 0. Procediendo de la misma forma hallamos los siguientes bits. En su
totalidad, el procedimiento es el siguiente:
.4375
2
0 .8750
2
1 .7500
2
1 .5000
2
1 .0000
Por lo tanto, la mantisa m es la cadena de bits
m = 0111000000000000000000000000000000000000000000000000
El valor de E es 3; por lo tanto,
(e)10 = 1026 = 210 + 2.
Ası́,
e = 10000000010
Colocando juntos s, m y e, obtenemos que el número decimal −11.6 está representado en
el formato doble IEEE por la cadena de bits:
1100000000100111000000000000000000000000000000000000000000000000
Note que el formato doble IEEE tiene una longitud finita de bits (o cualquier otro
formato de punto-flotante). Por ello, este formato puede representar un conjunto finito de
números decimales. El rango de números que pueden ser representados varı́a de 2−1023 ≈
10−308 a 21024 ≈ 10309 .
Existen números reales que no pueden ser representados de forma exacta en punto-
flotante, es decir, se incurre en un error al representar el número real como un número
punto-flotante. Dos medidas muy utilizadas del error cometido al aproximar un número
son el error absoluto y el error relativo.
Definición 1.38. Sea x un número distinto de cero y sea x̃ una aproximación de x.

|x − x̃|
Entonces, el error absoluto y el error relativo se definen como |x − x̃| y , respec-
|x|
tivamente.
En el siguiente ejemplo mostramos el error cometido en la representación en punto-

flotante del número 0.1.
Ejemplo 1.11. El número

x = 0.1 = (−1)0 × (1.6) × 2−4
puede ser representado en el formato doble IEEE mediante las siguientes cadenas de bits
c1 = 011111110111001100110011001100110011001100110011001100110011001
o
c2 = 011111110111001100110011001100110011001100110011001100110011010
Para la cadena de bits c1 , se observa que
s = 0,
e = 1111111011,
m1 = 1001100110011001100110011001100110011001100110011001;
por lo tanto, al realizar la conversión a decimal se obtiene:
(e)10 = 1019 ⇒ E = −4,
13
X
(m1 )10 = (2−4k+3 + 2−4k ),
k=1
ası́,
0.1 ≈ x1 = (−1)0 × (1.(m1 )10 ) × 2−4
13
X
= 2−4 + (2−4k−1 + 2−4k−4 )
k=1
13
X
−4 −1 −4
=2 + (2 +2 ) 2−4k
k=1

+123 −1 2−52
= 2−4 +
24 1 − 24
≈ 0.09999999999999999167332731531133
Para la cadena de bits c2 , se observa que
s = 0,
e = 1111111011,
m2 = 1001100110011001100110011001100110011001100110011010;
por lo tanto, al realizar la conversión a decimal se obtiene:
(e)10 = 1019 ⇒ E = −4,
−52
(m2 )10 = (m1 )10 − 2 + 2−51 ,
ası́,
0.1 ≈ x2 = (−1)0 × (1.(m2 )10 ) × 2−4
3 −52
−4 2 +1 2 −1
=2 + − 2−56 + 2−55
24 1 − 24
≈ 0.10000000000000000555111512312579
De esta forma, si se toma a x1 como la representación de punto-flotante del número 0.1,

se comete un error absoluto
|0.1 − x1 | ≈ 8.326673 × 10−18 ,
y un error relativo
|0.1 − x1 |
≈ 8.326673 × 10−17 .
0.1
Ahora, si se toma a x2 como la representación de punto-flotante del número 0.1, se comete
un error absoluto
|0.1 − x2 | ≈ 5.551115 × 10−18 ,
y un error relativo
|0.1 − x2 |
≈ 5.551115 × 10−17 .
0.1
Los números en punto-flotante x1 y x2 del Ejemplo 1.11 se obtuvieron mediante el pro-

ceso de truncamiento, que consiste en descartar los bits excedentes de la representación
binaria de la mantisa m. Observe que x1 se encuentra a la izquierda de 0.1 en la recta real;
en cambio, x2 queda a la derecha. El número x2 también se obtuvo por truncamiento,
pero a diferencia de x1 , se obtuvo por el proceso de redondeo por exceso. Esto significa
que se elimina el exceso de bits de la mantisa m, pero en esta ocasión se aumenta en una
unidad el último bit. Ahora bien, para que no exista ambigüedad en la representación en
punto-flotante de un número x, se escoge el número punto-flotante más cercano a x, el
cual denotaremos como f l(x). En el caso del Ejemplo 1.11, la representación en punto
flotante de 0.1 es f l(0.1) = x2 .
Construyamos una cota del error relativo de un número x ∈ [2−1023 , 21024 ] cuya repre-
sentación en el formato doble IEEE es f l(x). Se tiene que todo número real x 6= 0 se puede
representar como
x = ± 0.d1 d2 · · · d52 d53 · · · × 2α ,
donde 2−1023 ≤ α ≤ 21024 . Es claro que se satisface:
2α−1 ≤ |x| ≤ 2α , (1.2)
además, |⌊x⌋ − x| ≤ 1, donde ⌊x⌋ es el entero más cercano a x tal que ⌊x⌋ ≤ x. Luego,
⌊x253−α ⌋ = ±d1 d2 · · · d53 ,
y de esto se obtiene
f l(x) = ⌊x253−α ⌋2α−53 .
De esta forma, podemos escribir:
|f l(x) − x| = ⌊x252−α ⌋2α−53 − x 2α−53 253−α

= ⌊x253−α ⌋ − x 253−α 2α−53
≤ 2α−53 (1.3)
que es una cota del error absoluto. Por (1.2) se tiene 1/|x| ≤ 21−α . Entonces, por (1.3)
obtenemos
x − f l(x)
≤ 2α−53 21−α = 2−52 ,
x
que es una cota del error relativo. Si tomamos µ = (f l(x) − x)/x, la desigualdad anterior
adquiere la siguiente forma:
f l(x) = x(1 + µ), |µ| ≤ ǫ = 2−52 , (1.4)
para todo número real x.
Definición 1.39. El número ǫ que aparece en la ecuación (1.4) se denomina error de

redondeo unitario o, equivalentemente, epsilon o precisión de la máquina, y se denota
como eps. El epsilón de la máquina está directamente relacionado con el número de bits
asignados a la mantisa. En el formato doble IEEE el valor de eps es 2−52 ≈ 2.2 × 10−16 .
El épsilon de la máquina eps también puede definirse como el número positivo más
pequeño tal que (1 + eps) > 1. Con esta definición diseñaremos un pequeño programa que
nos permitirá calcular el eps de la computadora que estemos utilizando.
Ejemplo 1.12. Para encontrar la medida del eps, le sumaremos el valor 1 a una sucesión
decreciente de números positivos, hasta que los resultados no muestren cambios. De esta
forma, encontraremos el número positivo eps tal que (1 + eps) > 1. Seguidamente se
muestra un simple programa que calcula el epsilon de la máquina.
eps = 1 .
w h i l e 1 . + eps > 1 . :
eps /= 2 .
eps ∗= 2 .
p r i n t ( eps )
Cuando este programa corre en nuestra computadora (en Python), el valor que imprime es
eps = 2.2204 × 10−16 ,
exactamente lo que uno espera del formato doble IEEE.
En algunas situaciones, cuando realizamos cálculos en una computadora, se puede ge-

nerar un número de la forma ±p × 2E , donde E queda fuera del rango permitido por la
computadora. Si E es demasiado grande, se dice que tuvo lugar un overflow o desbor-
damiento por exceso. Si E es demasiado pequeño, se dice que ocurrió un underflow o
desbordamiento por defecto. En el formato doble IEEE esto significa que E > 1024 o que
E < −1023, respectivamente. En general, que ocurra un overflow se considera un error fa-
tal, que en muchas computadoras ocasionarı́a la interrupción de manera automática de la
ejecución del programa. Cuando ocurre un underflow, muchas computadoras simplemente
asignan a la variable el valor 0 y el cálculo continúa.
Bajo las especificaciones del formato estándar IEEE, se establecen las siguientes leyes
de la aritmética punto-flotante [15, 24].
Teorema 1.20 (Leyes de la aritmética punto flotante). Sean x e y dos números punto-
flotante y sean f l(x + y), f l(x − y), f l(xy), y f l(x/y) la suma, resta, producto y división,
respectivamente. Entonces:
i) f l(x ± y) = (x ± y)(1 + µ), donde |µ| < eps.
ii) f l(xy) = (xy)(1 + µ), donde |µ| < eps.
iii) Si y 6= 0, entonces f l(x/y) = (x/y)(1 + µ), donde |µ| < eps.
iv) Si x e y son números reales que no poseen representación exacta punto-flotante, en-
tonces se cumple la siguiente ley:
f l(x+y) = (x(1+µ1 )+y(1+µ2 ))(1+µ3 ), donde |µ1 | < eps, |µ2 | < eps, |µ3 | < eps.
Esta última propiedad se cumple de forma similar para la resta, el producto y la di-
visión.
El Ejemplo 1.11 es una evidencia que los errores de redondeo son inevitables y difı́ciles
de controlar. Existe otro tipo de errores que sı́ se pueden de alguna forma controlar. Un
error frecuente en análisis numérico es la pérdida de dı́gitos significativos que aparece
como consecuencia de descuidos durante el proceso de programación. El caso más em-
blemático es la sustracción de cantidades casi iguales. Como ejemplo veamos la magnitud
del error relativo que puede producirse al restar dos números muy cercanos.
Ejemplo 1.13. Sean x = 0.37214786932175 e y = 0.37214720021582. Se tiene que
x − y = 0.00000066910593
Si los cálculos se realizaran en una computadora decimal con mantisa de 8 dı́gitos, verı́amos
que
f l(x) = 0.37214786, f l(y) = 0.37214720, ⇒ f l(x) − f l(y) = 0.00000066
Por lo que, el error relativo es:
x − y − (f l(x) − f l(y)) 0.00000066910593 − 0.00000066

= ≈ 0, 013609 ⇒ 1.4%
x−y 0.00000066910593
En otras palabras, el error relativo es grande. Ahora, el número f l(x) − f l(y) se representa
en la computadora como 0.66000000 × 10−6 , aunque los ceros de la mantisa sólo sirven
para denotar el lugar decimal. Note que se perdieron 6 dı́gitos significativos.
La situación descrita en el ejemplo anterior también puede aparecer cuando se evalúan

expresiones algebraicas. Veamos el siguiente ejemplo.
Ejemplo 1.14. Supongamos que deseamos calcular

p
y = 1 + x2 − 1
cuando x es un número pequeño. Observamos que si x ≈ 0, entonces y ≈ 0. Este hecho

implica la cancelación por sustracción y pérdida de dı́gitos significativos para valores de x
muy cercanos a 0. ¿Cómo podemos evitar este problema?
Si consideramos una forma equivalente de calcular y:
!
p √ 1 + x2 + 1 x2
y= 1 + x2 − 1 √ =√ ,
1 + x2 + 1 1 + x2 + 1
entonces la dificultad se evita.
El siguiente teorema establece cotas del número de bits significativos que se pierden al
efectuar la sustracción x − y cuando x es cercano a y, en términos de la cantidad |1 − y/x|.
Teorema 1.21. Si x e y son números punto-flotantes, positivos y tales que x > y con
y
2−q ≤ 1 − ≤ 2−p ,
x
entonces en la resta x − y se pierden a lo sumo q y al menos p bits significativos.
1.3.3 Eficiencia, Estabilidad y Condicionamiento

En el cálculo numérico, además de medir el error en aproximaciones, también es nece-
sario medir la eficiencia y la estabilidad de un algoritmo.
Definición 1.40. Un algoritmo es un conjunto finito y ordenado de operaciones, lógicas

y aritméticas, las cuales al aplicarse a un problema computacional asociado a un con-
junto de datos, llamado entrada, produce una solución al problema. Esta solución com-
prende un conjunto de datos que se conocen como salida.
El trabajo computacional de los algoritmos o, equivalentemente, su eficiencia se mide

por el número de operaciones punto-flotante que son necesarias en su ejecución.
Definición 1.41. Una operación punto-flotante es toda operación de suma, resta, mul-
tiplicación o división.
Definición 1.42. Se dice que un algoritmo es de orden O(nk ), si el término dominante

en el número de operaciones punto-flotante requeridas en su ejecución es un múltiplo
de nk .
Ejemplo 1.15. Supongamos que deseamos calcular la siguiente suma
N
X 1
S= (−1)i+1 .
1 + i2
i=1
Desarrollando la suma se tiene

1 1 1
S = (−1)2 + (−1)3 + · · · + (−1)N +1
1 + 12 1 + 22 1 + N2
1 1 1
= (−1) · (−1) + (−1)2 · (−1) + · · · + (−1)N · (−1) .
1+1·1 1+2·2 1+N ·N
Luego, el número de operaciones punto-flotante necesarias para calcular suma S es: 2N −1
sumas, 3N multiplicaciones y N divisiones; por lo tanto, para calcular la suma S hacen
falta 6N − 1 operaciones punto-flotante. Ası́, un algoritmo para calcular la suma S es de
orden O(N ).
Por otro lado, la precisión de la solución de un problema obtenida por un algoritmo

depende de dos caracterı́sticas importantes ([4]):
• La estabilidad del algoritmo.
• La condición del problema considerado, es decir, que tan sensible es el problema a
pequeñas perturbaciones en el planteamiento del mismo.
Informalmente, diremos que un algoritmo es inestable si pequeños errores cometidos

en algún paso del algoritmo, se magnifican en los siguientes pasos y destruyen la precisión
de la solución obtenida. Una definición formal de estabilidad es la siguiente.
Definición 1.43. Un algoritmo es estable si produce la solución exacta de un problema

cuyos datos de entrada son cercanos a los datos originales.
Ejemplo 1.16. Un ejemplo simple de estabilidad son las operaciones aritméticas de punto-
flotante. Veamos que la suma de dos números punto-flotante es estable. Sean x e y dos
números punto-flotante. Por el Teorema 1.20 se tiene:
f l(x + y) = (x + y)(1 + µ),
donde |µ| < eps. De esta forma, podemos escribir
f l(x + y) − (x + y) = (x + y)(1 + µ) − (x + y) = µ(x + y).
Por tanto, como |µ| < eps, la suma de x e y obtenida por f l(x+y) es cercana a la suma x+y
de los datos originales x e y. En consecuencia, la operación de sumar dos números punto-
flotante es estable. Dde una manera similar se muestra la estabilidad de la multiplicación
y la división.
Ejemplo 1.17. Un ejemplo de un algoritmo numérico inestable es el proceso de cálculo de

la factorización LU de una matriz cuadrada, el cual lo estudiaremos con mayor detalle en
el Capı́tulo ??. Supongamos que se desea encontrar la factorización LU de la matriz

δ −1 1 0 u11 u12
A= = , 0 < δ ≪ 1,
1 1 l21 1 0 u22
es decir, se debe hallar l12 , u11 , u12 y u22 . Es claro que u11 = δ, u12 = −1, l21 = δ−1 y
u22 = 1 − l21 u12 = 1 + δ−1 . En aritmética punto-flotante, si δ es suficientemente pequeño
(δ < eps), entonces f l(u22 ) = δ−1 . Aun suponiendo que el resto de las cuentas se obtienen
de forma exacta, al multiplicar las dos matrices triangulares punto flotante, se obtiene la
matriz
δ −1
1 0
cuyo elemento (2, 2) difiere notablemente del elemento (2, 2) de la matriz A. Es decir, el
proceso de factorización realizado de esa manera reproduce una matriz que dista mucho
de la original, y por ende es un proceso inestable. Comprobemos esto numéricamente.
Inicialmente tomemos δ = 10−15 y calculemos el producto LU .
>>> d e l t a = 1 . 0 e−15
>>> L = np . a r r a y ( [ [ 1 , 0 ] , [ 1 / d e l t a , 1 ] ] )
>>> U = np . a r r a y ( [ [ d e l t a , −1] ,[0 ,1+(1/ d e l t a ) ] ] )
>>> L . dot (U)
a r r a y ( [ [ 1 . e −15, −1.e +00] ,
[ 1 . e+00, 1 . e +00]])
En este caso, el producto de las matrices L y U arroja la matriz A original. Ahora, tomemos
δ = 10−16 y calculemos nuevamente el producto LU .
>>> d e l t a = 1 . 0 e−16
>>> L = np . a r r a y ( [ [ 1 , 0 ] , [ 1 / d e l t a , 1 ] ] )
>>> U = np . a r r a y ( [ [ d e l t a , −1] ,[0 ,1+(1/ d e l t a ) ] ] )
>>> L . dot (U)
a r r a y ( [ [ 1 . e −16, −1.e +00] ,
[ 1 . e+00, 0 . e +00]])
La matriz obtenida difiere en el elemento (2, 2) de la matriz A original. Esta inestabilidad se

resolverá en el Capı́tulo ?? introduciendo diferentes estrategias de pivoteo para factorizar
matrices.
Hemos visto que la precisión de la solución de un problema depende de la estabilidad

de los algoritmos, pero exigir que el algoritmo sea estable no es suficiente para garantizar
la buena calidad de las soluciones obtenidas. Como indicamos arriba, el condicionamiento
también contribuye en la precisión de las soluciones esperadas. La condición de un pro-
blema es una propiedad intrı́nseca del problema y se refiere a como cambia la solución
cuando los datos contienen algún tipo de ruido o perturbación, independientemente del
algoritmo usado. El estudio teórico efectuado por analistas numéricos para investigar el
efecto que estas perturbaciones producen se conoce como análisis de perturbación. Este
tipo de análisis permite discernir si un problema dado es “bueno” o “malo”, es decir, si pe-
queñas perturbaciones en los datos producen pequeños o grandes cambios en la solución
del problema ([4]).
Definición 1.44. Un problema (con respecto a un conjunto de datos) se llama mal-

condicionado si un error relativo pequeño en los datos produce un error relativo grande
en la solución, independientemente del algoritmo usado. En caso contrario se llama
bien-condicionado.
Una interpretación matemática de esta definición es la siguiente. Supongamos que un

problema P se resuelve para un conjunto de datos d y denotamos con P (d) el valor de
la solución obtenida. Sea δd una perturbación de los datos d. Diremos que P está mal-
condicionado si el error relativo de la solución es mucho mayor que el error relativo de los
datos, en otras palabras, si se cumple la desigualdad
|P (d + δd ) − P (d)|
|P (d)|
κP = ≫ 1.
|δd |
|d|
El escalar κP > 0 se denomina número de condición del problema. Si κP ≪ 1, se dice que el
problema está bien-condicionado. Vale destacar que, en general, un problema puede estar
mal-condicionado para un conjunto de datos y bien-condicionado para otro conjunto de
datos. Veamos un ejemplo.
Ejemplo 1.18. Supongamos que se desea calcular un valor de la función
p
f (x) = x2 + 1 + ex ,
para x ∈ R. Calculemos el número de condición de este problema. Sea δ ∈ (0, 1) lo
suficientemente pequeño. Se tiene que el número de condición del problema está dado
por:
|f (x + δ) − f (x)|
|f (x)| f (x + δ) − f (x) |x| f ′ (x)
κP = = · ≃ x ,
|δ| δ |f (x)| f (x)
|x|
para todo x ∈ R. Calculemos el número de condición del problema para ciertos valores de
x. Se tiene que
x
f ′ (x) = √ + ex .
2
x +1
Para x = 100, √
f ′ (100) 100 e100 + 10000 10001
κP = 100 · = √ 10001 ≈ 100.
f (100) e100 + 10001
Para x = 10, √
f ′ (10) 10 e10 + 100101101
κP = 10 · = √ ≈ 10.
f (10) e10 + 101
Para x = 1, √
f ′ (1) e + 22
κP = 1 · = √ ≈ 0.82889109017393.
f (1) e+ 2
Para x = 10−4 ,
′ −4 1 √
−4 f (10 ) 100000001 e 10000 + 100000001
≈ 5.0 × 10−5 .
κP = 10 · = √
f (10−4 ) 1
100000001 10000 e 10000 + 100000001
De esta forma, calcular un valor de la función f (x) dada es un problema que está mal
o bien-condicionado dependiendo del valor de x considerado. Por ejemplo, si x = 100
el problema está mal-condicionado; por el contrario, si x = 10−5 el problema está bien-
condicionado.
Condicionamiento de Sistemas Lineales

En esta sección estudiaremos el efecto de pequeñas perturbaciones en sistemas lineales
de la forma Ax = b, donde A ∈ Rn×n y b ∈ Rn . Estos problemas aparecen en muchas
aplicaciones. En el Capı́tulo ?? los estudiaremos en detalle. En un sistema lineal, los datos
de entrada son una matriz cuadrada A de orden n y el vector b ∈ Rn , que pueden tener
perturbaciones. Veremos que un número especial, conocido como el número de condición
de la matriz A, jugará un papel fundamental en el condicionamiento de un sistema lineal.
Definición 1.45. El número de condición de la matriz A ∈ Rn×n no singular, denotado

por κp (A), se define por
κp (A) = kAkp kA−1 kp ,
donde k · kp es un norma matricial inducida.
Denotemos las perturbaciones de A, b y de la solución x de la siguiente manera:
A → A + △A (△A: perturbación de la matriz A),

b → b + δb (δb: perturbación del vector b),
x → x + δx (δx: cambio en la solución).
El siguiente teorema establece que el error relativo en la solución x está dominado por
el número de condición de la matriz A, siempre que las perturbaciones relativas en A y b
sean pequeñas. Para la demostración de este resultado general ver, por ejemplo, el libro
[9, p. 123].
Teorema 1.22. Sean △A y δb las perturbaciones de A y b 6= 0, respectivamente, y sea δx

el error absoluto en x al resolver Ax = b. Sea k · kp una norma matricial inducida. Si A
1
es no singular y k△Akp < −1
, entonces
kA kp

kδxkp κp (A) k△Akp kδbkp
≤ + .
kxkp (1 − k△Akp kA−1 kp ) kAkp kbkp
Del Teorema 1.22 se infiere:
• Si κp (A) es grande, entonces el sistema Ax = b está mal-condicionado.
• Si κp (A) es pequeño, entonces el sistema Ax = b está bien-condicionado.
El número de condición de una matriz depende de la norma utilizada, sin embargo,

recordando que en espacios de dimensión finita (como Rn o Rn×n ) todas las normas son
equivalentes (ver Ejercicio 4), podemos afirmar que si una matriz es mal-condicionada en
una norma, también lo es en cualquier otra norma. Algunas de las propiedades de κp (A)
que podemos destacar son las siguientes:
(i) κp (A) = κp (A−1 ).
(ii) κp (A) ≥ 1 para cualquier norma inducida.

(iii) κp (cA) = κp (A), para cualquier escalar c 6= 0.

(iv) κp (AT ) = κp (A).
(v) κp (AB) ≤ κp (A)κp (B), para cualquier norma inducida.
(vi) κp (AT A) = (κp (A))2 .
El comando cond del módulo linalg de la librerı́a Numpy, se utiliza para calcular el
número de condición κp (A), para valores de p = 1, 2, inf, ‘fro’, entre otros. Generalmente,
se emplea el número de condición κ2 (A) calculado con la norma-2 o norma Euclı́dea.
Veamos un ejemplo.
Ejemplo 1.19. Para la matriz A,
>>> A = np . a r r a y ( [ [ 1 , 4 , 6 ] , [ 3 , − 1 , 5 ] , [ 0 , 8 , 7 ] ] )
>>> A
array ([[ 1 , 4 , 6] ,
[ 3 , −1, 5 ] ,
[ 0 , 8 , 7]])
se tiene que el número de condición κp (A) se calcula como:

>>> from numpy import l i n a l g as LA
>>> LA . cond (A) # Número de c o n d i c i ó n con p=2
69.82769061243755
>>> LA . cond (A , 1) # Número de c o n d i c i ó n con p=1
127.38461538461519
>>> LA . cond (A , ’ f r o ’ ) # Número de c o n d i c i ó n con l a norma de F r o b e n i u s
75.18627754690965
>>> LA . cond (A , np . i n f ) # Número de c o n d i c i ó n con l a norma i n f i n i t o
107.30769230769216
Es claro que el número de condición κp (A) está definido para matrices no singulares.
Se puede establecer que κp (A) es el inverso de la distancia a la singularidad, es decir, si
κp (A) es grande entonces A es muy cercana a una matriz singular ([4]).
En el siguiente teorema se establece la versión simplificada del Teorema 1.22 para el
caso en que la matriz no sufre perturbaciones, esto es, cuando △A = 0. Este es un caso de
interés, ya que en muchas aplicaciones la matriz A representa el modelo matemático fijo y
el ruido o perturbación aparece en el vector b de mediciones experimentales.
Teorema 1.23. Sean δb y δx, respectivamente, la perturbación de b y x al resolver Ax =

b. Sea k · kp una norma matricial inducida. Si A es no singular y b 6= 0, entonces
kδxkp kδbkp
≤ κp (A) .
kxkp kbkp
El Teorema 1.23 dice que el error relativo en la solución puede llegar a ser tan grande
como la κp (A) multiplicada por el error relativo en el vector b. Por tanto, si κp (A) no es
muy grande, entonces una pequeña perturbación en b tendrá poco efecto en la solución
obtenida. Mientras que si κp (A) es grande, entonces aún una pequeña perturbación en b
puede ocasionar un cambio muy grande en la solución.
Ejemplo 1.20. Consideremos el sistema lineal Ax = b, donde

2 1 1
A= , b= .
2 1.0001 0
La solución del sistema es x = (5000.5, −10000). El efecto del mal condicionamiento

puede verificarse cambiando el elemento a22 de la matriz A; por ejemplo, si a22 = 1.0002,
entonces la solución del sistema de ecuaciones es x = (2500.5, −5000). Observe que un
cambio del 0.01% en el coeficiente a22 produce un cambio del 100% en la solución. Por
otro lado, si variamos un poco el vector del lado derecho, digamos b = (0.999, 0), entonces
la solución del sistema es x = (4995.5, −9990). Estos hechos nos permiten asegurar que la
matriz A está mal condicionada, cuyo número de condición es κ2 (A) = 50001.
Ejemplo 1.21. Consideremos un sistema lineal Ax = b, donde

   
2 1 1 4
A = 1 0.0003 0.0002 y b = 1.0005 .
1 0.0002 0.0003 1.0005
 
1
La solución exacta del sistema lineal es x = 1. Si cambiamos el vector b por
1
 
4
b = 1.0001 ,
b
1.0007
entonces la perturbación relativa de b es
b 2
kb − bk kδbk2
= = 1.054 × 10−4 (pequeña).
kbk2 kbk2
Si ahora resolvemos el sistema Ab b obtenemos la siguiente solución

x = b,
 
0.9999
b = x + δx = −1.9999 ,
x
4.0001
la cual es completamente diferente al vector x. Es de notar que en este caso el número de

condición de la matriz A es κ2 (A) = 2.7322 × 104 (grande).
1.3.4 Velocidad de Convergencia

En este libro trataremos con frecuencia la sucesión de vectores y su convergencia. Para
este propósito, recordemos que una sucesión de vectores {xk } en Rn , converge a un vector
x ∈ Rn , y escribimos lı́m xk = x, si
k→∞
lı́m xki = xi ,
k→∞
donde xki y xi son las componentes de los vectores correspondientes. Por otra parte, decir
que la sucesión {xk } converge x es equivalente a exigir que la sucesión {xk − x} converja a
0. Este hecho unido a que en Rn las normas son equivalentes, se establece el siguiente re-
sultado que permite definir convergencia de una sucesión de vectores utilizando cualquier
norma vectorial en Rn .
Teorema 1.24. Sea k · k una norma vectorial en Rn . Entonces,
lı́m xk = x ⇔ lı́m kxk − xk = 0,

k→∞ k→∞
donde x ∈ Rn y {xk } es una sucesión de vectores en Rn .
Un aspecto importante de un método iterativo es su velocidad de convergencia. A

continuación presentamos una definición de velocidad de convergencia que emplearemos
en distintos capı́tulos. Supongamos que la sucesión de vectores {xk } de Rn converge a un
vector x∗ ∈ Rn .
Definición 1.46. Si existe un número real p ≥ 1, un entero k0 > 0 y una constante

positiva C independiente de k0 > 0, tal que
kxk+1 − x∗ k ≤ Ckxk − x∗ kp para k ≥ K, (1.5)
se dice que la sucesión {xk } converge a x∗ con una Q-velocidad de convergencia de

orden p.
Se distinguen tres tipos de convergencias:

1. Cuando p = 1 y C ∈ (0, 1), se dice que la sucesión {xk } converge a x∗ Q-linealmente.
2. Cuando 1 < p < 2 y C > 0, se dice que la sucesión {xk } converge a x∗ Q-

superlinealmente; equivalentemente, una sucesión {xk } se dice que converge a x∗
Q-superlinealmente, si para alguna sucesión de números positivos {ck } que con-
verge a cero, se cumple
kxk+1 − x∗ k ≤ ck kxk − x∗ k.
3. Cuando p = 2, se dice que la sucesión {xk } tiene una convergencia Q-cuadrática.

Ejemplo 1.22. Demostrar que la sucesión { 21n }∞
n=0 converge Q-linealmente a 0, la sucesión
n n
{2−(3/2) }∞n=0 converge Q-superlinealmente a 0 y que la sucesión {3−2 }∞
n=0 converge Q-
cuadráticamente a 0.
n n
Solución. Es claro que las sucesiones { 21n }∞
n=0 , {2
−(3/2) }∞ y {3−2 }∞ convergen a 0.
n=0 n=0
Se tiene que
1 1 1
n+1
= · n , para todo n ≥ 0,
2 2 2
por tanto, la sucesión { 21n }∞
n=0 converge Q-linealmente a 0. Ahora,
n+1

n 3/2
2−(3/2) = 2−(3/2) , para todo n ≥ 0,
n
de lo cual se infiere que la sucesión {3−2 }∞n=0 converge Q-superlinealmente a 0. Por
último, se tiene que
n+1 n 2
3−2 = 3−2 , para todo n ≥ 0,
n
de lo cual se infiere que la sucesión {3−2 }∞
n=0 converge Q-cuadráticamente a 0. Numérica-
mente se puede mostrar la diferencia entre los órdenes de convergencia. En la siguiente
n
tabla se aprecian los siete primeros elementos de las sucesiones { 21n }∞
n=0 , {2−(3/2) }∞
n=0 y
n
{3−2 }∞n=0 . ¿Qué sucesión
se acerca más
{ 21n }∞
n n
n n=0 {2−(3/2) }∞
n=0 {3−2 }∞
n=0 rápido a 0 y
0 1.000000000000000000 0.500000000000000000 0.333333333333333315 en qué propor-
1 0.500000000000000000 0.353553390593273786 0.111111111111111105 ción?
2 0.250000000000000000 0.210224103813428626 0.012345679012345678
3 0.125000000000000000 0.096388176587996297 0.000152415790275873
4 0.062500000000000000 0.029925102521830428 0.000000023230573125
5 0.031250000000000000 0.005176705637342740 0.000000000000000540
6 0.015625000000000000 0.000372460486021616 0.000000000000000000
1.4 Ejercicios
1. Demuestre que Rn es un espacio vectorial de dimensión n.
2. Si I es la matriz identidad de orden n, demuestre que:
a) kIk2 = 1.
√
b) kIkF = n.
3. Si D = (dij ) ∈ Rn×n es una matriz diagonal, demuestre que:
a) kDk1 = kDk2 = kDk∞ = máx {|dii |}.

i=1,...,n
Pn
2 1/2 .
b) kDkF = i=1 dii
4. Las p-normas vectoriales y matriciales están relacionadas por varias desigualdades que
envuelven las dimensiones m y n. Estas desigualdades hacen que las normas sean
equivalentes. Verifique cada una de las siguientes desigualdades y de un ejemplo de un
vector o una matriz distintos de cero (para m, n generales) para el cual la desigualdad
se cumple. En este problema x ∈ Rm y A ∈ Rm×n .
a) kxk∞ ≤ kxk2 .
40 1.4. EJERCICIOS
√
b) kxk2 ≤ mkxk∞ .
√
c) kAk∞ ≤ nkAk2 .
√
d) kAk2 ≤ mkAk∞ .
e) kAk2 ≤ kAkF .
√
f) kAkF ≤ nkAk2 .
5. Encuentre en cada uno de los siguientes casos, una fórmula equivalente a la dada que
evite la pérdida de cifras significativas.
a) ln(x + 1) − ln(x), para x grande.

√
b) x2 + 1 − x, para x grande.
c) cos2 (x) − sen2 (x), para x ≈ π/4.
d) ex − x − 1, para valores negativos de x.
1 1
e) − , para x grande.
x x+1
f) x − sen(x), para valores de x cercanos a cero.
ex − 1
g) , para valores de x cercanos a cero.
x
6. Determine los errores absoluto y relativo, y el número de dı́gitos correctos, al aproxi-
mar:
a) 0.00001 con 0
1
b) con 0.333
3
c) 0.3000001 con 0.2999999
qP
n 2
7. Suponga que se desea calcular y = i=1 ai , para ai ∈ R. Determine una expresión
equivalente
Pn para calcular y de manera que se pueda evitar el overflow al calcular
2
i=1 ai , para valores grandes de los números ai .
8. El formato de simple precisión es la representación en punto-flotante que utiliza 32 bits:

un bit es usado para el signo del número, 23 bits son dedicados a la mantisa y los
restantes 8 bits representan el exponente. ¿Qué respuesta conseguirá usted si calcula
las siguientes expresiones en un computador con simple precisión?
√
a) 109 − 1
√
b) 10−21 − 1
c) 1016 − 5
9. En doble precisión, ¿cuál de los siguientes cálculos podrı́a resultar en una respuesta
distinta de cero?
a) (100 + 1/3) − 1/3 − 100.

b) (1/3 + 100) − 100 − 1/3.
c) (1/3 + 100) − (100 + 1/3).
Revise sus predicciones usando la computadora.

10. La Fórmula Mejorada para la Resolución de la Ecuación de Segundo Grado. Supongamos

que a 6= 0 y que b2 − 4ac > 0 y consideremos la ecuación ax2 + bx + c = 0. Sus raı́ces
pueden hallarse mediante la conocida fórmula
√ √
−b + b2 − 4ac −b − b2 − 4ac
x1 = y x2 = . (1.6)
2a 2a
Demuestre que estas raı́ces pueden calcularse mediante las fórmulas equivalentes:
−2c −2c
√
x1 = y x2 = √ . (1.7)
2
b + b − 4ac b − b2 − 4ac
√
Observación: Cuando |b| ≈ b2 − 4ac, se debe proceder con cuidado para evitar la
pérdida de precisión por cancelación. Si b > 0, entonces x1 deberı́a calcularse con (1.7)
y x2 deberı́a calcularse con (1.6); mientras que, si b < 0, entonces x1 deberı́a calcularse
con (1.6) y x2 deberı́a calcularse con (1.7).
11. Use la fórmula adecuada para calcular x1 y x2 , tal como se explica en el Ejercicio 9,
para hallar las raı́ces de las siguientes ecuaciones de segundo grado.
a) x2 − 106x + 1 = 0
b) 10−10 x2 − 1010 x + 1010 = 0
c) x2 − 1000000000001 x + 1 = 0
12. Determine el número de operaciones punto-flotante necesarias para realizar las siguien-
tes operaciones matriciales:
a) El producto AB, donde A ∈ Rm×n y B ∈ Rn×p .
b) El producto Av, donde A ∈ Rm×n y v ∈ Rn .
c) El producto escalar uT v, donde u, v ∈ Rn .
d) El producto uvT , donde u, v ∈ Rn con u un vector columna y vT un vector fila.
uvT
e) El cálculo de la matriz A = uT u
, donde u, v ∈ Rn .
f) El cálculo de la matriz B = A − uvT , donde A ∈ Rn×n y u, v ∈ Rn .
13. Para cada una de las siguientes funciones f (x), calcule el número de condición κP del
problema: calcular un valor de la función f (x), para los números x ∈ R indicados.
a) ln(x + 1) − ln(x), para x ∈ [1, 2].
√
b) x2 + 1 − x, para x ∈ [0, 1].
c) cos2 (x) − sen2 (x), para x ≈ π/4.
d) ex − x − 1, para x ∈ [0, 2π].
1 1
e) x − x+1 , para x ∈ [1000, 5000].
f) x − sen(x), para x ∈ [0, 1].
ex − 1
g) , para x ∈ [0, 1].
x
14. Sea U = (uij ) una matriz triangular superior no singular. Demuestre que que el número
de condición κ∞ (U ) satisface:
máxi=1,...,n {|uii |}
κ∞ (U ) ≥
mı́ni=1,...,n {|uii |}
42 1.4. EJERCICIOS

1 a
15. Dada la matriz A = , determine:
a 1
a) Los valores de a para que A sea mal-condicionada.
b) El número de condición de condición de A, κ2 (A).
16. Use los resultados de los Ejercicios 9 y 10 para construir un algoritmo y un programa
en Python, que calcule las raı́ces de una ecuación cuadrática
√ en todas las situaciones
posibles, incluyendo los casos problemáticos cuando |b| ≈ b2 − 4ac.
17. a) Escriba un programa para calcular la media x y la desviación estándar σ de una su-
cesión finita xi . Su programa debe aceptar un vector x = (x1 , . . . , xn ) de dimensión
n como entrada y producir la media y la desviación estándar de la sucesión como
salida. Para la desviación estándar, utilice la fórmula de dos pasos
" n
#1/2
1 X
σ= (xi − x)2
n−1
i=1
y la fórmula de un paso
" n
!#1/2
1 X
σ= x2i − nx2
n−1
i=1
y compare los resultados para una sucesión de entrada de su elección.

b) ¿Puede construir una sucesión de datos de entrada que ilustre dramáticamente la
diferencia numérica entre estas dos fórmulas matemáticamente equivalentes?
(Precaución: tenga cuidado de tomar la raı́z cuadrada de un número negativo).
18. a) Escriba un programa en Python para construir la matriz A = (aij ) ∈ Rn×n triangular
inferior ([4]): 

1, si i = j,
aij = −1, si i > j,


0, si i < j.
b) Realice un experimento para mostrar que la solución de Ax = b con A ∈ Rn×n como

en a) y b ∈ Rn tal que b = Ax, donde x = (1, 1, . . . , 1), es cada vez menos preciso
cuando n aumenta, a causa del mal-condicionamiento de A. Sea x∗ la solución
obtenida como:
x∗ = np.linalg.solve(A,b).
Presente sus resultados de la siguiente forma:
kx − x∗ k2 kb − Ax∗ k2
n κ2 (A) x∗
kxk2 kbk2
10
20
30
40
50
Bibliografı́a
[1] T. M. Apostol (1985). Calculus, Vol. I y II, Segunda Edición, Editorial Reverté, Barce-
lona.
[2] J. Butcher (1987). The Numerical Analysis of Ordinary Differential Equations: Runge-
Kutta and General Linear Methods, Wiley, Chichester.
[3] J.W. Cooley and J.W. Tukey (1965). An algorithm for the machine calculation of
complex Fourier series. Math. Comput. 19, 297—301.
[4] B.N. Datta, L.M. Hernández-Ramos y M. Raydan (2017). Análisis Numérico: Teorı́a y
Práctica, Editorial de la Universidad Nacional del Sur. Ediuns, Bahı́a Blanca.
[5] R. Fletcher and C.M. Reeves (1964). Function minimization by conjugate gradients.
Computer Journal 7, 149—154.
[6] W. Hackbusch (2016). Iterative Solution of Large Sparse Systems of Equations, Second
Edition, Springer.
[7] E. Hairer and G. Wanner (1980). Solving Ordinary Differential Equations I and II,
Springer-Verlag, Berlı́n.
[8] M. T. Heath 2018). Scientific Computing: An Introductory Survey, Second Edition,

SIAM, Philadelphia.
[9] N.J. Higham (2002). Accuracy and Stability of Numerical Algorithms, Second Edition,
SIAM, Philadelphia.
[10] D. Kincaid y W. Cheney (1994). Análisis Numérico. Las Matemáticas del Cálculo Ci-
entı́fico, Addison-Wesley Iberoamericana, Wilmintong, Delaware.
[11] C. T. Kelley (1995). Iterative Methods for Linear and Nonlinear Equations, SIAM, Phi-
ladelphia.
[12] B. Noble and J. W. Daniel (1977). Applied Linear Algebra, Prentice-Hall, Englewood
Cliffs, NJ.
[13] J. Nocedal Stephen and J. Wright (2006). Numerical Optimization, Second Edition,
Springer, New York.
[14] J. Ortega and W. Rheinboldt (1970). Iterative Solution of Nonlinear Equations in Se-
veral Variables, Academic Press, New York.
[15] M. L. Overton (2001). Computing with IEEE Floating Point Arithmetic, SIAM, Phila-
delphia.
43
44 BIBLIOGRAFÍA
[16] E. Polak and G. Ribière (1969). Note sur la convergence de méthodes de directions
conjuguées. Revue Française d’Informatique et de Recherche Opérationnelle 16, 35—43.
[17] A. Quarteroni, R. Sacco, and F. Saleri (2000). Numerical Mathematics, Springer-

Verlag, New York,
[18] W.F. Ramirez (1998). Computational Methods in Process Simulation, Elsevier Science
& Technology Books.
[19] S. Rosłoniec (2008). Fundamental Numerical Methods for Electrical Engineering,

pringer-Verlag, Berlı́n.
[20] W. Rudin (1980). Principios de Análisis Matemático, Tercera Edición, McGraw-Hill,

México.
[21] Y. Saad (2003). Iterative Methods for Sparse Linear Systems, Second Edition, SIAM,
Philadelphia.
[22] L. R. Scott (2011). Numerical Analysis, Princeton University Press, New Jersey.
[23] G. Strang (1988). Linear Algebra and its Applications, Academic Press, New York.
[24] L. N. Trefethen and D. Bau (1997). Numerical Linear Algebra. SIAM, Philadelphia.
[25] H. Zhang, C. Guo, X. Su, and C. Zhu. (2015). Measurement Data Fitting Based on
Moving Least Squares Method. Mathematical Problems in Engineering, V. 2015, Article
ID 195023, 10 pages.
[26] F. Zhang (1999). Matrix Theory: Basic Results and Techniques, Springer, New York.

Metodos Numericos

Cargado por

Copyright:

Formatos disponibles

Metodos Numericos

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Metodos Numericos

Cargado por

Copyright:

Formatos disponibles

Métodos Numéricos

Antes de leer y usar los programas de este libro iii

1 Herramientas Matemáticas Básicas 1

con la cual se cargan todos los programas del módulo.

En el libro también se utilizan las siguientes librerı́as cientı́ficas y gráficas de Python:

En este libro se usan sistemáticamente conceptos matemáticos básicos que el lector

1.1 Conceptos Básicos del Álgebra Lineal

1.1.1 Espacio Vectorial, Producto Escalar y Normas Vectoriales

Definición 1.1. Un conjunto no vacı́o V de elementos x, y, z, . . . se llama espacio vecto-

Definición 1.2. Decimos que un conjunto no vacı́o W de un espacio vectorial V es un

Definición 1.3. Un sistema de vectores {v1 , . . . , vm } de un espacio vectorial V es lla-

con α1 , . . . , αm ∈ R implica que α1 = · · · = αm = 0. De lo contrario, el sistema se dice

Definición 1.4. Si en un espacio vectorial V se pueden encontrar n vectores lineal-

Es fácil comprobar que Rn es de dimensión n, además, una base de Rn , denominada

Definición 1.5. El producto escalar de dos vectores x e y del espacio vectorial Rn , es un

1) hx, yi = hy, xi,

2) hx, y + zi = hx, yi + hx, zi, para todo z ∈ V ,

3) hλx, yi = λhx, yi, para todo escalar λ,

4) hx, xi ≥ 0 y hx, xi = 0 sólo si x = 0.

En adelante, consideraremos el espacio vectorial Rn . El producto escalar en Rn que

1) kxk ≥ 0, para x ∈ Rn y kxk = 0 si y sólo si x = 0.

2) kαxk = |α| kxk, para x ∈ Rn y α ∈ R.

3) kx + yk ≤ kxk + kyk, para x, y ∈ Rn .

En palabras, las condiciones de la Definición 1.6 requieren que: 1) la norma de un

Los casos particulares de p = 1, p = 2 y p = ∞, son las normas más importantes en la

kxk∞ = máx |xi |.

Definición 1.7. Dos vectores x, y ∈ Rn , se dicen ortogonales si

Un vector x para el cual kxk2 = 1 se dice normalizado. El conjunto de vectores

Ejemplo 1.2. La base canónica de Rn es un conjunto ortonormal. En efecto, se tiene que

eTi ej = (0, . . . , 0, 1 , 0, . . . , 0)T (0, . . . , 0, 1 , 0, . . . , 0) = 0,

siempre que i 6= j; además, es fácil ver que eTi ei = 1.

El ángulo entre dos vectores de Rn se define de la siguiente manera.

Definición 1.8. El águlo θ entre dos vectores x e y de Rn está dado por

xT y ≤ kxk2 kyk2 . (1.1)

En muchos métodos numéricos es necesario construir un sistema de vectores ortogo-

Teorema 1.1 (Proceso de Gram-Schmidt). Sea {v1 , . . . , vm } un sistema de vectores li-

Definición 1.9. Una matriz real A = (aij ) de orden m × n es un arreglo rectangular

Definición 1.10. La traspuesta de una matriz A ∈ Rm×n es la matriz AT ∈ Rn×m tal

En el siguiente ejemplo mostramos algunas matrices y sus traspuestas.

>>> import numpy as np

Definamos las principales operaciones con matrices.

Definición 1.11. Sean A = (aij ), B = (bij ) y C = (cij ) matrices.

• Suma: La suma de las matrices A ∈ Rm×n y B ∈ Rm×n , es una matriz C = A + B

• Multiplicación por un escalar: La multiplicación de un número real α por una

• Multiplicación por otra matriz: La multiplicación de dos matrices A ∈ Rm×p y

En general, el producto de matrices no es conmutativo.

En Python las operaciones: A + B, A − B, aC y AB, se realizan respectivamente como:

>>> np . a r r a y ([[ −1 ,3] ,[5 ,4] ,[0 , −2]])+ np . a r r a y ( [ [ 2 , 5 ] , [ − 6 , − 2 ] , [ 8 , 6 ] ] )

Una matriz cuadrada importante es la matriz identidad de orden n dada por:

en otras palabras, I = (δij ) es una matriz de orden n tal que

La matriz identidad satisface la igualdad IA = AI = A, para toda matriz cuadrada A de

Definición 1.13. Una matriz cuadrada A de orden n se dice invertible o no singular, si

Teorema 1.2. Si A y B son matrices no singulares de orden n, entonces:

ii) (AB)−1 = B −1 A−1 .

Ahora definiremos un número asociado con la matriz cuadrada A = (aij ) de orden n,

Definición 1.14. El determinante de una matrix de orden (n − 1) obtenida eliminando

Aij = (−1)i+j Mij ,

es llamado el cofactor de aij .

Definición 1.15. El determinante de una matriz cuadrada A de orden n está definido