Apuntes MIII

MATEMÁTICAS III
CÁLCULO NUMÉRICO
GRADO INGENIERÍA QUÍMICA
Departamento de Matemática Aplicada
Universidad de Salamanca
Mabel Asensio Sevilla
Julio-2011
Revisado Septiembre 2015
2
Índice general
0. Errores 7
0.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
0.2. Error absoluto y relativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
0.3. Errores de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1. Ecuaciones y sistemas de ecuaciones no lineales 9
1.1. Localización y separación de raı́ces de una ecuación. . . . . . . . . . . . . . . . . . . . . . 10
1.2. Ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.1. Método de bisección. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2.2. El método de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.3. El método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.2.4. Modificaciones del método de Newton. . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.2.5. Método de la secante. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.3. Sistemas de ecuaciones no lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.3.1. Método de punto fijo en varias variables. . . . . . . . . . . . . . . . . . . . . . . . . 22
1.3.2. Método de Newton en varias variables. . . . . . . . . . . . . . . . . . . . . . . . . . 23
2. Sistemas de ecuaciones lineales 27
2.1. Generalidades sobre matrices y vectores . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3
4 ÍNDICE GENERAL
2.2. Métodos directos de resolución de sistemas de ecuaciones lineales . . . . . . . . . . . . . . 32
2.2.1. Matrices triangulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.2. Eliminación gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.3. Técnicas de pivotaje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.4. Factorización LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.2.5. Matrices especiales: factorización LDLt , Cholesky . . . . . . . . . . . . . . . . . . 39
2.2.6. Aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.3. Métodos iterativos de resolución de sistemas de ecuaciones lineales . . . . . . . . . . . . . 40
2.3.1. Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.3.2. Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.3.3. Métodos de relajación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.3.4. Control de parada de las iteraciones . . . . . . . . . . . . . . . . . . . . . . . . . . 45
2.3.5. Resultados de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3. Interpolación 51
3.1. Interpolación polinómica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.1.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.1.2. Tipo de función interpoladora . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.1.3. Existencia y unicidad del polinomio interpolador . . . . . . . . . . . . . . . . . . . 52
3.1.4. Métodos de cálculo del polinomio interpolador. . . . . . . . . . . . . . . . . . . . . 53
3.1.5. Error de interpolación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.2. Interpolación de Hermite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.2.1. Ejemplo sencillo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2.2. Problema de Hermite generalizado . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2.3. Caso particular: el polinomio de Taylor . . . . . . . . . . . . . . . . . . . . . . . . 58

ÍNDICE GENERAL 5
3.2.4. Método de las diferencias divididas de Newton generalizado . . . . . . . . . . . . . 58
3.2.5. Ejemplo sencillo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4. Aproximación numérica. 61
4.1. Introducción. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.1.1. Conjunto de abscisas de aproximación . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.1.2. Funciones básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.1.3. Medida de la magnitud del error: normas funcionales . . . . . . . . . . . . . . . . . 62
4.2. Aproximación por mı́nimos cuadrados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.2.1. Definición del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.2.2. Productos escalares asociados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.2.3. Ecuaciones normales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.2.4. Un ejemplo sencillo: la recta de regresión . . . . . . . . . . . . . . . . . . . . . . . 66
4.3. Ortogonalización. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3.1. Ortogonalización de Gram-Schmidt . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5. Integración y derivación numéricas 69
5.1. Integración numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.1.1. Integración vı́a interpolación. Fórmulas de Newton-Cotes . . . . . . . . . . . . . . 70
5.1.2. Método de los coeficientes indeterminados . . . . . . . . . . . . . . . . . . . . . . . 72
5.1.3. Cambio de intervalo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.1.4. Cuadratura gaussiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.2. Derivación numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2.1. Derivadas primeras. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2.2. Derivadas de orden superior. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

6 ÍNDICE GENERAL
Capı́tulo 0
Errores
0.1. Introducción
Un método numérico es un método “aproximado” para la resolución de un problema matemático, éste,

a su vez, puede representar una modelización matemática de un problema fı́sico, quı́mico o del mundo
real. En la práctica, la solución al problema real que nosotros conoceremos será la que nos proporcione
el método numérico, que en general no va a coincidir con la solución exacta del problema real, ya que va
a estar afectada de diversos tipos de errores:
Experimentales: la presencia de errores puede comenzar en la misma formulación del problema

real, pues los datos se pueden haber obtenido de ciertas mediciones u otras observaciones experi-
mentales, siempre susceptibles de errores.
De modelización: debidos a la aproximación de la realidad del modelo matemático elegido.
De discretización o de truncamiento: debidos a la propia naturaleza del método numérico

elegido para resolver el problema matemático.
De redondeo: debidos a las restricciones aritméticas de los ordenadores y la limitada capacidad

humana, frente a la infinidad de cifras decimales de los números reales. Es necesario delimitar su
acumulación, ya que es habitual llevar a cabo un elevado número de operaciones en la resolución
de los métodos numéricos.
0.2. Error absoluto y relativo
Sea x el valor exacto de un número real y x0 un valor aproximado. Definimos:
Error absoluto de x: ϵ(x) = |x − x0 |.
7
8 Errores
ϵ(x)
Error relativo de x: e(x) = |x| .
El error absoluto da una referencia cuantitativa de la bondad de la aproximación, medida por la dis-
tancia que separa el valor exacto del aproximado. El error relativo proporciona una referencia cualitativa,
en tanto en cuanto refleja la proporción del error absoluto con respecto a la magnitud que se trata de
aproximar: en este sentido, no es lo mismo un error de una unidad cuando se aproxima el valor exacto
de π = 3.14159 . . . que cuando se aproxima el valor exacto del número de Avogadro (aproximadamente
igual a 6.022 · 1023 ).
Ejemplo 0.1 Comparar los errores absolutos y relativos en las aproximaciones 3.1 de 3 y 3099 de 3000.
Diremos que la aproximación x0 tiene p cifras decimales exactas si ϵ ≤ 10−p . Obsérvese que esto no
indica que hayan de coincidir las p primeras cifras decimales de x y x0 . Por ejemplo, si x = 1 y x0 = 0.9999
se tiene que ϵ ≤ 10−4 y, por tanto, 0.9999 aproxima a 1 con las cuatro cifras decimales exactas (aunque
no coincida ninguno de los decimales de ambas cifras).
0.3. Errores de redondeo
Dado un número real x expresado en su forma decimal
x = an an−1 . . . a0 .a−1 a−2 . . . a−k a−k−1 . . . , 0 ≤ ak ≤ 9, k ∈ Z
se llama parte decimal de x a la secuencia a−1 a−2 . . . a−k a−k−1 . . .. Por ejemplo, 17.352 tiene por
parte decimal a 352.
Si efectuamos los nuestros cálculos en una máquina que puede representar números con k cifras
decimales, esta representación se puede hacer de dos formas: por truncamiento, cortando la parte decimal
para dejarla en k cifras;
xt = an an−1 . . . a0 .a−1 a−2 . . . a−k
o por redondeo, si la cifra a−k−1 es menor que 5, entonces el resultado es el mismo que por truncamiento,
y si la cifra a−k−1 es igual o mayor que 5, entonces se añade 1 a la cifra k-ésima y se trunca el número
restante.
Observar que |x − xt | ≤ 10−k mientras que |x − xr | ≤ 12 10−k

Capı́tulo 1
Ecuaciones y sistemas de ecuaciones

no lineales
En este tema repasaremos uno de los problemas básicos del cálculo numérico:
Dada una función f real de variable real, hallar los valores de la variable x que satisfagan la ecuación
f (x) = 0.
La función f puede ser polinómica, transcendente o incluso puede que no dispongamos de una ex-
presión explı́cita de la misma, por ejemplo, si es la solución de una ecuación diferencial. Los valores que
buscamos son los valores x̄ que anulan dicha función. A estos valores se les denomina raı́ces o soluciones
de la ecuación, o también ceros de la función f (x). Geométricamente representan las abscisas de los
puntos de corte de la gráfica y = f (x) con el eje OX.
Definición 1.1 Multiplicidad de una raı́z.
Una raı́z x̄ de la ecuación f (x) = 0 se dice que tiene multiplicidad n si
f (x̄) = f ′ (x̄) = f ′′ (x̄) = . . . = f n−1) (x̄) = 0 y f n) (x̄) ̸= 0
Si la multiplicidad de una raı́z es 1, diremos que es una raı́z simple.
En general, las raı́ces de una ecuación no lineal no se pueden calcular de forma exacta, sino que se
recurre a métodos numéricos que permiten obtener aproximaciones numéricas de las mismas. El objetivo
de este capı́tulo es presentar algunos de estos métodos, pero antes veremos algunos resultados que nos
permitirán localizar y separar previamente las raı́ces de una ecuación. Posteriormente veremos algunos
de los métodos más clásicos para el cálculo de raı́ces de ecuaciones, como el método de la bisección, el
método de punto fijo y el método de Newton, y algunas de sus modificaciones. Por último veremos como
adaptar alguno de estos métodos al caso de sistemas no lineales, como el método de punto fijo y los
métodos Quasi-Newton.
9
10 Ecuaciones y sistemas de ecuaciones no lineales
El problema de hallar las raı́ces de una ecuación, f (x) = 0, aparece frecuentemente en ingenierı́a. Por
ejemplo, para calcular el volumen V de un gas de van der Waals como función de la temperatura absoluta
T , la presión P , el número de moles n y los parámetros de Van der Walls a, b, la ecuación de estado es,
( a )
P+ (V − b) = nRT,
V2
que conduce a una ecuación polinómica de grado 3 en V ,
P V 3 − (P b + nRT )V 2 + aV − ab = 0.
En teorı́a de la difracción de la luz necesitamos las raı́ces de la equación
x − tan x = 0.
En el cálculo de órbitas planetarias necesitamos las raı́ces de la ecuación de Kepler
x − a sin x = b, (1.1)
para varios valores de a y b. En teorı́a de la combustión
x = δeγx , (1.2)
para varios valores de γ y δ.
1.1. Localización y separación de raı́ces de una ecuación.
El proceso de localización y separación de raı́ces de una ecuación es una tarea previa a la aplicación de
un método numérico para el cálculo de estas raı́ces. Consiste en obtener información de las zonas donde
se encuentran las raı́ces reales de la ecuación, para posteriormente buscar intervalos [a1 , b1 ], [a2 , b2 ], . . .
que contengan una y sólo una de estas raı́ces.
Dada una ecuación no lineal f (x) = 0 con n raı́ces reales distintas, x̄1 , . . . , x̄n , se pretende hallar n
intervalos disjuntos Ii = [ai , bi ] para i = 1, . . . , n de modo que x̄i ∈ Ii , i = 1, . . . , n.
A veces puede obtenerse algún tipo de información gráfica si se transforma la ecuación f (x) = 0 en
otra del tipo g(x) = h(x) y se cotejan los puntos de corte de las gráficas de g(x) y h(x). Esto sólo da una
idea gráfica de donde están los ceros, pero no puede servir como prueba de localización y separación de
las raı́ces de una ecuación, ya que en algunos casos la información gráfica que proporciona un ordenador
puede no ajustarse a la realidad.
Es el estudio analı́tico de la función f (x) el que puede aportarnos la información necesaria, abordando
diversos aspectos:
1.1 Localización y separación de raı́ces de una ecuación. 11
(a) Crecimiento. Estudio de los intervalos de crecimiento de f (x). Si [a, b] es un intervalo de crecimiento
(resp. decrecimiento) monótono de f (x), entonces a lo más habrá una única raı́z de f (x) = 0 en ese
intervalo. El estudio de los intervalos de crecimiento de una función supone hallar los ceros de su
derivada, lo que en ocasiones puede ser tanto o más complejo que el problema de partida.
(b) Teorema de Bolzano. Se trata de aplicar el toerema de Bolzano a cada uno de los intervalos en los
que se sospecha que hay una raı́z. Esto require que se satisfagan las hipótesis de este teorema, lo
cual no siempre ocurre.
(c) Sucesiones de Sturm. No siempre se conoce como determinar una sucesión de Sturm para una
función dada. Estudiaremos el caso de las funciones polinómicas.
Ecuaciones polinómicas
Dado un polinomio
P (x) = a0 xn + a1 xn−1 + . . . + an−1 x + an
donde ai ∈ R para i = 0, 1, . . . , n diremos que P (x) = 0 es una ecuación polinómica.
El teorema fundamental del Álgebra nos dice que la ecuación polinómica P (x) = 0 con coeficientes
reales tiene n raı́ces reales y complejas contadas con sus multiplicidades. Las raı́ces complejas aparecen
en pares conjugados (si a + bi es raı́z, entonces a − bi también lo es).
Veamos algunos resultados que permiten localizar ceros de un polinomio.
Proposición 1.1 Si x̄ es una raı́z de P (x) = 0 entonces:
1 A
< |x̄| < 1 + siendo A = maxi≥1 |ai |
1 + |aAn | |a0 |
Proposición 1.2 (Regla de Laguerre) Dado c ∈ R+ podemos escribir P (x) = (x − c)C(x) + r con
C(x) = b0 xn−1 + . . . + bn−2 x + bn−1 y r ∈ R. Si r ≥ 0 y bi ≥ 0 para i = 0, 1, . . . , n − 1 ó r ≤ 0 y
bi ≤ 0 para i = 0, 1, . . . , n − 1, entonces el número real c es una cota superior de las raı́ces positivas de la
ecuación.
Proposición 1.3 Sea R(x) = an xn +an−1 xn−1 +. . .+a0 , es decir, R(x) = xn P ( x1 ) para x ̸= 0. Por tanto
P (x̄) = 0 ⇔ R( x̄1 ) = 0. Esto nos permite obtener una cota inferior de las raı́ces positivas de P (x) = 0
puesto que si c′ una cota superior de las raı́ces positivas de R(x) = 0 obtenida mendiante la regla de
1
Laguerre, entonces c′ es una cota inferior de las raı́ces positivas de P (x) = 0.
Proposición 1.4 Sea H(x) = P (−x), entonces P (x̄) = 0 ⇔ H(−x̄) = 0, esto es, si x̄ es una raı́z negativa
de P (x) = 0, entonces −x̄ es una raı́z positiva de H(x) = 0. Esto nos permite obtener cotas inferiores y
superiores de las raı́ces negativas de P (x) = 0: si c y c′ son cotas superior e inferior de las raı́ces positivas
de H(x) = 0, respectivamente, entonces −c′ es cota superior de las raı́ces negativas de P (x) = 0 y −c
cota inferior de las raı́ces negativas.
Ejemplo 1.1 Dado el polinomio P (x) = x4 + 2x3 − 3x2 − 4x − 1, acotar las raı́ces de P (x) = 0 tanto
como se pueda.
Definición 1.2 Una sucesión de Sturm para una función f (x) en [a, b] es un conjunto f0 (x) = f (x),
f1 (x), . . . , fn (x) de funciones continuas en dicho intervalo tales que:
(a) fn (x) ̸= 0∀x ∈ [a, b], es decir, el signo de fn (x) permanece constante en [a, b]
(b) Si fi (c) = 0 con c ∈ [a, b] entonces fi−1 (c) · fi+1 (c) < 0, es decir, tienen signos opuestos y no se
anulan en c.
f0 (x)
(c) Si f0 (c) = 0 con c ∈ [a, b] entonces f1 (x) pasa de negativa a positiva en c
La importancia de las sucesiones de Sturm radica en el resultado siguiente:
Teorema 1.1 (Teorema de Sturm) Sea f0 (x), f1 (x), . . . , fn (x) una sucesión de Sturm para f (x) = f0 (x)
en el intervalo [a, b]. Consideremos las siguientes sucesiones en las que sig(d) denota el signo de d (indis-
tintamente ± cuando d = 0)
sig(f0 (a)), sig(f1 (a)), . . . , sig(fn (a))
sig(f0 (b)), sig(f1 (b)), . . . , sig(fn (b)).
Denotemos por N1 el número de cambios de signo en la primera sucesión, y por N2 el número de cambios
de signo en la segunda (siempre ha de ser N1 ≥ N2 ). Entonces el número de raı́ces de la ecuación f0 (x) = 0
en el intervalo [a, b] viene dado por N1 − N2 .
Por tanto, si conocemos una sucesión de Sturm para una función f (x), podremos separar todos sus
ceros reales. Lamentablemente, no hay procedimientos sistemáticos para formar sucesiones de Sturm para
cualesquiera funciones dadas, salvo contadas excepciones, como es el caso de los polinomios, para los que
la sucesión de Sturm se construte de la siguiente forma:
f0 (x) = P (x), f1 (x) = P ′ (x), fi+1 (x) = −ri (x)

donde ri (x) es el resto de dividir fi−1 entre fi , es decir, fi−1 (x) = ci (x) · fi (x) + ri (x).
Ejemplo 1.2 Dado el polinomio P (x) = x4 + 2x3 − 3x2 − 4x − 1, construir una sucesión de Sturm para
este polinomio y separar las raı́ces de P (x) = 0.
1.2. Ecuaciones no lineales
1.2.1. Método de bisección.
Se basa en la aplicación reiterada del teorema de Bolzano: Si f es una función continua definida
sobre un intervalo cerrado [a, b] tal que f (a).f (b) < 0 entonces f debe tener un cero en (a, b).
1.2 Ecuaciones no lineales 13
El método de la bisección explota esta propiedad de la siguiente manera:
a+b
(a) Tomamos c = 2
(b) Si f (a).f (c) < 0, entonces f tiene un cero en (a, c) y hacemos b ←− c

Si f (a).f (c) > 0, entonces f (c).f (b) < 0 y f tiene un cero en (c, b) y hacemos a ←− c
Si f (a).f (c) = 0, está claro que f (c) = 0 y ya hemos encontrado un cero.
En las dos primeras situaciones del punto 2, hemos reducido el problema a la búsqueda de ceros en
un intervalo de longitud la mitad que la del intervalo original y repetimos el proceso.
La situación f (c) = 0 es poco probable que se dé en la práctica, debido a los errores de redondeo.
Ası́, el criterio para concluir no debe depender de que f (c) = 0, sino que permitiremos una tolerancia
razonable, tal como |f (c)| < ε, para cierto ε suficientemente pequeño.
Pseudo-código del algoritmo de la bisección
entrada a, b, M, δ, ε
u ←− f (a), v ←− f (b), e ←− b − a
si sign(u) = sign(v) entonces parar
para k = 1, ..., M hacer
• e ←− 2e , c ←− a + e, w ←− f (c)
• salida k, c, w, e
• si |e| < δ or |w| < ε entonces parar
• si sign(w) ̸= sign(u) entonces b ←− c, v ←− w
• sino a ←− c, u ←− w
• fin condicional
fin bucle
Varias de las partes de este pseudo-código necesitan una explicación adicional. En primer lugar, el
b−a a+b
punto medio c se calcula como c ←− a + 2 en lugar de c ←− 2 . Al hacerlo ası́ se sigue la estrategia
general de que, al efectuar cálculos numéricos, es mejor calcular una cantidad añadiendo un pequeño
término de corrección a una aproximación obtenida previamente. En segundo lugar, es mejor determinar
si la función cambia de signo en el intervalo recurriendo a que sign(w) ̸= sign(u) en lugar de utilizar
w.u < 0 ya que esta última requiere una multiplicación innecesaria. Por otra parte e corresponde al
cálculo de la cota del error que se establece más adelante.
En el algoritmo hay tres criterios que pueden detener la ejecución:

