Resumen Calculo
Resumen Calculo
Resumen Calculo
Estimación de Errores
Es importante poder distinguir entre el nuevo error producido durante cada calculo particular (error de
procesamiento) y el error heredado (propagado) desde los daros durante todos los cálculos.
Error en los datos iniciales: los datos iniciales pueden estar influenciados por errores sistemáticos o
perturbaciones temporales.
Error de redondeo: la limitación de los numero de PF en maquina lleva asociada una perdida frecuente de
información que puede o no según el caso ser importante.
Error de Truncamiento: este tipo de error ocurre un proceso límite es truncado antes de llegar al valor
límite.
Error absoluto
Sistemas Numéricos en PC
La magnitud con la que el digito a contribuye al valor del número depende de la posición de a en el
número.
Se pueden considerar otros sistemas posicionales con base distinta de 10. Todo número 𝛽 ≥ 2 puede
usasrse como base. Todo número real positivo puede expresarse de la forma:
1
Cuanto más pequeña es la base más simples son las reglas por eso la mayoría de las PC trabajan en base
binaria.
Sea a un entero dado en un sistema de numeración con base α. Es posible determinar su representación
en un sistema numérico de base β.
a = bn 𝛽 𝑛 + bn−1 𝛽 𝑛−1 + ⋯ + 𝑏1 𝛽 + 𝑏0 𝛽 0
Los cálculos de conversión se realizan en base α para que β este expresada en la representación. La
conversión se realiza por divisiones sucesivas de a con β.
Es la parte fraccional
Los dígitos 𝑏−𝑖 se obtienen de la parte entera de las sucesivas multiplicaciones de c con β.
Dado que una parte fraccionaria finita en una base α no implica necesariamente una parte fraccionaria
finita en una base β.
Una computadora es construida para manejar piezas de información de tamaño fijo llamado palabras. La
longitud de la palabra es definida por el número de dígitos de esto. Los enteros se pueden representar de
forma exacta siempre que la longitud de la palabra sea suficiente para su representación.
Punto Fijo: utiliza una cantidad fijas de dígitos después de la coma, pierde precisión en las operaciones
cuando los números se van de rango. Contra: los resultados intermedios de las operaciones pueden no
estar incluidos en el conjunto de número. Para que un sistema de punto fijo tenga una implementación
exitosa todos los resultados (incluidos los intermedios) deben estar incluidos en el conjunto.
𝐶𝑜𝑛 1 ≤ 𝑚 ≤ 𝛽
La parte fraccional m se denomina mantisa, q es el exponente y β es la base. Dentro de la maquina el
número de dígitos para m y q está limitado por la longitud de palabra del procesador.
̅𝛽 𝑒 donde 𝑚
Entonces solo es posible representar números de PF de la forma:𝑎̅ = ±𝑚 ̅ es la mantisa
redondeada a p dígitos y e está limitado a un rango finito.
Los números de este conjunto F caracterizado por la base β, la precisión p y los valores 𝑒𝑚𝑖𝑛 𝑦 𝑒𝑚𝑎𝑥 tiene
exactamente: 2 ∗ (𝛽 − 1) ∗ 𝛽 𝑝−1 (𝑒𝑚𝑎𝑥 − 𝑒𝑚𝑖𝑛 + 1) + 1 elementos.
2
Si el número a representar es mayor en magnitud que el número más grande de F a no puede ser
representado y ocurre un overflow y si a es menor que el menor número de F a no puede ser
representado y ocurre un underflow.
El espaciamiento entre los números se caracteriza por el épsilon de máquina, que es la distancia desde el 1
al siguiente mayor que él. La existencia de la épsilon de la máquina es una consecuencia de la precisión
finita de la aritmética en punto flotante.
Teorema 1
En un sistema de punto flotante F=F (β, p,𝑒𝑚𝑎𝑥 , 𝑒𝑚𝑖𝑛 ) cada número real dentro del rango de F puede ser
representado con un error relativo, que no excede la unidad de redondeo, μ, que se define como:
1
𝜇 1) ∗ 𝛽 −𝑝+1
2
2)𝛽 −𝑝+1
Para suma y resta una cota para el error absoluto en el resultado final es igual a la suma de las cotas de los
errores absolutos de los operandos:
𝑚 𝑚 𝑚 ∆𝑦 ∆𝑥𝑖
𝑦 = 𝑥1 1 . 𝑥2 2 … 𝑥𝑛 𝑛 → | | = ∑|𝑚𝑖 | | |
𝑦 𝑥𝑖
Para estudiar la propagación de errores en forma más general se centrara en el estudio de expresiones no
lineales.
Numero de condición
Nos da una idea de cuan sensible son los datos de salida con respecto a variaciones en los datos de
entrada. Si pequeñas variaciones en los datos de entrada resultan en grandes cambios en los datos de
salida el problema está mal condicionada. Se vio que |𝑓 ′ (𝜁)| puede ser interpretada como una medida de
la sensibilidad de f(x) con respecto a las perturbaciones sobre x. La proporción de le perturbación relativa
en f(x) y x proporciona más información.
3
|𝑓(𝑥 + ∆𝑥)|
|𝑓(𝑥)| 𝑥 ∗ 𝑓′(𝑥)
𝐾 = lim =| |
|∆𝑥|→0 |∆𝑥 | 𝑓(𝑥)
|𝑥 |
El número de condición es una propiedad de un problema numérico y no depende del algoritmo utilizado
en la resolución.
Capítulo 2
Método de Bisección
Es el método más simple de aplicar, con convergencia segura y se basa en el teorema del valor medio.
Por lo tanto si f es continua y definida en [a; b] y se verifica f(a)*f (b) <0 el teorema del valor medio
asegura la existencia de una raíz p en el intervalo.
Esquema iterativo
[A; b]=[𝑎𝑜 ; 𝑏𝑜 ]
𝑎𝑜 +𝑏𝑜
𝑝𝑜 =
2
1) 𝑓(𝑝𝑜 ) = 0 → 𝑝𝑜 𝑒𝑠 𝑙𝑎 𝑟𝑎𝑖𝑧
2) 𝑓(𝑎) ∗ 𝑓(𝑝𝑜 ) < 0 → 𝑒𝑙 𝑛𝑢𝑒𝑣𝑜 𝑖𝑛𝑡𝑒𝑟𝑣𝑎𝑙𝑜 [𝑎1 ; 𝑏1 ] = [𝑎; 𝑝𝑜 ]
3) 𝑓(𝑝𝑜 ) ∗ 𝑓(𝑏) < 0 → 𝑒𝑙 𝑛𝑢𝑒𝑣𝑜 𝑖𝑛𝑡𝑒𝑟𝑣𝑎𝑙𝑜 [𝑎1 ; 𝑏1 ] = [𝑝𝑜 ; 𝑏]
𝑎1 +𝑏1
En los casos 2 y 3 se sigue iterando: 𝑝1 =
2
lim 𝑃𝑛 = 𝑃
𝑛→∞
Sea {[a;b]} una secuencia de intervalos cerrados y encajados de forma tal que:
Además lim (𝑏𝑛 − 𝑎𝑛 ) = 0 entonces existe un único punto p ∈ [𝑎𝑛 ; 𝑏𝑛 ] para todo n que pertenece a los
𝑛→∞
naturales.
4
∩∞
𝑛=0 [𝑎𝑛 ; 𝑏𝑛 ] = {𝑝}
Para poder trabajar con el método, es necesario establecer un criterio de corte ya que debido a la
aritmética finita y errores de redondeo es casi imposible obtener un pn que cumpla exactamente f(pn)=0.
2)|𝑝𝑛 − 𝑝| < 𝜀
No conozco p.
3)𝑓(𝑝𝑛 ) < 𝜀
Que la imagen de pn sea muy cercana a cero no me asegura que estoy cerca de la raíz.
𝑝𝑛 −𝑝𝑛−1
4)| |<𝜀
𝑝𝑛
Es la opción usada.
El método de bisección tiene convergencia global no depende de la semilla utilizada siempre que f sea
continua en el intervalo y f(a)*f(b)<0.
𝑏𝑛 − 𝑎𝑛
𝜀𝑘 = |𝑝𝑘 − 𝑝| ≤
2𝑘
Lo que implica que el lim 𝜀𝑘 = 0. A partir de esta fórmula puede obtenerse el número de pasos para que
𝑘→∞
el error sea menor que un determinado valor.
Se basa en que dada una 𝑓: [𝑎; 𝑏] → ℜ, el problema de f(x)=0 puede transformarse en el problema
análogo de 𝑥 − ∅(𝑥) = 0, donde ∅: [𝑎; 𝑏] → ℜ debe cumplir que ∅(𝛼) = 𝛼 siempre que 𝑓(𝛼) = 0
Aproximar los ceros de una función se transforma en el problema de encontrar los puntos fijos del mapeo
∅.
Esquema iterativo: 𝑥𝑛+1 = ∅(𝑥𝑛 )
∅ No es único.
Sea ∅ ∈ 𝐶 1 [𝑎; 𝑏]
5
1) Si ∅(𝑥) ∈ [𝑎; 𝑏], ∀ 𝑥 ∈ [𝑎; 𝑏] → ∅ tiene al menos un punto fijo en el intervalo [𝑎; 𝑏].
2) Si ∅′(𝑥) esta definida en el intervalo (𝑎, 𝑏) y |∅′(𝑥)| < 1, ∀ 𝑥 ∈ (𝑎; 𝑏) → el punto fijo α en [a;b]
de ∅ es único.
Demostración
Si ∅(𝑎) = 0 𝑜 ∅(𝑏) = 0 el primer consecuente queda demostrado. Si esto no sucede, entonces se cumple:
∅(𝑎) ∈ (𝑎; 𝑏] 𝑦 ∅(𝑏) ∈ [𝑎; 𝑏). Por lo tanto la función 𝑓(𝑥) = 𝑥 − ∅(𝑥) tiene la propiedad
Al aplicar el teorema del valor medio, se asegura la existencia de un α∈ (𝑎; 𝑏) tal que 𝑓(𝛼) = 0entonces 𝛼 = ∅(𝛼)
con lo que se asegura la existencia del punto fijo.
Por el teorema del valor medio, es posible identificar un 𝑐 ∈ (𝛼1 , 𝛼2 ) tal que:
∅(𝛼1 ) − ∅(𝛼2 )
∅′ (𝑐) =
𝛼1 − 𝛼2
𝛼1 − 𝛼2
∅′ (𝑐) =
𝛼1 − 𝛼2
∅′ (𝑐) = 1
∅ Tiene un único punto fijo α en [a;b] y la secuencia {xn} converge a α ∀ 𝑥0 ∈ [𝑎; 𝑏]. Además se verifica que:
𝑥𝑛+1 − 𝛼
lim = ∅′(𝛼)
𝑛→∞ 𝑥𝑛 − 𝛼
Este teorema asegura la convergencia de la secuencia xn a la raíz α para cualquier punto que se elija como semilla
dentro del intervalo [a;b]. Asegura la convergencia global.
Teorema de OSTROWSKI
Sea α un punto fijo de la función ∅, la que es continua y derivable en una vecindad de α. Si |∅′(𝛼)| < 1entonces
existe 𝛿 > 0 tal que la secuencia xn converge a α para cualquier valor 𝑥0 que cumpla |𝑥0 − 𝛼 | < 𝛿. Este teorema
asegura la convergencia local del método.
6
Método de Newton-Rapshon
El método de Newton Rapshon es rápido y simple. Pero como utiliza f(x) y f’(x), se usa solo en casis donde f’(x) es
fácilmente calculable.
(𝑥𝑖+1 − 𝑥𝑖 )2 ′′
𝑓(𝑥𝑖+1 ) = 𝑓(𝑥𝑖 ) + (𝑥𝑖+1 − 𝑥𝑖 ) ∗ 𝑓 ′ (𝑥𝑖 ) + 𝑓 (𝑥𝑖 ) + ⋯
2
Si 𝑥𝑖+1 es raíz 𝑓(𝑥𝑖+1 ) = 0. Entonces:
(𝑥𝑖+1 − 𝑥𝑖 )2 ′′
0 = 𝑓(𝑥𝑖 ) + (𝑥𝑖+1 − 𝑥𝑖 ) ∗ 𝑓 ′ (𝑥𝑖 ) + 𝑓 (𝑥𝑖 ) + ⋯
2
Si 𝑥𝑖+1 es muy cercano a 𝑥𝑖 los términos de orden superior se vuelven despreciables.
𝑓(𝑥)
Para estudiar la convergencia: 𝑔(𝑥) = 𝑥 − . Recordando que la derivada de g proporciona la descripción de la
𝑓′(𝑥)
convergencia a través del método de punto fijo.
𝑓(𝑥)∗𝑓′′ (𝑥)
Derivando obtenemos: 𝑔′ (𝑥) = en la raíz 𝑥 = 𝛼 el valor de 𝑓(𝛼) = 0 por definición. Por lo tanto para el
[𝑓′ (𝑥)]2
proceso iterativo de Newton Rapshon 𝑔′ (𝛼) = 0 lo que implica convergencia óptima.
Cuando se analiza la convergencia a través del método de punto fijo es necesario calcular cuan alejado del valor α se
encuentra xn por lo tanto el análisis recaerá:
𝜀𝑛+1 = 𝑥𝑛+1 − 𝛼
𝜀𝑛+1 = 𝑔(𝑥𝑛 ) − 𝑔(𝛼)
𝜀𝑛 = 𝑥𝑛 − 𝛼 → 𝑥𝑛 = 𝜀𝑛 + 𝛼
𝜀𝑛+1 = 𝑔(𝑥𝑛 ) − 𝑔(𝛼)
𝜀𝑛+1 = 𝑔(𝜀𝑛 + 𝛼) − 𝑔(𝛼)
Como 𝑔′ (𝛼) = 0 se elimina ese término y en las cercanías de α se pueden despreciar los términos de orden superior
con lo que: 𝜀𝑛+1 𝛼 𝜀𝑛2 .
Método de la Secante
Una de las principales complicaciones del método de Newton Rapshon es la obtención analítica de la derivada. Este
proceso se simplifica si en vez de tomar la derivada se utiliza una aproximación numérica, es decir, la recta secante
en vez de la tangente.
7
𝑓(𝑥𝑛−1 ) − 𝑓(𝑥𝑛 )
𝑓 ′ (𝑥𝑛 ) = lim
𝑥𝑛−1 →𝑥𝑛 𝑥𝑛−1 − 𝑥𝑛
𝑓(𝑥𝑛−1 ) − 𝑓(𝑥𝑛 )
𝑓 ′ (𝑥𝑛 ) ≈
𝑥𝑛−1 − 𝑥𝑛
Introduciendo dicha aproximación de la derivada en el esquema iterativo de Newton Rapshon y sumando 1 a cada
uno de los subíndices llegamos a:
Global. Aprovecha la secuencia de rectas generadas por el método de la secante pero para la selección de los puntos
se basa en el método de bisección.
Dada una f(x) continúa y dados 𝑥0 𝑦 𝑥1 tales que 𝑓(𝑥0 ) ∗ 𝑓(𝑥1 ) < 0 se construye la recta secante que pasa por estos
puntos. Entonces la raíz de f se encuentra en el intervalo (𝑥0 , 𝑥2 ) siempre que 𝑓(𝑥0 ) ∗ 𝑓(𝑥2 ) < 0 si no se elige el
intervalo (𝑥2 , 𝑥1 ) para seguir iterando.
Análisis de Convergencia
Para caracterizar la eficiencia de un algoritmo iterativo, se toma en cuenta: la cantidad de flops y la velocidad de
convergencia.
Para conocer esta velocidad se comparan los distintos cocientes: ∆𝑥𝑘+2 = 𝑥𝑘+2 − 𝑥𝑘+1 ; ∆𝑥𝑘+1 = 𝑥𝑘+1 − 𝑥𝑘
Si los cocientes estabilizan en torno a 𝐶𝐿 , con 0 < 𝐶𝐿 < 1, el método tiene convergencia lineal, y la tasa de
convergencia es 𝐶𝐿 . Si 𝐶𝐿 > 1los iterados divergen.
∆𝑥𝑘+2
Si los cocientes tienden a cero y 𝑐𝑜𝑛 1 < 𝛼 < 2 estabilizan el método tiene convergencia superlineal.
∆𝑥𝑘+1 𝛼
∆𝑥𝑘+2
Si los cocientes tienden a cero y los cocientes estabilizan, el método tiene convergencia cuadrática.
∆𝑥𝑘+1 2
El tipo de convergencia indica el trabajo necesario para alcanzar determinada precisión. En convergencia lineal cada
cierta cantidad de iteraciones se obtiene un digito adicional exacto. En convergencia cuadrática, con la misma
cantidad de iteraciones se duplica el número de decimales exactos.
Bisección y Regula Falsi convergen linealmente. Punto Fijo suele tener convergencia lineal pero en algunos casos
particulares puede tener convergencia cuadrática (NR) o superior. Newton Rapshon, si se dan las condiciones
adecuadas (raíz simple), presenta convergencia cuadrática. El método de la secante suele presentar convergencia
superlineal.
8
Fallas en los métodos Iterativos.
El único método que asegura convergencia es bisección (se basa en la completitud de ℜ).
Punto fijo debe cumplir ciertas condiciones para asegurar la convergencia. Se debe verificar el teorema de Ostrowski
en cada paso. Pero lo mas difícil es la creación del mapeo.
Fallas de Newton-Rapshon
1) Finalización prematura o breakdown: para algún xn, se cumple que f’(x)=0, no puede continuar el algoritmo.
2) La sucesión {xn} se aleja cada vez más de la raíz.
3) La última falla es que genere una sucesión oscilante.
Capítulo 3
SEL-Métodos Directos
Para dar solución a SEL de n ecuaciones y n incógnitas, se desarrollan métodos directos de resolución.
Conceptos Básicos
.
.
El sistema tiene solución única si y solo si el determinante de A es distinto de 0, es decir, ninguna fila o columna es
combinación lineal de otra fila o columna.
En caso de que la matriz de coeficientes A sea singular (determinante de A=0) el sistema puede tener infinitas
solución o ninguna dependiendo del vector de constantes.
9
Existen distintas normas matriciales:
Sea Ѵ un espacio vectorial la función ‖. ‖ de valores no negativos se llama norma si cumple los siguientes axiomas:
‖𝑣 ‖ = 0 ↔ 𝑣⃗ = 0, 𝑣 ∈ Ѵ
‖ʎ𝑣 ‖ = |ʎ|‖𝑣 ‖, ∀ 𝑥 ∈ ℜ , ∀ 𝑣 ∈ Ѵ
‖𝑢 + 𝑣 ‖ ≤ ‖𝑢‖ + ‖𝑣 ‖, ∀ 𝑣, 𝑢 ∈ Ѵ
Ѵ con una norma asociada de denomina espacio vectorial normado. Cualquier norma de Ѵ= ℜ𝑛 se llama norma
vectorial. Tres normas son comunes:
La norma 2: ‖𝑣 ‖2 = √∑|𝑣𝑖 |2
𝑛 | |
La norma infinito: ‖𝑣 ‖∞ = 𝑚𝑎𝑥𝑖=1 𝑣𝑖
Una norma matricial es cualquier norma del espacio vectorial ℜ𝑛𝑥𝑛 . En particular se consideran normas matriciales
inducidas por normas vectoriales:
‖𝐴𝑣 ‖
‖𝐴‖ = 𝑚𝑎𝑥𝑣∈Ѵ
‖𝑣 ‖
Proceso: planteo un vector unitario. Encuentro ‖𝐴𝑣 ‖ y derivo para poder maximizar. Encuentro las raíces de 𝑑 ‖𝐴𝑣 ‖
(bisección) y calculo ‖𝐴𝑣 ‖ en las raíces. El máximo valor que toma ‖𝐴𝑣 ‖ en las raíces es la norma de A.
Se expresa como el mayor valor de las sumas de los valores absolutos de los elementos de las columnas.
Se expresa como el mayor valor de las sumas de los valores absolutos de los elementos por columnas.
Condicionamiento de matrices
El principal uso de las normas matriciales es el cálculo del número de condición de una matriz (que tan sensible es a
las perturbaciones).
10
Si 𝐴 ∈ ℜ𝑛𝑥𝑛 𝑦 ∃𝐴−1 : dada A.x=b
𝐴𝑒 = 𝐴. 𝑥 − 𝐴. 𝑥̅ = 𝑏 − 𝐴. 𝑥̅ = 𝑟
𝑒 = 𝐴−1 . 𝑟
Si el e=0 r=0. Pero si r es pequeño, e no es necesariamente pequeño. Un r pequeño no garantiza que 𝑥̅ sea cercano a
x. La ventaja es que r siempre puede ser calculado a diferencia del error ya que no conozco x. sin embargo relacionar
r con el número de condición permite verificar la cercanía de 𝑥̅ respecto de x. Ver demostración!
Definimos el número de condición: 𝐾𝑃 (𝐴) = ‖𝐴‖𝑃 ‖𝐴−1 ‖𝑝 como el número de condición de la matriz A asociado a
una norma p.
1 ‖𝑟‖𝑝 ‖𝑟‖𝑝
∗ ≤ 𝐸𝑟 ≤ 𝐾𝑃 (𝐴) ∗
𝐾𝑃 (𝐴) ‖𝑏‖𝑝 ‖𝑏‖𝑝
Una matriz con un 𝐾𝑃 (𝐴) >>> 1 esta mal condicionada y un mal condicionamiento no se soluciona al multiplicar a
la matriz por un escalar.
‖𝛿𝑥 ‖𝑝 ‖𝛿𝑏‖𝑝
≤ 𝐾𝑃 (𝐴).
‖𝑥 ‖𝑝 ‖𝑏‖𝑝
𝛿𝑥
Por efecto de redondeos, la solución de A.x=b NO SERA EXACTA. La solución se escribe 𝑥 + = 𝑑 + 𝛿𝑏 con 𝛿𝑏
𝐴(𝑥+𝛿𝑥)
muy pequeños. Si 𝐾𝑃 (𝐴) >>> 1 los 𝛿𝑥 pueden no ser pequeños.
Eliminación Gaussiana
El método ejecuta eliminaciones sucesivas hasta conseguir una matriz triangular superior (elimina los elementos
debajo de la diagonal principal). El sistema resultante se resuelve por sustitución hacia atrás.
11
𝐸1 = 𝐴11 𝐴12 … 𝐴1𝑛 𝑥1 𝑏1
𝐸2 = 𝐴21 𝐴22 … 𝐴2𝑛 ∗ 𝑥2 = 𝑏2
𝐸𝑛 = 𝐴𝑛1 𝐴𝑛2 … 𝐴𝑛𝑛 𝑥𝑛 𝑏𝑛
𝐸𝑖′ = 𝐸𝑗
𝐸𝐽′ = 𝐸𝑖
Desarrollado en forma algebraica, el método de Gauss siempre da resultados exactos. Numéricamente, la precisión
depende de la aritmética y del número de condición de la matriz de coeficientes principalmente.
Para utilizar el método es importante la elección de los pivotes. Conviene que el número con mayor valor absoluto
este en la posición 𝐴11 , o al menos que el 𝐴11 sea mayor que los otros elementos de la misma columna. Esto es
𝐴𝑗+1,𝑗
porque al eliminar los elementos de debajo de la diagonal principal se utiliza un factor: Si el valor de elemento
𝐴𝑗𝑗
superior (jj) es menor que el que se desea convertir en 0 (j+1,j) , el factor será muy grande.
Ejemplo:
Descomposición LU
Cualquier matriz 𝐴 ∈ ℜ𝑛𝑥𝑛 , bajo ciertas condiciones puede ser expresada como el producto de una matriz inferior L y
una matriz superior U.
A=L.U
Encontrar L y U se denomina factorización LU. La descomposición no es única. Las tres descomposiciones más
comunes son:
3 − 𝐿 = 𝑈 𝑇 → 𝐷𝑒𝑐𝑜𝑚𝑝𝑜𝑐𝑖𝑠𝑖𝑜𝑛 𝑑𝑒 𝐶ℎ𝑜𝑙𝑒𝑠𝑘𝑖
12
Luego de descomponer A en L y U se sigue así:
A.x=b
LU.x=b
La ventaja de LU respecto a Gauss es que una vez descompuesta A en L y U, se puede resolver el sistema con tanto
vectores b como se necesite. También autocompensa errores.
𝐷𝑒𝑠𝑐𝑜𝑚𝑝𝑜𝑐𝑖𝑠𝑖𝑜𝑛 𝐷𝑜𝑜𝑙𝑖𝑡𝑡𝑙𝑒 : En este método la matriz U se obtiene de manera similar a la eliminación Gaussiana,
pero sin la matriz aumentada. Los elementos de L son los factores utilizados para obtener U, y su diagonal principal
es de 1.
Consideraciones Generales
Una matriz de 𝐴 ∈ ℜ𝑛𝑥𝑛 es sparse si más de la mitad de sus 𝑛2 elementos son 0. Al realizar una sucesión de Caucha
de vectores 𝑥 𝑘 convergen a x, mediante métodos iterativos, se buscara solución cuando A es soparse o de gran
tamaño: los métodos directos necesitan movimiento considerable de información en las memorias de la PC haciendo
lento el proceso.
LU →autocompensa errores
Con los métodos iterativos se busca una sucesión de vectores 𝑥 𝑘 / lim 𝑥 𝑘 = 𝑥 , donde 𝑥 = {𝑥0 , 𝑥1 , … , 𝑥𝑛 }𝑇
𝑘→∞
Como la sucesión que se forma es de Cauchy para todo 𝜀 > 0 ∃ 𝑚 ∈ ℕ tal que ‖𝑥 𝑚 − 𝑥 ‖ ≤ 𝜀. El operador T se
define:
𝑥 𝑘+1 = 𝐵𝑥 𝑘 + 𝑓
𝑓 = (𝐼 − 𝐵)𝐴−1 𝑏
13
La convergencia solo ocurrirá si se elige bien B y si A cumple con ciertas condiciones:
Sea s(A) el conjunto de los autovalores de la matriz A ∈ ℜ𝑛𝑥𝑛 . Se define como el radio espectral como:
𝜌(𝐴) = 𝑚𝑎𝑥ʎ ∈𝑠(𝐴) |ʎ|
Teorema
Sea 𝐴 ∈ ℜ𝑛𝑥𝑛
1 1
Así, si 𝜌(𝐴) < 1, entonces la matriz (𝐼 − 𝐴) es invertible y también: ≤ ‖(𝐼 − 𝐴)−1 ‖ ≤ donde ‖. ‖ es una
1+‖𝐴‖ 1−‖𝐴‖
norma inducida de la matriz de forma tal que ‖𝐴‖ < 1.
Teorema
Sea 𝑓 ∈ ℜ𝑛 tal que 𝑓 = (𝐼 − 𝐵)𝐴−1 𝑏 𝑥 𝑘 converge a x satisfaciendo 𝐴. 𝑥 = 𝑏 para cualquier 𝑥 0 semilla si y solo si
𝜌(𝐵) < 1 . Este último teorema da una condición general para B talque la iteración converja.
𝜀 𝑘+1 = 𝑥 𝑘+1 − 𝑥
𝜀 𝑘+1 = 𝐵𝑥 𝑘 + 𝑓 − 𝑥
𝜀 𝑘+1 = 𝐵𝑥 𝑘 + (𝐼 − 𝐵)𝐴−1 𝑏 − 𝑥
𝜀 𝑘+1 = 𝐵𝜀 𝐾
Trabajando la expresión anterior 𝜀 𝐾 = 𝐵𝐾 . 𝜀 0 Sabemos que ‖𝐵‖ < 1pero no como hallar B.
Una forma general de construcción de matrices para métodos iterativos es la división aditiva de A.
𝐴 = 𝑃 − 𝑁; 𝑃 , 𝑁 ∈ ℜ𝑛𝑥𝑛 , ∃ 𝑃−1
𝑥 𝑘+1 = 𝑃−1 𝑁. 𝑥 𝑘 + 𝑃 −1 . 𝑏
14
Método de Jacobi o de Iteraciones Simultáneas.
Se agrega la condición de que todos los 𝐴𝑖𝑖 ≠ 0. Entonces 𝐴. 𝑥 = 𝑏 se puede expresar como:
𝑛−1
1
𝑥𝑖 = [𝑏 − ∑ 𝐴𝑖𝑗 𝑥𝑗 ] 𝑖 = 0,1,2, …
𝐴𝑖𝑖 𝑖
𝑗=0
𝑗≠𝑖
En secuencia iterativa:
𝑛−1
𝑘+1
1
𝑥𝑖 = [𝑏 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘 ] 𝑖 = 0,1,2, …
𝐴𝑖𝑖 𝑖
𝑗=0
𝑗≠𝑖
D es la matriz diagonal. La matriz L es una matriz triangular inferior tal que 𝐿𝑖𝑗 = −𝐴𝑖𝑗 𝑝𝑎𝑟𝑎 𝑖 > 𝑗 𝑦 𝐿𝑖𝑗 =
0 𝑝𝑎𝑟𝑎 𝑖 ≤ 𝑗. La matriz U es una matriz triangular superior tal que 𝑈𝑖𝑗 = −𝐴𝑖𝑗 𝑝𝑎𝑟𝑎 𝑖 < 𝑗 𝑦 𝑈𝑖𝑗 = 0 𝑝𝑎𝑟𝑎 𝑖 ≥ 𝑗.
𝐵 = 𝐵𝐽 = 𝐷 −1 ∗ (𝐿 + 𝑈)
Método de Gauss-Sede.
Es una alternativa al método de Jacobi. Se utilizan las versiones actualizadas de las variables.
𝑖−1 𝑛−1
𝑘+1
1
𝑥𝑖 = [𝑏 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘+1 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘 ] 𝑖 = 0,1,2, …
𝐴𝑖𝑖 𝑖
𝑗=0 𝑗=𝑖+1
En forma matricial se puede escribir:𝐷𝑥 𝑘 = 𝑏 + 𝐿𝑥 𝑘+1 + 𝑈𝑥 𝑘 donde D, L, U son las mismas matrices que las de
Jacobi. En Gauss-Seidel se implementa la siguiente división aditiva: 𝑃 = 𝐷 − 𝐿, 𝑁 = 𝑈.
𝐵𝐺 = (𝐷 − 𝐿)−1 . 𝑈
Tanto la iteración con B como con P o la versión lineal dan los mismos iterados y convergen al mismo x.
Pero él es que incluye a B es solo teórico, ya que para obtener f se requiere contar con la inversa de A, y teniendo la
inversa de A, simplemente 𝑥 = 𝐴−1 𝑏. B solo sirve para el análisis de convergencia. La forma iterativa a usar es la que
utiliza a la matriz P obtenida según la división aditiva del método.
Jacobi y Gauss Seidel se utilizan para resolver SEL mediante iteraciones sobre un vector inicial (semilla).
15
Método de Sobrerrelajacion de Jacobi (JOR)
𝑛−1
𝜔
𝑥𝑖 𝑘+1 = [𝑏 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘 ] + (1 − 𝜔)𝑥𝑖𝑘 𝑖 = 0,1,2, …
𝐴𝑖𝑖 𝑖
𝑗=0
𝑗≠𝑖
Gauss-Seidel
𝑖−1 𝑛−1
𝜔
𝑥𝑖 𝑘+1 = [𝑏 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘+1 − ∑ 𝐴𝑖𝑗 𝑥𝑗 𝑘 ] + (1 − 𝜔)𝑥𝑖𝑘 𝑖 = 0,1,2, …
𝐴𝑖𝑖 𝑖
𝑗=0 𝑗=𝑖+1
Para saber qué sistema será convergente con algún método iterativo, solo se necesita analizar la matriz de
coeficientes, independientemente de b:
A es diagonal dominante. Si A es diagonal dominante, los métodos de Jacobi y Gauss-Seidel convergen (Se demuestra
por teorema de punto fijo). Para terminar las secuencias iterativas se usa el error relativo. En este caso se intenta
estabilizar la convergencia de 𝐴𝑥 𝑘 = 𝑏 𝑘 con b.
‖𝑏 − 𝐴𝑥 𝑘 ‖
<𝜀
‖𝑏‖
Refinamiento Iterativo
Es uno de los métodos más básicos pero con mayor potencia de cálculo. Si 𝑥 0 es la solución aproximada de 𝐴. 𝑥 = 𝑏,
obtenido por algún método directo o iterativo lento sin llegar a la condición de convergencia.
Sistema inicial
𝐴. 𝑥 = 𝑏
Solución Aproximada
𝐴𝑥 0 = 𝑏 𝑜
16
Restando la solución aproximada menos el sistema inicial: 𝐴(𝑥 𝑜 − 𝑥) = 𝑏 𝑜 − 𝑏
𝐴𝑒 0 = 𝑟 0
Con este sistema se puede calcular el 𝑒 0 , y se obtiene una nueva solución (solución refinada) respecto a la solución
anterior.
𝐴𝑒 𝑘 = 𝑟 𝑘
𝑥 𝑘+1 = 𝑥 𝑘 + 𝑒 𝑘
𝑟 𝑘+1 = 𝑏 − 𝐴. 𝑥 𝑘+1
1) La norma de 𝑟 𝑘 es mayor que la norma de 𝑟 𝑘−1 , ósea no se puede satisfacer. Esto puede ocurrir por el
método en sí o por la aritmética utilizada.
2) El residuo relativo del paso k es menor a una tolerancia :
‖𝑏 − 𝐴𝑥 𝑘 ‖
<𝜀
‖𝑏‖
Capítulo 5.
Autovalores
Los autovalores dan información acerca del comportamiento de sistemas modelados por medio de una matriz u
operador.
Son útiles en cálculos de análisis de resonancia, inestabilidad y tasas de crecimiento o decrecimiento por ejemplo.
Introducción
En el estudio de matrices 𝐴 ∈ ℜ𝑛𝑥𝑛 son importantes ciertos vectores especiales cuyas direcciones no cambian
cuando se multiplican por A si no que solo sufren estiramientos o contracciones.
Si 𝐴 ∈ ℜ𝑛𝑥𝑛 , 𝑥 ≠ 0, ʎ ∈ ℜ
𝐴𝑥 = ʎ𝑥
17
𝑝(ʎ) = 𝑎1 ∗ (ʎ1 − ʎ) ∗ (ʎ2 − ʎ) … (ʎn − ʎ) Usando la relación entre raíces y coeficientes de una ecuación queda:
Teoremas de Gerschgorin
Con estos teoremas se obtienen regiones en las que se encuentran los autovalores. En general (sin condiciones):
Sea 𝐴 ∈ ₵𝑛𝑥𝑛 𝑦 𝑛 ≥ 2. Los discos de Gerschgorin 𝐷𝑖 , 𝑖 = 1,2, … , 𝑛, la matriz A se definen como las regiones circulares
cerradas:
Teorema 1
Sea 𝐴 ∈ ₵𝑛𝑥𝑛 𝑦 𝑛 ≥ 2. Todos los autovalores de la matriz A se ubican en la región 𝐷 = ⋃𝑛𝑖=1 𝐷𝑖 , 𝑖 = 1,2, … , 𝑛,son los
discos de Gerschgorin.
Demostración
Sea 𝑥𝑘 , 𝑐𝑜𝑛 𝑘 ∈ {1,2, … , 𝑛} la componente de x de mayor módulo entonces: |𝑥𝑗 | ≤ |𝑥𝑘 |, 𝑗 = 1,2, … , 𝑛
𝑛
| ʎ − Akk ||𝑥𝑘 | = || ∑ 𝐴𝐾𝑗 𝑥𝑗 ||
𝑗=1
𝐽≠𝐾
| ʎ − Akk ||𝑥𝑘 | = 𝑅𝐾 | 𝑥𝑗 | ≤ 𝑅𝑘 | 𝑥𝑘 |
| ʎ − Akk ||𝑥𝑘 | ≤ 𝑅𝑘 | 𝑥𝑘 |
| ʎ − Akk | ≤ 𝑅𝑘
Como conclusión ʎ se encuentra en el disco Dk , de radio R k , centrado en Akk por lo tanto ʎ ∈ 𝐷 = ⋃𝑛𝑖=1 𝐷𝑖 .
Teorema 2
Si la unión M de k discos es disjunta de los discos restantes, entonces M contiene exactamente k autovalores de A.
18
Método de la Potencia
Estos métodos sirven para aproximar los autovalores extremos ʎ1 (𝑚𝑎𝑦𝑜𝑟 𝑚𝑜𝑑𝑢𝑙𝑜) y ʎn ( 𝑚𝑒𝑛𝑜𝑟 𝑚𝑜𝑑𝑢𝑙𝑜).
𝐴 ∈ ℜ𝑛𝑥𝑛 Diagonalizable, y sea 𝑥 ∈ ℜ𝑛𝑥𝑛 la matriz de los autovectores 𝑥𝑖 , 𝑖 = 1,2, … , 𝑛. Suponiendo que:
Algoritmo Iterativo
𝑧 𝑘 = 𝐴𝑞 𝑘−1
𝑧𝑘
𝑞𝑘 =
‖𝑧 𝑘 ‖2
𝐿𝑘 = (𝑞 𝑘 )𝑇 𝑞 𝑘
𝑧 𝑘 = 𝐴𝑞 𝑘−1
𝑧 𝑘−1
𝑧𝑘 = 𝐴
‖𝑧 𝑘−1 ‖2
𝑧 𝑘−2
𝑧 𝑘 = 𝐴. 𝐴
‖𝑧 𝑘−2 ‖2
𝑧 𝑘−𝑘
𝑧 𝑘 = 𝐴𝑘
‖𝑧 𝑘−𝑘 ‖2
𝑧𝑜
𝑧 𝑘 = 𝐴𝑘
‖𝑧 𝑜 ‖2
𝑧 𝑘 = 𝐴𝑘 𝑞 0
𝐴𝑘 𝑞 0
𝑞𝑘 = 𝑘 0
‖𝐴 𝑞 ‖2
Como A es diagonalizable los autovectores 𝑥𝑖 forman una base de ℜ𝑛𝑥𝑛 . Por lo tanto puedo escribir 𝑞 0 como la
combinación lineal de los autovectores de A.
𝑛
𝑞0 = ∑ 𝛼𝑖 𝑥𝑖 𝑖 = 1,2, … , 𝑛 𝛼𝑖 ∈ ℜ
𝑖=1
Como 𝐴𝑛 𝑥𝑖 = ʎ𝑛𝑖 𝑥𝑖
19
𝑛
𝐴𝑘 𝑞 0 = 𝐴𝑘 ∑ 𝛼𝑖 𝑥𝑖
𝑖=1
𝐴𝑘 𝑞 0 = 𝐴𝑘 (α1 𝑥1 + α2 𝑥2 + ⋯ + α𝑛 𝑥𝑛 )
𝐴𝑘 𝑞 0 = ʎ𝑛𝑖 (α1 𝑥1 + α2 𝑥2 + ⋯ + α𝑛 𝑥𝑛 )
𝑛 𝛼𝑖 ʎ𝑖 𝑛
𝐴𝑘 𝑞 0 = ʎ1𝑛 α1 (𝑥1 + ∑ 𝑥𝑖 ( ) )
𝑖=2 𝛼1 ʎ1
𝐴𝑘 𝑞 0 = ʎ1𝑛 α1 (𝑥1 + 𝑦 𝑘 )
ʎ1𝑛 α1 ( 𝑥1 + 𝑦 𝑘 )
𝑞𝑘 =
‖ʎ1𝑛 α1 ( 𝑥1 + 𝑦 𝑘 )‖2
ʎ1𝑛 α1 ( 𝑥1 + 𝑦 𝑘 )
𝑞𝑘 =
‖( 𝑥1 + 𝑦 𝑘 )‖2 |ʎ1𝑛 α1 |
𝜇𝑘 ( 𝑥1 + 𝑦 𝑘 )
𝑞𝑘 =
‖( 𝑥1 + 𝑦 𝑘 )‖2
Donde 𝜇𝑘 es el símbolo de ʎ1𝑛 α1 . La sucesión {𝐿𝑘1 } converje a ʎ1 . Por lo tanto puede usarse esto como condición de
corte:
‖𝐴𝑘 𝑞 𝑘 − 𝐿𝑘1 𝑞 𝑘 ‖ ≤ 𝜀
Esquema Iterativo
20
(𝐴 − 𝜇𝐼) 𝑧 𝑘 = 𝑞 𝑘−1 → 𝑧 𝑘 = (𝐴 − 𝜇𝐼)−1 𝑞 𝑘−1
𝑧𝑘
𝑞𝑘 = 𝑘
‖𝑧 ‖2
𝜎 𝑘 = (𝑞 𝑘 )𝑇 𝑞 𝑘
Los autovectores de 𝑀𝜇 son los mismos que los de A. En cada paso debe resolverse un sistema de ecuaciones lineales
cuya matriz es 𝑀𝜇 . La factorización LU de 𝑀𝜇 es recomendada. El método es más pesado desde el punto de vista
computacional pero tiene una gran ventaja: converge a cualquier autovalor de A (el más cercano al desplazamiento
𝜇).
lim 𝑞 𝑘 = 𝑥𝑚 , lim 𝜎 𝑘 = ʎ𝑚
𝑘→∞ 𝑘→∞
Problemas de Implementación
La efectividad de estos métodos depende de la separación entre ʎ1 y los demás. Si hay dos autovalores dominantes
de igual modulo pueden presentarse tres situaciones:
1) ʎ1 = ʎ2 los dos autovalores dominantes son coincidentes. El método converge pero para un k suficientemente
grande.
𝐴𝑘 𝑞 0 ≈ ʎ1 (α1 𝑥1 + α2 𝑥2 )
Cuando k tiende a infinito 𝑞 0 converge a un vector que pertenece al subespacio generado por los autovectores
𝑥1 𝑦 𝑥2 mientas que 𝐿𝑘 converge a ʎ1 .
2)ʎ1 = −ʎ2 El autovalor dominante puede ser aproximado aplicando el método de la potencia a la matriz 𝐴2 ya
que para i=1,2,.., n [ʎ𝑖 (𝐴)]2 = ʎ𝑖 (𝐴2 ) → ʎ1 2 = ʎ2 2 . Se continua como en el caso anterior.
3) ʎ1 = ̅̅̅
ʎ2 surgen oscilaciones amortiguadas en la secuancia de vectores 𝑞 𝑘 y el método no converge.
Capítulo 7
En interpolación se construye una curva que pasa por todos los puntos. Se suponen que los puntos son distintos y
exactos. En ajuste de datos se busca que la curva pase lo más cercano de todos los puntos. Se consideran que los
datos pueden contener errores.
21
Interpolación Polinómica
Se usa dentro del rango de datos. Una de las formas más sencillas de construir una función interpoladora es elegir el
polinomio. Dado n+1 puntos de abscisas distintas existe un único polinomio p(x) de grado menor o igual a n que pasa
por todos los n+1 puntos.
Forma Normal
El planteamiento más directo se obtiene expresando el polinomio con respecto a la base canónica:
Se obtiene un sistema de ecuaciones lineales que es compatible determinado. La resolución del sistema está
afectada por el mal condicionamiento de la matriz de coeficientes (matriz de vandermonde). El condicionamiento
puede mejorarse mediante un desplazamiento del origen.
El sistema resulta ligeramente mejor condicionado, pero como la matriz sigue teniendo la misma estructura, al
aumentar el grado del polinomio volverá el mal condicionamiento.
Interpolación de Lagrange
Con Lagrange se obtiene en forma rápida una expresión polinómica de grado n-1 que interpola n puntos distintos.
𝑛
𝑃𝑛−1 (𝑥) = ∑ 𝑙𝑖 (𝑥) ∗ 𝑦𝑖
𝑖=1
𝑛 𝑥−𝑥𝑗
Donde 𝑙𝑖 (𝑥) = ∏𝑗=1 , 𝑖 = 1,2, … , 𝑛
𝑥𝑖 −𝑥𝑗
𝑗≠𝑖
𝑙𝑖 (𝑥𝑗 ) Son polinomios que tienen la propiedad de valer 1 cuando 𝑖 = 𝑗 y de tomar el valor cero cuando 𝑖 ≠ 𝑗.
Para probar que el polinomio pasa por los puntos se remplaza 𝑥 𝑝𝑜𝑟 𝑥𝑗 .
𝑛 𝑛
𝑃𝑛−1 (𝑥𝑗 ) = ∑ 𝑙𝑖 (𝑥𝑗 ) ∗ 𝑦𝑖 = ∑ 𝛿𝑖𝑗 ∗ 𝑦𝑖 = 𝑦𝑖
𝑖=1 𝑖=1
[(𝑥 − 𝑥1 ) ∗ (𝑥 − 𝑥2 ) … (𝑥 − 𝑥𝑛 )𝑓 𝑛 (ζ)]
𝑓(𝑥) − 𝑃𝑛−1 (𝑥) = , ζ ∈ [𝑥1 ; 𝑥𝑛 ]
𝑛!
Cuanto más lejos de x se toma el valor, mayor será su error.
Interpolación de Newton
Lagrange es conceptualmente simple pero no eficiente, ya que si se agregan o eliminan puntos debe
hacerse todo el procedimiento de nuevo. Newton ofrece una mejora computacional donde el polinomio
interpolante se escribe para n puntos:
22
Para cuatro puntos por ejemplo:
𝑃(𝑥1 ) = 𝑦1 = 𝑎0
𝑃(𝑥4 ) = 𝑦4 = 𝑎0 + (𝑥4 − 𝑥1 )𝑎1 + (𝑥4 − 𝑥1 )(𝑥4 − 𝑥2 )𝑎2 + (𝑥4 − 𝑥1 )(𝑥4 − 𝑥2 )(𝑥4 − 𝑥3 )𝑎4
.
∇𝑛−1 𝑦𝑛 − ∇𝑛−1 𝑦𝑛−1
∇𝑛 𝑦𝑛 =
𝑥𝑛 − 𝑥𝑛−1
𝑎0 = 𝑦1 , 𝑎1 = ∇𝑦2 , 𝑎3 = ∇2 𝑦3 , … , 𝑎𝑛 = ∇𝑛 𝑦𝑛
Generalmente es conveniente trabajarlo en tablas donde los datos ubicados en la diagonal son los
coeficientes del polinomio.
Si se copian en distintos orden los valores de la tabla cambian pero el polinomio será el mismo (el
polinomio de grado n-1 que interpola n puntos es único)
Interpolación de Neville
La interpolación polinómica involucra dos pasos: primero calcular los coeficientes y dos evaluar el
polinomio en el punto a interpolar. Si solo un punto debe ser interpolado, se utiliza Neville, que calcula
directamente el valor, en lugar de generar el polinomio.
Sea 𝑃𝑘 [𝑥𝑖 , 𝑥𝑖+1 , … , 𝑥𝑖+𝑘 ] el polinomio de grado k que pasa a través de k+1 puntos para un punto:
𝑃0 [𝑥𝑖 ] = 𝑦𝑖
23
[(𝑥 − 𝑥𝑖+𝑘 ) ∗ 𝑃𝑘−1 [𝑥𝑖 , 𝑥𝑖+1 , … , 𝑥𝑖+𝑘−1 ] − (𝑥𝑖 − 𝑥)𝑃𝑘−1 [𝑥𝑖 , 𝑥𝑖+1 , … , 𝑥𝑖+𝑘−1 , 𝑥𝑖+𝑘 ]
𝑃𝑘 [𝑥𝑖 , 𝑥𝑖+1 , … , 𝑥𝑖+𝑘 ] =
𝑥𝑖 − 𝑥𝑖+𝑘
La interpolación polinomial se debe realizar con la menor cantidad de puntos posibles: al utilizar muchos
puntos (6 o más), se obtienen polinomios de grado muy alto que presentan oscilaciones. Conocido como
el efecto Rungue.
Cuanto más cercano los datos mejor el resultado. Si un polinomio interpolante tiene bajo error y se
utilizaron más de seis puntos, debe revisarse con cuidado, ya que los puntos más alejados no aportan
demasiada información.
Hay que tener cuidado con la extrapolación polinómica. Si no puede evitarse, debe tenerse en cuenta lo
siguiente:
Graficar los datos y verificar que los valores extrapolados tengan sentido.
Realizar graficas de log(x) vs log (y), que suelen brindar curvas más suaves que x vs y.
Interpolación Segmentaria.
Si se divide el intervalo en un numero de subintervalos y se utiliza en cada uno de ellos una interpolación
polinómica de grado bajo, se pueden evitar oscilaciones polinomios de grado alto. Este tipo de
aproximación polinomiales a tramos se denominan splines, y los puntos que se utilizan se llaman nodos.
Si se requiere que la derivada de orden n fuera continua en todo punto, el spline sería un único polinomio.
Es imposible que dos polinomios distintos tengan el mismo valor en cada nodo y además sus derivadas
coincidan en los nodos.
Splines Lineales
Sea f una función real, definida y continua en [a; b]. Además, sea 𝑘 = {𝑥0 , 𝑥1 , … , 𝑥𝑚 }un subconjuto de [a;
b] con 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 = 𝑏, 𝑚 ≥ 2. El spline lineal 𝑆𝐿 que interpola a f en los puntos 𝑥𝑖 es:
𝑥𝑖 − 𝑥 (𝑥 − 𝑥𝑖−1 )
𝑆𝐿 (𝑥) = 𝑓(𝑥𝑖−1 ) + 𝑓(𝑥𝑖 ), 𝑥 ∈ [𝑥𝑖−1 , 𝑥𝑖 ]
𝑥𝑖 − 𝑥𝑖−1 𝑥𝑖 − 𝑥𝑖−1
𝑆𝐿 Interpola a f en los nodos por lo que 𝑆𝐿 (𝑥𝑖 ) = 𝑓(𝑥𝑖 ), y sobre cada intervalo [𝑥𝑖−1 , 𝑥𝑖 ], 𝑆𝐿 es una función
lineal. Por lo que 𝑆𝐿 es una función lineal a tramos sobre el intervalo [a;b].
24
Dada una serie de nodos 𝑘 = {𝑥0 , 𝑥1 , … , 𝑥𝑚 , ℎ𝑖 = 𝑥𝑖 − 𝑥𝑖−1 , 𝑦 ℎ = max{ℎ𝑖 }.
𝑖
𝑆𝐿 (𝑥) = 𝑆1 ∪ 𝑆2 ∪ …
Teorema
Sea 𝑓 ∈ 𝐶 2 [𝑎; 𝑏] 𝑦 𝑠𝑒𝑎 𝑆𝐿 el spline lineal que interpola a f en los nodos 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 = 𝑏
entonces la interpolación fue realizada con la siguiente cota de error:
1
‖𝑓 − 𝑆𝐿 ‖∞ ≤ ∗ ℎ2 ‖𝑓′′‖∞
8
Splines cúbicos
Sea ∈ 𝑓𝐶[𝑎; 𝑏] y 𝑘 = {𝑥0 , 𝑥1 , … , 𝑥𝑚 } un conjunto de m+1 nodo en [a;b], 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 = 𝑏,. Se
considera el conjunto S de todas las funciones S∈ 𝐶 2 [𝑎; 𝑏] tal que:
Cualquier elemento de S será un spline interpolador cubico. A diferencia de los splines lineales, existen
más de un spline S ∈ 𝐶 2 [𝑎; 𝑏] que satisface las condiciones antes mencionadas. Existen 4m coeficientes de
polinomio cúbicos (4 en cada subintervalo [𝑥𝑖−1 , 𝑥𝑖 ]), hay solo m+1 condiciones de interpolación y 3(m-1)
condiciones de continuidad.
Como S ∈ 𝐶 2 [𝑎; 𝑏], 𝑠, 𝑠′ 𝑦 𝑠′′son continuas en los nodos internos (𝑥1 , 𝑥2 , … , 𝑥𝑚−1 ) hay 4m-2 coeficientes
desconocidos. Dependiendo de la elección de las dos condiciones remanentes se pueden construir varios
tipos de splines cúbicos.
El spline natural cubico S2 es el elemento del conjunto S que satisface las condiciones terminales:
25
1 1
𝑓(𝑥𝑖−1 ) = ∗ 𝜎 ℎ2 + ℎ𝑖 𝛽𝑖 𝑓(𝑥𝑖 ) = ∗ 𝜎 ℎ2 + 𝛼𝑖 ℎ𝑖
6 𝑖−1 𝑖 6 𝑖 𝑖
Despejando 𝛼𝑖 𝑦 𝛽𝑖 , y aprovechando que 𝑠2′ (𝑥𝑖 +) = 𝑠′ 2 (𝑥𝑖 −) para i=1,2,….,m-1
𝑖 = 1,2, … , 𝑚-1
Teorema
Sea 𝑓 ∈ 𝐶 4 [𝑎; 𝑏] 𝑦 𝑠𝑒𝑎 𝑆 el spline cubico que interpola a f en los nodos 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 = 𝑏
entonces la interpolación fue realizada con la siguiente cota de error:
1
‖𝑓 − 𝑆‖∞ ≤ ∗ ℎ4 ‖𝑓 𝑖𝑣 ‖∞
24
Splines de Hermite
Se busca la suavidad y exactitud del spline interpolador respecto a la función a interpolar, exigiendo solo
que 𝑠 ∈ 𝐶 1 [𝑎; 𝑏]. Si y 𝑘 = {𝑥0 , 𝑥1 , … , 𝑥𝑚 } un conjunto de m+1 nodo en [a; b], 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 =
𝑏, 𝑚 ≥ 2, se define el spline cubico de Hermite como una función 𝑠 ∈ 𝐶 1 [𝑎; 𝑏] talque:
El spline cubico de Hermite tiene solo derivada continua en los nodos. No es un spline cubico en el sentido
de la definición original.
Los coeficientes del spline de Hermite se pueden calcular sin necesidad de resolver un sistema de
ecuaciones lineales diferencia de los splines cúbicos.
26
Teorema
Sea 𝑓 ∈ 𝐶 4 [𝑎; 𝑏] 𝑦 𝑠𝑒𝑎 𝑆 el spline de Hermite que interpola a f en los nodos 𝑎 = 𝑥𝑜 < 𝑥1 < ⋯ < 𝑥𝑚 = 𝑏
entonces la interpolación fue realizada con la siguiente cota de error:
1
‖𝑓 − 𝑆‖∞ ≤ ∗ ℎ4 ‖𝑓 𝑖𝑣 ‖∞
384
Hermite da una construcción casi perfecta de la función. De la derivada segunda no se da información. Es
un método local.
Ajustes de Datos
Para ajustar datos (generalmente obtenidos mediante mediciones) a una función que, teóricamente,
modela los datos, se puede utilizar mínimos cuadrados (clase general de métodos).
La función que modela los datos pertenece a una clase de funciones F. Se busca encontrar la mejor 𝑓 ∈ 𝐹
que ajuste los datos 𝑦 = 𝑓(𝑥).
Usualmente 𝐹 estará determinada por un número pequeño de parámetros, mucho más pequeño que la
cantidad de datos. Falta definir el concepto de mejor función y encontrar métodos para identificar los
parámetros que lo constituyen.
Mínimos cuadrados
Si 𝑥 = [𝑥0 , 𝑥1 , … 𝑥𝑛 ] , 𝑦 = [𝑦0 , 𝑦1 , … , 𝑦𝑛 ] son vectores de longitud n+1, el error o residuo de una función
dada con respecto a los datos es el vector: 𝑒 = 𝑓(𝑥) − 𝑦.
𝑒 = [𝑒0 , 𝑒1 , … , 𝑒𝑛 ]′ , 𝑒𝑖 = 𝑓(𝑥𝑖 ) − 𝑦𝑖
El método se llama método ordinario de los mínimos cuadrados, ya que se considera que no hay errores
en las mediciones de x.
𝐹 = {𝑓(𝑥) = 𝑎 ∗ 𝑥 + 𝑏, 𝑎, 𝑏 ∈ ℜ}
Debemos encontrar una f(x)=ax+b que minimice el error: ‖𝑒‖2 = ‖𝑓(𝑥) − 𝑦‖2 = √(∑𝑛𝑖=0(𝑓(𝑥𝑖 ) − 𝑦𝑖 )2 ).
27
2 𝑛 𝑛
𝑑 ‖𝑒‖2 𝑑
= ∑(𝑎𝑥𝑖 + 𝑏 − 𝑦𝑖 )2 = ∑ 2 ∗ (𝑎𝑥𝑖 + 𝑏 − 𝑦𝑖 )𝑥𝑖 = 0
𝑑𝑎 𝑑𝑎
𝑖=0 𝑖=0
2 𝑛 𝑛
𝑑 ‖𝑒‖2 𝑑
= ∑(𝑎𝑥𝑖 + 𝑏 − 𝑦𝑖 )2 = ∑ 2 ∗ (𝑎𝑥𝑖 + 𝑏 − 𝑦𝑖 )1 = 0
𝑑𝑏 𝑑𝑏
𝑖=0 𝑖=0
𝑛
𝑛 𝑛
𝑎 ∑ 𝑥𝑖 2 + 𝑏 ∑ 𝑥𝑖 = ∑ 𝑥𝑖 𝑦𝑖
𝑖=0 𝑖=0
𝑖=0
𝑛
𝑛 𝑛
𝑎 ∑ 𝑥𝑖 + 𝑏 ∑ 1=∑ 𝑦𝑖
𝑖=0 𝑖=0
𝑖=0
Nos queda un SEL. Si se observa que los datos pueden tomar una tendencia distinta a una recta, se puede
ajustar con una parábola.
Índice de Determinación
A medida que aumenta el grado del polinomio de ajuste, disminuye el error cuadrático del mismo. Dados
n+1 puntos, el mejor ajuste se logra con un polinomio de grado n, con ello el error cuadrático es 0. Pero a
medida es necesario un polinomio de grado reducido que ajuste una gran cantidad de datos.
2
El ‖𝑒‖2 no es un buen indicador de la aproximación polinómica ya que depende del número de puntos y
de las unidades de medida.
Se utiliza entonces el índice de determinación. Es una cantidad adimensional, va entre 0 y 1. Cuando más
cercano a uno mejore es el ajuste.
∑𝑛𝑖=0(𝑓(𝑥𝑖 ) − 𝑦̅)2
𝐼=
∑𝑛𝑖=0(𝑦𝑖 − 𝑦̅)2
Es sencillo ajustar cualquier conjunto de datos en dos dimensiones mediante una función arbitraria. Las
opciones más clásicas son:
Exponenciales: pueden dar errores altos debido a la linealizacion necesaria de la función de ajuste
Cuando las funciones de ajuste dependen de varios coeficientes el sistema resultante está mal
condicionado (el K(A) aumenta con el aumento del tamaño de la matriz). Cuando se linealita, se deben
28
ajustar los datos a la linealizacion, pero el error se calcula sobre la función de ajuste no sobre la
linealizacion.
Luego de obtener las constantes, se debe volver a transformar los valores hacia la función original.
Ejemplo
Sea 𝑦(𝑥) ∈ 𝐶[𝑎; 𝑏], y se busca un polinomio P(x) que ajuste lo mejor posible a dicha función. El error no
puede calcularse como lo veníamos haciendo por qué y no es una función discreta.
Se define el error como una sumatoria de términos infinitos, es decir, como una integral definida:
𝑏
𝐸 = ∫ [𝑃𝑛 (𝑥) − 𝑦(𝑥)]2 𝑑𝑥
2
𝑎
Por lo tanto hay que encontrar un P(x) tal que el 𝐸 2 sea mínimo.
𝐸2 = ∫ [∑ 𝑎𝑖 𝑥 𝑖 − 𝑦(𝑥)]2 𝑑𝑥
𝑎 𝑖=0
𝜕𝐸
Como 𝐸 = 𝐸(𝑎0 , 𝑎1 , … , 𝑎𝑛 ) para reducir E al mínimo, debe cumplirse: = 0.
𝜕𝑎𝑖
𝑏 𝑛
𝐸2 = ∫ ∑[𝑎𝑖 𝑥 𝑖 − 𝑦(𝑥)]2 𝑑𝑥
𝑎 𝑖=0
𝑏 𝑏
𝜕𝐸2
= ∫ 2 ∗ [𝑎𝑖 𝑥 ] 𝑥 𝑑𝑥 − 2 ∫ 𝑦(𝑥)𝑥 𝑗 𝑑𝑥 + 0 = 0
𝑖 𝑗
𝜕𝑎𝑗 𝑎 𝑎
Debe resolverse un sistema lineal de tamaño n+1. Las n+1 ecuaciones normales tienen solución única si
𝑓 ∈ [𝑎; 𝑏].
29
Si se quiere aumentar el grado del polinomio que ajuste a una función, solo se calculan 2 integrales más:
una para los coeficientes 𝑎𝑖 y otra para los términos libres. La complicación es la resolución del sistema
cada vez más grande ya que aumenta el número de condición.
Derivación Numérica
Si se desea calcular una de las derivadas de la función y=f(x) en el punto 𝑥 = 𝑥𝑘 , y se tiene un conjunto de
datos (𝑥𝑖 , 𝑦𝑖 ), 𝑖 = 1,2, … , 𝑛; solo se tiene acceso a un número finito de pares (x, y) con los que computar la
derivada.
Uno de los métodos más básico es aproximar la función a un polinomio y luego diferenciarlo. Otra opción
es utilizar una serie de Taylor alrededor del punto 𝑥𝑘 . La ventaja de este último método es que se puede
calcular una cota para el error.
La derivación numérica no es un proceso preciso existen errores de redondeo, truncamiento y los que son
inherentes a la interpolación.
La derivación por medio de diferencias finitas se basa en la expansión de las series de Taylor dadas como:
ℎ2 ′′ ℎ3 ℎ4
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ′′′ (𝑥) + 𝑓 ′′′′ (𝑥) + ⋯
2 3! 4!
ℎ2 ′′ ℎ3 ′′′ ℎ4 ′′′′
𝑓(𝑥 − ℎ) = 𝑓(𝑥) − ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) − 𝑓 (𝑥) + 𝑓 (𝑥) − ⋯
2 3! 4!
(2ℎ)2 ′′ (2ℎ)3 ′′′ (2ℎ)4 ′′′′
𝑓(𝑥 + 2ℎ) = 𝑓(𝑥) + 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥) + ⋯
2 3! 4!
(2ℎ)2 ′′ (2ℎ)3 ′′′ (2ℎ)4 ′′′′
𝑓(𝑥 − 2ℎ) = 𝑓(𝑥) − 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) − 𝑓 (𝑥) + 𝑓 (𝑥) − ⋯
2 3! 4!
También se pueden plantear las sumas y las restas de series.
ℎ4 ′′′′
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) = 2𝑓(𝑥) + ℎ2 𝑓 ′′ (𝑥) + 𝑓 (𝑥) + ⋯
12
ℎ3 ′′′
𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) = 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + ⋯
3
4ℎ4 ′′′′
𝑓(𝑥 + 2ℎ) + 𝑓(𝑥 − 2ℎ) = 2𝑓(𝑥) + 4ℎ2 𝑓 ′′ (𝑥) + 𝑓 (𝑥) + ⋯
3
8ℎ3 ′′′
𝑓(𝑥 + 2ℎ) − 𝑓(𝑥 − 2ℎ) = 4ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + ⋯
3
Las sumas contienen solos los términos pares (derivadas de orden par) y las restas solos los términos
impares (derivadas de orden impar).
Estas ecuaciones pueden ser vistas como ecuaciones simultáneas que deben ser resueltas para obtener las
derivadas de f(x). El número de ecuaciones y la cantidad de términos a utilizar depende del orden de la
derivada y del grado de precisión deseado.
30
Primera aproximación por diferencias centrales
Tomando la ecuación:
ℎ3 ′′′
𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ) = 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + ⋯
3
Se obtiene la solución para f’(x):
A partir de:
ℎ4 ′′′′
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) = 2𝑓(𝑥) + ℎ2 𝑓 ′′ (𝑥) + 𝑓 (𝑥) + ⋯
12
ℎ4 ′′′′
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) − 2𝑓(𝑥) = ℎ2 𝑓 ′′ (𝑥) + 𝑓 (𝑥) + ⋯
12
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) − 2𝑓(𝑥) ′′
ℎ4 ′′′′
=𝑓 (𝑥) + 𝑓 (𝑥) + ⋯
ℎ2 12ℎ2
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) − 2𝑓(𝑥) ′′ (𝑥) +
ℎ2 ′′′′
= 𝑓 𝑓 (𝑥) + ⋯
ℎ2 12
𝑓(𝑥 + ℎ) + 𝑓(𝑥 − ℎ) − 2𝑓(𝑥)
𝑓 ′′ (𝑥) = + O(ℎ2 )
ℎ2
Las aproximaciones por diferencias centrales de orden superior se obtienen de forma similar.
No siempre se pueden construir aproximaciones por diferencias centrales, por ejemplo cuando se
necesitan derivadas en los extremos.
Si solo se tiene información a derecha o a izquierda se utilizan las aproximaciones hacia adelante y hacia
atrás. Estas aproximaciones se obtienen con las ecuaciones planteadas al comienzo.
ℎ2 ′′ ℎ3 ℎ4
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ′′′ (𝑥) + 𝑓 ′′′′ (𝑥) + ⋯
2 3! 4!
ℎ2 ′′ ℎ3 ℎ4
𝑓(𝑥 + ℎ) − 𝑓(𝑥) = ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ′′′ (𝑥) + 𝑓 ′′′′ (𝑥) + ⋯
2 3! 4!
31
𝑓(𝑥 + ℎ) − 𝑓(𝑥) ℎ𝑓 ′ (𝑥) ℎ2 ′′ ℎ3 ′′′ ℎ4 ′′′′
= + 𝑓 (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥) + ⋯
ℎ ℎ 2ℎ 3! ℎ 4! ℎ
𝑓(𝑥 + ℎ) − 𝑓(𝑥) ℎ ℎ2 ℎ3
= 𝑓′(𝑥) + 𝑓 ′′ (𝑥) + 𝑓 ′′′ (𝑥) + 𝑓 ′′′′ (𝑥) + ⋯
ℎ 2 3! 4!
𝑓(𝑥 + ℎ) − 𝑓(𝑥)
𝑓 ′ (𝑥) = + 𝑂(ℎ)
ℎ
Primera aproximación por diferencias hacia adelante. Y de manera similar:
−𝑓(𝑥 − ℎ) + 𝑓(𝑥)
𝑓 ′ (𝑥) = + 𝑂(ℎ)
ℎ
Primera aproximación por diferencias hacia atrás.
Las aproximaciones anteriores, con error de orden h, no son muy utilizadas, ya que se prefieren
expresiones con error de orden ℎ2 . Para obtener fórmulas de aproximación no centrales de O (ℎ2 ) se
deben utilizar mayor cantidad de términos de la serie de Taylor.
ℎ2 ′′ ℎ3 ℎ4
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ′′′ (𝑥) + 𝑓 ′′′′ (𝑥) + ⋯
2 3! 4!
(2ℎ)2 ′′ (2ℎ)3 ′′′ (2ℎ)4 ′′′′
𝑓(𝑥 + 2ℎ) = 𝑓(𝑥) + 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥) + 𝑓 (𝑥) + ⋯
2 3! 4!
Se elimina 𝑓′′(𝑥) multiplicando la primera expresión por 4 y restándole la segunda expresión. Y
despejando 𝑓′(𝑥) queda:
Extrapolación de Richardson
1
𝑓 ′ (𝑥) = [𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ)] + 𝑎2 ℎ2 + 𝑎4 ℎ4 + ⋯
2ℎ
Donde 𝑎2 , 𝑎4 , … dependen de f y de x. cuando esta información está disponible en un proceso numérico→
se puede aplicar Extrapolacion de Richardson para conseguir mayor precisión: se utiliza para generar
resultados de gran precisión con el uso de fórmulas de bajo orden.
Se puede aplicar siempre que se conozca la expresión del error de una técnica de aproximación,
dependiente de un parámetro (generalmente la longitud de paso h).
32
1
∅(ℎ) = [𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ)]
2ℎ
Se ve que ∅(ℎ)es una aproximación a 𝑓 ′ (𝑥) con un error de orden ℎ2 . Se tratara de computar lim ∅(ℎ)
ℎ→0
para conseguir una aproximación de mejor calidad que una simple aproximación discreta. Pero no se
puede utilizar h=0 pero si se puede computar ∅(ℎ/2).
ℎ ℎ 4 ℎ 6
𝑓 ′ (𝑥) − 𝑎2 ( )2 − 𝑎4 ( ) − 𝑎6 ( ) … = ∅(ℎ/2)
2 2 2
El término dominante del error 𝑎2 se puede eliminar multiplicando la segunda ecuación por 4 y
restándosela a la primera ecuación:
ℎ 3 15
∅(ℎ) − 4∅ ( ) = −3𝑓 ′ (𝑥) − 𝑎4 ℎ4 − 𝑎6 ℎ6
2 4 16
Multiplicando por -1/3:
ℎ 1 ℎ 1 5
∅ ( ) + [∅ ( ) − ∅(ℎ)] = 𝑓 ′ (𝑥) + 𝑎4 ℎ4 + 𝑎6 ℎ6
2 3 2 4 16
1 ℎ
Al agregar [∅ ( ) − ∅(ℎ)] a ∅(ℎ/2) se mejora la aproximación inicial hasta un orden de 𝑂(ℎ4 ), y al ser h
3 2
pequeño, es una mejora de buena calidad. Este proceso teóricamente puede repetirse.
Este es el proceso de Extrapolación de Richardson, y puede ser aplicado en secuencia para ir eliminando
los términos de orden superior en la expresión del error.
En forma general, para un cálculo F que depende de un parámetro pequeño h, y se conoce la expresión
ℎ1
completa del error, se tiene que para ℎ1 > ℎ2 , 𝑐𝑜𝑛 ℎ2 = :
𝑘
𝛽 𝛾
∅(ℎ1 ) = 𝐹 − 𝑎α ℎ1𝛼 − 𝑎𝛽 ℎ1 − 𝑎𝛾 ℎ1 − ⋯
𝛽 𝛾
∅(ℎ2 ) = 𝐹 − 𝑎α ℎ2𝛼 − 𝑎𝛽 ℎ2 − 𝑎𝛾 ℎ2 − ⋯
ℎ1 ℎ1 𝛼 ℎ1 𝛽 ℎ1 𝛾
)
∅(ℎ2 = ∅ ( ) = 𝐹 − 𝑎α ( ) − 𝑎𝛽 ( ) − 𝑎𝛾 ( ) − ⋯
𝑘 𝑘 𝑘 𝑘
Multiplicando esta última expresión por 𝐾 𝛼 .
ℎ1 ℎ1 𝛽 ℎ1 𝛾
𝐾𝛼∅ ( ) = 𝐾𝛼 𝐹 − 𝑎α (ℎ1 )𝛼 − 𝑎𝛽 𝐾 𝛼 ( 𝛼
) − 𝑎𝛾 𝐾 ( ) − ⋯
𝑘 𝑘 𝑘
ℎ
Restando 𝐾 𝛼 ∅ ( 1 ) − ∅(ℎ1 ) queda:
𝑘
ℎ1 ℎ1 𝛽 ℎ1 𝛾
𝐾𝛼 ∅ ( ) − )
∅(ℎ1 = (𝐾 𝛼 − 1)𝐹 − 𝐶𝛽 ( ) − 𝐶𝛾 ( ) − ⋯
𝑘 𝑘 𝑘
33
ℎ 𝛽
La aproximación a F es de orden 𝑂 (( 1 ) ).
𝑘
Si bien el proceso puede repetirse indefinidamente, hay que tener cuidado con la aritmética, para evitar
que se agregue más error debido al redondeo.
Ambas son aproximaciones por centrales por lo que al eliminar el término de ℎ2 el siguiente término es
ℎ4.
2ℎ
Aplicamos la Extrapolación con 𝑘 = =2
ℎ
Las formulas vistas solo dependen de h. sin embargo, el error no solo involucre al truncamiento de la seria,
sino también el error generado por el redondeo aritmético.
𝑦1 = 𝑓(𝑥 + ℎ) + 𝑒1
𝑦0 = 𝑓(𝑥) + 𝑒0
34
𝑒1 − 𝑒0 |𝑘1 |
𝐷𝑓(ℎ, 𝑥) = 𝑓 ′ (𝑥) + + ℎ, 𝑘1 = 𝑓′′(𝑥)
ℎ 2
𝑒1 − 𝑒0 |𝑘1 | 2𝜀 |𝑘1 |
|𝐷𝑓(ℎ, 𝑥) − 𝑓 ′ (𝑥)| = | |+ ℎ≤ + ℎ
ℎ 2 ℎ 2
El primer término del error representa el error debido al redondeo aritmético y segundo término del error
representa el error por truncamiento de la serie de Taylor.
𝜕 2𝜀 |𝑘1 |
Para minimizar la expresión del error derivamos e igualamos a 0: ( + ℎ) = 0
𝜕ℎ ℎ 2
2𝜀 |𝑘1 |
− 2 =−
ℎ 2
1 |𝑘1 |
2 =
ℎ 4∗𝜀
4𝜀
ℎ2 =
|𝑘1 |
𝜀
ℎ𝑜𝑝𝑡𝑖𝑚𝑜 = 2 ∗ √
|𝑘1 |
A medida que h decrece, el error por redondeo se puede incrementar, mientras que el error por
truncamiento disminuye. Esto se llama dilema de la longitud de paso.
Las expresiones de h optimo son poco realistas, y solo de valor teórico que que f’(x), f’’(x), etc se
desconocen.
Integración Numérica
𝑏
Si se desea resolver: ∫𝑎 𝑓(𝑥)𝑑𝑥 si no se puede aplicar la regla de Barrow (𝐹(𝑏) − 𝐹(𝑎)) o es muy dificl
encontrar la primitiva 𝐹(𝑥), se utilizan los métodos numéricos que son:
1) Integrandos expansibles por series: se reemplaza f(x) por la expansión de una serie de funciones
elementales que se pueden integrar fácilmente. Este método tiene la ventaja de poder ajustar la
precisión del resultado final.
2) Métodos de paso finito: se subdivide el intervalo de integración [a; b] en n partes, y luego se
aproxima cada una de las partes en forma geométrica.
3) Métodos de Cuadratura: también llamados Gaussianos. Se reemplaza el integrando por
polinomios interpoladores y/o ortogonales. La principal complicación es determinar los nodos y
los valores de los pesos.
El integrando f(x) definido en [a;b] se puede representar por una serie fácilmente integrable. Las
expansiones más comunes son las series de Fourier y Taylor. El análisis de error es muy sencillo ya que la
componente principal es aportada por el error de truncamiento.
35
Métodos de Paso Finito
Dado un integrando f(x) definido en un [a; b], se puede dividir el intervalo en n segmentos distintos
(subintervalos) ∆𝑥𝑖 = 𝑥𝑖 − 𝑥𝑖−1 , 𝑝𝑎𝑟𝑎 𝑖 = 1,2, … , 𝑛. Los métodos siguientes se basan en la suma de
Riemann (teorema por el cual se define la integral definida). Su versión finita es:
𝑛 𝑛
𝑆 = 𝑆1 + 𝑆2 + ⋯ . +𝑆𝑛 = ∑ 𝑆𝑖 = ∫ 𝑓(𝜀𝑖 )∆𝑥𝑖
𝑖=1 𝑖=1
El intervalo [a; b] se divide en n subintervalos utilizando una partición de n+1 puntos. La base es
ℎ𝑖 = 𝑥𝑖 − 𝑥𝑖−1 , 𝑝𝑎𝑟𝑎 𝑖 = 1,2, … , 𝑛 Y la elección de la altura del rectángulo determinara cual forma de
aproximación de la integral definida se utilizara (se suman las áreas de los rectángulos).
𝑥𝑖 +𝑥𝑖−1
3) La altura del i-esimo rectángulo es 𝑓 ( ),𝑖 = 1,2, . . , 𝑛
2
𝑏 𝑏
𝑥𝑖 + 𝑥𝑖−1
∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝑓 ( ) ∗ (𝑥𝑖 − 𝑥𝑖−1 )𝑑𝑥 + 𝐸(ℎ)
𝑎 𝑎 2
El límite de las tres expresiones se da cuando max{𝑥𝑖 − 𝑥𝑖−1 } → 0 y es igual al valor de la integral resuelta
en forma analítica (valor exacto). Para acotar el error debe tenerse en cuenta que:
𝑑𝐹 𝑛 𝑑𝑓 𝑛−1
=
𝑑𝑥 𝑛 𝑑𝑥 𝑛−1
Esto es cierto porque F(x) es primitiva de f(x).
Como F y f son funciones continuas y derivables, se las expandir por series de Taylor.
ℎ2 ′′ ℎ3
𝐹(𝑥 + ℎ) = 𝐹(𝑥) + ℎ𝐹 ′ (𝑥) + 𝐹 (𝑥) + 𝐹 ′′′ (𝑥) + 𝑂(ℎ4 )
2 3!
ℎ2 ′′ ℎ3
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑓 ′′′ (𝑥) + 𝑂(ℎ4 )
2 3!
36
Lo que permite acotar el error de truncamiento. Partiendo de la aproximación central y analizando e área
del intervalo [𝑥𝑖 ; 𝑥𝑖+ℎ𝑖 ]:
𝐸(ℎ) = |𝐼𝐴 − 𝐼𝑅 |
𝑥𝑖 +ℎ𝑖
𝑥𝑖 + (𝑥𝑖 + ℎ𝑖 )
𝐸(ℎ) = |∫ 𝑓(𝑥)𝑑𝑥 − ℎ𝑖 𝑓 ( )|
𝑥𝑖 2
𝑥𝑖 + (𝑥𝑖 + ℎ𝑖 )
𝐸(ℎ) = |𝐹(𝑥𝑖 + ℎ𝑖 ) − 𝐹(𝑥𝑖 ) − ℎ𝑖 𝑓 ( )|
2
ℎ𝑖 )
𝐸(ℎ) = |𝐹(𝑥𝑖 + ℎ𝑖 ) − 𝐹(𝑥𝑖 ) − ℎ𝑖 𝑓 (𝑥𝑖 + )|
2
ℎ𝑖 )
Aplicando la relación entre F(x) y f(x) y expandiendo por serie de Taylor 𝐹(𝑥𝑖 + ℎ𝑖 ) 𝑦 𝑓 (𝑥𝑖 + ) y
2
operando llagamos a:
ℎ𝑖3 ′′
𝐸(ℎ𝑖 ) = 𝑓 (𝑥𝑖 ) + 𝑂(ℎ𝑖4 )
24
𝑂(ℎ𝑖4 ) Error cometido por truncamiento de la serie de Taylor.
Para acotar el error para todo el proceso de integración, se suma cada error:
𝑛
ℎ𝑖3 ′′ 𝑛ℎ𝑖3 ′′
𝐸(ℎ) = ∑ 𝑓 (𝑥𝑖 ) ≤ 𝑓 (𝜀𝑖 ), 𝜀𝑖 𝜖[𝑥𝑜 , 𝑥𝑛 ]
24 24
𝑖=1
𝑏 𝑛
(𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖−1 ))
∫ 𝑓(𝑥)𝑑𝑥 = ∑(𝑥𝑖 − 𝑥𝑖−1 ) ∗ + 𝐸(ℎ)
𝑎 2
𝑖=1
𝐸(ℎ) = |𝐼𝐴 − 𝐼𝑇 |
𝑥𝑖 +ℎ𝑖 (𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖 + ℎ𝑖 ))
𝐸(ℎ) = |∫ 𝑓(𝑥)𝑑𝑥 − ℎ𝑖 𝑓 |
𝑥𝑖 2
ℎ𝑖
𝐸(ℎ) = |𝐹(𝑥𝑖 + ℎ𝑖 ) − 𝐹(𝑥𝑖 ) − [𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖 + ℎ𝑖 )]|
2
Nuevamente aplicando la relación entre F(x) y f(x) y expandiendo por serie de Taylor 𝐹(𝑥𝑖 + ℎ𝑖 ) 𝑦 𝑓(𝑥𝑖 +
ℎ𝑖 ) y operando llagamos a:
ℎ𝑖3 ′′
𝐸(ℎ𝑖 ) = 𝑓 (𝑥𝑖 ) + 𝑂(ℎ𝑖4 )
12
Para acotar el error producido en el proceso de integración sumamos los errores de cada subintervalo:
37
𝑛
ℎ𝑖3 ′′ 𝑛ℎ𝑖3 ′′
𝐸(ℎ) = ∑ 𝑓 (𝑥𝑖 ) ≤ 𝑓 (𝜀𝑖 ), 𝜀𝑖 𝜖[𝑥𝑜 , 𝑥𝑛 ]
12 12
𝑖=1
El método de integración por trapecios tiene mayor error por truncamiento que el método de rectángulos
centrales.
Método de Simpson
El método de trapecios puede considerarse basado en interpolación, ya que por cada subintervalo se
construye un polinomio lineal. Para mejorar la precisión (reducir el error de truncamiento), se puede
considerar polinomios cuadráticos en vez de lineales.
Para integrar f(x) en [a; b] se debe generar una partición {𝑥0 , 𝑥1 , … , 𝑥𝑛 } con n par y cada tres puntos de la
partición construir un polinomio interpolador cuadrático.
𝑛/2
𝑏 𝑥2𝑖
∫ 𝑓(𝑥)𝑑𝑥 = ∑ ∫ 𝑃2𝑖−1 (𝑥)𝑑𝑥 + 𝐸(ℎ)
𝑎 𝑖=1 𝑥2𝑖−1
Para construir los 𝑃𝑖 (𝑥) de grado dos eficientemente, es importante centrarlos en el punto medio de cada
intervalo (simplifica el desarrollo).
𝑃𝑘 (𝑥) = 𝑎𝑘 (𝑥 − 𝑥𝑘 )2 + 𝑏𝑘 (𝑥 − 𝑥𝑘 ) + 𝑐𝑘
Polinomio cuadrático que pasa por los puntos (𝑥𝑘−1 , 𝑓(𝑥𝑘−1 ), (𝑥𝑘 , 𝑓(𝑥𝑘 )), (𝑥𝑘+1 , 𝑓(𝑥𝑘+1 )).
Para hallar 𝑎𝑘 , 𝑏𝑘 𝑦 𝑐𝑘 , se forma un sistema de ecuaciones lineales reemplazando los tres puntos a
interpolar.
38
𝑥𝑘+1 𝑥𝑘+1
∫ 𝑓(𝑥)𝑑𝑥 = ∫ (𝑎𝑘 (𝑥 − 𝑥𝑘 )2 + 𝑏𝑘 (𝑥 − 𝑥𝑘 ) + 𝑐𝑘 )𝑑𝑥 + 𝐸(ℎ)
𝑥𝑘−1 𝑥𝑘−1
ℎ
[𝑓(𝑥𝑘+1 + 4𝑓(𝑥𝑘 ) + 𝑓(𝑥𝑘−1 )] + 𝐸(ℎ)
3
Aplicando este resultado:
𝑏 𝑥2 𝑥4 𝑥𝑛
∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝑃1 (𝑥)𝑑𝑥 + ∫ 𝑃3 (𝑥)𝑑𝑥 + ⋯ + ∫ 𝑃𝑛 (𝑥)𝑑𝑥 + 𝐸(ℎ)
𝑎 𝑥0 𝑥2 𝑥𝑛−2
𝑛
2 −1 𝑛/2
ℎ
[𝑓(𝑥0 ) + 2 ∗ ∑ 𝑓(𝑥2𝑖 ) + 4 ∗ ∑ 𝑓(𝑥2𝑖−1 )𝑑𝑥 + 𝑓(𝑥𝑛 )] + 𝐸(ℎ)
3
𝑖=1 𝑖=1
Recordar: n debe ser par y mayor o igual a 2. Esta se conoce como la formula compuesta de Simpson.
Se puede demostrar que una cota para el error asociado a la aproximación por medio de la formula
compuesta de Simpson es:
𝑛ℎ4 𝑖𝑣
𝐸(ℎ) = 𝑓 (𝜀), 𝜀 ∈ [𝑥0 ; 𝑥𝑛 ]
180
Métodos de Cuadratura
Por lo tanto se utilizan las fórmulas de cuadratura, aproximando la integral mediante una suma ponderada
de valores de la función en los nodos (ciertos puntos que pertenecen al intervalo [a; b]), y utilizando
coeficientes llamados pesos.
Hay que elegir convenientemente los nodos (si es posible) y los pesos, para que el error de integración sea
lo más pequeño posible.
Formula de un punto
39
Se desea que la formula se exacta para polinomios de grado 0, es decir:
𝑏
∫ 𝑓(𝑥)𝑑𝑥 ≈ 𝑃0 𝑓(𝑥0 )
𝑎
𝑏 − 𝑎 = 𝑃0 + 𝑃1
𝑏 2 𝑎2
− = 𝑃0 𝑎 + 𝑃1 𝑏
2 2
𝑏−𝑎
Operando obtenemos que 𝑃0 = 𝑃1 =
2
Nota: la cota del error planteado por el teorema 26 para el método de Simpson (Newton Cotes de tres
puntos) difiere con el mostrado anteriormente. Esto es así por que el teorema 26 brinda una cota para el
40
error global, mientras que la expresión ya desarrollada muestra una cota para el error de la formula
compuesta de Simpson.
41