M , señala el máximo número de iteraciones, un algoritmo correctamente diseñado tiene que ser
finito.
Por otra parte la ejecución del programa se puede detener, ya sea cuando el error es suficientemente
pequeño o cuando lo es el valor de f (c). Los parámetros δ y ε controlan esta situación. Se pueden
dar ejemplos en los que se satisface uno de los dos criterios sin que el otro se satisfaga.
Teorema 1.2 : Análisis del error
Sea f continua en [a, b] = [a0 , b0 ] con f (a).f (b) < 0. Sean [a0 , b0 ], [a1 , b1 ], ..., [an , bn ] los intervalos
sucesivos generados por el método de la bisección. Entonces los lı́mites lı́mn→∞ an , lı́mn→∞ bn existen,
an +bn
son iguales y representan un cero de f . Si r = lı́mn→∞ cn con cn = 2 , entonces
b0 − a0
|r − cn | ≤
2n+1
Demostración:
Por la propia construcción del algoritmo, tenemos,
a0 ≤ a1 ≤ ... ≤ b0
b0 ≥ b1 ≥ ... ≥ a0
bn − a n
bn+1 − an+1 = , n≥0
2
La sucesión {an } converge debido a que es creciente y está acotada superiormente.
La sucesión {bn } converge por ser decreciente y estar acotada inferiormente.
Además, se tiene,
bn−1 − an−1 b0 − a0
bn − a n = = ... =
2 2n
Ası́
b0 − a 0
lı́m bn − lı́m an = lı́m =0
n→∞ n→∞ n→∞ 2n
Si escribimos r = lı́m an = lı́m bn , tomando lı́mites en la desigualdad f (an ).f (bn ) < 0, resulta f (r)2 =
f (r).f (r) ≤ 0, es decir f (r) = 0.
Finalmente, en la etapa en la que se ha construido el intervalo [an , bn ], si se detiene en este momento

el algoritmo, sabemos que la raı́z de la ecuación se encuentra en ese intervalo. La mejor estimación para
an +bn
esa raı́z será el punto medio cn = 2 y el error cometido verificará
bn − an 1 b0 − a 0 b0 − a 0
|r − cn | ≤ ≤ n
= n+1
2 2 2 2
1.2.2. El método de punto fijo
Utilizaremos este método para resolver ecuaciones de la forma x = g(x). Observemos que si queremos
hallar las raı́ces de una ecuación f (x) = 0, podemos ponerla de la forma anterior, por ejemplo, haciendo
g(x) = x − f (x) o más generalmente g(x) = x − ρ(x)f (x) donde ρ(x) ̸= 0, es una función adecuadamente
elegida, que puede ser constante o no.
De manera más precisa el problema planteado es el siguiente:
Dada g : [a, b] −→ [a, b] función continua, hallar x ∈ [a, b] tal que x = g(x).
Teorema 1.3 : Existencia del punto fijo.
Sea g : [a, b] −→ [a, b] continua, entonces existe al menos un x ∈ [a, b] tal que x = g(x).
Demostración:
Si a = g(a) o b = g(b) entonces a o b es una solución. Supongamos pues que a ̸= g(a) y que b ̸= g(b).
Pongamos f (x) = x − g(x), tendremos, f (a) = a − g(a) < 0 y f (b) = b − g(b) > 0. Por el teorema de
Bolzano existe al menos x̄ ∈ (a, b) tal que f (x̄) = 0, es decir, x̄ = g(x̄).
Teorema 1.4 : Unicidad del punto fijo.
Sea g : [a, b] −→ [a, b] continua y contractiva, es decir, existe k < 1 tal que |g(x) − g(y)| ≤ k|x −
y|, ∀x, y ∈ [a, b], entonces el punto fijo x̄ es único.
Demostración:
Sean x¯1 y x¯2 dos puntos fijos de g, x¯1 ̸= x¯2 , es decir, x¯1 , x¯2 ∈ [a, b], x¯1 = g(x¯1 ) y x¯2 = g(x¯2 ).
|x¯1 − x¯2 | = |g(x¯1 ) − g(x¯2 )| ≤ k|x¯1 − x¯2 | < |x¯1 − x¯2 |

Observación: Si g es diferenciable y existe un número k < 1 tal que |g ′ (x)| < k para todo x ∈ [a, b],
entonces para ξ ∈ [a, b], resulta |g(x) − g(y)| = |g ′ (ξ)||x − y| ≤ k|x − y|.
El algoritmo de punto fijo o iteración funcional es:
Dado un x0 ∈ [a, b],
calculado xn , obtenemos xn+1 = g(xn )
Pseudo-código del algoritmo de punto fijo
entrada x0 , M, ε
x ← x0
Para k = 1, ..., M hacer
• x1 ← x, x ← g(x), e ← |x − x1 |
• salida k, x, e
• si e < ε entonces parar
fin bucle
Teorema 1.5 : Teorema de convergencia y análisis del error.
Sea g : [a, b] −→ [a, b] continua y contractiva, es decir, tal que
|g(x) − g(y)| < k|x − y| ∀x, y ∈ [a, b], k<1
entonces la sucesión xn generada por el algoritmo de punto fijo verifica
lı́m xn = x̄
n→∞
siendo x̄ el único punto fijo de g, y además,
kn
|xn − x̄| ≤ |x1 − x0 |
1−k
Demostración:
|xn+1 − x̄| = |g(xn ) − g(x̄)| ≤ k|xn − x̄| ≤ ... ≤ k n |x0 − x̄|
de donde
lı́m |xn − x̄| ≤ |x0 − x̄| lı́m k n = 0
n→∞ n→∞
pues k < 1.
Por otro lado,
|xn+1 − xn | = |g(xn ) − g(xn−1 )| ≤ k|xn − xn−1 | ≤ ... ≤ k n |x1 − x0 |
Para m > n ≥ 1 tendremos,
|xm − xn | = |xm − xm−1 + xm−1 − xm−2 + xm−2 − ... + xn+1 − xn |

≤ |xm − xm−1 | + |xm−1 − xm−2 | + ... + |xn+1 − xn |
≤ (k m−1 + k m−2 + ... + k n )|x1 − x0 |
≤ k n (1 + k + ... + k m−n−1 )|x1 − x0 |
Pasando al lı́mite cuando m → ∞ se obtiene
kn
|xn − x̄| ≤ |x1 − x0 |
1−k
Definición 1.3 Orden de convergencia, convergencia lineal, cuadrática y orden α.
Supongamos que {xn }∞ n=1 es una sucesión convergente cuyo lı́mite es p. Sea en = xn − p. Si existen
dos constantes λ > 0 y α > 0 tales que
|en+1 |
limn→∞ =λ
|en |α
diremos que {xn } converge hacia p, con orden α. En particular:
Si α = 1, diremos que la convergencia es lineal.
Si α = 2, diremos que la convergencia es cuadrática.
Si 1 < α < 2, diremos que la convergencia es superlineal.
Orden de convergencia del método de punto fijo

El método de punto fijo tiene convergencia lineal si g ′ es continua y g ′ (x̄) ̸= 0 siendo x̄ el punto fijo de g.
En efecto,
en+1 = xn+1 − x̄ = g(xn ) − g(x̄) = g ′ (ξn )(xn − x̄) = g ′ (ξn )en
donde ξn ∈ [xn , x̄], finalmente
|en+1 |
limn→∞ = limn→∞ |g ′ (ξn )| = |g ′ (x̄)| = λ > 0
|en |
1.2.3. El método de Newton
Consideremos de nuevo el problema de buscar las raı́ces de una ecuación del tipo f (x) = 0. Si
f (x), f ′ (x) y f ′′ (x) son continuas cerca de una raı́z x̄, esta información adicional sobre la naturaleza
de f (x) puede usarse para desarrollar algoritmos que produzcan sucesiones {xk } que converjan a x̄ más
rápidamente que el método de bisección o de punto fijo. El método de Newton-Raphson, o simplemente
de Newton, que descansa en la continuidad de f ′ (x) y f ′′ (x), es uno de los algoritmos más útiles y mejor
conocidos.
Supongamos que x̄ es una raı́z de la ecuación anterior y supongamos además que f es dos veces
derivable con continuidad. Si x es una aproximación de x̄, usando el desarrollo de Taylor, podemos
escribir,
1
0 = f (x̄) = f (x) + f ′ (x)(x̄ − x) + f ′′ (ξ)(x̄ − x)2
2
Si x está cerca de x̄, (x̄ − x)2 es un número pequeño y podremos despreciar el último término frente a
los otros, y x̄ vendrá dado aproximadamente por
f (x)
x̄ ≈ x −
f ′ (x)
Como hemos despreciado el término cuadrático este valor no será exactamente x̄, pero es de esperar que
será una mejor aproximación que el valor x de partida. De aquı́ se obtiene el algoritmo de Newton:
x0 , valor cercano a x̄.
Calculado xn , obtenemos xn+1 ,
f (xn )
xn+1 = xn − Fórmula de Newton-Raphson
f ′ (xn )
El método de Newton también es conocido como método de la tangente, ya que si trazamos la tangente
a la curva y = f (x) en el punto (xn , f (xn )) obtenemos la recta y = f (xn ) + f ′ (xn )(x − xn ) que corta al
f (xn )
eje y = 0 en el punto de abscisa x = xn − f ′ (xn ) , que es precisamente el valor de xn+1 en la fórmula de
Newton-Raphson.
El método de Newton se puede interpretar como un método de punto fijo, pues buscamos el punto
fijo de la función x − f (x)/f ′ (x).
Pseudocódigo del algoritmo de Newton
entrada x0 , M, δ, ε
v ← f (x0 )
salida 0, x0 , v
si |v| < ε entonces parar

v
• x1 ← x0 − f ′ (x0 )
• v ← f (x1 )
• salida k, x1 , v
• si |x1 − x0 | < δ o |v| < ε entonces parar
• x0 ← x1
fin bucle
Un inconveniente de la sucesión generada por la fórmula de Newton-Raphson es que no siempre se tiene

asegurada la convergencia hacia x̄, incluso tomando x0 próximo a x̄. Nos preguntamos que condiciones
hay que exigir a x0 y a f para que la sucesión {xn } generada por la fórmula de Newton-Raphson sea
convergente a x̄.
El resultado más general de convergencia del método de Newton es el siguiente:
Teorema 1.6 : Convergencia del método de Newton
Supongamos que f ∈ C 2 [a, b] y que x ∈ [a, b] es una raı́z simple de f , es decir, f (x) = 0 y f ′ (x) ̸= 0.
Entonces existe una constante δ > 0 tal que la sucesión {xn }∞
0 generada por el método de Newton
converge a x cualquiera que sea la aproximación inicial x0 ∈ [x − δ, x + δ], y además la convergencia es
cuadrática, es decir, existe una contante C > 0 tal que
|xn+1 − x| ≤ C|xn − x|2 , n ≥ 0.
Existe un resultado que, partiendo de un intervalo inicial adecuado, cuando este existe, nos permite
asegurar la convergencia del método de Newton y nos indica el valor incical con el que comenzar la
iteración.
Teorema 1.7 (Regla de Fourier) Sea f (x) : [a, b] → R, continua y dos veces diferenciable con continuidad
en [a, b] y tal que verifica:
f (a) · f (b) < 0, es decir sigf (a) ̸= sigf (b),
f ′ (x) ̸= 0, ∀x ∈ [a, b],
f ′′ (x) ̸= 0, ∀x ∈ [a, b].

Entonces el método de Newton converge si tomamos x0 = a o x0 = b de tal forma que f (x0 ) · f ′′ (x0 ) > 0,
es decir, tomando como valor inicial x0 el extremo del intervalo en el que la función y su segunda derivada
tienen el mismo signo.
Ejemplo 1.3 Aproximar la raı́z cuadrada de 3 con el método de Newton partiendo del intervalo [1, 2].
1.2.4. Modificaciones del método de Newton.
El método de Newton presenta problemas cuando x̄, la raı́z de f (x) = 0 que se busca es múltiple.
Esta situación se detecta porque la convergencia del método se hace especialmente lenta. La fórmula de
Newton-Raphson puede modificarse para adpatarse a este caso:
f (xn )
xn+1 = xn − k
f ′ (xn )
donde k representa la multiplicidad de x̄.
En la práctica, el problema es que no conocemos k, pero a esto nos puede ayudar el comportamiento
de f y sus derivadas al aplicar el método.
Ejemplo 1.4 La ecuación x − sin x = 0 tiene una raı́z triple en x = 0. Aplicar el método de Newton y
su modificación a este ejemplo partiendo del intervalo [−1, 1].
1.2.5. Método de la secante.
En muchas aplicaciones, f (x) no viene dada por una fórmula explı́cita, por ejemplo si f (x) es el resul-
tado de algún algoritmo numérico o de un proceso experimental. Como f ′ (x) no estará en consecuencia
disponible, el método de Newton deberá modificarse de modo que únicamente requiera valores de f (x).
Cuando f ′ (x) no está disponible, podemos reemplazarlo por una aproximación suya, por ejemplo,
tomando la pendiente de la secante formada a partir de dos puntos sobre la gráfica de la función, es decir,
aproximamos la derivada en un punto xn mediante
f (xn + hn ) − f (xn )
f ′ (xn ) ≈ an =
hn
Una manera de aproximar f ′ (xn ) es utilizar los valores de f en xn y xn−1 , es decir,
f (xn ) − f (xn−1 )
f ′ (xn ) ≈ an =
xn − xn−1
obtenemos ası́ el llamado método de la secante,

1.3 Sistemas de ecuaciones no lineales. 21
x0 , x1
xn+1 = xn − f (xn ) f (xxnn)−f

−xn−1
(xn−1 )
y que necesita de una sola evaluación de la función en cada iteración.
Pseudocódigo del algoritmo de la secante
entrada x0 , x1 , M, δ, ε
v0 ← f (x0 ), v1 ← f (x1 )
salida 0, x0 , v0 , 1, x1 , v1
si |v1 | < ε entonces parar
• x2 ← x1 − v1 xv11 −x
−v0
0
• v2 ← f (x2 )
• salida k, x2 , v2
• si |x2 − x1 | < δ o |v2 | < ε entonces parar
• x 0 ← x 1 , v0 ← v1
• x 1 ← x 2 , v1 ← v2
fin bucle
Si hay convergencia en el método de la secante, esta es superlineal. El orden de convergencia es el

número áureo.
√
1+ 5
α= ≈ 1.62
2
1.3. Sistemas de ecuaciones no lineales.
Nos ocupamos ahora del problema del cálculo numérico de los ceros de funciones vectoriales de varias
variables reales que tienen la forma general F (x) = 0, donde F : ℜn → ℜn viene definida por sus n
componentes fi : ℜn → ℜ para i = 0, . . . , n, esto es, un sistema de n ecuaciones no lineales con n
incógnitas:


 f1 (x1 , . . . , xn ) = 0

f2 (x1 , . . . , xn ) = 0

 ...... ...

fn (x1 , . . . , xn ) = 0
El salto de una a varias variables conlleva la introducción de nuevos conceptos que generalicen los
habituales en ℜ. El concepto de norma como distancia, compatible con las operaciones de la estructura
de espacio vectorial, generaliza el concepto de valor absoluto en ℜ, y con él es fácil expresar el análisis
en varias variables manteniendo una semejanza casi total con el caso de una variable. La dificultad es
mayor, pero las herramientas son las mismas generalizando lo que hemos estudiado en una variable.
Recordemos algunos resultados del análisis en varias variables que necesitaremos.
Sea D un conjunto cerrado de ℜn y f : D ⊂ ℜn → ℜ, entonces,
• f tiene lı́mite l en x0 (lı́mx→x0 f (x) = l), si ∀ϵ > 0 ∃δ > 0 tal que |f (x) − l| < ϵ ∀x ∈ D con
0 < ∥x − x0 ∥ < δ.
• f es continua en x0 ∈ D si ∃ lı́mx→x0 f (x) y lı́mx→x0 f (x) = f (x0 ).

• f es continua en D si lo es en cada punto de D.
• Sea x0 ∈ D, si ∃δ > 0, K > 0 con | ∂f∂x(x)

j
| ≤ K para cada j = 1, . . . , n siempre que ∥x − x0 ∥ < δ
y x ∈ D ⇒ f es continua en x0 .
Sea F : D ⊂ ℜn → ℜn , con F = (f1 , . . . , fn )t entonces,
• lı́mx→x0 F (x) = L = (l1 , . . . , ln )t ⇔ lı́mx→x0 fi (x) = li para cada i = 1, . . . , n.
• F es continua en x0 ∈ D si ∃ lı́mx→x0 F (x) y lı́mx→x0 F (x) = F (x0 ).

• F es continua en D si lo es en cada punto de D.
1.3.1. Método de punto fijo en varias variables.
Sea G : D ⊂ ℜn → ℜn , decimos que tiene un punto fijo en p ∈ D si G(p) = p.
Como en el caso de una variable, cuando tenemos un sistema de ecuaciones no lineales, F (x) = 0,
podemos escribirlo en la forma G(x) = x, de varias formas, por ejemplo haciendo G(x) = x − F (x), y
transformar el problema de calcular una raı́z de F en calcular un punto fijo de G.
Pseudo-código del algoritmo de punto fijo en varias variables
entrada x0 , M, ε
x ← x0
Para k = 1, ..., M hacer
• x1 ← x, x ← G(x), e ← ∥x − x1 ∥
• salida k, x, e
• si e < ε entonces parar
fin bucle
Tenemos el siguiente resultado que nos da la convergencia del método de punto fijo en varias variables.
Observar que la convergencia depende de las propiedades de G, por tanto, la elección de esta función a
la hora de escribir el sistema que queremos resolver F (x) = 0, en la forma G(x) = x, es crucial.
Teorema 1.8 :
Sea D = {(x1 , . . . , xn )t /ai ≤ xi ≤ bi , i = 1, . . . , n} y G : D → ℜn tal que sea continua y G(x) ∈ D,

∀x ∈ D. Entonces G tiene al menos un punto fijo p ∈ D. Si además, G tiene derivadas parciales primeras
continuas y ∃K < 1 tal que
∂gi (x) K
≤ ∀x ∈ D, i, j = 1, . . . , n
∂xj n
entonces la sucesión {x(m) }∞

m=0 definida por la iteración funcional x
(m)
= G(x(m−1) ) para m ≥ 1,
partiendo de un x(0) ∈ D arbitrario, converge a dicho punto fijo p ∈ D y
Km
∥x(m) − p∥∞ ≤ ∥x(1) − x(0) ∥∞
1−K
1.3.2. Método de Newton en varias variables.
Supongamos que x(0) es un valor próximo a x solución de F (x) = 0, es decir, x(0) = x + h, haciendo
el desarrollo de Taylor,
0 = F (x) = F (x(0) − h) = F (x(0) ) − DF (x(0) )(h) + Resto
Despreciando el resto,
DF (x(0) )(h) ≈ F (x(0) )
es decir,
h ≈ DF (x(0) )−1 F (x(0) )
de donde,
x(1) = x(0) − h = x(0) − DF (x(0) )−1 F (x(0) )
será una mejor aproximación de x que x(0) , donde DF (x(0) ) viene dado por la matriz Jacobiana,
 ∂f1 (x) ∂f1 (x) ∂f1 (x)

∂x1 ∂x2 ... ∂xn
 ∂f1 (x) ∂f1 (x) ∂f1 (x) 
 ... 
J(x) =  ∂x1 ∂x2 ∂xn 
 ... ... ... ... 
∂f1 (x) ∂f1 (x) ∂f1 (x)
∂x1 ∂x2 ... ∂xn
De aquı́ se obtiene el algoritmo de Newton:
x(0) , valor cercano a x.
Calculado x(n) , obtenemos x(n+1) ,
x(n+1) = x(n) − J(x(n) )−1 F (x(n) )
Este algoritmo plantea la dificultad de tener que calcular la inversa de la matriz Jacobiana en cada
iteración, en la práctica el método se realiza en don pasos,
se resuelve el sistema J(x(n) )y = −F (x(n) )
se calcula x(n+1) = x(n) + y,
Pseudocódigo del algoritmo de Newton
entrada x(0) , M, ε
• calcular F (x) y J(x)

• resolver el sistema lineal n × n J(x)y = −F (x)
• hacer x = x + y
• salida k, x
• si ∥y∥ < ε entonces parar
• sino hacer k = k + 1
fin bucle
En el método de Newton en una variable se podı́a interpretar como un método de punto fijo donde
tratábamos de encontrar una función ρ(x) tal que la iteración funcional de g(x) = x − ρ(x)f (x) diera
convergencia cuadrática al punto fijo p de g(x), escogiendo ρ(x) = 1/f ′ (x). Para el caso n-dimensional,
buscamos una matriz n × n, A(x) = (aij (x)), donde cada componente es una función aij (x) : ℜn → ℜ tal
que,
G(x) = x − A(x)−1 F (x)
de convergencia cuadrática a la solución de F (x) = 0, siempre que, desde luego, A(x) sea no singular en
el punto fijo de G(x). Esta matriz va a ser precisamente la matriz Jacobiana.
Tenemos el siguiente resultado que nos da la convergencia del método de Newton en varias variables,
de nuevo, el valor incial debe ser próximo a la solución buscada.
Teorema 1.9 :
Sea F = (f1 , . . . , fn )t : D ∈ ℜn → ℜn con D un abierto convexo y F diferenciable con continuidad en

D. Supongamos que existe un punto x̄ ∈ D, dos constantes r, β > 0 y otra γ ≥ 0 tales que:
F(x̄) = 0
∥J(x̄)−1 ∥ ≤ β
∥J(x) − J(y)∥ ≤ γ∥x − y∥ ∀x, y ∈ B(x̄, r)
Entonces para ϵ = min{r, 1/2γβ}, y ∀x(0) ∈ B(x̄, ϵ) la sucesión generada por el método de Newton,
converge a x̄ y además
∥x(m+1) − x̄∥ ≤ βγ∥x(m) − x̄∥2
.
Capı́tulo 2
Sistemas de ecuaciones lineales
2.1. Generalidades sobre matrices y vectores
Tipos de matrices. Notación
Denotaremos Mn,n el espacio vectorial de matrices cuadradas (n filas, n columnas) con coeficientes
reales. Sea A = (aij ) ∈ Mn,n , At = (aji ) la matriz traspuesta, A−1 la matriz inversa.
A es singular si detA = 0, A es regular si detA ̸= 0.
A es simétrica si At = A, A es ortogonal si A−1 = At .
A es definido positiva si es simétrica y xt Ax > 0∀x ̸= 0.
A es semidefinido positiva si es simétrica y xt Ax ≥ 0∀x ̸= 0.

∑
A es diagonal dominante si |aii | ≥ j̸=i |aij |, i = 1, . . . , n
∑
A es estrictamente diagonal dominante si |aii | > j̸=i |aij |, i = 1, . . . , n
A y B son matrices semejantes si existe una matriz regular C tal que B = C −1 AC.
A es digonalizable si es semejante a una matriz diagonal.
A es banda p, q si aij = 0 para i ≥ j + p y j ≥ i + q (triangular superior:p=1,q=n; triangular

inferior:p=n,q=1; Hessemberg superior:p=2,q=n; Hessemberg inferior:p=n,q=2; diagonal:p=q=1;
tridiagonal:p=q=2; etc.)
∏n
Sea L una matriz triangular, entonces detL = i=1 lij .
27
28 Sistemas de ecuaciones lineales
El producto de dos matrices triangulares superiores (resp. inferiores) es una matriz triangular supe-
rior (resp. inferior), y los elementos de la diagonal son el producto de los elementos de las diagonales.
La inversa de una matriz triangular superior (resp. inferior) es una matriz triangular superior (resp.
inferior), y los elementos de la diagonal son los inversos de los elementos de la diagonal.
Propiedades
Sea A una matriz cuadrada n×n con coeficientes reales. Las siguientes afirmaciones son equivalentes.
• A−1 existe.
• detA ̸= 0 .
• El sistema lineal Ax = 0 tiene solamente la solución x = 0.
• Para cualquier vector b, el sistema lineal Ax = b tiene solución única.
• Las filas y las columnas de A son linealmente independientes.
• El rango de la matriz A es n.
(At )t = A, (A + B)t = At + B t , (αA)t = αAt .
(At )−1 = (A−1 )t , (AB)t = B t At .
detId = 1, detAt = detA, det(αA) = αn detA.
det(AB) = detAdetB, detA−1 = 1/detA si A−1 existe.
detB = det(C −1 AC) = detC −1 detAdetC = detA
Una matriz cuadrada tiene inversa ⇔ es regular.
Lema de Schur: Toda matriz es semejante a una matriz triangular superior mediante una trans-
formación de semejanza con una matriz ortogonal.
Criterio de Sylvester: Una matriz simétrica es definido positiva ⇔ todos los determinantes
principales son estrictamente positivos.
Las submatrices principales de una matriz definido positiva son definido positivas.
Los elementos diagonales de una matriz definido positiva son estrictamente positivos.
Si Q1 y Q2 son ortogonales ⇒ Q1 Q2 es ortogonal.
Si Q es ortogonal ⇒ |detQ| = 1.
Si Q es ortogonal ⇒ ∥Qx∥2 = ∥x∥2 .

2.1 Generalidades sobre matrices y vectores 29
Valores y vectores propios
Dada una matriz cuadrada A, diremos que un vector v ̸= 0 es un vector propio de A de valor propio
λ cuando Av = λv.
Los valores propios de A son las raı́ces del polinomio caracterı́stico pA (λ) = det(A − λId).
El espectro de A, Esp(A) = conjunto de valores propios de A.
El radio espectral de A, ρ(A) = máxλi ∈Esp(A) |λi |.
Los vectores propios de A y At se llaman vectores propios por la derecha y vectores propios por la
izquierda de A (respect.).
Las matrices A y At tienen los mismos valores propios. Los vectores propios por la derecha son
ortogonales a los vectores propios por la izquierda de valores propios diferentes.
A es regular ⇔ tiene todos los valores propios diferentes de cero. Los valores propios de A−1 son
los inversos de los valores propios de A. v es un vector propio de valor propio λ de A ⇔ v es un
vector propio de valor propio 1/λ de A−1 .
Los valores propios de una matriz simétrica son reales.
Una matriz simétrica es definido positiva ⇔todos los valores propios son estrictamente positivos.
Los espectros de dos matrices semejantes son iguales. Si v es un vector propio de valor propio λ de
A y B = CAC −1 , entonces C −1 v es vector propio de valor propio λ de B.
Se llaman valores singulares µ de A a las raı́ces cuadradas positivas de los valores propios de At A.
Normas vectoriales
Sea E un espacio vectorial, una norma en E es una aplicación

∥ ∥: E → ℜ+
x → ∥x∥
que cumple:
• ∥x∥ = 0 ⇔ x = 0,
• ∥cx∥ = |c|∥x∥∀ escalar c, ∀x ∈ E,
• ∥x + y∥ ≤ ∥x∥ + ∥y∥∀x, y ∈ E.
Las normas vectoriales son las que están definidas sobre espacios vectoriales de la forma E = ℜn o
E = Cn .
Normas Hölder:
n
∑
∥x∥p = ( |xi |p )1/p , p ≥ 1
i=1
Norma suma de módulos:

n
∑
∥x∥1 = |xi |
i=1
Norma euclı́dea:
n
∑
∥x∥2 = ( |xi |2 )1/2
i=1
Norma infinito:
∥x∥∞ = máx |xi |
i=1÷n
Normas matriciales
Una norma matricial es una norma en el espacio de las matrices cuadradas Mn,n que sea multiplicativa,
es decir, que cumpla
∥AB∥ ≤ ∥A∥∥B∥ ∀A, B ∈ Mn,n .
Una norma matricial ∥ ∥ es consistente con una norma vectorial (que denotaremos igual) si y sólo si
∥Ax∥ ≤ ∥A∥∥x∥ ∀A ∈ Mn,n , ∀x ∈ Rn .
Dada una norma vectorial ∥ ∥, siempre se puede definir una norma matricial consistente con ella, llamada
norma matricial subordinada, mediante
∥Ax∥
∥A∥ = máx = máx ∥Ax∥
x̸=0 ∥x∥ ∥x∥=1
por ejemplo, las normas subordinadas a las normas Hölder,
∥A∥p = máx ∥Ax∥p

∥x∥p =1
Se verifican las siguientes propiedades,
Norma 1:
n
∑
∥A∥1 = máx |aij |
1≤j≤n
i=1
Norma euclı́dea:
√
∥A∥2 = ρ(At A) = µmax
Norma infinito:
n
∑
∥A∥∞ = máx |aij |
1≤i≤n
j=1
Si A es simétrica ⇒ ρ(A) = ∥A∥2

2.1 Generalidades sobre matrices y vectores 31
Propiedades
ρ(A) ≤ ∥A∥ para una norma matricial cualquiera.
∥A∥ ≤ ρ(A) + ε para al menos una norma matricial inducida.
Las siguientes condiciones son equivalentes:
• lı́mk→∞ Ak = 0
• lı́mk→∞ Ak x = 0 ∀x ∈ Rn
• ρ(A) < 1
• ∥A∥ < 1 para alguna norma matricial inducida.
Condicionamiento de una matriz
Sea A una matriz cuadrada regular, y ∥ ∥ una norma matricial, se define el condicionamiento de A
asociado a la norma ∥ ∥, como
cond(A) = ∥A∥∥A−1 ∥
cond(αA) = cond(A), ∀α ̸= 0, ∀A ∈ Mn,n regular.
cond(A) ≥ 1 si se calcula con una norma inducida.

µmax
cond2 (A) = µmin , donde µmax , µmin son los valores singulares de A más grande y más pequeño.
Una matriz está bien condicionada si su número de condicionamiento es próximo a 1.
Las matrices mejor condicionadas son las matrices ortogonales.
En este capı́tulo consideraremos el problema de resolución de sistemas de ecuaciones lineales, unos de

los problemas numéricos que en la práctica aparece con mayor frecuencia, por sı́ mismo o como parte de
un problema mayor.
Los métodos que utilizaremos se clasifican en dos grupos: métodos directos y métodos iterativos. En
teorı́a, los métodos directos permiten calcular la solución exacta con un número finito de operaciones
aritméticas, en la práctica debido a la finitud de las cifras con que podemos representar los números
reales en un ordenador, la acumulación de errores de redondeo produce soluciones aproximadas. En los
métodos iterativos la solución se define como lı́mite de una sucesión infinita de vectores, en la práctica se
puede obtener una solución aproximada fijando un número finito de iteraciones.
Un sistema de n ecuaciones lineales con n incógnitas se puede escribir,



 a11 x1 + a12 x2 + . . . + a1n xn = b1

a21 x1 + a22 x2 + . . . + a2n xn = b2

 ...

an1 x1 + an2 x2 + . . . + ann xn = bn
donde los coeficientes aij (i, j = 1, 2, . . . , n) y los términos independientes bi (i = 1, 2, . . . , n) son

constantes dadas. El problema es determinar las incógnitas del sistema x1 , x2 , . . . , xn de forma que se
satisfagan las n ecuaciones simultaneamente.
El sistema en forma matricial se escribe,
Ax = b,
siendo A = (aij ) la matriz de coeficientes, b = (bi ) el vector de términos de independientes, y x = (xi ) el

vector de incógnitas. A partir de ahora, nos centraremos en el caso de sistemas determinados, es decir,
con solución única.
2.2. Métodos directos de resolución de sistemas de ecuaciones

lineales
2.2.1. Matrices triangulares
Comenzaremos presentando dos algoritmos de resolución muy sencillos para sistemas con matriz de
coeficientes triangular. Tengamos en cuenta que para que un sistema de este tipo sea determinado los
coeficientes de la diagonal principal deben ser no nulos.
Consideremos el siguiente sistema de n ecuaciones con n incógnitas, cuya matriz de coeficientes es

triangular superior invertible,


 a11 x1 + a12 x2 + . . . . . . + a1n xn = b1


 a22 x2 + . . . . . . + a2n xn = b2
......



 a n−1,n−1 x n−1 + an−1,n xn = bn−1

ann xn = bn
Teniendo en cuenta que aii ̸= 0, para i = 1, . . . , n, de la última ecuación podemos despejar,
xn = bn /ann ,
sustituyendo este valor en la penúltima ecuación, tenemos
xn−1 = (bn−1 − an−1,n xn )/an−1,n−1 ,
y ası́ sucesivamente hasta llegar a
x1 = (b1 − a1,2 x2 − . . . − a1n xn )/a1,1 .
Este método se llama método de sustitución inversa.

2.2 Métodos directos de resolución de sistemas de ecuaciones lineales 33
Algoritmo 2.1 : Método de sustitución inversa.
Para resolver Ax = b siendo A triangular superior invertible
xn = bn /ann ,
para i desde ∑n − 1 hasta 1
n
xi = (bi − j=i+1 aij xj )/aii ,
Algoritmo 2.2 : Método de sustitución directa.
Para resolver Ax = b siendo A triangular inferior invertible
x1 = b1 /a11 ,
para i desde 2 hasta n
∑i−1
xi = (bi − j=1 aij xj )/aii ,
Es importante tener una idea del costo computacional de estos algoritmos para poder hacer compa-
raciones entre distintos métodos. Para algoritmos de álgebra lineal, una forma de estimar este coste es
contar el número de operaciones algebraicas del algoritmo. En casi todos los ordenadores las sumas y
restas llevan aproximadamente el mismo tiempo de ejecución, al igual que productos y divisiones, ası́,
normalmente para contar el número de operaciones de un algoritmo se cuentan el número de sumas/restas
y el número de productos/divisiones.
En los algoritmos anteriores, el número de operaciones realizadas es,
n−1
∑ (n − 1)n n2 n
i= = − sumas/restas
i=1
2 2 2
n
∑ n(n + 1) n2 n
i= = + productos/divisiones
i=1
2 2 2
2.2.2. Eliminación gaussiana
Entre los métodos directos de resolución de sistemas de ecuaciones lineales es el más popular. También
se usa para calcular determinantes e inversas de matrices.
La idea básica consiste en transforma (reducir) el sistema inicial en uno equivalente (misma solución)
cuya matriz de coeficientes sea triangular superior. El nuevo sistema se puede resolver fácilmente mediante
el método de sustitución inversa. La reducción se realiza mediante transformaciones elementales sobre
las ecuaciones del sistema: permutar dos ecuaciones y sustituir una ecuación por su suma con otra
multiplicada por una constante.
Veamos como es posible transformar un sistema determinado Ax = b de n ecuaciones con n incógnitas,

realizando un número finito de transformaciones elementales, en un sistema equivalente con matriz de
coeficientes triangular superior. La reducción se realiza en n − 1 pasos.
e donde cada fila representa

Para simplificar la notación usaremos la matriz ampliada del sistema A,
una ecuación.
   (1) (1) (1) (1)

a11 a12 ... a1n b1 a11 a12 ... a1n b1
 a21 a22 ... a2n b2   (1)
a21
(1)
a22 ...
(1)
a2n
(1)
b2

e1 = 
e=A   
A  .. .. .. .. .. = .. .. .. .. .. 
 . . . . .  
 . . . . .


an1 an2 ... ann bn an1
(1) (1)
an2 ...
(1)
ann
(1)
bn
(1)
Primer paso de la reducción: si a11 ̸= 0 cada fila i, (i = 2, 3, . . . , n) se sustituye por ella misma menos
(1) (1)
la primera fila multiplicada por li1 = ai1 /a11 , obteniéndose la siguiente matriz ampliada,
 (1) (1) (1) (1)

a11 a12 ... a1n b1
 (2) (2) (2) 
 0 a22 ... a2n b2 
e
A2 =  .. .. .. .. 
 .. 
 . . . . . 
(2) (2) (2)
0 an2 ... ann bn
siendo
(2) (1) (1)
aij = aij − li1 a1j (i = 2, . . . , n; j = 2, . . . , n),
(2) (1) (1)
bi = bi − li1 b1 (i = 2, . . . , n).
(2)
Segundo paso de la reducción: si a22 ̸= 0 cada fila i, (i = 3, 4, . . . , n) se sustituye por ella misma menos
(2) (2)
la segunda fila multiplicada por li2 = ai2 /a22 , obteniéndose la siguiente matriz ampliada,
 (1) (1) (1) (1) (1) 
a11 a12 a13 ... a1n b1
 (2) (2) (2) (2) 
 0 a22 a23 ... a2n b2 
 (3) (3) (3) 
e
A3 =  0 0 a33 ... a3n b3 
 
 .. .. .. .. .. 
 . . . . . 
(3) (3) (3)
0 0 an2 ... ann bn
siendo
(3) (2) (2)
aij = aij − li2 a2j (i = 3, . . . , n; j = 3, . . . , n),
(3) (2) (2)
bi = bi − li2 b2 (i = 3, . . . , n).
y ası́ sucesivamente hasta llegar al paso n − 1 en el que se obtiene la matriz ampliada,

 (1) (1) (1) (1)

a11 a12 ... a1n b1
 (2) (2) (2) 
 0 a22 ... a2n b2 
en = 
A .. .. .. .. 
 .. 
 . . . . . 
(n) (n)
0 0 ... ann bn
de un sistema equivalente al primero triangular superior que se resuelve por sustitución inversa.
(1) (2) (n)

Los números a11 , a22 , . . . , ann se llaman pivotes. Los números lij se llaman multiplicadores.
(k)
Si en el paso correspondiente se tiene que akk = 0, por ser el sistema determinado, siempre podemos
(k)
encontrar ajk ̸= 0, (j ≥ k), permutar las filas j y k y continuar el proceso, esto se denomina pivotaje.
El proceso de eliminación gaussiana se puede realizar sin pivotaje, si y sólo si todos los determinantes
principales de la matriz de coeficientes A son no nulos.
Ejemplo 2.1 Resolver el siguiente sistema de ecuaciones lineales mediante el método de eliminación
gaussiana,

 x1


− x2 + 2x3 − x4 = −8
2x1 − 2x2 + 3x3 − 3x4 = −20
 x1
 + x2 + x3 = −2

x1 − x2 + 4x3 + 3x4 = 4
Algoritmo 2.3 Método de eliminación gaussiana
Para resolver el sistema determinado Ax = b siendo A de orden n
para k desde 1 hasta n

para i desde k + 1 hasta n
l = aik /akk
para j desde k + 1 hasta n
aij = aij − lakj ,
bi = bi − lbk ,
xn = bn /ann ,
para i desde ∑n − 1 hasta 1
n
xi = (bi − j=i+1 aij xj )/aii ,
El número de operaciones aritméticas realizadas en el proceso de reducción es,
n−1
∑ n−1
∑ 1 1
(n − k)(n − k + 1) = m(m + 1) = (n − 1)n(2n − 1) + n(n − 1) sumas/restas
m=1
6 2
k=1
n−1
∑ n−1
∑ 1
(n − k)(n − k + 2) = m(m + 2) = (n − 1)n(2n − 1) + n(n − 1) productos/divisiones
m=1
6
k=1
Junto con el coste operacional del proceso de sustitución inversa, tenemos un total de,
n3 n2 5n
+ − sumas/restas
3 2 6
n3 n
+ n2 − productos/divisiones
3 3
En la siguiente tabla podemos ver como aumenta el coste operacional con el tamaño del sistema.
n +/− ∗/÷
3 17 11
10 430 375
50 44150 42875
100 343300 338250
Una variante del método de eliminación gaussiana es el Método de Gauss-Jordan. Consiste en

reducir a 0 no sólo los elementos por debajo de la diagonal, sino también por encima, es decir, en el paso
k de la reducción, se transforman en 0 todos los elementos de la columna k fuera de la diagonal, obteniendo
al final un sistema equivalente con matriz de coeficientes diagonal, que se resuelva fácilmente dividiendo
los términos independientes por el correspondiente elemento de la diagonal. El coste operacional de este
método es superior al de eliminación gaussiana.
n3 n
− sumas/restas
2 2
n3 n
+ n2 − productos/divisiones
2 2
2.2.3. Técnicas de pivotaje
Si en algún paso del proceso de eliminación gaussiana el pivote es no nulo pero muy pequeño, aunque
podemos aplicar el método, éste es muy inestable numéricamente en el sentido de que los errores de la
solución, propagados a partir de los errores de los datos, pueden aumentar notablemente. Conviene en
estos casos modificar el método de eliminación gaussiana con alguna técnica de pivotaje, por ejemplo:
Pivotaje parcial o maximal por columnas
(k)
En el paso k-ésimo, en lugar de tomar como pivote akk , se toma el elemento de máximo valor absoluto
(k) (k) (k)
entre los elementos aik para i = k, . . . , n, es decir, si |apk | = máxk≤i≤n |aik |, se permutan las filas p y
k.
También es recomendable utilizar alguna técnica de pivotaje cuando las magnitudes entre las distintas
filas (distintas ecuaciones) es muy diferente, en este caso se recomienda:
Pivotaje escalado
:
Calculamos para cada fila un factor escalar si = máx1≤j≤n |aij |, es decir, elegimos en cada fila el
elemento de máximo valor absoluto. Si detA ̸= 0 entonces si ̸= 0 ∀i = 1, . . . , n. En cada paso k,
permutamos la fila k por la p si
(k) (k)
|apk | |aik |
= máx
sp k≤i≤n si
De este modo conseguimos eliminar la diferencia de magnitud relativa entre los elementos de distintas
filas. Los factores escalares se eligen una sola vez al principio del proceso y no en cada paso, pues esto
elevarı́a en exceso el coste operacional.
2.2.4. Factorización LU
El método de factorización LU produce la factorización de la matriz del sistema A en el producto de

una matriz triangular inferior con unos en la diagonal L, y otra triangular superior U . Esta factorización
A = LU se puede realizar de manera única siempre que todos los determinantes principales de A sean no
nulos, la misma condición que nos permite realizar el proceso de eliminación gaussiana sin pivotaje. En
caso de que no se cumpla esta condicón y siempre que A sea regular, será posible permutar las ecuaciones
de manera que la nueva matriz P A admita una tal factorización, P A = LU .
Una vez realizada la factorización LU , la solución del sistema Ax = b, (o bien P Ax = P b si ha sido

necesario hacer alguna permutación), se realiza en dos pasos, resolviendo sucesivamente dos sistemas
triangulares,
{
1) Ly = b ⇒ y
Ax = LU x = b ⇒
2) Ux = y ⇒ x
{
1) Ly = P b ⇒ y
P Ax = LU x = P b ⇒
2) Ux = y ⇒ x
El modo de calcular los coeficientes de las matrices L y U nos lo da el proceso de eliminación gaussiana
sin más que interpretarlo en forma matricial.
Permutar las filas p y k de una matriz A consiste en multiplicar dicha matriz por la matriz de
permutación P correspondiente, que no es mas que la matriz identidad en la que se han permutado
las filas p y k. Observar que las matrices de permutación verifican P −1 = P . Observar también que el
determinate de una matriz a la que se han permutado dos filas es el mismo cambiado de signo.
Ası́, si en el primer paso del proceso de eliminación gaussiana tenemos que permutar las filas 1 y p,
estamos multiplicando el sistema Ax = b por la matriz de permutación P1 , que es la matriz identidad
donde se han permutado las filas 1 y p. Después, al hacer ceros por debajo del pivote estamos multiplicando
por una matriz L1 que se obtiene a partir de los multiplicadores, a saber,
P1 A1 = L1 A2
donde,
 
1 0 0 ... 0
 l21 1 0 ... 0 
 
 l31 0 1 ... 0 
L1 =  
 .. .. .. .. .. 
 . . . . . 
l31 0 0 ... 1
Entonces, multiplicando por P1 , tenemos,
A = A1 = P1 L1 A2
En el segundo paso,
P2 A2 = L2 A3
donde,
 
1 0 0 ... 0
 0 1 0 ... 0 
 
 0 l32 1 ... 0 
L2 =  
 .. .. .. . . .. 
 . . . . . 
0 ln2 0 ... 1
de modo que,
A = A1 = P1 L1 A2 = P1 L1 P2 L2 A3
Al llegar al último paso,

A = P1 L1 P2 L2 . . . Pn−1 Ln−1 An
donde An = U es la matriz triangular superior buscada. Si no hiciera falta ninguna permutación
tendrı́amos A = L1 L2 . . . Ln−1 U = LU donde L = L1 L2 . . . Ln−1 es la matriz triangular inferior con
unos en la diagonal construida con los multiplicadores lij ,
 
1 0 0 ... 0
 l21 1 0 ... 0 
 
 l31 l32 1 . . . 0 
L=  
 .. .. .. .. .. 
 . . . . . 
ln1 ln2 ln3 . . . 1
En el caso de ser necesarias las permutaciones, como el producto de matrices no es conmutativo,

tenemos,
A = P1 L1 P2 L2 . . . Pn−1 Ln−1 An = P1 P2 . . . Pn−1 L′1 L′2 . . . L′n−1 U = P LU
donde,
P = P1 P2 . . . Pn−1
L = L′1 L′2 . . . L′n−1
L′k = Pn−1 Pn−2 . . . Pk+1 Lk Pk+1 . . . Pn−2 Pn−1
de modo que multiplicando por la matriz de permutación P tenemos, P A = LU .
Ejemplo 2.2 Resolver el sistema de ecuaciones lineales del ejemplo anterior mediante el método de
factorización LU.
2.2.5. Matrices especiales: factorización LDLt , Cholesky
Si la matriz A es simétrica y admite una factorización A = LU , podemos escribir U = DR con

D diagonal y R triangular superior con unos en la diagonal. Al ser A simétrica, A = At , por lo que
LDR = Rt DLt , y como la factorización es única, entonces R = Lt , de donde se puede escribir,
A = LDLt
Si además A es definido positiva, entonces los elementos diagonales de D son positivos y podemos
calcular su raı́z cuadrada. Sea D1/2 la matriz diagonal cuyos coeficientes son las raı́ces cuadradas de los
de D. Entonces,
A = LDLt = LD1/2 D1/2 Lt = LLt
donde L = LD1/2 que es la factorización de Cholesky para matrices simétricas definido positivas.
Ejemplo 2.3 Calcular la factorización LDLt y la de Cholesky para la siguiente matriz,

 
13 11 11
 11 13 11 
11 11 13
2.2.6. Aplicaciones
Cálculo del determinante de una matriz.
Si en una matriz se permutan dos filas o dos columnas, el valor absoluto del determinante no varı́a,
pero si cambia de signo. Si en una matriz a los elementos de una fila o columna se les suman los de otra
multiplicados por un escalar, el determinante no varı́a. Estas son las transformaciones que se realizan en
el proceso de eliminación gaussiana, por tanto,
n
∏
detA = εdetU = ε uii
i=1
ya que U es triangular, siendo ε = 1 si el número de permutaciones de filas realizadas es par, y ε = −1

si es impar.
3
El número de operaciones realizadas al calcular el determinate de este modo es del orden de ϑ( n3 ),
mucho menor que el número de operaciones realizadas al usar la regla de Laplace para el cálculo del
determinante,
∑
detA = (−1)signσ a1,σ(1) . . . an,σ(n)
σ
que son n! − 1 sumas/restas, y n!(n − 1) productos/divisiones.

Cálculo de la inversa de una matriz.
Si A es una matriz regular, y X es su inversa, entonces su producto es la matriz identidad AX = Id.

Separando las columnas de esta ecuación, si x(k) es la k-ésima columna de X, y e(k) la k-ésima columna
de la matriz identidad, esto es, el k-ésimo vector de la base canónica, tenemos
Ax(k) = e(k)
Entonces el cálculo de X se reduce a la resolución de n sistemas de ecuaciones lineales de dimensión
n con la misma matriz de coeficientes A y distintos términos independientes, las columnas de la matriz
identidad. Entonces basta con aplicar el método de eliminación gaussiana la la siguiente matriz ampliada,
 
a11 ... a1n 1 ... 0
 .. .. .. .. . . .. 
(A|Id) =  . . . . . . 
an1 ... ann 0 ... 1
Ejemplo 2.4 Invertir las siguiente matrices,

 
1 2 −1
 2 1 0 
−1 1 2
 
1 2 3 4
 1 4 9 16 
 
 1 8 27 64 
1 16 81 256
2.3. Métodos iterativos de resolución de sistemas de ecuaciones

lineales
Los métodos directos son eficaces para sistemas de tamaño moderado, por ejemplo n ≈ 1000 o en el
caso de matrices huecas n ≈ 5000 , 10000. Para valores significativamente mayores los métodos directos
pierden eficacia, no sólo porque el número de operaciones necesario crece desmesuradamente sino también
porque la acumulación de errores de redondeo puede desvirtuar el resultado.
En el caso de grandes sistemas de ecuaciones, los llamados métodos iterativos, resultan más conve-
nientes. De forma genérica: Para resolver un sistema Ax = b, se transforma en otro equivalente (es decir,
con la misma solución) que tenga la forma
x = Bx + c
expresión que sugiere el siguiente método iterativo
{
x(0) , arbitrario
x(k+1) = Bx(k) + c
2.3 Métodos iterativos de resolución de sistemas de ecuaciones lineales 41
Diremos que el método es convergente, si
lı́m x(k) = x
k→∞
cualquiera que sea el valor inicial x(0) elegido.
Teorema 2.1 El método iterativo anterior es convergente si ρ(B) < 1, o de forma equivalente si ∥B∥ < 1
para al menos una norma matricial (que podemos elegir subordinada).
Demostración:
x = Bx + c
(k+1)
x = Bx(k) + c
restando
x(k+1) − x = B(x(k) − x)
Llamando e(0) = x(0) − x al error inicial y e(k) = x(k) − x al error en la iteración k, resulta e(k+1) = Be(k)
y también
e(k) = Be(k−1) = ... = B k e(0)
de donde
∥e(k) ∥ = ∥B k e(0) ∥ ≤ ||B k ∥∥e(0) ∥ ≤ ∥B∥k ∥e(0) ∥
Si ∥B∥ < 1 entonces

lı́m ∥e(k) ∥ = 0
k→∞
es decir,
lı́m x(k) = x
k→∞
Hemos obtenido,
∥e(k) ∥ ≤ ∥B k ∥∥e(0) ∥
esto se puede interpretar como que en cada una de las k primeras iteraciones el error se ha reducido en
un factor de ∥B k ∥1/k y, en consecuencia, se puede estimar que para que el error se reduzca en un factor
de 10−m se deben realizar N iteraciones cumpliéndose,
m
(∥B k ∥1/k )N ≤ 10−m , o bien N ≥ .
−log10 (∥B k ∥1/k ))
Al número −log10 (∥B k ∥1/k )) se le llama velocidad media de convergencia en k iteraciones.
Se puede demostrar que ρ(B) = lı́mk→∞ ∥B k ∥1/k . Al número −log10 (ρ(B)) se le llama velocidad
asintótica de convergencia.
Teorema 2.2 El método iterativo anterior se verifica,
∥B∥k
∥e(k) ∥ = ∥x(k) − x∥ ≤ ∥x(1) − x(0) ∥
1 − ∥B∥
Demostración: En efecto, tomando normas en la siguiente expresión,
x(k) − x = x(k) − x(k+1) + x(k+1) − x

resulta
∥x(k) − x∥ ≤ ∥x(k) − x(k+1) ∥ + ∥x(k+1) − x∥ ≤ ∥x(k+1) − x(k) ∥ + ∥B∥∥x(k) − x∥
es decir,
(1 − ∥B∥)∥x(k) − x∥ ≤ ∥x(k+1) − x(k) ∥ ≤ ∥B∥∥x(k) − x(k−1) ∥ ≤ . . . ≤ ∥B∥k ∥x(1) − x(0) ∥
de donde se deduce la desigualdad buscada.
2.3.1. Método de Jacobi
Supongamos que queremos resolver el sistema
a11 x1 + ... +a1n xn = b1

a21 x1 + ... +a2n xn = b2
...
an1 x1 + ... +ann xn = bn
que podemos escribir

∑
a11 x1 = b1 − a1j xj
j̸=1
∑
a22 x2 = b2 − a2j xj
j̸=2
...
∑
ann xn = bn − anj xj
j̸=n
El algoritmo de Jacobi se escribe: Dado x(0) ∈ Rn arbitrario, una vez calculada una aproximación x(k) ,
calculamos x(k+1) de la manera siguiente,
∑ (k)
(k+1) b1 − j̸=1 a1j xj
x1 =
a11
∑ (k)
(k+1) b2 − j̸=2 a2j xj
x2 =
a22
...
∑ (k)
bn − j̸=n anj xj
x(k+1)
n =
ann
Este método está definido sólo si aii ̸= 0 para i = 1, . . . , n. La ecuación i-ésima
(k+1) 1 ( ∑ (k)
)
xi = bi − aij xj
aii
j̸=i
(k)
se puede escribir también, restando en los dos miembros xi ası́:
1 ( ) r(k)
∑n
(k+1) (k) (k)
xi − xi = bi − aij xj = i
aii j=1
aii
donde hemos designado mediante r(k) al vector residuo
r(k) = b − Ax(k)
correspondiente al valor x(k) .
Vamos a escribir el método anterior en forma matricial. Pondremos
A=D−E−F
donde
D es la parte diagonal de A, Dii = aii , i = 1, . . . , n
−E es la parte estrictamente triangular inferior

{
(−E)ij = aij i > j
(−E)ij = 0 i≤j
−F es la parte estrictamente triangular superior

{
(−F )ij = aij i < j
(−F )ij = 0 i≥j
Entonces la iteración de Jacobi se escribe
x(k+1) = D−1 (E + F )x(k) + D−1 b
o bien,
x(k+1) = (I − D−1 A)x(k) + D−1 b
es pues de la forma general x = Bx + c, con B = J = I − D−1 A y c = D−1 b.

2.3.2. Método de Gauss-Seidel
Si observamos con atención la expresión general de una iteración del algoritmo de Jacobi, podemos ver
(k+1) (k+1) (k+1) (k+1)
que si procedemos en el orden natural, i = 1, 2, . . . , n, al calcular xi , los valores x1 , x2 , . . . , xi−1
ya los hemos obtenido. Si el método es convergente, tenemos la esperanza que estos i − 1 valores estén
(k) (k) (k)
más cerca de la solución que los anteriores x1 , x2 , . . . , xi−1 . Por lo tanto podemos utilizarlos en lugar
(k+1)
de estos en la expresión que sirve para calcular xi , quedando
1 ( )
i−1
∑ ∑n
(k+1) (k+1) (k)
xi = bi − aij xj − aij xj
aii j=1 j=i+1
Obtenemos ası́ el llamado método de Gauss-Seidel, que podemos escribir de la forma
i
∑ n
∑
(k+1) (k)
aij xj = bi − aij xj
j=1 j=i+1
o en forma matricial
(D − E)x(k+1) = b + F x(k)
y también
x(k+1) = (D − E)−1 F x(k) + (D − E)−1 b
que es de la forma general x = Bx + c, con B = L1 = (D − E)−1 F y c = (D − E)−1 b.
En el método de Gauss-Seidel no aparece el residuo explı́citamente, sino
i−1
∑ n
∑
(k+1) (k)
rei = bi − aij xj − aij xj
j=1 j=i
entonces,
(k+1) (k) rei
xi = xi +
aii
2.3.3. Métodos de relajación
Se pueden generalizar los dos métodos anteriores de Jacobi y Gauss-Seidel, introduciendo un parámetro
(k) (k+1) (k)
ω > 0. Sea xi ya calculado y x̂i obtenido a partir de xi por uno de los dos métodos anteriores. Se
define entonces la combinación lineal
(k+1) (k+1) (k)

xi = ωx̂i + (1 − ω)xi
Si el método de partida es el de Jacobi, obtenemos para i = 1, . . . , n
(k+1) ω ∑ (k) (k)

xi = (bi − aij xj ) + (1 − ω)xi
aii
j̸=i
o bien multiplicando por aii
(k+1)
∑ (k) (k)
aii xi = ω(bi − aij xj ) + (1 − ω)aii xi
j̸=i
y con notación matricial

Dx(k+1) = (1 − ω)Dx(k) + ωb + ω(E + F )x(k)
y también
x(k+1) = (I − ωD−1 A)x(k) + ωD−1 b
En el caso del método de Gauss-Seidel, el correspondiente método de relajación se llama S.O.R.

(Succesive Over Relaxation) y se escribe
i−1
∑ ∑n
(k+1) ω (k+1) (k) (k)
xi = (bi − aij xj − aij xj ) + (1 − ω)xi
aii j=1 j=i+1
y con notación matricial

(D − ωE)x(k+1) = ωb + ((1 − ω)D + ωF )x(k)
es decir
x(k+1) = (D − ωE)−1 ωb + (D − ωE)−1 ((1 − ω)D + ωF )x(k)
que es de la forma general x = Bx + c, con B = Lω = (D − ωE)−1 ((1 − ω)D + ωF ) y c = (D − ωE)−1 ωb
2.3.4. Control de parada de las iteraciones
Designemos mediante r(k) al vector residuo correspondiente a la iteración k-ésima, es decir,
r(k) = b − Ax(k)
Un posible control de parada consiste en parar en la k-ésima iteración si
∥r(k) ∥
≤ε
∥b∥
para ε elegido convenientemente pequeño.

Esta relación implica que el error e(k) = x − x(k) verifica
∥e(k) ∥
≤ εcond(A)
∥x∥
siendo x la solución exacta de Ax = b. En efecto, como
∥e(k) ∥ = ∥A−1 r(k) ∥ ≤ ∥A−1 ∥∥r(k) ∥
entonces
∥e(k) ∥ ≤ ε∥A−1 ∥.∥b∥ ≤ ε∥A−1 ∥.∥Ax∥ ≤ εcond(A)∥x∥
y de ahı́ la afirmación realizada.
En los métodos de Gauss-Seidel y S.O.R. no aparece el residuo explı́citamente, sino
i−1
∑ n
∑
(k+1) (k)
rek = bi − aij xj − aij xj ,
j=1 j=i
entonces el criterio de podrı́a ser,

r(k) ||
||e
≤ε
||b||
lo que evita cálculos suplementarios.
Otro posible criterio de parada consiste en interrumpir las iteraciones cuando,
∥x(k) − x(k−1) ∥
≤ε
∥x(k) ∥
que es un control cómodo desde el punto de vista del cálculo. Presenta sin embargo el inconveniente que
podrı́a darse en ciertos casos en los que se verificase el control sin que x(k) estuviese cerca de la solución
x. Por ejemplo si para algún k resulta x(k) = 0 sin ser ésta la solución buscada.
2.3.5. Resultados de convergencia
Los métodos anteriores son de la forma general
x(k+1) = Bx(k) + c
La condición necesaria y suficiente de convergencia es
ρ(B) < 1
Para el método de Jacobi

B = J = D−1 (E + F ) = Id − D−1 A
Para el método de Gauss-Seidel

B = L1 = (D − E)−1 F
Para el método S.O.R.
B = Lω = (D − ωE)−1 ((1 − ω)D + ωF )
Estas matrices se pueden expresar en función de L = D−1 E y de U = D−1 F que son respectivamente
dos matrices triangulares inferior y superior con diagonal nula
   
0 0 ... 0 0 − aa12
11
... − aa1n
11
 − aa21 0 ... 0   0 0 ... − aa2n 
L=

22 


U = 22 

... ... ... ...
− aann
n1
− aann
n2
... 0 0 0 ... 0
Teniendo en cuenta que D−1 A = D−1 (D − E − F ) = I − L − U , podemos escribir fácilmente, para el

método de Jacobi
J = D−1 (E + F ) = D−1 E + D−1 F = L + U
para el método de Gauss-Seidel
L1 = (D − E)−1 F = (I − D−1 E)D−1 F = (I − L)−1 U
y para el método S.O.R.
Lω = (D − ωE)−1 ((1 − ω)D + ωF ) = (I − ωD−1 E)−1 ((1 − ω)I + ωD−1 F )

= (I − ωL)−1 ((1 − ω)I + ωU )
Vamos a ver una condición necesaria para que el radio espectral de la matriz del método S.O.R. sea
menor que la unidad.
Teorema 2.3 Para toda matriz A, el radio espectral de la matriz del método de relajación S.O.R. es
superior o igual a |ω − 1| en consecuencia una condición necesaria para que el método sea convergente es
0 < ω < 2.
Demostración: Los valores propios de la matriz Lω del método de relajación verifican la relación
det( 1−ω
ω D + F) ( 1−ω n
ω ) Πaii
Πni=1 λi (Lω ) = det(Lω ) = = 1 n = (1 − ω)n
det( D
ω − E) ( ω ) Πaii
y como por otra parte

ρ(Lω ) ≥ |λi |
lo que implica
ρn (Lω ) ≥ Πni=1 |λi | = |ω − 1|n
resulta finalmente
ρ(Lω ) ≥ |ω − 1|
Corolario 2.1 Para toda matriz A, una condición necesaria de convergencia del método de S.O.R. es
0<ω<2
Matrices diagonal dominantes
Definición 2.1 Una matriz A cuadrada de orden n se dice que es estrictamente diagonal dominante si
n
∑
|aii | > |aij | para i = 1, . . . , n
j=1
j̸=i
Teorema 2.4 Si A es una matriz de orden n estrictamente diagonal dominante entonces es no singular.
Demostración: Consideremos el sistema de ecuaciones
Ax = 0
y veamos que tiene como única solución x = 0.
Por reducción al absurdo, supongamos que x = [x1 , . . . , xn ]t es una solución distinta de cero. En este
caso para algún k, 0 < |xk | = máx1≤j≤n |xj |
∑n
Como j=1 aij xj = 0 para todo i = 1, . . . , n, tomando i = k resulta
n
∑
akk xk = − akj xj
j=1
j̸=k
de donde
n
∑
|akk ||xk | ≤ |akj ||xj |
j=1
j̸=k
es decir
n
∑ n
|xj | ∑
|akk | ≤ |akj | ≤ |akj |
j=1
|xk | j=1
j̸=k j̸=k
en contradición con la propiedad de A de ser estrictamente diagonal dominante.
Teorema 2.5 Sea A, matriz cuadrada de orden n estrictamente diagonal dominante. Entonces el método
de Jacobi para resolver un sistema de ecuaciones lineales asociado a dicha matriz es convergente.
Demostración: La matriz de iteración para el método de Jacobi es J = D−1 (E +F ) = L+U . Vamos

a demostrar que ||J||∞ < 1. En efecto,
 
0 − aa11
12
... ... − aa1n
11
 − a21 0 ... ... − aa2n 
J =L+U =

a22 22 

... ... ...
− aann
n1
− aann
n2
... ... 0
de donde
n ∑
∑ aij j̸=i |aij |
||J||∞ = máx | | = máx <1
1≤i≤n
j=1
aii 1≤i≤n |aii |
j̸=i
pues A es estrictamente diagonal dominante.
Teorema 2.6 Sea A una matriz estrictamente diagonal dominante, entonces el método de Gauss-Seidel
para resolver un sistema de ecuaciones lineales asociado a dicha matriz es convergente.
Demostración: La matriz asociada a la iteración de Gauss-Seidel es
L1 = (D − E)−1 F = (I − L)−1 U
Para determinar el radio espectral de L1 , calcularemos primero los valores propios, es decir, las raı́ces del
polinomio caracterı́stico
p(λ) = det(λI − L1 ) = 0
Observando que det(I − L) = 1 resulta
p(λ) = det(I − L)det(λI − L1 )
= det(I − L)det(λI − (I − L)−1 U )

= det(λ(I − L) − U )
U
= det(λ(I − L − ))
λ
U
= λn det(I − L − )
λ
U
de donde p(λ) = 0 si λ = 0 o bien si det(I − L − λ) = 0.
Queremos demostrar que todas las raı́ces de p(λ) = 0 verifican |λ| < 1. Supongamos por reducción
U
al absurdo que existe al menos una raiz λ tal que |λ| ≥ 1. Entonces por una parte det(I − l − λ) =0
y por otra parte como A = D − E − F es estrictamente diagonal dominante, también lo es I − L − U
U U
y lo será también I − L − λ si |λ| ≥ 1. Por lo tanto I − L − λ es no singular en contradicción con
U
det(I − L − λ) = 0.
Matrices simétricas y definidas postivas
En esta sección vamos a ver algunos resultados interesantes que relacionan la convergencia de los méto-
dos iterativos cuando la matriz del sistema es definido positiva, caso que aparece en muchas aplicaciones
prácticas, pero no nos detendremos en los detalles de las demostraciones.
Teorema 2.7 Sea A una matriz simétrica no singular descompuesta en la forma A = M − N donde M
es no singular. Sea B = M −1 N = Id − M −1 A la matriz de iteración. Supongamos que M t + N (que es
simétrica) es definido positiva. Entonces si A es definido positiva ⇒ ρ(B) < 1.
Teorema 2.8 Si A es simétrica y definido positiva ⇒ el método SOR es convergente si 0 < w < 2. En
particular, Si A es simétrica y definido positiva ⇒ el método de Gauss-Seidel es convergente (w = 1).
Teorema 2.9 Si A es simétrica, definido positiva y 2D − A es definido positiva ⇒ el método de Jacobi

es convergente.
Comparación de los métodos de Jacobi y Gauss-Seidel. Búsqueda del parámetro de relaja-

ción óptimo en el método S.O.R.
Teorema 2.10 Si A es una matriz tridiagonal ⇒ ρ(L1 ) = (ρ(J))2 . Entonces los dos métodos Jacobi y
Gauss-Seidel convergen o divergen simultáneamente, y si convergen, Gauss-Seidel lo hace más rápida-
mente.
Teorema 2.11 Si A es definido positiva y tridiagonal, y 0 < w < 2 ⇒ los tres métodos Jacobi, Gauss-
Seidel y SOR convergen y la elección óptima del parámetro es
2
wop = √
1 + 1 − ρ(J)2
siendo ρ(Lwop ) = w − 1.
Capı́tulo 3
Interpolación
En este tema y el siguiente intentaremos dar respuesta a una situación bastante habitual en el ámbito
cientı́fico: investigamos un fenómeno fı́sico/quı́mico que se está desarrollando ante nuestros ojos, podemos
tomar muestras experimentales y a partir de estas mediciones obtener más información. Para ello podemos
intentar recrear/reconstruir el fenómeno en su totalidad (en un dominio continuo del espacio/tiempo o
cualquier otra magnitud) con una función que represente lo mejor posible esos datos.
Las técnicas que utilizan funciones continuas y que vamos a estudiar en este y el próximo capı́tulo
son de dos tipos:
Interpolación: cálculo de funciones que pasan (interpolan es el término matemático) exactamente

por los puntos señalados por los datos.
Aproximación: cálculo de funciones que aproximan los datos en un cierto sentido (para una deter-
minada forma de medir el error).
En este capı́tulo trataremos el problema de la interpolación, que además tiene mucha utilidad al tratar
la derivación y la integración numérica.
Problema de interpolación: sean (xi , yi ) para i = 0, . . . , n, pares de valores reales (puntos del plano)
tales que xi ̸= xj para i ̸= j, buscamos una función p(x) de un determinado tipo, tal que p(xi ) = yi para
i = 0, . . . , n.
Los datos a interpolar pueden proceder de mediciones experimentales como hemos mencionado an-
tes: conocida experimentalmente la respuesta yi obtenida bajo condiciones xi , nos interesa encontrar el
resultado y que obtendrı́amos al tomar condiciones x no experimentales. Pero también podemos pensar
que los puntos dados forman parte de la gráfica de una función f que queremos conocer al menos apro-
ximadamente y de la que únicamente sabemos su valor en ciertos puntos xi . A partir de ahora, para una
mayor generalidad, hablaremos del problema de interpolación de funciones.
51
52 Interpolación
Al plantearse un problema de interpolación, uno debe contestar a tres preguntas:
¿De qué tipo debe ser la función p(x) buscada? Polinomial, trigonométrica, racional, exponencial,
etc. El comportamiento de los datos a interpolar nos puede orientar sobre el tipo de función in-
terpoladora a elegir: si f tiene un comportamiento periódico, elegiremos funciones trigonométricas;
si sospechamos que f puede tener ası́ntotas, convendrá que p sea racional; si f responde a un
comportamiento polinómico, buscaremos p entre las funciones polinómicas. En este capı́tulo sólo
trataremos este último caso, el de la interpolación polinómica.
Una vez elegido el conjunto de funciones en el que debemos buscar p, ¿existe la función buscada?,
y si existe, ¿es única?.
¿Es la función p una buena aproximación de la función f fuera de los puntos de interpolación?.
3.1. Interpolación polinómica.
3.1.1. Planteamiento del problema
Dados n + 1 puntos de interpolación (xi , yi ) para i = 0, . . . , n, con xi ̸= xj para i ̸= j, llamamos

interpolación polinomial a la determinación de un polinomio p(x) de grado ≤ N tal que
p(xi ) = yi i = 0, . . . , n
Si yi es el valor de una función f en xi para i = 0, . . . , n, hablaremos de la interpolación polinomial de la

función f en las abscisas de interpolación o nodos xi , i = 0, . . . , n.
3.1.2. Tipo de función interpoladora
La función p buscada formará parte del conjunto de polinomios de grado ≤ N , para un cierto N que
determinaremos más adelante, es decir, será de la forma,
p(x) = aN xN + aN −1 xN −1 + . . . + a1 x + a0
y para determinarla habrá que encontrar los N + 1 coeficientes a0 , a1 . . . . , aN .
3.1.3. Existencia y unicidad del polinomio interpolador
Teorema 3.1 Dados x0 , x1 , . . . , xn n + 1 valores reales distintos, para cada conjunto de n + 1 valores
arbitrarios y0 , y1 , y2 , . . . , yn existe un único polinomio pn (x) de grado a lo más n tal que p(xi ) = yi para
i = 0, 1, . . . , n.
3.1 Interpolación polinómica. 53
Demostración:
Demostremos en primer lugar la unicidad : supongamos que hubiera dos polinomios pn (x) y qn (x)
verificando las condiciones del teorema. Por tanto, pn (x) − qn (x) es un polinomio de grado a lo más n
verificando (pn − qn )(xi ) = 0 para i = 0, 1, . . . , n, es decir, tiene n + 1 raı́ces pero es de grado n, por lo
tanto, pn − qn ≡ 0, es decir, pn ≡ qn .
La existencia la demostraremos por inducción sobre n: para n = 0, si p0 (x0 ) = y0 , se trata de la

función constante p0 (x) = y0 . Supongamos que el teorema es cierto para n ≤ k − 1, demostrémoslo para
n = k. Por hipótesis de inducción, existe un polinomio pk−1 de grado a lo más k − 1 tal que pk−1 (xi ) = yi
para i = 0, 1, . . . , k − 1. Tratemos de construir pk de la siguiente forma,
pk (x) = pk−1 (x) + c(x − x0 )(x − x1 ) . . . (x − xk−1 )
que es un polinomio de grado a lo más k verificando pk (xi ) = yi para i = 0, 1, . . . , k − 1. Para determinar

pk basta calcular el valor de c despejando de
pk−1 (xk ) + c(xk − x0 )(xk − x1 ) . . . (xk − xk−1 ) = yk
posible puesto que todos los nodos xi son distintos.
3.1.4. Métodos de cálculo del polinomio interpolador.
Podemos dar otra demostración del Teorema 4.1. que nos permite calcular el polinomio interpolador
p(x) = a0 + a1 x + . . . + an xn simplemente imponiendo las n + 1 condiciones que debe cumplir,

a0 + a1 x0 + . . . + an xn0 = y0 


a0 + a1 x1 + . . . + an xn1 = y1
... 


a0 + a1 xn + . . . + an xnn = yn
que es un sistema lineal de n + 1 ecuaciones con n + 1 incógnitas que son los coeficientes a0 , a1 , . . . , an
del polinomio pn . El determinante de la matriz de coeficientes es el determinante de Vandermonde que
tiene la forma,
1 x0 ... xn0
1 x1 ... xn1 ∏
= (xk − xi )
. . ... .
k>i
1 xn ... xnn
y es no nulo ya que xk ̸= xi si k ̸= i. Por tanto el sistema es compatible determinado y tiene solución

única.
Este proceso para calcular pn es excesivamente laborioso cuando n no es pequeño. Hay que entender
que el polinomio interpolador es único pero se puede expresar de muy diversas formas y llegar hasta él a
través de diferentes algoritmos. Estudiaremos dos métodos para calcular el polinomio interpolador.
54 Interpolación
Método de Lagrange
Se toma como expresión del polinomio interpolador la fórmula de interpolación de Lagrange,
n ∏
∑ j̸=i (x − xj )
pn (x) = yi li (x), li (x) = ∏ , i = 0, 1, . . . , n
i=0 j̸=i (xi − xj )
donde li (x) son los llamados polinomios de Lagrange, que son de grado n y verifican li (xj ) = δij para
i, j = 0, 1, . . . , n, por tanto pn (xi ) = yi para i = 0, 1, . . . , n como se deseaba.
Ejemplo 3.1 Encontrar el polinomio interpolador de la siguiente tabla de datos, mediante el método de
Lagrange.
x 1 2 4 5
y 0 2 12 21
Ejemplo 3.2 Encontrar el polinomio interpolador de la siguiente tabla de datos, mediante el método de
Lagrange. Observar que es la ecuación de la recta que pasa por los puntos del plano (x0 , y0 ) y (x1 , y1 ).
x x0 x1
y y0 y1
Método de diferencias divididas de Newton
Expresamos el polinomio interpolador de la siguiente forma,
pn (x) = c0 + c1 (x − x0 ) + c2 (x − x0 )(x − x1 ) + . . . + cn (x − x0 )(x − x1 ) . . . (x − xn−1 )
El método de las diferencias divididas de Newton permite calcular los coeficientes cj para j = 0, 1, . . . , n,
mediante la construcción de las llamadas diferencias divididas:
f [xi ] = yi , (i = 0, . . . , n)
f [xi+1 , . . . , xi+j+1 ] − f [xi , . . . , xi+j ]
f [xi , xi+1 , . . . , xi+j , xi+j+1 ] = ,
xi+j+1 − xi
(i = 0, . . . , n − j), (j = 0, . . . , n − 1)
de forma que cj = f [x0 , x1 , . . . , xj ], (j = 0, . . . , n), es decir, el polinomio interpolador viene dado por
la siguiente fórmula de interpolación de Newton:
pn (x) = f [x0 ] + f [x0 , x1 ](x − x0 ) + f [x0 , x1 , x2 ](x − x0 )(x − x1 )+

+ . . . + f [x0 , . . . , xn ](x − x0 )(x − x1 ) . . . (x − xn−1 )
3.1 Interpolación polinómica. 55
Ejemplo 3.3 Veamos el esquema de construcción de las diferencias divididas de Newton para n = 2,
x0 f [x0 ] = y0
↘
f [x1 ]−f [x0 ]
f [x0 , x1 ] = x1 −x0
↗ ↘
f [x1 ,x2 ]−f [x0 ,x1 ]
x1 f [x1 ] = y1 f [x0 , x1 , x2 ] = x2 −x0
↘ ↗
f [x2 ]−f [x1 ]
f [x1 , x2 ] = x2 −x1
↗
x2 f [x2 ] = y2
El método de diferencias divididas de Newton tiene la ventaja de que cuando se añaden puntos de
interpolación puede aprovecharse todo el trabajo hecho, basta con continuar el esquema de construcción de
diferencias divididas y calcular las nuevos coeficientes cn+1 , cn+2 , . . ., aprovechando ası́ todos los cálculos
previos.
Ejemplo 3.4 Con la siguiente tabla de datos,
x 1 2 4 5
f (x) 0 2 12 21
mediante el método de diferencias divididas de Newton, aproximar el valor de f (3), usando
el polinomio interpolador de grado 2 calculado usando los tres primeros nodos de la tabla,
el polinomio interpolador de grado 2 calculado usando los tres últimos nodos de la tabla,
el polinomio interpolador de grado 3 calculado usando todos los datos de la tabla.
3.1.5. Error de interpolación
Nos interesa tener un criterio para medir la proximidad del polinomio pn a la función f fuera de los
puntos de interpolación xk .
Teorema 3.2 Sea f ∈ C n+1 (a, b) y sea pn un polinomio de grado a lo más n que interpola a f en n + 1
puntos distintos x0 , x1 , . . . , xn del intervalo (a, b). Entonces para cada x ∈ (a, b) existe un ξx ∈ (a, b) tal
que
∏n
1
f (x) − pn (x) = f n+1) (ξx ) (x − xi )
(n + 1)! i=0
Si analizamos el error podemos observar tres términos diferentes:

56 Interpolación
1
(n+1)! −
−−−→
n→∞
0,
f n+1) (ξx ) que depende de si la derivada n + 1-ésima de la función a interpolar está acotada,
∏n
i=0 (x − xi ) que depende de la colocación de los nodos de interpolación.
Una pregunta natural en este contexto es la siguiente: Supongamos que dado un intervalo (a, b) lo
vamos subdividiendo en más puntos, concretamente, xj = a+jh para j = 0, 1, 2, ..., n, donde h = (b−a)/n
y supongamos que construimos con estos puntos el polinomio de interpolación pn (x) para una función
dada f , esto es, que pn (xi ) = f (xi ), para estos n puntos. La pregunta es, ¿tenderá a 0 el error a medida que
crece en número de nodos de interpolación, es decir, el grado del polinomio interpolador?. La respuesta
es NO.
Ejemplo 3.5 Comparar lo que sucede con el error de interpolación al aumentar el número de nodos de
interpolación para las siguientes funciones:
Sen(πx) en el intervalo (0, 1.5),

1
1+25x2 en el intervalo (−1, 1).
1
Lo que ocurre al aproximar la función 1+25x2 en el intervalo (−1, 1) con polinomios de grado alto es
lo que se conoce como el efecto Runge. La aproximación es mala en los extremos del intervalo, ası́ que
una idea para mejorar dicha aproximación es la de olvidarnos de tomar nodos igualmente espaciados y
tomar nodos que se concentren más cerca de los extremos. De este modo al obligar al polinomio pn (x) a
pasar por estos puntos quizás se mejore la aproximación. Por supuesto que tiene que haber un equilibrio
en la disposición de los nodos xi , pues si ponemos pocos puntos en la región central del intervalo quizás
perderı́amos allı́. Estas ideas son las que llevan a una teorı́a de aproximación muy bonita, donde resulta
que los nodos a usar son los ceros de los llamados polinomios de Chebyshev Tn (x), dados por,
(2k + 1)π
xk = cos , k = 0, 1, . . . , n
2(n + 1)
Otra posibilidad para reducir el error de interpolación debido al uso de polinomios interpoladores de
grado alto es el uso de la interpolación polinómica a trozos o splines. Consiste en trazar una serie de
puntos que uniremos por pedazos de curvas cúbicas. Esto es, tomamos un polinomio de grado 3 distinto
que une cada par de puntos consecutivos a interpolar. Los coeficientes de cada polinomio se tienen que
tomar adecuadamente para que hasta las segundas derivadas coincidan en los puntos de enganche. El
resultado es una curva suave agradable a la vista.
3.2. Interpolación de Hermite.
El término interpolación de Hermite hace referencia a la interpolación de una función y de algunas

de sus derivadas en un conjunto de nodos.
3.2 Interpolación de Hermite. 57
3.2.1. Ejemplo sencillo
Sean x0 , x1 dos puntos donde conocemos el valor de una función f y también de su primera derivada
′
f . Buscamos el polinomio p de menor grado que verifique,
p(x0 ) = f (x0 ), p(x1 ) = f (x1 ),

p′ (x0 ) = f ′ (x0 ), p′ (x1 ) = f ′ (x1 ).
En vista de que hay cuatro condiciones, parece lógico buscar p en el espacio de polinomios de grado ≤ 3,
escribámoslo de la siguiente forma,
p(x) = a + b(x − x0 ) + c(x − x0 )2 + d(x − x0 )2 (x − x1 ),
cuya derivada se escribe,
p′ (x) = b + 2c(x − x0 ) + 2d(x − x0 )(x − x1 ) + d(x − x0 )2 .
Imponiendo las condiciones y denotando h = x1 − x0 , obtenemos,
f (x0 ) = a f (x1 ) = a + bh + ch2 ,

f (x0 ) = b f (x1 ) = b + 2ch + dh2 .
′ ′
Por tanto, despejando,

a = f (x0 ) c = (f (x1 ) − a − bh)/h2 ,
b = f ′ (x0 ) d = (f ′ (x1 ) − b − 2ch)/h2 .
3.2.2. Problema de Hermite generalizado
La interpolación de Hermite puede generalizarse al caso en que conocemos la función f en una serie
de nodos xi para i = 0, 1, . . . , n y sus respectivas derivadas hasta un cierto orden que puede ser distinto
en cada nodo.
Teorema 3.3 Dados x0 , x1 , . . . , xn n + 1 nodos distintos dos a dos, y los valores de la función f y
derivadas sucesivas en esos nodos,
f j) (xi ), j = 0, 1, . . . , ki − 1, i = 0, 1, . . . , n,
entonces existe un único polinomio pN de grado a lo más N con N + 1 = k0 + k1 + . . . + kn verificando

las condiciones de interpolación
j)
pN (xi ) = f j) (xi ), j = 0, 1, . . . , ki − 1, i = 0, 1, . . . , n.
Teorema 3.4 Buscamos un polinomio de grado a lo más N , que tiene N + 1 coeficientes, e imponemos
N + 1 condiciones. Por tanto, tenemos que resolver un sistema lineal de N + 1 ecuaciones con N + 1
incógnitas y deseamos asegurarnos de que la matriz de coeficientes es no singular para que exista una
58 Interpolación
solución única. Para demostrar que una matriz cuadrada es no singular basta con demostrar que el
correspondiente sistema homogéneo tiene como única solución la idénticamente nula. En nuestro caso
esto se corresponderı́a con encontrar un polinomio q de grado a lo más N verificando,
q j) (xi ) = 0, j = 0, 1, . . . , ki − 1, i = 0, 1, . . . , n.
es decir, buscamos un polinomio q de grado a lo más N que tiene un cero con multiplicidad ki en xi para
∏n
i = 0, 1, . . . , n, y por tanto debe ser múltiplo de i=0 (x − xi )ki que es de grado N + 1, imposible a no
ser que q ≡ 0 como deseábamos.
En lo que se refiere al error en este caso, puede decirse que si f ∈ C N +1 (a, b) y xi , ∈ (a, b) para
i = 0, 1, . . . , n, entonces para cada x ∈ (a, b) existe un ξx ∈ (a, b) tal que
1
f (x) − pN (x) = f N +1) (ξx )(x − x0 )k0 (x − x1 )k1 . . . (x − xn )kn
(N + 1)!
3.2.3. Caso particular: el polinomio de Taylor
El caso de interpolación de Hermite en un solo nodo se trata del conocido polinomio interpolador de
Taylor: sea f ∈ C n+1 (a, b), para cada x0 ∈ (a, b), existe un único polinomio pn de grado a lo más n tal
j)
que pn (x0 ) = f j) (x0 ) para j = 0, 1, . . . , n, que es el polinomio de Taylor,
f ′′ (x0 ) f n) (x0 )
pn (x) = f (x0 ) + f ′ (x0 )(x − x0 ) + (x − x0 )2 + . . . + (x − x0 )n ,
2! n!
junto con la fórmula del error de la interpolación de Taylor,
1
f (x) − pn (x) = f n+1) (ξx )(x − x0 )n+1 ,
(n + 1)!
para cierto ξx ∈ (a, b).
3.2.4. Método de las diferencias divididas de Newton generalizado
Para calcular el polinomio interpolador de Hermite se usa una generalización del método de diferencias
divididas de Newton en la que el esquema triangular se construye de la siguiente manera: en la primera
columna se coloca cada nodo repetido tantas veces como condiciones haya sobre él; en la segunda columna
los respectivos valores de la función a interpolar en los nodos correspondientes, es decir, f [xi ] = f (xi )
tantas veces como condiciones sobre el nodo i tengamos, para i = 0, 1, . . . , n; en la tercera columna cuando
aparezcan dos nodos iguales, tendremos en cuanta que,
f (x) − f (xi )
f ′ (xi ) = lı́m = lı́m f [xi , x] = f [xi , xi ],
x→xi x − xi x→xi
3.2 Interpolación de Hermite. 59
y en general,
f j) (xi )
f [xi , xi , . . . , xi ] = ,
| {z } j!
j
de modo que el polinomio interpolador serı́a,
pN (x) = f [x0 ] + f [x0 , x0 ](x − x0 ) + . . . + f [x0 , x0 , . . . , x0 ](x − x0 )k0 −1 + . . .

| {z }
k0
+f [x0 , x0 , . . . , x0 , x1 ](x − x0 )k0 + . . .
| {z }
k0
+f [x0 , . . . , x0 , x1 , . . . , x1 ](x − x0 )k0 (x − x1 )k1 −1 + . . .
| {z } | {z }
k0 k1
+...
f [x0 , . . . , x0 , . . . , xn , . . . , xn ](x − x0 )k0 . . . (x − xn )kn −1
| {z } | {z }
k0 kn
3.2.5. Ejemplo sencillo
El triángulo de diferencias divididas de Newton que deberı́amos construir para el ejemplo sencillo
propuesto antes, en el que conocemos el valor de la función y su primera derivada en dos nodos, serı́a,
x0 f [x0 ] = f (x0 )
↘
f [x0 , x0 ] = f ′ (x0 )
↗ ↘
x0 f [x0 ] = f (x0 ) f [x0 , x0 , x1 ]
↘ ↗ ↘
f [x0 , x1 ] f [x0 , x0 , x1 , x1 ]
↗ ↘ ↗
x1 f [x1 ] = f (x1 ) f [x0 , x1 , x1 ]
↘ ↗
f [x1 , x1 ] = f ′ (x1 )
↗
x1 f [x1 ] = f (x1 )
y el correspondiente polinomio interpolador,
p(x) = f [x0 ] + f [x0 , x0 ](x − x0 ) + f [x0 , x0 , x1 ](x − x0 )2 +

+f [x0 , x0 , x1 , x1 ](x − x0 )2 (x − x1 ).
Ejemplo 3.6 Siendo f (x) = x12 , hallar el polinomio p11 (x) verificando,
i)
p11 (−1) = f i) (−1), i = 0, 1, 2, 3,
i)
p11 (0) = f i) (0), i = 0, 1, 2,
i)
p11 (1) = f i) (1), i = 0, 1, 2, 3, 4.
60 Interpolación
Capı́tulo 4
Aproximación numérica.
4.1. Introducción.
En el capı́tulo anterior hemos hablado de aproximación de funciones mediante interpolación que tiene
muchas ventajas: el polinomio interpolador es fácil de calcular y se dispone de una fórmula explı́cita para
el error de interpolación; la interpolación es muy útil para generar fórmulas de derivación e integración
numérica; la interpolación es especialmente apropiada para el cálculo de funciones dadas por tablas, es
decir, para funciones bien conocidas sobre conjuntos discretos de abscisas donde el error de redondeo de
los valores es menor que el error propio de interpolación.
Sin embargo, la interpolación presenta ciertos problemas en otros casos, por ejemplo: si tenemos un
conjunto discreto de valores (xk , yk ) (k = 0, 1, . . . , m) que tienen errores de redondeo apreciables, no es
conveniente utilizar el polinomio interpolador que interpole exactamente esos datos ya que su carácter
oscilante puede provocar que el error fuera de los puntos de interpolación sea muy grande; otra situación en
la que tampoco es conveniente la interpolación es cuando conocemos una función f en todo un intervalo
I, generar una tabla de valores y buscar el polinomio interpolador no es la manera más eficiente de
aproximar dicha función.
Si nos encontramos con el caso discreto, es decir, un conjunto discreto de valores (xk , yk ) (k =
0, 1, . . . , m), podemos pensar en buscar un polinomio pn de grado n ≤ m tal que los errores, ek =
yk − pn (xk ) (k = 0, 1, . . . , m) sean lo más pequeños posibles en un sentido que determinaremos más
adelante. Este es el proceso de aproximación polinomial. Podemos aproximar estos datos con otro tipo de
función, fn (x) = a0 φ0 (x) + . . . + an φn (x), donde debemos encontrar los parámetros a0 , . . . , an de forma
que los errores ek = yk − fn (xk ) (k = 0, 1, . . . , m) sean lo más pequeños posibles en un cierto sentido.
Este es el problema de aproximación discreta.
En el caso continuo, cuando conocemos la función a aproximar f en todo un intervalo I, buscamos

una función fn de forma que la función de error de aproximación en (x) = f (x)−fn (x) sea lo más pequeña
61
62 Aproximación numérica.
posible sobre el intervalo I en algún sentido que se determina a priori. La función aproximadora puede ser
polinómica o de forma general, fn (x) = a0 φ0 (x) + . . . + an φn (x), donde φ0 , . . . , φn son funciones dadas,
fácilmente calculables, y el problema se reduce a calcular los parámetros a0 , . . . , an .
Por tanto, el problema general de aproximación es: dado un conjunto I de abscisas de aproximación,
y unas funciones básicas φj (j = 0, 1, . . . , n) definidas sobre I, para cada función f definida sobre I,
buscamos los coeficientes a0 , . . . , an , de forma que fn (x) = a0 φ0 (x) + . . . + an φn (x) haga que la magnitud
del error de aproximación en (x) = f (x) − fn (x) sea lo más pequeña posible.
Por tanto, para determinar totalmente un problema de aproximación es necesario especificar: el con-
junto I de abscisas de aproximación, las funciones básicas y la forma de medir la magnitud del error.
4.1.1. Conjunto de abscisas de aproximación
Si I es finito (I = x0 , . . . , xm ) hablaremos de aproximación discreta; si I es un intervalo de extremos

a y b (a < b), hablaremos de aproximación continua. Dar una función f sobre un conjunto finito I =
x0 , . . . , xm equivale a dar yk = f (xk ) (k = 0, . . . , m).
4.1.2. Funciones básicas
Las funciones φ0 , . . . , φn , definidas sobre I, pueden escogerse de diversas formas dependiendo del
comportamiento de la función f a aproximar. Si f es periódica, elegiremos las funciones básicas entre
las funciones trigonométricas, por ejemplo: φ0 (x) = 1, φ1 (x) = senx, φ2 (x) = cosx, . . . , φ2s−1 (x) =
sen2sx, φ2s (x) = cos2sx, donde n = 2s, que llamaremos aproximación trigonométrica; Si f responde a un
comportamiento polinómico, elegiremos cada φj (x) = pj (x) entre los polinomios de grado j (j = 0, . . . , n
(por ejemplo, φj (x) = xj aunque esta no será siempre la elección más adecuada) y hablaremos de
aproximación polinomial.
Observar que nos estamos limitando al caso de aproximación lineal, es decir, buscamos la función
aproximadora fn (x) en el espacio vectorial generado por las funciones básicas,
En = ⟨φ0 , . . . , φn ⟩,
esto es,
n
∑
fn (x) = aj φj (x), x ∈ I.
j=0
4.1.3. Medida de la magnitud del error: normas funcionales
La magnitud del error de aproximación se puede medir de diferentes maneras según sea el caso discreto
o el caso continuo, y según la norma que utilicemos.
4.1 Introducción. 63
Caso discreto
El error de aproximación es un vector de m + 1 valores,
ek = f (xk ) − fn (xk ), (k = 0, . . . , m)
por tanto, para medirlo utilizaremos una norma vectorial. Las dos normas más usadas son:
la norma euclı́dea
(∑
m ) 21
∥e∥2 = |ek |2
k=0
la norma del máximo

∥e∥∞ = máx |ek |
k=0÷m
Cuando quiere darse una importancia diferente a los distintos términos del error se usan normas
ponderadas introduciendo coeficientes positivos llamados pesos w = {wk }k=0÷m ,
(∑
m ) 21
∥e∥2,w = wk |ek |2 , ∥e∥∞,w = máx wk |ek |.
k=0÷n
k=0
Caso continuo
El error de aproximación es una función definida en el intervalo I = [a, b], definimos:
la norma euclı́dea
(∫ b ) 21
∥e∥2 = |e(x)|2 dx
a
la norma del máximo

∥e∥∞ = máx |e(x)|
x∈I
Se puede probar que estas definiciones cumplen las propiedades de norma sobre el conjunto C([a, b])
de funciones continuas sobre el interval [a, b].
Como en el caso discreto, se pueden definir las correspondientes normas ponderadas introduciendo
una función peso w ∈ C([a, b]) positiva (w(x) > 0 sobre I),
(∫ b ) 21
∥e∥2,w = w(x)|e(x)|2 dx , ∥e∥∞ = máx |e(x)|w(x).
a x∈I
Tanto en el caso discreto como en el continuo, si se elige la norma euclı́dea hablaremos de aproximación
por mı́nimos cuadrados, si se elige la norma del máximo, hablaremos de aproximación minimax. En este
capı́tulo nos centraremos en la aproximación por mı́nimos cuadrados.
4.2. Aproximación por mı́nimos cuadrados.
4.2.1. Definición del problema
Consideremos un conjunto de abscisas I, ya sea continuo o discreto, unas funciones básicas φj (j =

0 ÷ n), y el espacio vectorial que generan En . Para cada función f definida sobre I, buscamos una función
fn∗ ∈ En tal que ∥f − fn∗ ∥2 sea mı́nima en En , es decir,
∥f − fn∗ ∥2 = mı́n ∥f − fn ∥2 ,
fn ∈En
donde ∥ ∥2 representa aquı́ cualquiera de las normas euclı́deas, ponderada o no, tanto en el caso continuo
como discreto.
en el caso discreto, I = x0 , . . . , xm y si e = (e0 , . . . , em ),

(∑
m ) 12
∥e∥2,w = wk |ek |2
k=0
donde w = w0 , . . . , wm es una colección de pesos positivos;
en el caso continuo, I es un intervalo de la recta real de extremos a y b,

(∫ b ) 21
∥e∥2,w = w(x)|e(x)|2 dx
a
donde w(x) > 0 es una función peso sobre I.
4.2.2. Productos escalares asociados
La propiedad fundamental de las normas euclı́deas es que provienen de sendos productos escalares:
en el caso discreto,
m
∑
(u, v) = w k uk v k ,
k=0
en el caso continuo,
∫ b
(u, v) = w(x)u(x)v(x)dx
a
en el sentido que se cumple, en ambos casos,
∥e∥22 = (e, e).
Estos productos escalares cumplen las propiedades de definición de producto escalar:

4.2 Aproximación por mı́nimos cuadrados. 65
(u, u) ≥ 0 y (u, u) = 0 si y sólo si u = 0,
(u, v) = (v, u),
(a1 u1 + a2 u2 , v) = a1 (u1 , v) + a2 (u2 , v), para funciones u1 , u2 , v sobre I y números reales a1 , a2

cualesquiera.
4.2.3. Ecuaciones normales.
Sea fn∗ una función sobre I tal que,
(f − fn∗ , fn ) = 0, ∀fn ∈ En ,
entonces tenemos,
∥f − fn ∥22 = (f − fn , f − fn ) = (f − fn∗ + fn∗ − fn , f − fn∗ + fn∗ − fn ) =

= (f − fn∗ , f − fn∗ ) + 2(fn − fn∗ , fn∗ − f ) + (fn∗ − fn , fn∗ − fn ) =
= ∥f − fn∗ ∥22 + ∥fn∗ − fn ∥22 ,
por tanto,
∥f − fn ∥22 = ∥f − fn∗ ∥22 + ∥fn∗ − fn ∥22 , ∀fn ∈ En .
En particular,
∥f − fn ∥2 ≥ ∥f − fn∗ ∥2 , ∀fn ̸= fn∗ ∈ En ,
es decir, fn∗ es la única función de En que satisface la condición de aproximación por mı́nimos cuadrados,
∥f − fn∗ ∥2 = mı́n ∥f − fn ∥2 .
fn ∈En
Dado que En está generado por las funciones básicas φi (i = 0 ÷ n), y fn∗ ∈ En , podemos escribir,
n
∑
fn∗ (x) = a∗j φj (x),
j=0
y la condición anterior equivale a encontrar los coeficientes a∗j (j = 0÷n), tales que satisfagan las llamadas
ecuaciones normales,
n
∑
(φi , φj )a∗j = (φi , f ) (i = ÷n).
j=0
Este sistema puede escribirse en forma matricial,
Aa∗ = b,
donde A = ((φi , φj ))i,j=0÷n , a∗ = (a∗j )j=0÷n y b = ((φi , f ))i=0÷n .

La matriz A es semidefinido positiva, es decir, simétrica y para cualquier vector a = (a0 , a1 , . . . , an )t ,

se tiene,
(∑n n
∑ ) ∑n
at Aa = aj φj , ai φi = ∥ aj φj ∥22 = ∥fn ∥22 ≥ 0
j=0 i=0 j=0
∑n
donde fn viene dada por, fn (x) = j=0 aj φj (x).
Esta relación nos muestra además que las funciones básicas φj son linealmente independientes si y
sólo si detA ̸= 0, y que las ecuaciones normales tienen solución única para cualquier f si y sólo si las
funciones básicas son linealmente independientes.
4.2.4. Un ejemplo sencillo: la recta de regresión
Tenemos un conjunto de puntos del plano (xk , yk ) (k = 0 ÷ m), con m > 2, y buscamos una recta
∑m
y = a0 +a1 x que los aproxime de modo que minimice k=0 d2k la suma de los cuadrados de las desviaciones
dk = yk − a0 − a1 x (k = 0 ÷ m).
Este no es más que un problema de aproximación discreta por mı́nimos cuadrados con I = x0 , x1 , . . . , xm ,
∑m
φ0 (x) = 1, φ1 (x) = x, todos los pesos iguales a 1, y el producto escalar (u, v) = k=0 uk vk .
Las correspondientes ecuaciones normales,

( )( ∗ ) ( )
(φ0 , φ0 ) (φ0 , φ1 ) a0 (φ0 , f )
=
(φ1 , φ0 ) (φ1 , φ1 ) a1 ∗ (φ1 , f )
forman el siguiente sistema lineal de dos ecuaciones con dos incógnitas,

( ∑m ) ( ∗ ) ( ∑m )
m
∑m +1 ∑m k=0 xk a0 ∑ k=0 yk
2 = m
k=0 xk k=0 xk a1 ∗ k=0 xk yk
cuya solución es:

a∗1 = xy−x·y ,
x2 −x2
∗ ∗
a0 = y − a1 x,
donde la barra indica la media, es decir,
1
∑m 1
∑m
x = m+1 xk y = m+1 k=0 yk
1
∑k=0
m 1
∑ m
x = m+1 k=0 x2k
2 xy = m+1 k=0 xk yk
4.3. Ortogonalización.
Una vez reducido el problema de aproximación por mı́nimos cuadrados, es necesario resolver el sistema
de ecuaciones normales asociadas,
Aa∗ = b,
4.3 Ortogonalización. 67
donde A = ((φi , φj ))i,j=0÷n , a∗ = (a∗j )j=0÷n y b = ((φi , f ))i=0÷n .
Como la matriz A es semidefinido positiva siempre que las funciones básicas sean linealmente inde-
n3
pendientes, un método especialmente adecuado es el de Cholesky que requiere 6 + ϑ(n2 ) operaciones.
El trabajo preliminar de construcción de las ecuaciones normales requiere generalmente 12 p(n+1)(n+4)

operaciones, donde p es el número de operaciones necesarias para cada producto escalar, que normalmente
será mayor que n + 1: en el caso discreto, p = m + 1 ≥ n + 1, si I consta de m + 1 elementos; en
el caso continuo hay que calcular las correspondientes integrales. Por tanto, la mayor parte del cálculo
corresponde a la formación de las ecuaciones normales. Ahora bien, esta parte del cálculo está fuertemente
condicionada por la elección de las funciones básicas de En . En general, debemos calcular todos los
productos escalares, es decir, todos los coeficientes de la matriz A. También debemos tener en cuenta
que la matriz A puede estar mal condicionada si las funciones básicas son ”poco independientes”desde
el punto de vista numérico. Todo esto nos lleva a pensar que una buena elección de de una base de
funciones de En reduce considerablemente los cálculos, por ejemplo con una base de funciones ortogonales
ψj (j = 0 ÷ n) respecto al producto escalar, es decir (ψi , ψj ) = 0, ∀i ̸= j y (ψi , ψi ) > 0 (i = 0 ÷ n), el
sistema es diagonal y la solución es inmediata:
n
∑ (ψj , f )
fn∗ (x) = c∗j ψj (x), c∗j = .
j=0
(ψj , ψj )
Dada la simplicidad de estas expresiones, los métodos estándar de resolución de las ecuaciones nor-
males, están basados en la ortogonalización de las funciones básicas, es decir, en la expresión de las
ecuaciones normales en una base de funciones ortogonales.
4.3.1. Ortogonalización de Gram-Schmidt
Consideramos una base de funciones φi (x) (i = 0, ÷n) de nuestro espacio vectorial En que está dotado
del correspondiente producto escalar (·, ·). Buscamos otra base de funciones de En , ψi (x) (i = 0, ÷n), que
sean ortogonales respecto a ese producto escalar, es decir,
(ψi (x), ψj (x)) = 0 ∀i ̸= j

(ψi (x), ψi (x)) > 0 ∀i = 0, . . . , n
El proceso es
ψ0 (x) = φ0 (x)
∑i−1
ψi (x) = φi (x) − j=1 αij ψj (x) , i = 1, . . . , n
(φi (x), ψj (x))
con αij =
(ψj (x), ψj (x))
Cuando el espacio de funciones es el espacio de polinomios de grados ≤ n, partiendo de la base de

polinomios φi (x) = xi (i = 0, ÷n), tenemos la siguiente recurrencia para calcular una base de polinomios
ortogonales respecto a un producto escalar (·, ·) determinado.
p0 (x) = 1
p1 (x) = x − a1
pj (x) = (x − aj )pj−1 (x) − bj pj−2 (x) j≥2
(xpj−1 (x), pj−1 (x))
con aj =
(pj−1 (x), pj−1 (x))
(xpj−1 (x), pj−2 (x))
bj =
(pj−2 (x), pj−2 (x))
Capı́tulo 5
Integración y derivación numéricas
5.1. Integración numérica.
La integración numérica es el proceso por medio del cual se genera un valor numérico que aproxima
el valor de la integral definida de una función que no posee una primitiva fácil de calcular. Para calcular,
∫ b
f (x)dx,
a
buscamos primero una primitiva, es decir, una función F tal que F ′ = f , y entonces
∫ b
f (x)dx = F (b) − F (a).
a
2
Pero existen muchas funciones elementales que no poseen primitivas sencillas, por ejemplo, f (x) = ex .
Una primitiva de esta función es,
∞
∑ x2k+1
F (x) = .
(2k + 1)k!
k=0
∫b
Una estrategia muy poderosa para calcular el valor numérico de la integral a
f (x)dx, consiste en reem-
plazar f por otra función g que aproxime f de manera adecuada en el intervalo de integración y que sea
fácil de integrar. Entonces,
∫ b ∫ b
f ≈g⇒ f (x)dx ≈ g(x)dx.
a a
Por ejemplo, g puede ser un polinomio que interpole a f en un conjunto de nodos o una serie de Taylor.
En el ejemplo anterior,
∫ 1 ∫ 1∑ n ∑n
x2 x2k 1
e dx ≈ dx ≈
0 0 k! (2k + 1)k!
k=0 k=0
69
70 Integración y derivación numéricas
5.1.1. Integración vı́a interpolación. Fórmulas de Newton-Cotes
Deseamos calcular,
∫ b
f (x)dx.
a
Elegimos los nodos x1 , x1 , . . . , xn en [a, b], e iniciamos un proceso de interpolación polinómica de Lagrange,
n ∏
∑ j̸=i (x − xj )
p(x) = f (xi )li (x), li (x) = ∏ , i = 0, 1, . . . , n.
i=0 j̸=i (xi − xj )
Aproximamos,
∫ b ∫ b n
∑ ∫ b
f (x)dx ≈ p(x)dx = f (xi ) li (x)dx.
a a k=0 a
Entonces, para cualquier función f (x) tenemos,
∫ b n
∑ ∫ b
f (x)dx ≈ Ai f (xi ), donde Ai = li (x)dx.
a k=0 a
llamadas Fórmulas de Newton-Cotes.
Regla del trapecio: n = 1
El ejemplo más sencillo de una fórmula de Newton-Cotes es la regla del trapecio que se obtiene para
n = 1, es decir, dos nodos que son los extremos del intervalo de integración, x0 = a, x1 = b. Por tanto,
los correspondientes polinomios de Lagrange son,
b−x x−a
l0 (x) = , l1 (x) = ,
b−a b−a
e integrando,
∫ b ∫ b
1
A0 = l0 (x)dx = (b − a) = l1 (x)dx = A1 ,
a 2 a
obteniéndose la conocida regla del trapecio,
∫ b
b − a( )
f (x)dx ≈ f (a) + f (b)
a 2
Si en el intervalo [a, b] se hace una partición como la siguiente,
a = x0 < x1 < . . . < xn = b,

5.1 Integración numérica. 71
aplicando la regla del trapecio en cada uno de los subintervalos, obtenemos la regla del trapecio compuesta,
∫ n ∫
∑ n
b xi
1∑
f (x)dx = f (x)dx ≈ (xi − xi−1 )(f (xi−1 ) + f (xi )).
a i=1 xi−1 2 i=1
b−a
Los nodos no tiene porque estar espaciados uniformemente, pero si lo están, es decir, tomando h = n ,
definiendo xi = a + hi para i = 0, 1, . . . , n, la regla del trapecio compuesta se escribe,
∫
h[ ]
b n−1
∑
f (x)dx ≈ f (a) + 2 f (a + hi) + f (b)
a 2 i=1
Regla de Simpson: n = 2
Un ejemplo más complicado de fórmula de Newton-Cotes es la regla de Simpson que se obtiene para
n = 2, es decir, con tres nodos, que son los extremos del intervalo de integración y el punto medio,
a+b
x0 = a, x1 = 2 , x2 = b. Procediendo como en el caso anterior,
∫ b ∫ b a+b
(x − 2 )(x − b) b−a
A0 = l0 (x)dx = a+b
dx = ,
a a (a − 2 )(a − b)
6
∫ b ∫ b
(x − a)(x − b) b−a
A1 = l1 (x)dx = dx = 4 ,
a a ( a+b
2 − a)( a+b
2 − b) 6
∫ b ∫ b a+b
(x − a)(x − 2 ) b−a
A2 = l2 (x)dx = a+b
dx = ,
a a (b − a)(b − 2 )
6
obtenemos la regla de Simpson
∫ b
b − a( a+b )
f (x)dx ≈ f (a) + 4f ( ) + f (b)
a 6 2
Si en el intervalo [a, b] se hace una partición con un número par de intervalos, es decir, eligiendo n un
b−a
número par, definimos xi = a + hi para i = 0, 1, . . . , n con h = n , y aplicamos la regla de Simpson a
cada par de intervalos, obtenemos la regla de Simpson compuesta,
∫ n/2
b
h ∑[ ]
f (x)dx ≈ f (x2i−2 ) + 4f (x2i−1 ) + f (x2i ) ,
a 3 i=1
reordenando,
∫ n/2 n/2
b
h[ ∑ ∑ ]
f (x)dx ≈ f (x0 ) + 2 f (x2i−2 ) + 4 f (x2i−1 ) + f (xn ) .
a 3 i=2 i=1
Las fórmulas de Newton-Cotes nos llevan a otras fórmulas de integración más generales del tipo,
∫ b n
∑ ∫ b
f (x)w(x)dx ≈ Ai f (xi ), donde Ai = li (x)w(x)dx,
a k=0 a
donde w(x) es cualquier función de peso.
5.1.2. Método de los coeficientes indeterminados
A medida que elegimos más nodos en las fórmulas de Newton-Cotes, las integrales que tenemos que
calcular para obtener los coeficientes Ai se complican. Hay otro procedimiento para calcular el valor de
estos coeficientes, el llamado método de los coeficientes indeterminados, que consiste en imponer en la
fórmula de integración correspondiente, las condiciones que debe cumplir, es decir, que sea exacta para
polinomios del grado correspondiente.
Veamos cómo se obtiene la regla de Simpson con este método. La regla de Simpson es una expresión
del tipo,
∫ b
a+b
f (x)dx ≈ A0 f (a) + A1 f ( ) + A2 f (b),
a 2
donde tenemos que calcular los coeficientes A0 , A1 y A2 de forma que dicha fórmula sea exacta para
polinomios de grado más alto posible, como tenemos tres grados de libertad, podemos imponer que esta
fórmula sea exacta para polinomios de grado ≤ 2. Basta imponer que la fórmula de integración sea exacta
para f (x) = 1, x, x2 , obteniendo el siguiente sistema de ecuaciones,
∫ b
b−a = 1dx = A0 + A1 + A2
a
∫ b
b2 − a2 a+b
= xdx = A0 a + A1 + A2 b
2 a 2
3 3 ∫ b
b −a ( a + b )2
= x2 dx = A0 a2 + A1 + A2 b2
3 a 2
b−a
de donde podemos despejar A0 = A2 = 6 y A1 = 4 b−a
6 , como corresponde.
5.1.3. Cambio de intervalo
A partir de una fórmula de integración numérica en un intervalo de integración determinado, podemos

deducir la correspondiente fórmula de integración numérica para cualquier otro intervalo de integración
mediante un cambio de variable lineal. Si la primera fórmula es exacta para polinomios de un cierto grado,
lo mismo será cierto para la segunda fórmula. Veamos cómo se lleva a cabo.
5.1 Integración numérica. 73
Supongamos que contamos con la siguiente fórmula de integración numérica,

∫ d n
∑
f (t)dt ≈ Ai f (ti )
c i=0
No nos importa el origen de esta fórmula, sin embargo, supongamos que es exacta para todos los polino-
mios de grado ≤ m. Si necesitamos esta fórmula para algún otro intervaalo, [a, b], definimos primero una
función lineal λ(t) tal que, si t recorre [c, d], entonces λ(t) recorre [a, b]. La expresión explı́cita de λ(t) es,
b−a ad − bc
λ(t) = t+ .
d−c d−c
b−a
Por tanto, para el cambio de variable en la integral tenemos x = λ(t) ⇒ dx = λ′ (t)dt = d−c dt, de donde,
∫ ∫ λ−1 (b)=d n
b
b−a b−a∑
f (x)dx = f (λ(t))dt ≈ Ai f (λ(ti )),
a d−c λ−1 (a)=c d − c i=0
fórmula que seguirá siendo exacta para todos los polinomios de grado ≤ m.
Ejemplo 5.1 Deducir la fórmula de Simpson para el intervalo [0, 1] por el método de los coeficientes
indeterminados, es muy sencillo, después utilizando este cambio de variable podremos tener la fórmula
de Simpson para cualquier otro intervalo.
5.1.4. Cuadratura gaussiana.
En la sección anterior hemos visto cómo generar fórmulas de integración numérica, también llamadas
fórmulas de cuadratura, del tipo,
∫ b ∑ n
f (x)dx ≈ Ai f (xi ),
a k=0
que son exactas para polinomios de grado ≤ n. En estas fórmulas, la elección de los nodos x0 , x1 , . . . , xn se
hace a priori, y una vez fijados los nodos, los coeficientes Ai se determinan de manera unı́voca imponiendo
la igualdad en la fórmula de cuadratura para todos los polinomios de grado ≤ n. Nos preguntamos ahora
si una elección de nodos puede ser mejor que otra, por ejemplo, nos preguntamos si podrı́a haber un
conjunto particular de nodos para los que todos los coeficientes Ai fueran todos iguales, simplificando ası́
la fórmula de cuadratura.
Partiendo de las fórmulas de cuadratura más generales, a saber,

∫ b n
∑
f (x)w(x)dx ≈ Ai f (xi ),
a k=0
donde w(x) es una función de peso positiva, sabemos que esta fórmula es exacta para polinomios de grado
≤ n si y sólo si,
∫ b j=n
∏ x − xj
Ai = w(x) dx.
a x i − xj
j=0,j̸=i
En vista de que se cuenta con n + 1 coeficientes Ai y n + 1 nodos xi , sin que exista ninguna restricción a
priori sobre estos últimos, sospechamos que se pueden encontrar fórmulas de cuadratura que sean exactas
para polinomios de grado ≤ 2n + 1. El siguiente resultado nos indica dónde colocar los nodos para que
esto sea posible, obteniendo las llamadas fórmulas de cuadratura gaussianas.
Teorema 5.1 Dada una función de peso positiva w, y un polinomio q no nulo de grado n + 1 que sea
w-ortogonal a πn , espacio de polinomios de grado ≤ n, en el sentido de que para cualquier p ∈ πn se
tiene,
∫ b
q(x)p(x)w(x)dx = 0,
a
entonces, si x0 , x1 , . . . , xn son las raı́ces de q, la fórmula de cuadratura,
∫ b n
∑ ∫ b j=n
∏ x − xj
f (x)w(x)dx ≈ Ai f (xi ), Ai = w(x) dx,
a a xi − x j
k=0 j=0,j̸=i
será exacta para todo polinomio de grado ≤ 2n + 1.
Demostración:
Sea f ∈ π2n+1 , dividimos f entre q obteniendo un cociente p y un resto r, f = pq + r. Por tanto

p, r ∈ πn y f (xi ) = r(xi ), para i = 0, 1, . . . , n. Integrando,
∫ b ∫ b ∫ b n
∑ n
∑
f (x)w(x)dx = q(x)p(x)w(x)dx + r(x)w(x)dx = Ai r(xi ) = Ai f (xi ).
a
|a {z } a k=0 k=0
=0
como querı́amos demostrar.
Para poder aplicar la fórmula de integración en ese conjunto de nodos que son las raı́ces de q, es
necesario que éstas sean reales y simples. Esto se deduce de forma inmediata del siguiente resultado.
Teorema 5.2 Sea w una función de peso positiva en C[a, b]. Sea q un elemento no nulo de C[a, b] que sea
w-ortogonal a πn . Entonces q cambia de signo en (a, b) al menos n + 1 veces.
Demostración:
∫b
Como 1 ∈ πn , entonces a
q(x)w(x)dx = 0, mostrando que q cambia de signo al menos una vez en
(a, b) ya que la función de peso w es positiva.
Supongamos que q cambia de signo en sólo r ocasiones, con r ≤ n. Escogemos puntos ti de manera
que a = t0 < t1 < . . . < tr < tr+1 = b, y tal que q sólo tiene un signo en cada intervalo definido por
∏r
estos puntos. Entonces el polinomio p(x) = i=1 (x − ti ) tiene la misma propiedad respecto al signo que
∫b
q y por lo tanto a q(x)p(x)w(x)dx ̸= 0, pero esto es una contradicción puesto que p ∈ πn , a no ser que
r = n + 1 como querı́amos demostrar.
5.2 Derivación numérica. 75
El cálculo de los coeficientes Ai en las fórmulas de cuadratura gaussianas, se realiza del mismo modo
que en el caso de las fórmulas anteriores no gaussianas, una vez determinados los nodos xi . Podemos
calcular directamente su valor mediante las integrales de los correspondientes polinomios de Lagrange, o
mediante el método de los coeficientes indeterminados.
A su vez, los nodos son las raı́ces de un cierto polinomio qn+1 que queda univocamente determinado
mediante dos condiciones:
qn+1 es un polinomio mónico de grado n + 1, es decir, el coeficiente de xn+1 es la unidad.

∫b
qn+1 es w-ortogonal a πn , es decir, a
qn+1 (x)w(x)p(x)dx = 0, ∀p ∈ πn
Estos son los llamados polinomios ortogonales que podemos calcular con la fórmula recurrente vista en
el Tema 5,
p0 (x) = 1
p1 (x) = x − a1
pj (x) = (x − aj )pj−1 (x) − bj pj−2 (x) j≥2
(xpj−1 (x), pj−1 (x))
con aj =
(pj−1 (x), pj−1 (x))
(xpj−1 (x), pj−2 (x))
bj =
(pj−2 (x), pj−2 (x))
de forma que el polinomio pn+1 calculado con esta fórmula será ortogonal a πn en el sentido del producto
escalar usado en la misma, que en nuestro caso debe ser,
∫ b
(p, q) = p(x)w(x)q(x)dx.
a
Ejemplo 5.2 Encontrar la fórmula de cuadratura gaussiana para [a, b] = [−1, 1], w(x) = 1 y n = 1.
5.2. Derivación numérica.
Aunque haya reglas bien conocidas para derivar las funciones más usuales, no siempre pueden ser
utilizadas (por ejemplo, en funciones dadas por tablas de valores), o no es conveniente (por ejemplo,
en funciones con expresiones analı́ticas muy complicadas). En estos casos debemos recurrir a técnicas
numéricas que, partiendo de los valores de la función en diversas abscisas, nos permitirá calcular una
aproximación al valor de alguna de sus derivadas en una abscisa próxima.
5.2.1. Derivadas primeras.
La derivada de una función f en un punto x0 es,
f (x0 + h) − f (x0 )
f ′ (x0 ) = lı́m
h→0 h
lo que nos da una forma obvia de generar una aproximación de f ′ (x0 ),
f (x0 + h) − f (x0 )
f ′ (x0 ) ≈
h
para valores pequeños de h. Aunque esto parezca muy evidente no es demasiado útil debido a los errores
de redondeo, pero es un buen punto de partida.
Fórmulas de derivación interpolatoria.
Para conocer mejor el error que se comete con este tipo de aproximaciones vamos a utilizar las fórmulas
de interpolación polinómica de las que conocemos el error.
Sean x0 , x1 , . . . , xn , n + 1 puntos distintos de un intervalo I en el que f ∈ C n+1 (I), por las fórmulas
de interpolación polinómica sabemos que, para algún ξx ∈ I,
n
∑ ∏n
n+1) k=0 (x
− xk )
f (x) = f (xk )Lk (x) + f (ξx )
(n + 1)!
k=0
donde Lk (x) es el k-ésimo polinomio de Lagrange de los nodos x0 , x1 , . . . , xn , es decir,

∏
i̸=k (x − xi )
Lk (x) = ∏
i̸=k (xk − xi )
Si derivamos esta expresión, obtenemos,

∏n ( ∏n
− xk ) )
n
∑
′ k=0 (x
− xk ) k=0 (x
f (x) = f (xk )L′k (x) + Dx (f n+1)
(ξx )) + f n+1) (ξx )Dx
(n + 1)! (n + 1)!
k=0
que en el caso en que x sea una de los nodos xj , se reduce a,

n
∑ n
f n+1) (ξj ) ∏
f ′ (xj ) = f (xk )L′k (xj ) + (xj − xk )
(n + 1)!
k=0 k=0
k̸=j
llamada fórmula de derivación de n + 1 puntos para aproximar f ′ (xj ).
En términos generales, la utilización de más puntos produce una mayor exactitud aunque no es
conveniente dada la cantidad de evaluaciones funcionales y el aumento del error de redondeo. Las fórmulas
más comunes son las de 2, 3 y 5 puntos, que veremos con más detenimiento.
Fórmulas de 2 puntos: n = 1
Supongamos que x0 ∈ (a, b), donde f ∈ C 2 [a, b], y que x1 = x0 + h para algún h ̸= 0 suficiente-
mente pequeño para asegurarnos que x1 ∈ [a, b]. Construimos el primer polinomio de Lagrange para f
determinado por x0 y x1 con su término de error,
x − x1 x − x0 (x − x0 )(x − x1 )
f (x) = f (x0 ) + f (x1 ) + f ′′ (ξx )
x0 − x1 x1 − x 0 2
para cierto ξx ∈ [a, b]. Sustituyendo x1 = x0 + h,
x − x0 − h x − x0 (x − x0 )(x − x0 − h)
f (x) = f (x0 ) + f (x0 + h) + f ′′ (ξx )
−h h 2
Al diferenciar, obtenemos,
−1 1 ( ) (x − x0 )(x − x0 − h) 2(x − x0 ) − h
f ′ (x) = f (x0 ) + f (x0 + h) + Dx f ′′ (ξx ) + f ′′ (ξx )
h h 2 2
de donde, tomando x = x0 tenemos,
f (x0 + h) − f (x0 ) h ′′
f ′ (x0 ) = − f (ξx0 )
h 2
Para valores pequeños de h podemos utilizar (f (x0 + h) − f (x0 ))/h para aproximar f ′ (x0 ) con un
error acotado por M |h|/2 donde M es una cota de |f ′′ (x)| en [a, b]. Esta fórmula se llama fórmula de la
diferencia progresiva si h > 0, y fórmula de la diferencia regresiva si h < 0.
Supongamos que x0 ∈ (a, b), donde f ∈ C 3 [a, b], x1 = x0 + h y x2 = x0 + 2h para algún h ̸=

0 suficientemente pequeño para asegurarnos que x1 , x2 ∈ [a, b]. Construimos el primer polinomio de
Lagrange para f determinada por x0 , x1 y x2 con su término de error,
(x − x0 )(x − x1 )(x − x2 )
f (x) = f (x0 )L0 (x) + f (x1 )L1 (x) + f (x2 )L2 (x) + f ′′′ (ξx )
6
para cierto ξx ∈ [a, b], donde,
(x − x1 )(x − x2 ) 2x − x1 − x2
L0 (x) = ⇒ L′0 (x) =
(x0 − x1 )(x0 − x2 ) (x0 − x1 )(x0 − x2 )
(x − x0 )(x − x2 ) ′ 2x − x0 − x2
L1 (x) = ⇒ L1 (x) =
(x1 − x0 )(x1 − x2 ) (x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 ) ′ 2x − x0 − x1
L2 (x) = ⇒ L2 (x) =
(x2 − x0 )(x2 − x1 ) (x2 − x0 )(x2 − x1 )
de modo que para x = xj para j = 0, 1, 2, tenemos,
2x − x1 − x2 2x − x0 − x2 2x − x0 − x1
f ′ (xj ) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )
2
∏
1
+ f ′′′ (ξx ) (xj − xk )
6
k=0
k̸=j
Tomando ahora x1 = x0 + h y x2 = x0 + 2h, la fórmula anterior queda:
para xj = x0 ,
1[ 3 1 ] h2
f ′ (x0 ) = − f (x0 ) + 2f (x0 + h) − f (x0 + 2h) + f ′′′ (ξ0 )
h 2 2 3
para xj = x1 = x0 + h,
1[ 1 1 ] h2
f ′ (x0 + h) = − f (x0 ) + f (x0 + 2h) − f ′′′ (ξ1 )
h 2 2 6
para xj = x2 = x0 + 2h,
1 [1 3 ] h2
f ′ (x0 + 2h) = f (x0 ) − 2f (x0 + h) + f (x0 + 2h) + f ′′′ (ξ2 )
h 2 2 3
Por razones de comodidad, podemos sustituir en la segunda fórmula x0 por x0 + h y en la tercera

fórmula x0 por x0 + 2h, obteniendo,
1 [ ] h2
f ′ (x0 ) = − 3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) + f ′′′ (ξ0 )
2h [ ] h2 3
1
f ′ (x0 ) = − f (x0 − h) + f (x0 + h) − f (ξ1 )′′′
2h [ 6 ]
1 h2
f ′ (x0 ) = f (x0 − 2h) − 4f (x0 − h) + 3f (x0 ) + f ′′′ (ξ2 )
2h 3
donde la primera y última fórmula son iguales sin más que sustituir h por −h. Por tanto, en realidad hay
dos fórmulas de 3 puntos, la fórmula de diferencias finitas centrada:
1 [ ] h2
f ′ (x0 ) = − f (x0 − h) + f (x0 + h) − f ′′′ (ξ0 ) ,
2h 6
que emplea datos a ambos lados de x0 y por ello tiene un error aproximadamente la mitad que la otra
fórmula, ya sea para h > 0 o para h < 0, que emplea únicamente datos a un lado de x0 :
1 [ ] h2
f ′ (x0 ) = − 3f (x0 ) + 4f (x0 + h) − f (x0 + 2h) + f ′′′ (ξ1 )
2h 3
Para obtener las fórmulas de 5 puntos se evalúa la función en otros dos puntos más, por ejemplo,
x0 − 2h, x0 − h, x0 , x0 + h y x0 + 2h, pero cuyo término de error tiene la forma θ(4). Una de estas
fórmulas es,
1 [ ] h4
f ′ (x0 ) = f (x0 − 2h) − 8f (x0 − h) + 8f (x0 + h) − f (x0 + 2h) + f v (ξ)
12h 30
donde ξ está entre x0 − 2h y x0 + 2h.
5.2.2. Derivadas de orden superior.
Se pueden obtener fórmulas para aproximar derivadas de orden superior de una función en un punto
x0 utilizando exclusivamente los valores de la función en varios puntos. La obtención de estas fórmulas
por el procedimiento anterior es muy laboriosa, pero usando desarrollos de Taylor alrededor de un punto
se pueden obtener dichas fórmulas de modo más sencillo.
Veamos un ejemplo: hagamos el desarrollo de Taylor de grado 3 de una funcióm f en un entorno de

x0 y evaluemos en x0 − h y x0 + h.
1 1 1
f (x0 + h) = f (x0 ) + f ′ (x0 )h + f ′′ (x0 )h2 + f ′′′ (x0 )h3 + f iv (ξ1 )h4
2 6 24
1 1 1
f (x0 − h) = f (x0 ) − f ′ (x0 )h + f ′′ (x0 )h2 − f ′′′ (x0 )h3 + f iv (ξ−1 )h4
2 6 24
donde x0 − h < ξ−1 < x0 < ξ1 < x0 + h.
Sumando,
h4 ( iv )
f (x0 + h) + f (x0 − h) = 2f (x0 ) + f ′′ (x0 )h2 + f (ξ1 ) + f iv (ξ−1 )
24
y despejando f ′′ (x0 ),
1[ ] h2 ( iv )
f ′′ (x0 ) = 2
f (x0 − h) − 2f (x0 ) + f (x0 + h) − f (ξ1 ) + f iv (ξ−1 )
h 24
Suponiendo que f iv es continua en [x0 − h, x0 + h], por el teorema del valor intermedio, existe un ξ entre
ξ−1 y ξ1 con
1 iv
f iv (ξ) = (f (ξ−1 ) + f iv (ξ1 ))
2
por tanto,
1[ ] h2 iv
f ′′ (x0 ) = f (x 0 − h) − 2f (x 0 ) + f (x 0 + h) − f (ξ)
h2 12
para ξ entre ξ−1 y ξ1 .

Apuntes MIII

Cargado por

Copyright:

Formatos disponibles

Apuntes MIII

Cargado por

Información del documento

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

Apuntes MIII

Cargado por

Copyright:

Formatos disponibles

MATEMÁTICAS III

GRADO INGENIERÍA QUÍMICA

Departamento de Matemática Aplicada

Mabel Asensio Sevilla

0.2. Error absoluto y relativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

0.3. Errores de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1. Ecuaciones y sistemas de ecuaciones no lineales 9

1.1. Localización y separación de raı́ces de una ecuación. . . . . . . . . . . . . . . . . . . . . . 10

1.2. Ecuaciones no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.2.1. Método de bisección. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.2.2. El método de punto fijo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.2.3. El método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.2.4. Modificaciones del método de Newton. . . . . . . . . . . . . . . . . . . . . . . . . . 20

1.2.5. Método de la secante. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

1.3. Sistemas de ecuaciones no lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

1.3.1. Método de punto fijo en varias variables. . . . . . . . . . . . . . . . . . . . . . . . . 22

1.3.2. Método de Newton en varias variables. . . . . . . . . . . . . . . . . . . . . . . . . . 23

2. Sistemas de ecuaciones lineales 27

2.1. Generalidades sobre matrices y vectores . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.2. Métodos directos de resolución de sistemas de ecuaciones lineales . . . . . . . . . . . . . . 32

2.2.1. Matrices triangulares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.2.2. Eliminación gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.2.3. Técnicas de pivotaje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.2.5. Matrices especiales: factorización LDLt , Cholesky . . . . . . . . . . . . . . . . . . 39

2.3. Métodos iterativos de resolución de sistemas de ecuaciones lineales . . . . . . . . . . . . . 40

2.3.1. Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.3.2. Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

2.3.3. Métodos de relajación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

2.3.4. Control de parada de las iteraciones . . . . . . . . . . . . . . . . . . . . . . . . . . 45

2.3.5. Resultados de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.1. Interpolación polinómica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1.2. Tipo de función interpoladora . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1.3. Existencia y unicidad del polinomio interpolador . . . . . . . . . . . . . . . . . . . 52

3.1.4. Métodos de cálculo del polinomio interpolador. . . . . . . . . . . . . . . . . . . . . 53

3.1.5. Error de interpolación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.2. Interpolación de Hermite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.2.1. Ejemplo sencillo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.2.2. Problema de Hermite generalizado . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.2.3. Caso particular: el polinomio de Taylor . . . . . . . . . . . . . . . . . . . . . . . . 58

3.2.4. Método de las diferencias divididas de Newton generalizado . . . . . . . . . . . . . 58

3.2.5. Ejemplo sencillo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.1.1. Conjunto de abscisas de aproximación . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.1.2. Funciones básicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.1.3. Medida de la magnitud del error: normas funcionales . . . . . . . . . . . . . . . . . 62

4.2. Aproximación por mı́nimos cuadrados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.2.1. Definición del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.2.2. Productos escalares asociados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.2.3. Ecuaciones normales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.2.4. Un ejemplo sencillo: la recta de regresión . . . . . . . . . . . . . . . . . . . . . . . 66

4.3.1. Ortogonalización de Gram-Schmidt . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

5. Integración y derivación numéricas 69

5.1. Integración numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

5.1.1. Integración vı́a interpolación. Fórmulas de Newton-Cotes . . . . . . . . . . . . . . 70

5.1.2. Método de los coeficientes indeterminados . . . . . . . . . . . . . . . . . . . . . . . 72

5.1.3. Cambio de intervalo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

5.1.4. Cuadratura gaussiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.2. Derivación numérica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.2.1. Derivadas primeras. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.2.2. Derivadas de orden superior. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79