1

Métodos numéricos
Unidad 1
Cifras significativas: Las cifras significativas de un número son aquellas que pueden usarse de
manera confiable. Las cifras significativas pueden usarse como criterio para especificar qué tan
confiables son los resultados que nos da un método numérico.

Exactitud: Se refiere a qué tan cercano está el valor calculado o medido del valor verdadero.

Precisión: Se refiere a qué tan cercanos se encuentran, unos de otros. Diversos valores calculados
o medidos.

Inexactitud o Sesgo: Se define como una desviación sistemática del valor verdadero.

Imprecisión o incertidumbre: Es la magnitud de la dispersión de los valores calculados.

Definiciones de Error

Errores de truncamiento: Resultan del uso de aproximaciones como un procedimiento

matemático exacto.

Errores de redondeo: Se producen cuando se usan números que tienen un límite de cifras
significativas para representar números exactos.

Error verdadero: 𝐸𝑡 = 𝑣𝑎𝑙𝑜𝑟 𝑣𝑒𝑟𝑑𝑎𝑑𝑒𝑟𝑜 − 𝑣𝑎𝑙𝑜𝑟 𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑑𝑜

 2

(𝑒𝑟𝑟𝑜𝑟 𝑣𝑒𝑟𝑑𝑎𝑑𝑒𝑟𝑜)
Error relativo porcentual verdadero: 𝜀𝑡 = ∗ 100
𝑣𝑎𝑙𝑜𝑟 𝑣𝑒𝑟𝑑𝑎𝑑𝑒𝑟𝑜

Error relativo porcentual aproximado:

(𝑒𝑟𝑟𝑜𝑟 𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑑𝑜)
 𝜀𝑎 = ∗ 100
𝑣𝑎𝑙𝑜𝑟 𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑑𝑜
(𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖ó𝑛 𝑎𝑐𝑡𝑢𝑎𝑙−𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖ó𝑛 𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟)
 𝜀𝑎 = ∗ 100, para métodos
𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖𝑜𝑛 𝑎𝑐𝑡𝑢𝑎𝑙
iterativos.

En el caso de los métodos iterativos un criterio para detenerlo cuando:

|𝜀𝑎 | < 𝜀𝑠

*donde 𝜺𝒔 es una tolerancia porcentual prefijada.

Relación entre el número de cifras significativas (𝒏) y 𝜺𝒔 :

𝜀𝑠 = (0.5 ∗ 102−𝑛 )%

Binario

*Se sabe dejar el primer bit para el signo

Punto Flotante

Forma en la que representa la computadora las cantidades fraccionarias, está lo hace con una
parte fraccionaria llamada mantisa (𝒎) y una parte entera llamada exponente (𝒆).

𝑚 ∗ 𝑏𝑒

*Con b = base del sistema numérico

 3

*La mantisa es normalizada si tiene primero cero dígitos (esto permite guardar más cifras
significativas). Por Ej: 0,1234 ∗ 10−1 en vez de 0,0123 ∗ 100, usando base 10.

Esto nos permite que tanto fracciones como números muy grandes puedan ser representados. El
uso de representación de punto flotante introduce error de redondeo, ya que la mantisa conserva
un número finito de cifras significativas (que depende de la cantidad de bits que se disponga).

Ej: número positivo de punto flotante más pequeño posible con 7 bits de base 2 utilizando los primeros dos
bits para los signos de la mantisa y el exponente, los siguientes dos para la magnitud del exponente y los
últimos tres para la magnitud de la mantisa normalizada.

Aspectos importantes de los puntos flotantes:

 El rango de cantidades que se puede representar es limitado. Hay números más grandes o
más chicos de los que pueden ser representados con la cantidad de bits que tenemos a mano,
overflow y underflow.
 Existe un número finito de cantidades que puede representarse dentro de un rango. Cuando
se quiera guardar en un punto flotante un número que no puede ser representado
exactamente por el punto flotante en cuestión estos se guardaran como la cantidad más
cercana a esta que si pueda representarse, empleando redondeo o corte . Este error es
conocido como error de cuantificación.
*Corte: Para cualquier cantidad que esté dentro de un intervalo de longitud ∆𝑥, esta se
representara como la cantidad en el extremo inferior del intervalo.
*Redondeo: Cualquier cantidad en un intervalo ∆𝑥, se representara como la cantidad más
cercana permitida.

 El intervalo entre dos cantidades que se pueden representarse en un punto flotante (∆𝒙)
aumenta a medida que estas crecen en magnitud.

Error de cuantificación: Este error es directamente proporcional a ∆𝑥. Tenemos dos casos:

Cuando se emplea corte:

 4

|∆𝑥|
≤ ℰ
|𝑥|

Cuando se emplea redondeo:

|∆𝑥| ℰ
≤
|𝑥| 2

A “𝓔” se lo denomina é𝒑𝒔𝒊𝒍𝒐𝒏 𝒅𝒆 𝒍𝒂 𝒎á𝒒𝒖𝒊𝒏𝒂 y se calcula como: 𝓔 = 𝒃𝟏−𝒕 , donde b=base y t

es el número de cifras significativas de la mantisa.

Errores de truncamiento

Serie de Taylor: Nos permite predecir el valor de una función (continua y diferenciable) en un
punto en términos del valor de la función y sus derivadas en otro punto.

Con residuo:

*Con ℎ = 𝑥𝑖+1 − 𝑥𝑖 , que recibe el nombre de incremento.

*En el texto uso ℰ en vez de

La expansión de la serie de Taylor de n-ésimo orden es exacta para un polinomio del mismo
orden. Para otro tipo de funciones se necesitaría un número infinitos de términos para qué la
expansión sea exacta.

Inconvenientes de la serie de Taylor: No se conoce ℰ con exactitud, solo que se encuentra entre
𝑥𝑖 y 𝑥𝑖+1 (por teorema de valor medio) y se requiere de la n-ésima primera derivada de 𝑓(𝑥). Por
lo tanto lo único que sabemos es que 𝑹𝒏 = 𝑶(𝒉𝒏+𝟏 )

*La “O” indica que una proporcionalidad directa con lo que hay en su argumento.

Error numérico total: Es la suma de los errores de truncamiento y de redondeo. La relación que
hay entre estos es que a medida que disminuye uno aumenta el otro, a estos los determina el
tamaño del paso (ℎ).
 5

*Notar que los ejes no son lineales sino logarítmicos.

Unidad 2

Métodos cerrados
Son métodos que para encontrar una raíz se valen del cambio de signo que tiene una función en
la vecindad de esta. Para estos tipos de métodos se necesitan de dos valores iniciales qué deben
“encerrar” la raíz. Siempre generan aproximaciones cada vez más cercanas a la raíz, son métodos
convergentes.

Criterio: Si 𝑓(𝑥) es una función real y continúa en el intervalo delimitado por 𝑥𝑙 y 𝑥𝑢 , entonces sí:

𝑓(𝑥𝑙 )𝑓(𝑥𝑢 ) < 0

Entonces existe al menos una raíz real entre 𝑥𝑙 y 𝑥𝑢 .

Criterio de paro: Se utiliza como criterio de paro para detener los cálculos, cuando 𝜀𝑎 < 𝜀𝑠 ; donde
𝜀𝑠 es un valor prefijado. Siempre el 𝜺𝒂 es mayor a 𝜺𝒕 (error relativo porcentual verdadero). Por lo
tanto lo que estamos asegurando con este criterio es 𝜀𝑡 < 𝜀𝑠 que es lo que realmente buscamos.

Estimaciones de errores: tenemos como error relativo porcentual aproximado:

(𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖ó𝑛 𝑎𝑐𝑡𝑢𝑎𝑙 − 𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖ó𝑛 𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟) 𝑥𝑟𝑛𝑢𝑒𝑣𝑜 − 𝑥𝑟𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟

𝜀𝑎 = | ∗ 100| = | | ∗ 100
𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑐𝑖𝑜𝑛 𝑎𝑐𝑡𝑢𝑎𝑙 𝑥𝑟𝑛𝑢𝑒𝑣𝑜
 6

Método de bisección

Es un método de búsqueda incremental, es decir el método repite el proceso hasta obtener una
mejor aproximación. La posición de la raíz (𝑥𝑟 ) se determina situándola en el punto medio del
subintervalo, dentro del cual ocurre un cambio de signo. Los pasos a seguir son:

Estimaciones de errores: para este método podemos el error relativo porcentual aproximado nos
queda:

𝑥𝑟𝑛𝑢𝑒𝑣𝑜 − 𝑥𝑟𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 𝑥𝑢 − 𝑥𝑙
𝜀𝑎 = | 𝑛𝑢𝑒𝑣𝑜 | ∗ 100 = | | ∗ 100
𝑥𝑟 𝑥𝑢 + 𝑥𝑙

Numero de iteraciones: La cantidad de iteraciones se puede calcular a priori, debido a que cada
iteración reduce el error asociado a la iteración cero (antes de ejecutar el método), ∆𝑥 0 , a la
mitad. La ecuación que relaciona el error y el número de iteraciones, n, es:

∆𝒙𝟎 ∆𝒙𝟎
𝑬𝒏𝒂 = 𝟐𝒏
, 𝑑𝑒 𝑑𝑜𝑛𝑑𝑒 𝑠𝑒 𝑑𝑒𝑠𝑝𝑒𝑗𝑎 𝒏 = 𝒍𝒐𝒈𝟐 ( 𝑬𝒏 )
𝒂

Pseudocodigo:
 7

Método de la falsa posición

A este método se lo conoce también como “regula falsi” o como “método de interpolación lineal”.
Consiste en unir 𝑓(𝑥𝑙 ) y 𝑓(𝑥𝑢 ) con una recta

Por triángulos semejantes:

𝑓(𝑥𝑙 ) 𝑓(𝑥𝑢 )
=
𝑥𝑟 − 𝑥𝑙 𝑥𝑢 − 𝑥𝑟

Despejando:

𝑓(𝑥𝑢 )(𝑥𝑙 − 𝑥𝑢 )
𝑥𝑟 = 𝑥𝑢 −
𝑓(𝑥𝑙 ) − 𝑓(𝑥𝑢 )

Luego para decidir en que subintervalo se

encuentra la raíz usamos es criterio de los
métodos cerrados.

*Nota: no se puede calcular el 𝜺𝒂 en la iteración cero

Pseudocodigo:
 8

Desventajas del método: La desventaja que presenta Comparación de los errores

este método es que tiende a mostrar problemas de relativos de los métodos
convergencia a cuando unos de los límites del intervalo busqueda de raices
permanece fijo o “estancado” durante muchas
iteraciones, a esto se conoce como unilateralidad.

Método falsa posición modificada

Resuelve el problema de unilateralidad del método de

falsa posición. Detecta cuando se produce un
estancamiento, toma como estancamiento cuando
alguno de los límites del intervalo permanece fijo por
más de dos iteraciones. Luego si esto ocurre el valor de
la función en este punto de divide a la mitad.

Pseudocodigo método falsa posición modificada:

 9

Para el LibreOfficeCalc: Supongamos que el método se ha estancado en el límite inferior del

intervalo, 𝑥𝑙 , el cual se encuentra cargado en la columna A y la función valuada en este punto,
𝑓(𝑥𝑙 ), cargada en la columna B. Aplicamos el método de falsa posición modificada a partir de la
segunda iteración. (Análogo para 𝑥𝑢 y 𝑓(𝑥𝑢 ))

Iteración 𝑥𝑙 𝑓(𝑥𝑙 ) Para poder generar B3 y los siguientes a él, deberemos emplear la
0 A1 B1 función “SI (condición; caso verdadero; caso falso)” de la siguiente
1 A2 B2 manera:
2 A3 B3

B3 = SI (A1 =A3; SI (B1=B2; B2/2; B2); 𝒇(A3))

Donde:

 A1=A2=A3, asegura el estancamiento durante dos iteraciones previas.

 En el caso verdadero: B1=B2, prueba si hace dos iteraciones que está
estancado con el mismo valor de 𝑓(𝑥𝑙 ). Si esto se cumple se divide el valor de
𝑓(𝑥𝑙 ) a la mitad, en caso contrario se mantiene el valor de 𝑓(𝑥𝑙 ).
 En caso falso: se calcula B3 a partir de 𝑓(𝑥).
 10

Métodos abiertos
Son métodos que se basan en fórmulas que requieren de un valor inicial o dos que no
necesariamente encierran una raíz. En algunos casos divergen pero cuando convergen suelen
hacerlo mucho más rápido que los métodos cerrados.

Iteración simple de punto fijo:

Es un método que solo necesita de un valor inicial y se basa de despejar la variable 𝒙 de

𝒇(𝒙) = 𝟎, para obtener:

𝑥 = 𝑔(𝑥)

Este método nos permite obtener una nueva aproximación a la raíz"𝑥", a partir de una
aproximación anterior:

𝒙𝒊+𝟏 = 𝒈(𝒙𝒊 )

𝒙𝒊+𝟏 −𝒙𝒊
Error relativo porcentual aproximado: 𝜺𝒂 = | 𝒙𝒊+𝟏
| ∗ 𝟏𝟎𝟎

Convergencia: El error relativo porcentual disminuye por un factor de 0,5 – 0,6 por cada iteración,
esto se lo conoce como convergencia lineal. Un método grafico para observar la convergencia o
divergencia consiste en graficar las curvas 𝑥 y 𝑔(𝑥). Los valores de 𝑥 correspondientes a las
intersecciones de estas curvas son las raíces de 𝑓(𝑥) = 0. La convergencia y divergencia pueden
ser:

𝑎) Convergencia
monótona

𝑏) Convergencia
oscilatoria

𝑐) Divergencia monótona

𝑑) Divergencia oscilatoria
 11

*Notar que la convergencia ocurre cuando |𝒈´(𝒙)| < 𝟏 nos garantiza de que 𝜀𝑎 disminuya con
cada iteración. (Se puede llegar a esto restando, 𝑥𝑖+1 = 𝑔(𝑥𝑖+1 ) a 𝑥𝑟 = 𝑔(𝑥𝑟 ) y usando el TVM)

Pseudocodigo para el método de la falsa posición:

Método Newton-Raphson
Utiliza un valor inicial 𝑥𝑖 que utiliza para trazar una recta tangente desde el puntos [𝒙𝒊 , 𝒇(𝒙𝒊 )]
sobre la curva. Donde corta la recta tangente con el eje de abscisas (x), representa una
aproximación mejorada de la raíz.

Usando la serie de Taylor podemos

calcular la función como:

𝑓(𝑥) = 𝑓(𝑥𝑖 ) + 𝑓´(𝑥𝑖 )(𝑥 − 𝑥𝑖 ) +

𝑓´´(ℰ)
(𝑥 − 𝑥𝑖 )2
2!

Truncando la expresión exacta se obtiene

una aproximación (una recta):

𝑓(𝑥) ≅ 𝑓(𝑥𝑖 ) + 𝑓´(𝑥𝑖 )(𝑥 − 𝑥𝑖 )

Luego la raíz de la recta puede pensarse

como una aproximación a la raíz de la función 𝑓(𝑥), denotando como 𝑥𝑖+1 a la raíz de la
recta y despejando:
𝒇(𝒙 )
𝑓(𝑥𝑖+1 ) = 𝑓(𝑥𝑖 ) + 𝑓´(𝑥𝑖 )(𝑥𝑖+1 − 𝑥𝑖 ) = 0 → 𝒙𝒊+𝟏 = 𝒙𝒊 − 𝒇´(𝒙𝒊 )
𝒊

Errores: Para calcular el error verdadero del método restamos la serie de Taylor de la función
valuada en su raíz 𝑓(𝑥𝑟 ), con su la aproximación (recta), 𝑓(𝑥𝑖+1 ). Queda:
 12

𝑓´´(ℰ)
0 = 𝑓´(𝑥𝑖 )(𝑥𝑟 − 𝑥𝑖+1 ) + (𝑥𝑟 − 𝑥𝑖 )2
2!

Tomando 𝐸𝑡,𝑖+1 = 𝑥𝑟 − 𝑥𝑖+1 , 𝐸𝑡,𝑖 = 𝑥𝑟 − 𝑥𝑖 y despejando:

𝒇´´(𝓔)
𝑬𝒕,𝒊+𝟏 = (𝑬 )𝟐
𝟐! 𝒇´(𝒙𝒊 ) 𝒕,𝒊

De esto se observa que el error es proporcional al cuadrado del error anterior, este tipo de error

Convergencia: el método presenta convergencia cuadrática, ya su error de cada iteración es

proporcional al cuadrado del error de la anterior. Por el teorema de Newton Raphson se
demuestra que se asegura la convergencia si |𝑔(𝑥𝑖 )| < 1 , ∀𝑥𝑖 𝜖(𝑥0 , 𝑥𝑟 ). Donde 𝑥𝑖+1 = 𝑔(𝑥𝑖 ).

Desventajas del método: si la convergencia no se asegura el método tiende a divergir cuando la

raíz se encuentra cerca de un punto de inflexión, tiende a oscilar alrededor de máximos o mínimos
locales lo cual puede persistir o no y no tolera cuando 𝑓´(𝑥𝑖 ) = 0.

Pseudocodigo del método Newton Raphson:

Método de la secante

Surge de evitar el cálculo de la derivada del método

Newton-Raphson y usar una aproximación

𝑓(𝑥𝑖−1 ) − 𝑓(𝑥𝑖 )
𝑓´(𝑥𝑖 ) ≅
𝑥𝑖−1 − 𝑥𝑖

Reemplazando en la aproximación de la raíz del

método Newton- Raphson queda:

𝒇(𝒙𝒊 )(𝒙𝒊−𝟏 − 𝒙𝒊 )
𝒙𝒊+𝟏 = 𝒙𝒊 −
𝒇(𝒙𝒊−𝟏 ) − 𝒇(𝒙𝒊 )
 13

Este método requiere de dos valores iniciales.

Diferencias con el método de la falsa posición: La diferencia que hay entre estos dos métodos es
que el método de la secante reemplaza las aproximaciones, 𝑥𝑖 𝑦 𝑥𝑖−1 , en una secuencia estricta.
Cuando el método de la falsa posición emplea el criterio de los métodos cerrados.

Error y convergencia: Se puede demostrar que 𝑬𝒕,𝒊+𝟏 = 𝒄𝒕𝒆 ∗ 𝑬𝟏,𝟔𝟏𝟖

𝒕,𝒊 , se dice entonces que tiene
convergencia superlineal.

Pseudocodigo del método de la secante:

Método de la secante modificado

Utiliza un solo valor inicial para aproximar la derivada haciendo uso de un pequeño fraccionario
𝛿 en la variable independiente.

𝑓(𝑥𝑖 + 𝛿𝑥𝑖 ) − 𝑓(𝑥𝑖 ) 𝒇(𝒙𝒊 )(𝜹𝒙𝒊 )

𝑓´(𝑥𝑖 ) ≅ → 𝒙𝒊+𝟏 = 𝒙𝒊 −
𝛿𝑥𝑖 𝒇(𝒙𝒊 + 𝜹𝒙𝒊 ) − 𝒇(𝒙𝒊 )

Raíces múltiples

Para el método de la secante y Newton-Raphson modificados en el caso de raíces múltiples, se

define una nueva función 𝑢(𝑥) que tiene raíces en las mismas posiciones que la función original:

𝑓(𝑥)
𝑢(𝑥) =
𝑓´(𝑥)

El uso 𝑢(𝑥) viene a resolver en el caso de raíces múltiples, donde 𝑓´(𝑥) tiende a aproximarse a
cero cerca de la raíz. Esto que genera una división por cero en los métodos de la secante y
𝑓(𝑥)
Newton-Raphson. Cómo 𝑓(𝑥) tiende a cero más rápido que 𝑓´(𝑥) el cociente 𝑓´(𝑥)
se vuelve cero
antes que se provoque la división por cero.

Ambos métodos tienen convergencia cuadrática.

 14

Método Newton-Raphson modificado para raíces múltiples:

𝑓(𝑥) 𝒇(𝒙𝒊 )𝒇´(𝒙𝒊 )

𝑥𝑖+1 = 𝑥𝑖 − → 𝒙𝒊+𝟏 = 𝒙𝒊 −
𝑓´(𝑥) [𝒇´(𝒙𝒊 )]𝟐 − 𝒇(𝒙𝒊 )𝒇´´(𝒙𝒊 )

Método de la secante modificado para raíces múltiples:

𝑓(𝑥𝑖 )(𝑥𝑖−1 − 𝑥𝑖 ) 𝒖(𝒙𝒊 )(𝒙𝒊−𝟏 − 𝒙𝒊 )

𝑥𝑖+1 = 𝑥𝑖 − → 𝒙𝒊+𝟏 = 𝒙𝒊 −
𝑓(𝑥𝑖−1 ) − 𝑓(𝑥𝑖 ) 𝒖(𝒙𝒊−𝟏 ) − 𝒖(𝒙𝒊 )

Pseudocodigo método Newton-Raphson modificado:

Raíces de Polinomios
Deflación polinomial

Es el proceso de eliminar una raíz ya encontrada de un polinomio lo que nos permite buscar las
demás raíces de este. Una forma de lograr esto es dividir nuestro polinomio de n-ésimo grado
entre un factor monomial 𝑥 − 𝑡, donde t es la raíz encontrada. Como solo se conoce la raíz de
forma aproximada, la deflación polinomial es sensible al error de redondeo.

Método de Müller

Construye una parábola haciendo uso de tres puntos, el método consiste en obtener los
coeficientes 𝑎, 𝑏 𝑦 𝑐 de la parábola. La parábola debe de pasar por los puntos
𝑓(𝑥0 ), 𝑓(𝑥1 ) 𝑦 𝑓(𝑥2 ), entonces:
 15

𝑓(𝑥0 ) = 𝑎(𝑥0 − 𝑥2 )2 + 𝑏(𝑥0 − 𝑥2 ) + 𝑐

𝑓(𝑥1 ) = 𝑎(𝑥1 − 𝑥2 )2 + 𝑏(𝑥1 − 𝑥2 ) + 𝑐

𝑓(𝑥2 ) = 𝑎(𝑥2 − 𝑥2 )2 + 𝑏(𝑥2 − 𝑥2 ) + 𝑐

Despejando:

𝛿 −𝛿
𝑎 = 𝑥1 −𝑥0 𝑏 = 𝛿0 − 𝑎(𝑥0 − 𝑥2 ) 𝑐 = 𝑓(𝑥2 )
1 0

𝑓(𝑥0 )−𝑓(𝑥2 ) 𝑓(𝑥1 )−𝑓(𝑥2 )

Donde: 𝛿0 = ; 𝛿1 = ;
𝑥0 −𝑥2 𝑥1 −𝑥2

Para hallar la raíz se utilizara una forma alternativa a la fórmula cuadrática (Bascara):

De las dos raíces se elige la raíz que coincida con el signo de 𝒃. Esta elección proporciona el
denominador más grande y por lo tanto una raíz estimada más cercana a 𝑥2 , lo cual disminuye el
error.

Una vez calculado 𝑥3 para la siguiente iteración se hace lo siguiente:

 Si se buscan raíces reales se descarta el valor más alejado de la raíz calculada 𝑥3 .

 Si se buscan raíces reales y complejas 𝑥1 , 𝑥2 𝑦 𝑥3 toman el lugar de 𝑥0 , 𝑥1 𝑦 𝑥2 .

𝒙𝟑 −𝒙𝟐
Error relativo porcentual aproximado: 𝜺𝒂 = | 𝒙𝟑
|∗ 𝟏𝟎𝟎

Pseudocodigo método Müller

 16

Unidad 3: Sistemas de ecuaciones lineales

Métodos de eliminación
Eliminación de Gauss: Es una técnica que implica combinación de ecuaciones para eliminar
incógnitas. Se analizaran sistemas de ecuaciones del tipo:

Solución para pequeños sistemas de ecuaciones (𝒏 ≤ 𝟑): Las formas de resolver estos sistemas
de ecuaciones sin el uso de una computadora son las siguientes:

 Método gráfico: para 𝑛 = 2 (dos incógnitas), se grafican unas coordenadas cartesianas

donde un eje corresponda a 𝑥1 y el otro a 𝑥2 . Se despeja en ambas ecuaciones 𝑥2
obteniéndose así dos rectas, en la intersección de estas los valores de 𝑥1 𝑦 𝑥2 .
 Regla de Cramer: Establece que cada incógnita de un S.E.L puede expresar como una
fracción de dos determinantes.

Donde D es el determinante de la matriz de coeficientes y la matriz del numerador se obtiene

de reemplazar en la matriz de coeficientes la columna correspondiente a la incógnita que se
quiera expresar por la matriz de términos independientes (que es una columna).

 Eliminación de incógnitas: consiste en multiplicar ecuaciones por constantes, de manera que

se eliminen incógnitas cuando se combinen las ecuaciones y así poder despejar la incógnita
restante. Esta se luego se reemplaza en cualquier otra de las ecuaciones con el fin de
despejar el resto de las incógnitas. (Intro. a la matemática)

Eliminación de Gauss simple

Es una técnica de eliminación de incógnitas para sistemas de ecuaciones simultáneas lineales. El

procedimiento consiste de dos pasos: la eliminación hacia adelante y la sustitución hacia atrás
 17

Eliminación hacia adelante: Se multiplican las ecuaciones por escalares y se combinan hasta
obtener una matriz triangular superior o escalón reducida.

Por lo que para cada incógnita, 𝑥𝑖 , se la elimina desde la i-ésima primera hasta la n-ésima
ecuación; con 𝑖 = 1,2, … , 𝑛. Para lograr esto se resta la i-ésima ecuación multiplicada por 𝑎𝑗𝑖 /𝑎𝑖𝑖 a
la j-ésima ecuación, para 𝑗 = 𝑖 + 1, 𝑖 + 2, … , 𝑛. Obteniendo como resultado:

La i-ésima ecuación recibe el nombre de ecuación pivote y 𝑎𝑖𝑖 coeficiente o elemento pivote.

Sustitución hacia atrás: Se empiezan a despejar las incógnitas a partir de la n-ésima ecuación,
sustituyéndolas hacia atrás para despejar las demás. El proceso que se repite para evaluar las
incógnitas restantes es representado mediante la fórmula.

(𝑖−1) (𝑖−1)
𝑏𝑖 − ∑𝑛𝑗=𝑖+1 𝑎𝑖𝑗 𝑥𝑗
𝑥𝑖 = (𝑖−1)
; 𝑝𝑎𝑟𝑎 𝑖 = 𝑛 − 1, 𝑛 − 2, … , 1
𝑎𝑖𝑖

Donde el superíndice (i-1) indica la cantidad de modificaciones que se le han hecho al coeficiente
o término independiente.

Pseudocodigo de Gauss simple:

Eliminación hacia adelante: Donde “A” es la matriz de coeficientes,” b” es la matriz de términos

independiente, “x” un vector de tamaño n que inicia en cero y “a” es la matriz aumentada [A|b]

Versión 1:

Versión 2
 18

Eliminación hacia atrás

Versión 1

Versión 2

𝟐𝒏𝟑 𝒏𝟐 𝟐𝒏𝟑
Costo computacional: 𝐥𝐢𝐦𝒏→∞ [ + 𝑶(𝒏𝟐 ) + + 𝑶(𝒏)] = + 𝑶(𝒏𝟐 ) ; donde
𝟑 𝟐 𝟑
2𝑛3 𝑛2
3
+ 𝑂(𝑛2 ) corresponde al costo de la eliminación hacia adelante y 2
+ 𝑂(𝑛) a la sustitución
hacia atrás. Esto se calcula contando la cantidad de operaciones que se efectúan para llevar a cabo
la técnica.

*Nota: La cantidad de FLOP (operaciones de punto flotante) se utilizan para determinar el costo
computacional de la técnica. Estas comprenden la suma multiplicación, división, suma y resta.

Problemas de la eliminación de Gauss simple

1) Errores de redondeo: El método de eliminación de Gauss es sensible al error de redondeo,

ya que cada resultado depende de los anteriores. Una regla generalizada para considerar
importante el error de redondeo es cuando el sistema es de 100 o más ecuaciones.
2) División por cero durante la eliminación hacia adelante: esto ocurre cuando el elemento
pivote es cero.
3) Sistemas mal condicionados y sistemas singulares: Los sistemas mal condicionados son
aquellos en donde pequeños cambios en uno o más coeficientes generan grandes cambios en
la solución. Se caracterizan por tener un determinante cercano a cero. Los sistemas
singulares son aquellos en los cuales el determinante es cero por lo que la única solución es la
trivial.
 19

*Nota: El determinante es igual al producto de los elementos de la diagonal principal de la matriz

´´ ´´´ 𝑛−1
triangular superior. 𝐷 = 𝑎11 𝑎22 𝑎33 … 𝑎𝑛𝑛

Modificaciones a la eliminación de Gauss simple (pivoteo y escalamiento)

1) El uso de más cifras significativas reduce el error de redondeo.

2) Para solucionar el problema de la división por cero es usar el pivoteo parcial. Este consiste
en que antes de normalizar cada reglón se determine el coeficiente más grande disponible
en la columna debajo del elemento pivote y se intercambien reglones tal que el coeficiente
más grande sea el elemento pivote. El pivoteo parcial también reduce el error de redondeo,
ya que en la eliminación hacia adelante se evita tener una división con un elemento pivote
igual – cercano a cero.
3) El uso del escalamiento permite la estandarización del determinante (𝐷 = 1), así como
también disminuye el error de redondeo.

Pseudocodigo eliminación de Gauss con pivoteo y escalamiento: Donde “A” es la matriz de

coeficientes,” b” es la matriz de términos independiente, “x” y “s” son vectores de tamaño n que
se inician en cero y “a” es la matriz aumentada [A|b]
 20

Costo computacional de la eliminación de Gauss con pivoteo y escalamiento: El costo original

2𝑛3 𝑛2 3
era: + 𝑂(𝑛2 ) + + 𝑂(𝑛), con el pivoteo y el escalamiento se agregan 𝑛(𝑛 − 1)
3 2 2
𝑛(𝑛+1)
comparaciones y − 1 divsiones. Luego para grandes sistemas tenemos
2
2𝑛3 𝑛2 3 𝑛(𝑛+1) 2𝑛3
lim𝑛→∞ [ + 𝑂(𝑛2 ) + + 𝑂(𝑛) + 2 𝑛(𝑛 − 1) + − 1] = + 𝑂(𝑛2 ), se puede
3 2 2 3
concluir que para grandes sistemas de ecuaciones el pivoteo y el escalamiento no
aumentan significativamente el costo computacional.

Descomposición LU

Se usa para resolver sistemas de ecuaciones lineales, es una alternativa más eficiente que la
eliminación de Gauss cuando deben resolverse ecuaciones con los mismos coeficientes [A] pero
diferentes términos independientes, {B}. Separa el tiempo utilizado para en las eliminaciones de la
matriz [A] de las manipulaciones del lado derecho {B}, esto permite evaluar los múltiples vectores
{B} de manera eficiente una vez que se ha descompuesto [A].

En la descomposición LU se parte de:

[𝐴]{𝑋} = {𝐵}

Y se busca descomponer [A] y {B} en:

[𝐴] = [𝐿][𝑈] 𝑦 {𝐵} = [𝐿]{𝐷}

𝑐𝑜𝑛 [𝑈]{𝑋} = {𝐷}

Donde [𝑈] es una matriz triangular superior y [L] es una matriz triangular inferior de la forma:

𝟏 𝟎 𝟎
[𝑳] = 𝒍𝟐𝟏 𝟏 𝟎 ; Ej: para un sistema de 3 ecuaciones
𝒍𝟑𝟏 𝒍𝟑𝟐 𝟏
 21

Entonces la descomposición LU se puede reducir a dos pasos:

 Paso de la descomposición: se descompone

[A] en las matrices triangulares superior e
inferior [L] y [U]
 Paso de la sustitución: Una vez obtenidas las
matrices [L] y [U] se utilizan para determinar
una solución {X} para un {B} dado. Esto se
puede descomponer en dos pasos:
 De la relación [𝐿]{𝐷} = {𝐵} se debe
obtener {𝐷} mediante la sustitución
hacia adelante (por ser [𝐿] una matriz
triangular inferior).
 Se reemplaza {𝐷} en [𝑈]{𝑋} = {𝐷} y
se obtiene {𝑋} mediante sustitución
hacia atrás (por ser [𝑈] una matriz
triangular superior).

Eliminación de Gauss utilizando descomposición LU

El resultado de la eliminación hacia delante de la eliminación de Gauss es la matriz

[U]. Usando como ejemplo sistema de tres ecuaciones la matriz obtenida será de la
forma:
𝑎11 𝑎12 𝑎13
, ,
𝑈= 0 𝑎22 𝑎23
,,
0 0 𝑎33

Donde la matriz [L] se obtiene en el proceso de eliminación y se puede demostrar que

es de la forma:

1 0 0
𝑓
𝐿 = 21 1 0
𝑓31 𝑓32 1
𝑎 𝑎 𝑎,
Donde 𝑓21 = 𝑎21 , 𝑓31 = 𝑎31 𝑦 𝑓32 = 𝑎32
, . De manera que [𝐿][𝑈] = [𝐴], esta matriz es
11 11 22
ligeramente diferente a la original debido a los errores de redondeo.

Una vez terminada la descomposición para generar una solución para un vector {B}
dado, primero debemos obtener la matriz {D} de [𝐿]{𝐷} = {𝐵} mediante la sustitución
hacia adelante usado:
𝑖−1

𝑑𝑖 = 𝑑𝑖 − ∑ 𝑎𝑖𝑗 𝑏𝑗 , 𝑝𝑎𝑟𝑎 𝑖 = 1,2,3, … , 𝑛

𝑗=1
 22

Finalmente se obtiene la matriz {𝑋} de la relación [𝑈]{𝑋} = {𝐷} mediante la sustitución hacia
atrás, usando:
(𝑖−1) (𝑖−1)
𝑏𝑖 − ∑𝑛𝑗=𝑖+1 𝑎𝑖𝑗 𝑥𝑗
𝑥𝑖 = (𝑖−1)
; 𝑝𝑎𝑟𝑎 𝑖 = 𝑛 − 1, 𝑛 − 2, … , 1
𝑎𝑖𝑖

Pseudocodigo eliminación de Gauss por descomposición LU

Costo computacional: El costo computacional de la eliminación de Gauss por descomposición LU

𝑛3 𝑛 𝒏𝟑
es: 3
− 3 + 𝑛2 , 𝑐𝑜𝑛𝑓𝑜𝑟𝑚𝑒 𝑛 𝑎𝑢𝑚𝑒𝑛𝑡𝑎 → 𝟑
+ 𝑶(𝒏𝟐 ). Se puede observar que el costo total es
𝑛3 𝑛
idéntico al de la eliminación de Gauss. Donde 3 − 3 son los FLOPS de la descomposición de la
2
matriz A en las matrices L y U; y 𝑛 corresponde a los FLOPS para la sustitución hacia adelante y
hacia atrás.
 23

Descomposición LU con pivoteo

Si en la descomposición LU se pivotean filas ya no se cumple que [𝐿][𝑈] = [𝐴], sino que

[𝐿][𝑈] = [𝑃][𝐴]. Donde a la matriz [𝑃] se la conoce como matriz de permutación. Al inicio de la
descomposición [𝑃] = [𝐼], cada vez que se pivotean filas durante la descomposición, también se
pivotean las filas en [P].

Matriz inversa

La matriz inversa se puede calcular mediante la descomposición LU, usando como {B} a las
columnas de la matriz identidad. Donde {X} serán columnas de la matriz inversa. Para calcular la
primera columna de la inversa se utiliza la primera columna de la identidad el procedimiento es
exactamente igual para el resto de estas. EJ: una matriz 3x3, para calcular la primer columna:
−1
𝑎11 𝑎12 𝑎13 𝑎11 1
𝑎21 𝑎22 𝑎23 ∗ 𝑎−| = 0
21
𝑎31 𝑎32 𝑎33
𝑎31 0
−1

Pseudocodigo para la obtención de la matriz inversa

 24

*Nota: Este pseudocodigo se construyó a partir del pseudocodigo de la eliminación de Gauss con
descomposición LU. Se resaltó en negrita el código agregado.

Utilidad de la matriz inversa (condicionamiento del sistema)

La matriz inversa puede utilizarse para determinar si los sistemas de ecuaciones simultáneas están
mal condicionados. Hay tres formas de detectar el mal funcionamiento:

 Escalar la matriz A e invertirla. Si la inversa tiene elementos mucho mayores a 1 es

posible que el sistema este mal condicionado.
 Hacer el producto [𝐴] ∗ [𝐴−1 ] y estimar si el resultado es lo suficientemente
parecido a la matriz identidad.
 Invertir la matriz inversa, si el resultado no está lo suficientemente cercano a la
matriz original, el sistema está mal condicionado.

Normas vectoriales, matriciales y número de condición de una matriz

Una norma es una función que toma valores reales, que nos da una medida del tamaño o
“longitud” de entidades formadas por múltiples elementos como los vectores y las matrices.
 25

 ||𝑋||𝑒 y ||𝐴||𝑒 son normas euclidianas las cuales son la raíz cuadrada de la suma de los
cuadrados de los componentes de la matriz o vector respectivamente.
 ||𝑋||1 es la suma del valor absoluto de los elementos del vector.
 ||𝐴||1 conocida como norma columna-suma realiza sumatorias de los valores absolutos de los
elementos de cada columna de la matriz y la mayor de estas sumatorias se toma como norma.
 ||𝑋||∞ conocida como norma máxima la cual es igual al elemento de mayor valor absoluto.
 ||𝐴||∞ conocida como norma reglón-suma o matriz uniforme realiza sumatorias de los valores
absolutos de los elementos de cada reglón y la mayor de estas sumas se toma como norma.

Número de condición de una matriz: Este es igual a 𝑪𝒐𝒏𝒅[𝑨] = ||𝑨|| . ||𝑨−𝟏 || este número será
mayor o igual a uno. Es un producto punto entre los valores absolutos de los elementos de la
matriz y su inversa.

*Nota: En LibreOffice usar la función [Link] para hacer el producto punto.

Matrices bandeadas

Son matrices cuadradas en la que todos sus elementos son

cero a excepción de una banda centrada sobre la diagonal
principal. Sus dimensiones se cuantifican mediante el ancho
de banda, 𝐵𝑊, y el ancho de media banda,
𝐻𝐵𝑊, con 𝐵𝑊 = 2𝐻𝐵𝑊 + 1. En una matriz bandeada
𝒂𝒊𝒋 = 𝟎 si |𝒊 − 𝒋| > 𝑯𝑩𝑾. Si bien se pueden utilizar el
método de eliminación de Gauss o descomposición LU para
resolver estos sistemas resultan ineficientes. Esto se debe a
que si no fuera a ser necesario el pivoteo ninguno de los
elementos de la fuera de la banda dejara de ser cero, por lo
que desperdiciaria tiempo y espacio de almacenamiento
para almacenar estos ceros inútiles.

Algoritmo de Thomas

Es un algoritmo para solucionar sistemas tridiagonales, es

decir para resolver sistemas de ancho de banda tres (BW=3).
Para evitar guardar ceros inútiles utiliza solo los vectores e,
f, g y r. luego aplica descomposición LU, obteniendo las
matrices [L] y [U] que también son tridiagonales.
 26

Pseudocodigo algoritmo de Thomas

Descomposición de Cholesky

Utilizado para hallar una matriz triangular superior [𝐿] tal que [𝐴] = [𝐿] ∗ [𝐿]𝑇 , si 𝐴 es una
matriz simétrica. El resultado se expresa de forma simple mediante relaciones de recurrencia,
para el reglón k-ésimo:

𝑎𝑘𝑖 − ∑𝑖−1
𝑗=1 𝑙𝑖𝑗 𝑙𝑘𝑗
𝑙𝑘𝑖 = 𝑝𝑎𝑟𝑎 𝑖 = 1,2, … , 𝑘 − 1
𝑙𝑖𝑖

Para los elementos de la diagonal principal:

𝑘−1
2
𝑙𝑘𝑘 = √𝑎𝑘𝑘 − ∑ 𝑙𝑘𝑗
𝑗=1

*Nota: Esta matriz L no necesariamente tiene sus elementos de su diagonal principal igual a uno,
es de la forma:

𝑙11 0 0
𝑙21 𝑙22 0
𝑙31 𝑙32 𝑙33
 27

Desventajas de la descomposición de Cholesky: Puede presentar problemas de calcular la raíz a

la hora de evaluar los elementos de la diagonal principal de la matriz A. Se puede tener la raíz
cuadrada de un número negativo si la matriz A no es definida positiva.

Pseudocodigo descomposición de Cholesky

Métodos iterativos
Son métodos que se basan en suponer un conjunto de valores que satisfagan un conjunto de
ecuaciones y luego usar un método sistemático para obtener una mejor aproximación.

Método de Gauss-Seidel

Sea un sistema de ecuaciones [𝐴]{𝑋} = {𝐵}, si fuese un conjunto de ecuaciones de 3x3. Si los
elementos de la diagonal no son todos cero, se obtiene:

𝑏1 − 𝑎12 𝑥2 − 𝑎13 𝑥3 𝑏2 − 𝑎21 𝑥1 − 𝑎23 𝑥3 𝑏1 − 𝑎31 𝑥1 − 𝑎32 𝑥2

𝑥1 = , 𝑥2 = , 𝑥3 =
𝑎11 𝑎22 𝑎33

Una forma simple de obtener los valores iniciales es suponer que todos son cero 𝑥1 , 𝑥2 , 𝑥3 = 0.
Se calcula 𝒙𝟏 con 𝑥2 y 𝑥3 iniciales, 𝒙𝟐 se calcula con el 𝑥1 calculado previamente y 𝑥3 inicial, 𝒙𝟑 se
calcula con 𝑥1 y 𝑥2 calculados previamente. Luego se vuelven a hacer todos los cálculos pero
tomando como valores iniciales los previamente calculados. Este proceso se repite hasta que la
solución converja lo suficientemente cerca de los valores verdaderos.

Las incógnitas se pueden calcular de forma general para un sistema de n ecuaciones:

𝑏𝑖 − ∑𝑖−1 𝑛
𝑗=1 𝑎𝑖𝑗 𝑥𝑗 − ∑𝑗=𝑖+1 𝑎𝑖𝑗 𝑥𝑗
𝑥𝑖 = , 𝑐𝑜𝑛 𝑖 = 1,2, … , 𝑛
𝑎𝑖𝑖
 28

𝑗 𝑗−1
𝑥𝑖 −𝑥𝑖
Error del método: El error relativo porcentual es |𝜀𝑎,𝑖 | = | 𝑗 | ∗ 100, la 𝑖 denota la incógnita y
𝑥𝑖

𝑗, 𝑗 − 1 son las iteraciones actuales y previas, respectivamente.

Desventajas del método: En algunas ocasiones no es convergente y cuando converge, con

frecuencia lo hace de forma muy lenta

Criterio de convergencia: Para asegurar la convergencia, un criterio suficiente pero no necesario

es: los elementos de la diagonal deben ser mayores que los elementos fuera de la diagonal para
cada reglón, es decir |𝒂𝒊𝒊 | > ∑𝒏𝒋=𝟏|𝒂𝒊,𝒋 |. Esto se debe a que para evitar que el error aumente con
𝒋≠𝒊
respecto a la iteración anterior para cualquiera de las variables. La sumatoria de las pendientes de
𝜕𝑥
las funciones 𝑥𝑖 debe ser menor a uno, ∑𝑛𝑗=1 | 𝜕𝑥 𝑖 | < 1, para 𝑖 = 1,2, … , 𝑛; con 𝑥𝑗 como variables y
𝑗

𝑥𝑖 como función. (Ver iteración de punto fijo). Ejemplos de convergencia y divergencia para un
sistema 2x2, donde la función 𝑢 𝑦 𝑣 son las funciones correspondientes a los despejes de 𝑥1 𝑦 𝑥2 :

*Nota: Los sistemas que cumplen el criterio de convergencia anterior reciben el nombre de
sistemas diagonalmente dominantes.

Pseudocodigo del método de Gauss-Seidel

 29

Método de iteración de Jacobi

El método de Jacobi realiza aproximaciones de los valores que satisfacen el sistema (las X), solo
con los valores de la iteración anterior a la actual. A diferencia del método de Gauss-Seidel que usa
los valores de la iteración anterior y los de la actual.

Pseudocodigo método de Jacobi

*Nota: Aunque el método de Jacobi puede llegar a ser útil, el método de Gauss-Seidel da la mejor
aproximación.

Método de Gauss-Seidel con relajación

El uso de la relajación en el método de Gauss-Seidel mejora la convergencia. Este consiste en

después de calcular cada valor de 𝑥, modificar este mediante un promedio ponderado de la
iteración actual y anterior de la siguiente manera:

𝑥𝑖𝑛𝑢𝑒𝑣𝑜 = 𝜆𝑥𝑖𝑛𝑢𝑒𝑣𝑜 + (1 − 𝜆)𝑥𝑖𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟

Donde 𝝀 es un factor ponderado que tiene un valor entre 0 y 2.

Tipos de relajaciones:

 Subrelajación: Si 𝝀 se encuentra entre 0 y 1, se obtiene un promedio ponderado de los

resultados actuales y anteriores. Se utiliza para que un sistema que no converge, converja y
para apresurar la convergencia al amortiguar sus oscilaciones.
 Sobrerelajación: Si 𝝀 se encuentra entre 1 y 2, se utiliza en sistemas convergentes permite
acelerar la convergencia.
 30

La elección de un 𝝀 adecuado se determina de forma empírica para cada problema, su utilización

resulta útil cuando si el sistema bajo estudio se va a resolver de manera repetida.

Pseudocodigo método de Gauss-Seidel con relajación

Unidad 4: Optimización
La optimización es la búsqueda de un óptimo, el cual puede ser unidimensional o
multidimensional, dependiendo de si se involucran funciones de una variable o de más de una
variable respectivamente. Para los problemas unidimensionales y multidimensionales la
búsqueda del óptimo consiste en la búsqueda del mínimo o máximo o máximo.

Se utilizan para la determinar la solución óptima de un problema, equilibran el funcionamiento y

las limitaciones. Los problemas de optimización tenemos:
 31

 Función objetivo: función cuyo resultado se quiera optimizar, por ej una

función que calcule costos
 Variables de diseño: son variables que pueden ser números reales o enteros,
por ej: el radio r de la una semiesfera que representa un paracaídas.
 Restricciones: limitaciones debajo las cuales se trabaja.

Clasificación problemas de optimización

Generalmente los problemas de optimización se establecer como la determinación de un 𝑥 que

minimiza o maximiza 𝑓(𝑥), sujeto a:

𝑑𝑖 (𝑥) ≤ 𝑎𝑖 ; 𝑐𝑜𝑛 𝑖 = 1,2, … , 𝑚 𝑦 𝑒𝑖 (𝑥) = 𝑏𝑖 ; 𝑐𝑜𝑛 𝑖 = 1,2, … , 𝑝

Donde 𝑥 es un vector de diseño n-dimensional; 𝑓(𝑥) la función objetivo; 𝑑𝑖 (𝑥) son las
restricciones de desigualdad; 𝑒𝑖 (𝑥) son las restricciones de igualdad y 𝑎𝑖 , 𝑏𝑖 son constantes. Si
𝑝 + 𝑚 > 𝑛 se dice que el problema está sobrerestringido. Los problemas con estas restricciones
se conocen como problemas de optimización restringidos, de otra forma, se trata de un problema
de optimización no restringido.

Según la forma de 𝑓(𝑥) se clasifican los problemas de optimización como problemas de:

 Programación lineal: si 𝑓(𝑥) y las restricciones son lineales.

 Programación cuadrática: si 𝑓(𝑥) es cuadrática y las restricciones lineales.
 Programación no lineal: si 𝑓(𝑥) no es lineal o cuadrática y/o las restricciones no son
lineales.

Optimización unidimensional no restringida

Son técnicas para encontrar el mínimo o máximo de funciones de una variable. En general se está
interesado en la búsqueda de óptimo global, no del óptimo local. Una dificultad a la que se
enfrentan los problemas de optimización es la distinguir entre estos.

Método de búsqueda de la sección dorada (método cerrado)

Es un que utiliza cuatro puntos dentro de un intervalo determinado por dos valores iniciales (los
primeros dos puntos) 𝒙𝒍 y 𝒙𝒖 , que contengan un máximo o mínimo. Se necesitan tres puntos y
sus valores en la función para detectar si hay máximo o mínimo en un intervalo. Una vez que se
determinan el tercer y cuarto punto, se hace la prueba para detectar el máximo o mínimo está
entre los primeros o últimos tres puntos.

El método minimiza las evaluaciones de la función reemplazando los valores anteriores con los
nuevos, esto se logra si se satisfacen las siguientes dos condiciones:

𝑙1 𝑙2
𝑙0 = 𝑙1 + 𝑙2 (1) 𝑦 = (2)
𝑙0 𝑙1
 32

Reemplazando (1) en (2) y tomando 𝑅 = 𝑙2 /𝑙1 y

reordenando se obtiene:

𝟏
=𝟏+𝑹
𝑹

Despejando 𝑅 de 0 = 𝑅 2 + 𝑅 − 1 se obtiene:

√5 − 1
𝑅 = 0,6183 … =
2

Donde 𝑅 es conocida como razón dorada o razón áurea.

El tercer y cuarto punto se determinan usando la razón dorada de la siguiente manera:

√𝟓 − 𝟏
𝑻𝒆𝒓𝒄𝒆𝒓 𝒑𝒖𝒏𝒕𝒐: 𝒙𝟏 = 𝒙𝒍 + 𝒅; 𝑪𝒖𝒂𝒓𝒕𝒐 𝒑𝒖𝒏𝒕𝒐: 𝒙𝟐 = 𝒙𝒖 − 𝒅; 𝒄𝒐𝒏 𝒅 = ∗ (𝒙𝒖 − 𝒙𝒍 )
𝟐

Criterio para determinar en qué intervalo se encuentra el máximo:

 Si 𝑓(𝑥1 ) > 𝑓(𝑥2 ), entonces el intervalo la izquierda de 𝑥2

se puede eliminar y en la siguiente iteración 𝑥𝑙 = 𝑥2 y
𝑥2 = 𝑥1𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 .

 Si 𝑓(𝑥1 ) < 𝑓(𝑥2 ), entonces el intervalo a la derecha de 𝑥1

se puede eliminar y en la siguiente iteración 𝑥𝑢 = 𝑥1 y
𝑥1 = 𝑥2𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 .

Criterio para determinar en qué intervalo se encuentra el mínimo

 Si 𝑓(𝑥1 ) > 𝑓(𝑥2 ), entonces el intervalo la derecha de 𝑥1 se puede eliminar y en la

siguiente iteración 𝑥𝑢 = 𝑥1 y 𝑥1 = 𝑥2𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 .

 Si 𝑓(𝑥1 ) < 𝑓(𝑥2 ), el intervalo a la izquierda de 𝑥2 se puede eliminar y en la siguiente

iteración 𝑥𝑙 = 𝑥2 y 𝑥2 = 𝑥1𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟 .

Después de cada iteración el intervalo que encierra el máximo o mínimo se reduce en un factor de
la razón dorada. Por lo tanto después de n iteraciones el intervalo se acorta a 𝑹𝒏 % de su valor
inicial.

Error: Se puede calcular un límite superior para el error de cada iteración. Al final de cada
iteración se obtienen un nuevo intervalo y un punto dentro de este (el óptimo), el máximo error
posible es igual al subintervalo más grande de este nuevo intervalo. El error relativo aproximado
queda:
 33

𝒙𝒖 − 𝒙𝒍
𝓔𝒂 = (𝟏 − 𝑹) | | ∗ 𝟏𝟎𝟎
𝒙𝒐𝒑𝒕

Desventajas del método de la sección dorada: El método puede consumir mucho tiempo en la
evaluación de la función, si esta es compleja, ya que hace muchas evaluaciones por iteración.
Converge muy lentamente al extremo, necesita muchas iteraciones.

Pseudocodigo método de la sección dorada

Interpolación cuadrática (método cerrado)

La interpolación cuadrática hace uso de un polinomio de segundo grado el cual es una buena
aproximación de 𝑓(𝑥) en las cercanías de un óptimo. Si se tienen tres puntos que contienen un
máximo o mínimo se ajusta una parábola a estos. Después se puede derivar e igualar el
resultado a cero para encontrar el óptimo quedaría:

Donde 𝑥0 , 𝑥1 𝑦 𝑥2 son los valores iniciales y 𝑥3 es el valor estimado del óptimo.

−𝑏
*Nota: 𝑥3 se obtiene de donde 𝑎 y 𝑏 se obtienen despejando de: 𝑓(𝑥0 ) = 𝑎. 𝑥02 + 𝑏. 𝑥0 + 𝑐,
2𝑎
𝑓(𝑥1 ) = 𝑎. 𝑥12 + 𝑏. 𝑥1 + 𝑐 𝑦 𝑓(𝑥2 ) = 𝑎. 𝑥22 + 𝑏. 𝑥2 + 𝑐.

𝑥3𝑎𝑐𝑡𝑢𝑎𝑙 −𝑥3𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟
Error: El error relativo porcentual se calcula como ℰ𝑎 = | | ∗ 100
𝑥3𝑎𝑐𝑡𝑢𝑎𝑙
 34

Desventajas del método: Su convergencia puede ser lenta, ya que puede ser que durante las
iteraciones solo se retenga un extremo del intervalo.

Pseudocodigo método de la interpolación cuadrática

Método de Newton (método abierto)

Es un método abierto. Este método define una función 𝑔(𝑥) = 𝑓´(𝑥), tal que para el 𝑥𝑜𝑝𝑡 :

𝑓´(𝑥𝑜𝑝𝑡 ) = 𝑔(𝑥𝑜𝑝𝑡 ) = 0

De esta relación para estimar el óptimo se emplea lo siguiente:

𝒇´(𝒙𝒂𝒏𝒕𝒆𝒓𝒊𝒐𝒓
𝒐𝒑𝒕 )
𝒙𝒂𝒄𝒕𝒖𝒂𝒍
𝒐𝒑𝒕 = 𝒙𝒂𝒏𝒕𝒆𝒓𝒊𝒐𝒓
𝒐𝒑𝒕 −
𝒇´´(𝒙𝒂𝒏𝒕𝒆𝒓𝒊𝒐𝒓
𝒐𝒑𝒕 )

*Nota: La función 𝑔(𝑥) se obtiene de la serie de Taylor de segundo orden para 𝑓(𝑥) igualando la
serie a cero.

Desventajas del método: Al ser un método abierto el método de Newton puede divergir,
dependiendo del comportamiento de la función y del valor inicial.

𝑥3𝑎𝑐𝑡𝑢𝑎𝑙 −𝑥3𝑎𝑛𝑡𝑒𝑟𝑖𝑜𝑟
Error: El error relativo porcentual se calcula como ℰ𝑎 = | | ∗ 100
𝑥3𝑎𝑐𝑡𝑢𝑎𝑙
 35

Pseudocodigo método de Newton

Optimización multidimensional no restringida

Son técnicas para encontrar el mínimo o el máximo de funciones de varias variables. Estas
técnicas se pueden clasificar como métodos sin gradiente o directos son métodos que no necesitan
la evaluación de la derivada y métodos de gradientes o de descenso (o ascenso).

Método de la búsqueda aleatoria (método directo)

Este método evalúa la función seleccionando aleatoriamente valores de la variable

independiente. Si el método se lleva a cabo con un número suficiente de muestras el óptimo
eventualmente se localizara. Este método funciona para cualquier tipo de funciones.

Los generadores de números aleatorios dan números entre cero y uno, designando tal número
como 𝑟. Luego una fórmula para generar números para cualquier variable dentro de un intervalo
es: 𝒙 = 𝒙𝒍 + (𝒙𝒖 − 𝒙𝒍 ). 𝒓; con 𝑥 que es una variable cualesquiera y 𝑥𝑙 , 𝑥𝑢 son los extremos del
intervalo para esta variable.

Desventajas del método: Es costoso computacionalmente ya que al no tomar en cuenta la forma

de la función, el método es muy ineficiente.

Pseudocodigo método de la búsqueda aleatoria

 36

Búsquedas univariadas y direcciones patrón

Búsqueda univariada: Consiste en trabajar con solo una

variable a la vez mientras las otras permanecen constantes,
por lo que el método resulta a una secuencia de búsquedas
en una dimensión utilizando. En el caso que se tengan dos
variables 𝑥 e 𝑦, se parte de un punto inicial 1 y se deja fija
una de las variables, en este caso 𝑦 luego se mueve a lo largo
del eje 𝑥 hasta encontrar el máximo (hasta encontrar la
función) en el punto 2. Luego se hace lo mismo solo con 𝑥 fijo
hasta encontrar el próximo máximo, 3. Este proceso se repite
hasta que se considere que se ha aproximado lo suficiente al
óptimo.

Direcciones patrón: las direcciones patrón son líneas que

unen puntos alternados como 1-3, 3-5 o 2-4, 4-6 que apuntan
al óptimo. Estas trayectorias nos llevan directamente a el máximo.

*Nota: [1-2 (o 3-4; 5-6) y 2-4]; [2-3 (o 4-5) y 3-5] son direcciones conjugadas

Método Powell (método directo)

Es un método de convergencia cuadrática. Se basa en que si dos puntos

se obtienen por búsqueda en una dimensión en la misma dirección pero
con diferentes puntos de partida (como 1 y 2), la línea que los une está
dirigida al máximo. Tales líneas se llaman direcciones conjugadas.

El método inicia la búsqueda en un punto 0 con direcciones

iniciales ℎ1 y ℎ2 las cuales no son necesariamente conjugadas. Se
mueve por la dirección ℎ1 hasta que encuentra un máximo en un
punto 1. Después se mueve desde 1 por la dirección ℎ2 hasta
encontrar el máximo en un punto 2. Ahora se define la dirección
ℎ3 que se obtiene de la unión de los puntos 0 y 2, se mueve por
esta hasta encontrar el máximo en el punto 3. Desde 3 se mueve
en la dirección ℎ2 hasta llegar al máximo en 4. Desde 4 se mueve
en la dirección ℎ3 hasta llegar al máximo en 5. Se define la
 37

dirección ℎ4 que se obtiene de la unión de los puntos 3 y 5. Se puede ver que ℎ3 y ℎ4 son
direcciones conjugadas. Así buscando desde el punto 5 por ℎ4 nos llevara al máximo.

Derivada direccional

Nos da la pendiente en una dirección arbitraria ℎ que forma un ángulo 𝜃con el eje 𝑥:

𝜕𝑓 𝜕𝑓
𝑔´(0) = cos(𝜃) + 𝑠𝑒𝑛(𝜃)
𝜕𝑥 𝜕𝑦

El gradiente

Nos indica de manera local nos indica la dirección de la máxima pendiente:

𝜕𝑓 𝜕𝑓
∇𝑓 = 𝒊+ 𝒋
𝜕𝑥 𝜕𝑦

El Hessiano

Es el equivalente a la segunda derivada en funciones univariadas, pero para funciones

multivariadas. El Hessiano, |𝐻|, es igual al determinante de una matriz Hessiana, 𝐻, formada con
las segundas derivadas:

𝜕𝑓 𝜕𝑓
𝜕𝑥 2 𝜕𝑥𝜕𝑦
|𝐻| =
𝜕𝑓 𝜕
𝜕𝑦𝜕𝑥 𝜕𝑦 2

Método de máxima inclinación

Este método comienza desde un punto inicial y determina la

dirección ℎ0 de la pendiente máxima utilizando el gradiente, para
luego moverse por esta dirección hasta encontrar un máximo o
mínimo. Una vez encontrado, se vuelve a evaluar el gradiente para
obtener una nueva dirección ℎ1 por la que moverse hasta encontrar
un nuevo máximo o mínimo. Este proceso se repite hasta encontrar
el óptimo. El método tiene convergencia lineal.

Comenzando en un punto 𝑥𝑖 y 𝑦𝑖 cualquier punto en la dirección del

gradiente, ℎ𝑖 , se expresa como

𝜕𝑓 𝜕𝑓
𝑥 = 𝑥𝑖 + ℎ; 𝑦 = 𝑦𝑖 + ℎ
𝜕𝑥 𝑖 𝜕𝑦 𝑖

Encontrar la localización del máximo o mínimo ℎ∗ en la dirección del gradiente equivale a

encontrar el máximo de una función univariada con única variable ℎ.
 38

𝜕2 𝑓
Si |𝐻|(𝐻𝑒𝑠𝑠𝑖𝑎𝑛𝑜) > 0 y 𝜕𝑥 2 < 0, entonces en 𝑓(𝑥𝑜𝑝𝑡 , 𝑦𝑜𝑝𝑡 ) se tiene un máximo.
𝑜𝑝𝑡

𝜕2 𝑓
Si |𝐻|(𝐻𝑒𝑠𝑠𝑖𝑎𝑛𝑜) > 0 y 𝜕𝑥 2 > 0, entonces en 𝑓(𝑥𝑜𝑝𝑡 , 𝑦𝑜𝑝𝑡 ) se tiene un mínimo.
𝑜𝑝𝑡

Si |𝐻|(𝐻𝑒𝑠𝑠𝑖𝑎𝑛𝑜) < 0, se tiene un punto silla.

Desventajas del método de máxima inclinación: Tiende a converger muy lentamente cuando la
función presenta crestas largas y angostas, lo que le lleva a dar muchos pasos hasta llegar al
óptimo.

Unidad 5: Ajuste de Curvas

En esta unidad se ven técnicas para ajustar curvas a un conjunto de datos discretos para obtener
estimaciones intermedias (continuas). También se ven técnicas para obtener una versión
simplificada de una función complicada.

Existen dos métodos generales para el ajuste de curvas:

- Si los datos presentan mucho error se utiliza la regresión por mínimos cuadrados, esta
construye una curva que siga la tendencia de los puntos tomados como un grupo.
- Si los datos son muy precisos se coloca una curva o una serie de curvas que pasen por
cada uno de los puntos. La estimación de valores entre este tipo de datos se conoce como
interpolación.

Aplicaciones:

- Análisis de tendencia: Sirve para predecir o pronosticar valores de la variable

dependiente, dentro del intervalo de los datos o más allá.
- Prueba de hipótesis: Sirve para validar un modelo matemático existente con los
resultados experimentales o adecuar el modelo a los datos.

Regresión por mínimos cuadrados

Regresión lineal

Ajusta una línea recta a un conjunto de datos:

𝑦 = 𝑎0 + 𝑎1 𝑥 + 𝑒

Con 𝑎0 y 𝑎1 representan la intersección con el eje 𝑥 y la pendiente y 𝑒

es la discrepancia entre el valor verdadero y el valor aproximado.

𝑒 = 𝑦 − 𝑎0 − 𝑎1 𝑥
 39

Criterios inadecuados para un “mejor “ajuste

Minimizar la suma de residuos: (para 𝑛 puntos)

∑𝑛𝑖=1 𝑒𝑖 = ∑𝑛𝑖=1 𝑦𝑖 − 𝑎0 − 𝑎1 𝑥

Problema: Con minimizar los residuos es que cualquier otra

línea recta que pase por el punto medio de la línea que ajusta
a los datos, hace que el error sea cero. Ya que los errores se
cancelan.

Minimizar la suma de los valores absolutos de las

discrepancias:

∑𝑛𝑖=1 |𝑒𝑖 | = ∑𝑛𝑖=1 |𝑦𝑖 − 𝑎0 − 𝑎1 𝑥|

Problema: Cualquier línea entre las líneas punteadas

minimizara el valor absoluto de las sumas de las discrepancias.

Criterio minimax: Minimiza la máxima distancia a la que un

punto se encuentra de la línea que ajusta los datos.

Problema: Es sensible a puntos con mucho error.

Mejor criterio para el ajuste: Minimiza la suma de los cuadrados de los residuos entre la
𝑦𝑚𝑒𝑑𝑖𝑑𝑎 y 𝑦𝑚𝑜𝑑𝑒𝑙𝑜
𝑛 𝑛

𝑆𝑟 = ∑ 𝑒𝑖2 = ∑(𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 )2
𝑖=1 𝑖=1

Ajuste de una línea por mínimos cuadrados

Para determinar 𝑎0 y 𝑎1 derivamos 𝑆𝑟 con respecto a estos:

𝑛 𝑛
𝜕𝑆𝑟 𝜕𝑆𝑟
= −2 ∑(𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 ) ; = −2 ∑[𝑥𝑖 (𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 )] ;
𝜕𝑎0 𝜕𝑎1
𝑖=1 𝑖=1

Igualando ambas ecuaciones a 0, correspondiente al 𝑆𝑟 mínimo y despejando se obtiene:

 40

𝑛 ∑ 𝑥𝑖 𝑦𝑖 −∑ 𝑦𝑖 ∑ 𝑥𝑖
𝑎0 = (𝑦̅ − 𝑎1 𝑥̅ ) 𝑦 𝑎1 = 𝑛 ∑ 𝑥𝑖2 −(∑ 𝑥𝑖 )2

Error en la regresión lineal

Si la dispersión de los puntos alrededor de la línea de ajuste es de

magnitud similar en todos los datos y la distribución de los puntos
cerca de la línea es normal, la regresión por mínimos cuadrados
proporcionara la mejor estimación de 𝑎0 y 𝑎1 . Esto se conoce como
principio de máxima verosimilitud.

𝑟 𝑆
Error estándar del estimado: 𝑆𝑦/𝑥 = √𝑛−2 es el error para un valor
predicho de 𝑦 correspondiente a un valor particular de 𝑥. Este error
cuantifica la dispersión de los datos alrededor de la línea de regresión.

𝑺𝒕 −𝑺𝒓
Coeficientes de determinación y correlación: 𝒓𝟐 = 𝑦 𝑟 son los coeficientes de
𝑺𝒕
determinación y correlación, respectivamente. El coeficiente 𝑟 2 representa cuanto explica la línea
la variabilidad de los datos. Con 𝑺𝒕 = (𝒚𝒊 − 𝒚 ̅)𝟐 . La diferencia 𝑆𝑡 − 𝑆𝑟 cuantifica la mejora o
reducción del error y al dividir por 𝑆𝑡 esta se normaliza. El coeficiente de determinación 𝑟 2 toma
valores entre 0 y 1 representando la explicación del 100% y 0% de los datos respectivamente. El
coeficiente de correlación se puede expresar como:

Linealización de relaciones no lineales

Consiste en usar manipulaciones matemáticas para transformar funciones en las que la relación
entre su variable independiente y dependiente no es lineal, para que si lo sea. Después se utiliza la
regresión lineal para ajustar los datos.
 41

Por ejemplo:

𝑦 = 𝛼1 𝑒 𝛽1𝑥 → ln(𝑦) = ln(𝛼1 ) + 𝛽1 . 𝑥. 𝑙𝑛(𝑒)

𝑦 = 𝛼2 𝑥 𝛽2 → log(𝑦) = log(𝛼2 ) + 𝛽2 . log(𝑥)

𝑥 1 𝛽3 1 1
𝑦 = 𝛼3 → = +
𝛽3 +𝑥 𝑦 𝛼3 𝑥 𝛼3

Regresión polinomial

La regresión polinomial consiste en ajustar polinomios a los datos. Se utiliza en los casos en el que
las variables independientes y dependientes de una función no presentan una relación lineal, al
igual que en el caso de las transformaciones.

Los coeficientes del polinomio utilizado para el ajuste se pueden obtener de la misma manera los
coeficientes de la regresión lineal. Ya que la regresión lineal puede considerarse como un ajuste
mediante un polinomio de grado uno.

Criterio para el ajuste: Minimiza la suma de los cuadrados de los residuos entre la
𝑦𝑚𝑒𝑑𝑖𝑑𝑎 y 𝑦𝑚𝑜𝑑𝑒𝑙𝑜
 42

𝑛 𝑛

𝑆𝑟 = ∑ 𝑒𝑖2 = ∑(𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 − 𝑎2 𝑥𝑖2 … − 𝑎𝑛 𝑥𝑖𝑛 )2

𝑖=1 𝑖=1

Se hacen las derivadas parciales de 𝑆𝑟 con respecto a cada uno de los coeficientes y se igualan a
cero, obteniéndose un sistema de ecuaciones, de donde se despejan los coeficientes.

Para el caso de una regresión con un polinomio de segundo grado:

𝑛 𝑛
𝜕𝑆𝑟 𝜕𝑆𝑟
= −2 ∑(𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 − 𝑎2 𝑥𝑖2 ) ; = −2 ∑[𝑥𝑖 (𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 − 𝑎2 𝑥𝑖2 )] ;
𝜕𝑎0 𝜕𝑎1
𝑖=1 𝑖=1

𝜕𝑆𝑟
= −2 ∑[𝑥𝑖2 (𝑦𝑖 − 𝑎0 − 𝑎1 𝑥𝑖 − 𝑎2 𝑥𝑖2 )]
𝜕𝑎2

Se iguala a cero y se arma el siguiente sistema de ecuaciones, del cual se obtienen los
coeficientes:

𝑛 ∑ 𝑥𝑖 ∑ 𝑥𝑖2 ∑ 𝑦𝑖
𝑎0
∑ 𝑥𝑖 ∑ 𝑥𝑖2 ∑ 𝑥𝑖3 ∗ 𝑎1 = ∑ 𝑥𝑖 𝑦𝑖
𝑎2
∑ 𝑥𝑖2 ∑ 𝑥𝑖3 ∑ 𝑥𝑖4 ∑ 𝑥𝑖2 𝑦𝑖

𝒓 𝑺
Error estándar: 𝑆𝒚/𝒙 = √𝒏−(𝒎+𝟏) , donde 𝑚 es el grado del polinomio utilizado en el ajuste.

𝑆𝑡 −𝑆𝑟
Coeficiente de correlación: 𝑟 2 = 𝑆𝑡
, con 𝑆𝑡 = ∑(𝑦𝑖 − 𝑦̅)2 .

Pseudocodigo regresión polinomial

 43

Regresión lineal múltiple

Se utiliza en el caso de funciones con más de una variable

independiente, por ejemplo 𝑦 = 𝑓(𝑥1 , 𝑥2 ). En este caso la línea
de regresión se vuelve un plano. El criterio sigue siendo el mismo
que en la regresión lineal y polinomial:

Los coeficientes se obtienen de la matriz obtenida de las

ecuaciones de las derivadas parciales de 𝑆𝑟 con respecto a sus
coeficientes 𝑎0 , 𝑎1 y 𝑎2 , igualadas a cero.

𝒓𝑺
Error estándar: 𝑆𝒚/𝒙 = √𝒏−(𝒎+𝟏)

𝑺𝒕 −𝑺𝒓
Coeficiente de determinación: 𝒓𝟐 = 𝑺𝒕
, con 𝑆𝑡 = (𝑦𝑖 − 𝑦̅)2 .

Interpolación
Se utiliza para la estimación entre datos intermedios entre datos que suelen ser datos definidos
por puntos. La interpolación polinomial, consiste en la determinación de un único polinomio de n-
ésimo grado que se ajusta a 𝑛 + 1 puntos. Este polinomio proporciona una forma de estimar los
valores intermedios. La interpolación polinomial es sensible al error de redondeo y a los puntos
lejanos.

Interpolación lineal de Newton en diferencias divididas

Interpolación lineal: Consiste en unir dos puntos con una línea recta. Por triángulos semejantes se
obtiene:
 44

Donde el subíndice de 𝑓1 (𝑥) indica que designa que este es

un polinomio de primer grado. Conforme el intervalo entre
los datos disminuya, mejor será la aproximación. Cabe notar
que 𝑏1 es similar a la aproximación en diferencias divididas de
la primera derivada.

Interpolación cuadrática: Mejora la estimación en comparación a la interpolación lineal. Si se

tienen tres puntos el método consiste en ajustar un polinomio de segundo grado a estos. Su
fórmula es:

𝑓2 (𝑥) = 𝑏0 + 𝑏1 (𝑥 − 𝑥0 ) + 𝑏2 (𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) (1)

Aprovechando el hecho de que 𝑓2 (𝑥0 ) = 𝑓(𝑥0 ), 𝑓2 (𝑥1 ) = 𝑓(𝑥1 ) y 𝑓2 (𝑥2 ) = 𝑓(𝑥2 ) se despejan los
coeficientes de la siguiente manera:

Evaluando 𝑥0 , 𝑓2 (𝑥0 ), queda: 𝑏0 = 𝑓(𝑥0 )

𝑓(𝑥1 )−𝑓(𝑥0 )
Evaluando 𝑥1 , 𝑓2 (𝑥1 ), reemplazando 𝑏0 y despejando, queda: 𝑏1 =
𝑥1 −𝑥0

𝑓(𝑥2 )−𝑓(𝑥1 ) 𝑓(𝑥1 )−𝑓(𝑥0 )

−
𝑥2 −𝑥1 𝑥1 −𝑥0
Evaluando 𝑥2 , 𝑓2 (𝑥2 ), reemplazando 𝑏0 𝑦 𝑏1 y despejando, queda: 𝑏2 =
𝑥2 −𝑥0

Se puede observar que los dos primeros términos de (1) son semejantes a la interpolación lineal y
que 𝑏2 es similar a la aproximación en diferencias divididas de la segunda derivada. Al igual que en
la interpolación lineal conforme menor sea el intervalo entre los datos mejor será la estimación.

Interpolación de Newton forma general: Generalizando lo anterior. Para 𝑛 + 1 datos el

polinomio de n-ésimo grado es:

𝑓𝑛 (𝑥) = 𝑏0 + 𝑏1 (𝑥 − 𝑥0 ) + 𝑏2 (𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) + ⋯ + 𝑏𝑛 (𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) … (𝑥 − 𝑥𝑛 )

Los coeficientes se evalúan como:

 45

Donde las evaluaciones de la función puesta entre paréntesis son diferencias divididas finitas.

La primera, segunda y n-ésima diferencia dividida finita es de la forma:

𝑓(𝑥𝑖 )−𝑓(𝑥𝑗 ) 𝑓[𝑥𝑖 ,𝑥𝑗 ]−𝑓[𝑥𝑗 ,𝑥𝑘 ]

𝑓[𝑥𝑖 , 𝑥𝑗 ] = 𝑥𝑖 −𝑥𝑗
𝑓[𝑥𝑖 , 𝑥𝑗 , 𝑥𝑘 ] = 𝑥𝑖 −𝑥𝑘
𝑓[𝑥𝑛 ,𝑥𝑛−1 ,…,𝑥1 ]−𝑓[𝑥𝑛−1 ,𝑥𝑛−2 ,…,𝑥0 ]
𝑓[𝑥𝑛 , 𝑥𝑛−1 , … , 𝑥1 , 𝑥0 ] = 𝑥𝑛 −𝑥0

Luego la expresión:

𝒇𝒏 (𝒙) = 𝒇(𝒙𝟎 ) + 𝒇[𝒙𝟏 , 𝒙𝟎 ](𝒙 − 𝒙𝟎 ) + 𝒇[𝒙𝟐 , 𝒙𝟏 , 𝒙𝟎 ](𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 ) + ⋯

+ 𝒇[𝒙𝒏 , 𝒙𝒏−𝟏 , … , 𝒙𝟎 ](𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 ) … (𝒙 − 𝒙𝒏 )

Es conocida como polinomio de interpolación de Newton en diferencias divididas, las ecuaciones

de polinomios de orden ≥ 2 son recursivas.

*Nota: Si la función original es un polinomio de n-ésimo grado el polinomio de interpolación de n-

ésimo grado basado en n+1 puntos da resultados exactos.

Errores: El polinomio de interpolación es similar a la serie de Taylor, en el aspecto de que cada

término aproxima el comportamiento del polinomio a la función original, ya que las diferencias
divididas de cada término aproximan las derivadas de los distintos órdenes. De esta manera el
error de truncamiento del polinomio de interpolación es análogo al de la serie de Taylor, este es
igual a:

𝑓 𝑛+1 (ℰ)
𝑅𝑛 = (𝑥 − 𝑥0 )(𝑥 − 𝑥1 ) … (𝑥 − 𝑥𝑛 )
(𝑛 + 1)!

Una fórmula alternativa que requiere conocer un punto más (𝑥𝑛+1 , 𝑓(𝑥𝑛+1 )), pero no necesita
conocer la función original es:

𝑹𝒏 ≅ 𝒇[𝒙𝒏+𝟏 , 𝒙𝒏 , 𝒙𝒏−𝟏 , … , 𝒙𝟎 ](𝒙 − 𝒙𝟎 )(𝒙 − 𝒙𝟏 ) … (𝒙 − 𝒙𝒏 )

 46

Polinomios de interpolación de Lagrange

Es una reformulación del polinomio de interpolación de Newton, este evita el cálculo de las
diferencias divididas, se presenta como:
𝒏 𝒏
𝒙 − 𝒙𝒋
𝒇𝒏 (𝒙) = ∑ 𝑳𝒊 (𝒙)𝒇(𝒙𝒊 ) , 𝑐𝑜𝑛 𝑳𝒊 (𝒙) = ∏
𝒙𝒊 − 𝒙𝒋
𝒊=𝟎 𝒋=𝟎
𝒋≠𝒊

Para n=1 y n=2 se obtiene:

El razonamiento detrás de la formulación está en que

𝐿𝑖 (𝑥) = 1, cuando 𝑥 = 𝑥𝑖 y 0 en cualquier otro punto. Por lo
tanto 𝑓𝑛 (𝑥) = 𝑓(𝑥), cuando 𝑥 = 𝑥𝑖 . En consecuencia la
sumatoria de todos estos productos nos da el único polinomio
de n-ésimo que pasa por los n+1 datos que se tienen como
datos.

Error: Como solo hay un único polinomio de n-ésimo grado

que pase por n+1 puntos, el polinomio de interpolación de
Newton y el de Lagrange son el mismo y tienen el mismo error:
𝒏

𝑹𝒏 ≅ 𝒇[𝒙𝒏+𝟏 , 𝒙𝒏 , 𝒙𝒏−𝟏 , … , 𝒙𝟎 ] ∏( 𝒙 − 𝒙𝒏 )
𝒊=𝟎

Polinomio de Newton vs polinomio de Lagrange

El polinomio de Newton es ventajoso cuando se desconoce el grado del polinomio, es decir para
cálculos exploratorios. Ya que al tener formulas recursivas, a la hora de generar un polinomio
puede reutilizar cálculos realizados para polinomios de menor grado. Esto permite hacer pruebas
con polinomios de diferentes grados con menor costo computacional.

El polinomio de Lagrange generalmente se usa cuando el grado del polinomio se conoce a priori.
Ya que aunque tiene un costo computacional semejante al polinomio de Newton y no necesita
almacenar diferencias divididas.
 47

Interpolación mediante trazadores (Splines)

La interpolación mediante trazadores consiste en colocar polinomios de grado inferior en

subconjuntos de los datos, tales polinomios se denominan trazadores o splines. Esta técnica no es
tan sensible a los errores de redondeo ni a la influencia de puntos lejanos.

Los trazadores tienen mejor comportamiento que un polinomio de grado superior en funciones
que presentan un cambio abrupto cerca de la región de interés. Se puede observar desde la
figura 𝑎) a 𝑐) que a medida de que el grado del polinomio aumenta, se forman oscilaciones
bruscas cerca del cambio abrupto.

Trazadores lineales: Los trazadores de primer grado o lineales, se definen como un conjunto de
ecuaciones lineales:

𝑓(𝑥) = 𝑓(𝑥0 ) + 𝑚0 (𝑥 − 𝑥0 ), 𝑥0 ≤ 𝑥 ≤ 𝑥1

𝑓(𝑥) = 𝑓(𝑥1 ) + 𝑚1 (𝑥 − 𝑥1 ), 𝑥1 ≤ 𝑥 ≤ 𝑥2

∶ ∶ ∶ ∶

𝑓(𝑥) = 𝑓(𝑥𝑛−1 ) + 𝑚𝑛−1 (𝑥 − 𝑥𝑛−1 ), 𝑥𝑛−1 ≤ 𝑥 ≤ 𝑥𝑛

Con:

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

De esta manera se puede estimar cualquier punto en el intervalo (𝑥0 , 𝑥𝑛 ), localizando el

subintervalo en cual se encuentra y luego utilizando la ecuación correspondiente.
 48

Trazadores cuadráticos: Los trazadores cuadráticos son polinomios de segundo grado. Estos
tienen primeras derivadas continuas en todos los nodos (los extremos de su intervalo). No
aseguran derivadas iguales en los nodos.

Para 𝑛 + 1 puntos, existen 𝑛 intervalos con

un polinomio de segundo grado por
intervalo y en consecuencia 3𝑛 coeficientes
(𝑎, 𝑏 y 𝑐) desconocidos. Por lo tanto se
requieren de 3𝑛 ecuaciones o condiciones
para evaluar las incógnitas:

- Los valores de los polinomios adyacentes deben ser iguales en los nodos interiores,
esto nos da 𝟐𝒏 − 𝟐 condiciones.
2
𝑎𝑖−1 𝑥𝑖−1 + 𝑏𝑖−1 𝑥𝑖−1 + 𝑐𝑖−1 = 𝑓(𝑥𝑖−1 )
2
𝑎𝑖 𝑥𝑖−1 + 𝑏𝑖 𝑥𝑖−1 + 𝑐𝑖 = 𝑓(𝑥𝑖−1 )
- La primera (i=1) y última función (i=n-1) deben pasar a través de los puntos extremos,
esto nos da 2 condiciones.
𝑎1 𝑥02 + 𝑏1 𝑥0 + 𝑐1 = 𝑓(𝑥0 )
𝑎𝑛−1 𝑥𝑛2 + 𝑏𝑛−1 𝑥𝑛 + 𝑐𝑛−1 = 𝑓(𝑥𝑛 )

- Las primeras derivadas en los nodos deben ser iguales, esto nos da 𝒏 − 𝟏 condiciones.
2𝑎𝑖−1 𝑥𝑖−1 + 𝑏𝑖−1 = 2𝑎𝑖 𝑥𝑖−1 + 𝑏𝑖
- Se supone que en el primer punto la segunda derivada es cero, última condición, de
donde queda: (en consecuencia los dos primeros puntos se unen con una línea recta)
2𝑎1 = 0 → 𝑎1 = 0

Trazadores cúbicos: Se basa en tener un polinomio de tercer grado por cada intervalo entre
nodos. Por lo tanto para 𝑛 + 1 datos se tienen 𝟒𝒏 incógnitas (𝑎, 𝑏, 𝑐 y 𝑑) y por lo tanto se
requieren la misma cantidad de ecuaciones o condiciones. Las primeras tres condiciones son
iguales a la de los trazadores cuadráticos y las restantes son:

 Las segundas derivadas de los nodos deben ser iguales, esto nos da 𝒏 − 𝟏 condiciones.
6𝑎𝑖−1 𝑥𝑖−1 + 2𝑏𝑖−1 = 6𝑎𝑖 𝑥𝑖−1 + 2𝑏𝑖
 Las segundas derivadas en los nodos extremos son cero, esto nos da las últimas 𝟐
condiciones. En consecuencia la función se vuelve una línea recta nodos extremos.
6𝑎0 𝑥0 + 2𝑏0 = 0

6𝑎𝑛−1 𝑥𝑛 + 2𝑏𝑛−1 = 0

En la última condición si las segundas derivadas no son cero en los nodos extremos se puede
utilizar como alternativa las siguientes dos condiciones:

6𝑎0 𝑥0 + 2𝑏0 ≠ 0

6𝑎𝑛−1 𝑥𝑛 + 2𝑏𝑛−1 ≠ 0
 49

Trazadores cúbicos forma alternativa (𝒏 − 𝟏 incógnitas): Como cada par de nodos está unido por
un polinomio de tercer grado, la segunda derivada de este es una recta. Ahora para un intervalo
entre los nodos 𝑥𝑖−1 y 𝑥𝑖 , representando la recta con un polinomio de interpolación de Lagrange:

Está expresión une las derivadas segundas del primer nodo con los del segundo, luego integrando
dos veces obtenemos la función original:

La ventaja de representar el polinomio de segundo grado de esta manera, es que solo se deben
determinar dos coeficientes 𝑓𝑖´´ (𝑥𝑖−1 ) y 𝑓𝑖´´ (𝑥𝑖 ). Diferente de determinar 𝑎, 𝑏, 𝑐, y 𝑑 de un
polinomio convencional.

Derivando la última expresión para los intervalos (𝑖 − 1; 𝑖) y (𝑖; 𝑖 + 1), utilizando la condición
𝒇´𝒊−𝟏 (𝒙𝒊 ) = 𝒇´𝒊 (𝒙𝒊 ) y reordenando se llega a la siguiente expresión:

Esta ecuación es válida para todos los nodos interiores lo que nos da 𝑛 − 1 ecuaciones y 𝑛 + 1
coeficientes a determinar. Utilizando la condición de que las segundas derivadas en los nodos
extremos deben ser cero, 𝑓´´(𝑥0 ) = 𝑓´´(𝑥𝑛 ) = 0. Por lo que el problema se reduce a un sistema de
𝑛 − 1 ecuaciones y 𝑛 − 1 incognitas.
 50

Unidad 6: Diferenciación e integración numéricas

Métodos sin computadoras:

Diferenciación grafica por áreas iguales: Se

utiliza para determinar la derivada a partir de
datos conocidos. Se tabulan los datos 𝑥 e 𝑦,
luego hace la diferencia dividida ∆𝑦/∆𝑥 para
estimar la derivada de los datos conocidos.
Después se dibuja una curva escalonada que
se suaviza tratando de balancear las áreas
negativas y positivas.

Métodos de cuadraturas: Se utilizan para estimar las

derivadas. Las alturas de la función se multiplican por el
ancho de las barras y se suman, obteniendo una
aproximación del área debajo de la curva.

Integración a través de datos tabulados

Las técnicas que se muestran a continuación son para uso de información tabulada, por lo que se
está limitado a la cantidad de puntos que se tengan.

Fórmulas de integración de Newton-Cotes (formulas cerradas)

Se basa en reemplazar funciones complicadas o datos tabulados por un polinomio fácil de

integrar;
𝑏 𝑏
𝐼 = ∫𝑎 𝑓(𝑥)𝑑𝑥 ≅ ∫𝑎 𝑓𝑛 (𝑥)𝑑𝑥

Donde 𝑓𝑛 (𝑥) es un polinomio de n-ésimo grado.

Existen formas abiertas y cerradas de este método.

Las formas cerradas son aquellas en las que se
conocen los datos al inicio y al final de los límites de
integración. Las formas abiertas tienen límites de
integración que se encuentran más allá del intervalo
de los datos.

Generalmente las formas abiertas se utilizan para evaluar integrales definidas, sino integrales
impropias y E.D.O.
 51

Regla del trapecio (formula cerrada): Caso en el que la función reemplazada es un

𝑏
polinomio de primer grado, 𝐼 ≅ ∫𝑎 𝑓1 (𝑥)𝑑𝑥 .

La integral de 𝑓(𝑥) se aproxima con el área debajo de la línea recta que une la función en
los dos extremos del intervalo.
𝑏
𝑓(𝑏) − 𝑓(𝑎) 𝒇(𝒂) + 𝒇(𝒃)
𝐼 ≅ ∫ 𝑓(𝑎) + (𝑥 − 𝑎)𝑑𝑥 → 𝑰 ≅ (𝒃 − 𝒂)
𝑎 𝑏−𝑎 𝟐

Esta fórmula se denomina regla del trapecio.

Error: El error de truncamiento local verdadero para una sola aplicación es:

1
𝐸𝑡 = − 𝑓´´(ℰ)(𝑏 − 𝑎)3 , 𝑐𝑜𝑛 ℰ𝜖(𝑎, 𝑏)
12
Y el error de truncamiento local aproximado es:

1
𝐸𝑎 = − ̅
𝑓´´(𝑥)(𝑏 − 𝑎)3
12
𝑏
∫𝑎 𝑓´´(𝑥)𝑑𝑥
̅
𝐶𝑜𝑛 𝑓 ´´(𝑥) = 𝑏−𝑎

Si la función que se quiere integrar es exacta la regla del trapecio es exacta

Recordatorio
Primera diferencia dividida hacia adelante: Es una forma de estimar la derivada
𝑓(𝑥𝑖+1 ) − 𝑓(𝑥𝑖 ) ∆𝒇𝒊
𝑓´(𝑥𝑖 ) = + 𝑂(𝑥𝑖+1 − 𝑥𝑖 ) = + 𝑶(𝒉)
𝑥𝑖+1 − 𝑥𝑖 𝒉
Donde ℎ es el tamaño del paso. Se dice que es hacia adelante porque se usan los datos 𝑖
e 𝑖 + 1.
∆𝒇𝒊 = 𝑵𝒖𝒎𝒆𝒓𝒂𝒅𝒐𝒓 𝒅𝒆 𝒇[𝒙𝒊 , 𝒙𝒊+𝟏 ]
∆´´ 𝒇𝒊 = 𝑵𝒖𝒎𝒆𝒓𝒂𝒅𝒐𝒓 𝒅𝒆 𝒇[𝒙𝒊 , 𝒙𝒊+𝟏 , 𝒙𝒊+𝟐 ]
:
𝒏
∆ 𝒇𝒊 = 𝑵𝒖𝒎𝒆𝒓𝒂𝒅𝒐𝒓 𝒅𝒆 𝒇[𝒙𝒊 , 𝒙𝒊+𝟏 , 𝒙𝒊+𝟐 , 𝒙𝒊+𝟑 ]

Regla del trapecio (fórmula alternativa y ERROR): La regla del trapecio se puede
obtener integrando un polinomio de primer grado Newton-Gregory de interpolación hacia
adelante:

Forma general de un polinomio de interpolación de Newton-Gregory:

 52

∆´´ 𝑓(𝑎) ∆𝑛 𝑓(𝑎)

𝑓𝑛 (𝑎) = 𝑓(𝑎) + ∆𝑓(𝑎)𝜶 + 𝜶(𝜶 − 1) + ⋯ + 𝜶(𝜶 − 1) … (𝜶 − 𝑛 + 1) + 𝑅𝑛
2! 𝑛!
𝒙−𝒂 𝒙−𝒂
𝑐𝑜𝑛: 𝜶 = = , 𝑑𝑜𝑛𝑑𝑒 ℎ 𝑒𝑠 𝑒𝑙 𝑡𝑎𝑚𝑎ñ𝑜 𝑑𝑒𝑙 𝑝𝑎𝑠𝑜
𝒉 𝒃−𝒂

𝑓 𝑛+1 (ℰ) 𝑛+1

𝑒𝑙 𝑟𝑒𝑠𝑖𝑑𝑢𝑜: 𝑅𝑛 = ℎ 𝛼(𝛼 − 1)(𝛼 − 2) … (𝛼 − 𝑛)
(𝑛 + 1)!

Luego la integral del polinomio de primer grado de Newton-Gregory en un intervalo comprendido

entre 𝑎 y 𝑏 es:

𝑏
𝒇´´ (𝓔) 𝟐
𝐼 = ∫ [ 𝑓(𝑎) + ∆𝑓(𝑎)𝜶 + 𝒉 𝛼(𝛼 − 1)]𝑑𝑥
𝑎 𝟐!
𝑑𝑥
Realizando un cambio de variable de 𝑥 por 𝜶, 𝑑𝜶 = ℎ
→ 𝑑𝑥 = 𝒉. 𝒅𝜶, entonces queda:

1
𝒇´´ (𝓔) 𝟐
𝐼 = ℎ ∫ [ 𝑓(𝑎) + ∆𝑓(𝑎)𝜶 + 𝒉 𝛼(𝛼 − 1)]𝑑𝜶
0 𝟐!

Integrando y aplicando Barrow:

𝒇(𝒂) + 𝒇(𝒃) 𝒇´´(𝓔) 𝟑

𝑰= 𝒉− 𝒉
𝟐 𝟏𝟐

Donde:

𝒇(𝒂)+𝒇(𝒃)
 𝟐
𝒉 es la regla del trapecio
𝒇´´(𝓔) 𝟑
 − 𝒉 es el error de truncamiento verdadero
𝟏𝟐

Regla del trapecio de aplicación múltiple (Regla

compuesta): Consiste en separar un intervalo (𝑎, 𝑏) que se
quiere integrar en variaos subsegmentos y aplicar la regla del
trapecio en cada uno.

Si se tienen 𝑛 + 1 puntos igualmente espaciados se tienen 𝑛

segmentos del mismo ancho:

𝑏−𝑎
ℎ=
𝑛
La integral de la función 𝑓 en el intervalo (𝑎, 𝑏) es:
 53

𝑥1 𝑥2 𝑥𝑛
𝐼 = ∫ 𝑓(𝑥)𝑑𝑥 + ∫ 𝑓(𝑥)𝑑𝑥 + ⋯ + ∫ 𝑓(𝑥)𝑑𝑥
𝑥0 𝑥1 𝑥𝑛−1

Después de aplicar la regla del trapecio queda:

𝑓(𝑥0 ) + 𝑓(𝑥1 ) 𝑓(𝑥1 ) + 𝑓(𝑥2 ) 𝑓(𝑥𝑛−1 ) + 𝑓(𝑥𝑛 )

𝐼 ≅ ℎ( + + ⋯+ )
2 2 2
𝑛−1
ℎ
𝐼 ≅ [𝑓(𝑥0 ) + 𝑓(𝑥𝑛 ) + 2 ∑ 𝑓(𝑥𝑖 )]
2
𝑖=1

𝑓(𝑥0 ) + 𝑓(𝑥𝑛 ) + 2 ∑𝑛−1

𝑖=1 𝑓(𝑥𝑖 )
𝐼≅⏟
(𝑏 − 𝑎)
⏟ 2𝑛
𝐴𝑛𝑐ℎ𝑜
𝐴𝑙𝑡𝑢𝑟𝑎 𝑝𝑟𝑜𝑚𝑒𝑑𝑖𝑜

A los puntos interiores se le da el doble de peso que a los puntos extremos 𝑥0 y 𝑥1 .

Error de la regla del trapecio de aplicación múltiple: Error de truncamiento verdadero:

𝒏
(𝒃 − 𝒂)𝟑
𝑬𝒕 = − ∑ 𝒇´´(𝓔𝒊 )
𝟏𝟐𝒏𝟑
𝒊=𝟏

𝑐𝑜𝑛 ℰ𝑖 𝑢𝑛 𝑝𝑢𝑛𝑡𝑜 𝑙𝑜𝑐𝑎𝑙𝑖𝑧𝑎𝑑𝑜 𝑒𝑛 𝑒𝑙 𝑠𝑒𝑔𝑚𝑒𝑛𝑡𝑜 𝑖

𝒏
∑ 𝒇´´(𝓔 )
Ahora utilizando 𝒇̅´´ ≅ 𝒊=𝟏 𝑛 𝒊 , se obtiene el error de truncamiento aproximado:

(𝒃 − 𝒂)𝟑
𝑬𝒂 = − 𝒇̅´´
𝟏𝟐𝒏𝟐

Regla del trapecio de aplicación múltiple (Regla compuesta):

 54

Reglas de Simpson

Obtiene una estimación más exacta que la regla del trapecio haciendo uso de polinomios de grado
superior. Donde 𝑛 + 1 puntos se unirán con un polinomio de grado 𝑛. Las formulas resultantes de
tomar las integrales de estos polinomios se conocen como reglas de Simpson.

Regla de Simpson 1/3: El 1/3 viene de que en la fórmula de la

regla ℎ está dividido por tres. La regla se obtiene cuando se utiliza un
polinomio de segundo grado para estimar la integral:

𝑏
𝐼 ≅ ∫ 𝑓2 (𝑥)𝑑𝑥
𝑎

Designando los extremos del intervalo en el que se quiere integrar

como 𝑥0 y 𝑥2. Con un punto intermedio 𝑥1 , se puede representar𝑓2 (𝑥)
mediante un polinomio de interpolación de Lagrange. Entonces la
 55

integral queda:
𝑥2 (𝑥 − 𝑥1 )(𝑥 − 𝑥2 ) (𝑥 − 𝑥0 )(𝑥 − 𝑥2 )
𝐼≅∫ [ 𝑓(𝑥0 ) + 𝑓(𝑥1 )
𝑥0 (𝑥0 − 𝑥1 )(𝑥0 − 𝑥2 ) (𝑥1 − 𝑥0 )(𝑥1 − 𝑥2 )
(𝑥 − 𝑥0 )(𝑥 − 𝑥1 )
+ 𝑓(𝑥2 )]𝑑𝑥
(𝑥2 − 𝑥0 )(𝑥2 − 𝑥1 )

De donde se obtiene la regla de Simpson 1/3:

ℎ 𝑏−𝑎
𝐼≅ [𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )]; 𝑐𝑜𝑛 ℎ =
3 2
Otra forma de expresar la regla:

𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )

𝐼 ≅ (𝑥
⏟ 2 − 𝑥0 ) ⏟
6
𝐴𝑛𝑐ℎ𝑜
𝐴𝑙𝑡𝑢𝑟𝑎 𝑝𝑟𝑜𝑚𝑒𝑑𝑖𝑜

En esta última expresión el punto medio esta ponderado por dos tercios y los extremos están
ponderados por un sexto. La regla de Simpson 𝟏/𝟑 da resultados exactos para polinomios de
tercer grado o menores.

Error de la regla Simpson 1/3: El error se puede obtener con el mismo procedimiento que
aplicado para la regla del trapecio, utilizando un polinomio de interpolación de Newton-Gregory
hacia adelante de tercer grado con su respectivo residuo. Cambiando los límites de integración,
integrando, aplicando Barrow y reordenando se obtiene:

ℎ 1
𝐼= [𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )] − 𝑓´´(ℰ)ℎ5
⏟
3 ⏟ 90
𝑅𝑒𝑔𝑙𝑎 𝑆𝑖𝑚𝑝𝑠𝑜𝑛 1/3 𝑒𝑟𝑟𝑜𝑟 𝑑𝑒 𝑡𝑟𝑢𝑛𝑐𝑎𝑚𝑖𝑒𝑛𝑡𝑜

Entonces el error de truncamiento verdadero de una sola aplicación es:

1
𝐸𝑡 = − 𝑓´´(ℰ)ℎ5
90
Y el error de truncamiento aproximado es:

2 𝐼𝑉 𝑥
(𝑥2 − 𝑥0 )5 ∫𝑥0 𝑓 (𝑥)𝑑𝑥 (𝑥2 − 𝑥0 )5 𝐼𝑉
𝐸𝑎 = − 5 ( )=− 𝑓 ̅ (ℰ)
2 ∗ 90 𝑥2 − 𝑥0 2880

Regla de Simpson 1/3 de aplicación múltiple: Divide el intervalo de integración en

varios segmentos de un mismo tamaño:
 56

𝑏−𝑎
ℎ= 𝑛

La integral queda:
𝑥2 𝑥4
𝐼 = ∫ 𝑓(𝑥)𝑑𝑥 + ∫ 𝑓(𝑥)𝑑𝑥 + ⋯
𝑥0 𝑥2
𝑥𝑛
+∫ 𝑓(𝑥)𝑑𝑥
𝑥𝑛−2

Aplicando la regla de Simpson 1/3 queda:

2ℎ 2ℎ
𝐼≅ [𝑓(𝑥0 ) + 4𝑓(𝑥1 ) + 𝑓(𝑥2 )] + [𝑓(𝑥2 ) + 4𝑓(𝑥3 ) + 𝑓(𝑥4 )] + ⋯
3 3
2ℎ
+ [𝑓(𝑥𝑛−2 ) + 4𝑓(𝑥𝑛−1 ) + 𝑓(𝑥𝑛 )]
3
Como cada subintervalo está formado por tres puntos los cuales están igualmente espaciados, su
ancho es de 2ℎ.

Reordenando se obtiene la regla de Simpson 1/3 de aplicación múltiple:

𝑓(𝑥0 ) + 4 ∑𝑛−1 𝑛−2

𝑖=1,3,5[𝑓(𝑥𝑖 )] + 2 ∑𝑗=2,4,6[𝑓(𝑥𝑗 )] + 𝑓(𝑥𝑛 )
𝐼 ≅ (𝑏 − 𝑎)
3𝑛

Error de la regla Simpson 1/3 de aplicación múltiple: El error de truncamiento es:

(𝒙𝟐 − 𝒙𝟎 )𝟓 ̅𝑰𝑽
𝑬𝒂 = − 𝒇 (𝓔)
𝟏𝟖𝟎𝒏𝟒
𝑥 𝐼𝑉
𝐼𝑉 ∫𝑥 2 𝑓 (𝑥)𝑑𝑥
̅ (ℰ) =
𝑐𝑜𝑛 𝑓 0
𝑥𝑛 − 𝑥0

Limitaciones de la regla Simpson 1/3 de aplicación múltiple: La primera limitación es que esté
método solo es válido para cuando se tienen valores equidistantes. La segunda limitación es solo
es válida cuando se tiene números impar de segmentos y un número impar de puntos.

Pseudocodigo de la regla Simpson 1/3 de aplicación múltiple (regla compuesta)

57
 58

Regla Simpson 3/8: Consiste en ajustar un polinomio de Lagrange de tercer grado a cuatro
puntos. Da resultado exactos para polinomios de grado 𝑛 ≤ 3. Es una regla útil para cuando se
tienen un número de intervalos impar.

𝑏
𝐼 ≅ ∫ 𝑓3 (𝑥)𝑑𝑥
𝑎

𝑥3 (𝑥
− 𝑥1 )(𝑥 − 𝑥2 )(𝑥 − 𝑥3 ) (𝑥 − 𝑥0 )(𝑥 − 𝑥2 )(𝑥 − 𝑥3 )
𝐼≅∫ [ 𝑓(𝑥0 ) + 𝑓(𝑥1 )
𝑥0 (𝑥0 − 𝑥1 )(𝑥0 − 𝑥2 )(𝑥0 − 𝑥3 ) (𝑥1 − 𝑥0 )(𝑥1 − 𝑥2 )(𝑥1 − 𝑥3 )
(𝑥 − 𝑥0 )(𝑥 − 𝑥1 )(𝑥 − 𝑥3 ) (𝑥 − 𝑥0 )(𝑥 − 𝑥1 )(𝑥 − 𝑥2 )
+ 𝑓(𝑥2 ) + 𝑓(𝑥2 )]𝑑𝑥
(𝑥2 − 𝑥0 )(𝑥2 − 𝑥1 )(𝑥2 − 𝑥3 ) (𝑥3 − 𝑥0 )(𝑥3 − 𝑥1 )(𝑥3 − 𝑥2 )

Para obtener:

𝟑𝒉 (𝑏 − 𝑎)
𝑰≅ [𝒇(𝒙𝟎 ) + 𝟑𝒇(𝒙𝟏 ) + 𝟑𝒇(𝒙𝟐 ) + 𝒇(𝒙𝟑 )], 𝑐𝑜𝑛 ℎ =
𝟖 3
Otra forma de expresar la ecuación como:

𝒇(𝒙𝟎 ) + 𝟑𝒇(𝒙𝟏 ) + 𝟑𝒇(𝒙𝟐 ) + 𝒇(𝒙𝟑 )

𝑰≅⏟
(𝒃 − 𝒂)
⏟ 𝟖
𝒂𝒏𝒄𝒉𝒐
𝒂𝒍𝒕𝒖𝒓𝒂 𝒑𝒓𝒐𝒎𝒆𝒅𝒊𝒐

Los puntos interiores tienen un peso de tres octavos mientras que los puntos exteriores tienen un
peso de un octavo.

Error de la Regla Simpson 3/8: El error verdadero es:

3 5 𝐼𝑉 (𝑏 − 𝑎)5 𝐼𝑉
𝐸𝑡 = − ℎ 𝑓 (ℰ) = − 𝑓 (ℰ)
80 6480

Regla Simpson 3/8 de aplicación múltiple: Se divide el intervalo a integrar [𝑎, 𝑏] en 𝑛

(𝑏−𝑎)
segmentos iguales ℎ = 𝑛
, con 𝑛 múltiplo de tres. Para luego aplicar la regla de Simpson 3/8 en
cada uno de ellos obteniendo:

3ℎ 3ℎ
𝐼≅ [𝑓(𝑥0 ) + 3𝑓(𝑥1 ) + 3𝑓(𝑥2 ) + 𝑓(𝑥3 )] + [𝑓(𝑥3 ) + 3𝑓(𝑥4 ) + 3𝑓(𝑥5 ) + 𝑓(𝑥6 )] + ⋯
8 8
3ℎ
+ [𝑓(𝑥𝑛−3 ) + 3𝑓(𝑥𝑛−2 ) + 3𝑓(𝑥𝑛−1 ) + 𝑓(𝑥𝑛 )],
8

Reordenando queda:

𝒏−𝟏 𝒏−𝟑
𝟑𝒉
𝑰≅ [𝒇(𝒙𝟎 ) + 𝟑 ∑ 𝒇(𝒙𝒊 ) + 𝟐 ∑ 𝒇(𝒙𝒋 ) + 𝒇(𝒙𝒏 )]
𝟖
𝒊=𝟏,𝟐,𝟒,𝟓 𝒋=𝟑,𝟔,𝟗

Reemplazando ℎ:
 59

𝒇(𝒙𝟎 ) + 𝟑 ∑𝒏−𝟏 𝒏−𝟑

𝒊=𝟏,𝟐,𝟒,𝟓 𝒇(𝒙𝒊 ) + 𝟐 ∑𝒋=𝟑,𝟔,𝟗 𝒇(𝒙𝒋 ) + 𝒇(𝒙𝒏 )
𝑰 ≅ (𝒃 − 𝒂) [𝟑 ∗ ]
𝟖∗𝒏

Pseudocodigo regla Simpson 3/8 de aplicación múltiple:

 60

Comparación de los errores de la regla de trapecio y la regla Simpson 1/3

*Nota: Para cuando se tiene una gran cantidad de segmentos predomina el error de redondeo.

Integración de ecuaciones
En la integración de ecuaciones a diferencia de la integración a través de datos tabulados, se
conoce la función, por lo que se pueden generar tantos valores de 𝑓(𝑥) como se requieran para
llegar alcanzar una exactitud aceptable.

Extrapolación de Richardson

Este método usa dos estimaciones de una integral para calcular una tercera más exacta. La
fórmula de la extrapolación de Richardson se obtiene de la regla del trapecio.
 61

Se parte de la regla del trapecio de aplicación múltiple con su error correspondiente, expresado
en su fórmula general:

𝐼 = 𝐼(ℎ) + 𝐸𝑡 (ℎ)

Con: 𝐼(ℎ) la aproximación de la integral en un intervalo [𝑎, 𝑏] obtenida de una aplicación de la

regla del trapecio sobre 𝑛 segmentos de ancho ℎ.

Ahora si se hacen dos aproximaciones sobre un intervalo con tamaños de paso ℎ1 y ℎ2 se

obtiene:

𝐼(ℎ1 ) + 𝐸𝑡 (ℎ1 ) = 𝐼(ℎ2 ) + 𝐸𝑡 (ℎ2 ) (1)

El error aproximado de la regla del trapecio es:

(𝑏 − 𝑎) 2
𝐸𝑎 (ℎ) = − ̅
ℎ 𝑓´´; 𝑑𝑜𝑛𝑑𝑒 𝑓 ̅´´ 𝑒𝑠 𝑢𝑛𝑎 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒
12

Del cual se puede obtener la razón entre los dos errores:

𝑓 ̅´´ 𝐸𝑎 (ℎ1 ) 𝑓 ̅´´ 𝐸𝑎 (ℎ2 )

− (𝑏 − 𝑎) = 𝑦 − (𝑏 − 𝑎) =
12 ℎ 21 12 ℎ2 2

𝐸𝑎 (ℎ2 ) ℎ2 2
=( )
𝐸𝑎 (ℎ1 ) ℎ1

Como 𝐸𝑎 ≅ 𝐸𝑡 , entonces:

𝐸𝑡 (ℎ2 ) ℎ2 2
≅( )
𝐸𝑡 (ℎ1 ) ℎ1

De donde se despeja:

ℎ1 2
𝐸𝑡 (ℎ1 ) ≅ ( ) 𝐸𝑡 (ℎ2 )
ℎ2

Reemplazando en (1) y despejando 𝐸𝑡 (ℎ2 ):

ℎ1 2 𝑰(𝒉𝟐 ) − 𝑰(𝒉𝟏 )
𝐼(ℎ1 ) + ( ) 𝐸𝑡 (ℎ2 ) ≅ 𝐼(ℎ2 ) + 𝐸𝑡 (ℎ2 ) → 𝑬𝒕 (𝒉𝟐 ) ≅
ℎ2 𝒉 𝟐
𝟏 − ( 𝟏)
𝒉𝟐

Luego

𝐼(ℎ2 ) − 𝐼(ℎ1 )
𝐼 ≅ 𝐼(ℎ2 ) +
ℎ 2
1 − ( 1)
ℎ2

Para el caso en que ℎ2 = ℎ1 /2 queda:

 62

𝟒 𝟏
𝑰≅ 𝑰(𝒉𝟐 ) − 𝑰(𝒉𝟏 )
𝟑 𝟑

Esta última expresión es la fórmula del método de extrapolación de Richardson.

Error de la extrapolación de Richardson: Este método combina dos estimaciones de la integral

con la regla del trapecio de error 𝑂(ℎ2 ), para obtener una aproximación con un error 𝑶(𝒉𝟒 ).
Luego si se tienen tres estimaciones originales con la regla del trapecio basadas en la división
sucesiva del tamaño del paso. De estas aproximaciones se obtienen dos nuevas aproximaciones de
error 𝑂(ℎ4 ) con las cuales se obtiene una aproximación de 𝑶(𝒉𝟔 ) a través de la siguiente
ecuación:

𝟏𝟔 𝟏
𝑰≅ 𝑰𝒎 − 𝑰
𝟏𝟓 𝟏𝟓 𝒍

Con 𝐼𝑚 estimación de mayor valor y 𝐼𝑙 la estimación de menor valor.

De manera similar con dos estimaciones de error 𝑂(ℎ6 ) se puede obtener una nueva estimación
de error 𝑶(𝒉𝟖 ), utilizando:

𝟔𝟒 𝟏
𝑰≅ 𝑰𝒎 − 𝑰
𝟔𝟑 𝟔𝟑 𝒍

Algoritmo de integración de Romberg

El algoritmo de Romberg nos da una forma general de la ecuación de una aproximación por
extrapolación de Richardson de la integral, con un error 𝑂(ℎ2∗𝑘 ):

4𝑘−1 𝐼𝑗+1,𝑘−1 − 𝐼𝑗,𝑘−1

𝐼𝑗,𝑘 =
4𝑘−1 − 1

Donde:

 𝑘 indica la aproximación, siendo 𝑘 = 0 la aproximación de la regla del

trapecio
 𝑗 distingue entre la estimación más exacta (𝑗 + 1) y menos exacta (𝑗)

Esta aplicación sistemática recibe el nombre de integración de Romberg.

 63

Error aproximado y criterio de paro: El error aproximado se obtiene mediante:

𝐼1,𝑘 − 𝐼1,𝑘−1
|ℰ𝑎 | = | | ∗ 100
𝐼1,𝑘

Se toma como criterio de paro

|ℰ𝑎 | < ℰ𝑠

Integrales impropias

En esta sección se estudia integrales con límite inferior −∞ y/o límite superior +∞.

Este tipo de integrales se resuelve con un cambio de variable, que transforma el límite infinito en
1 1 1
uno finito. El cambio de variable es 𝑥 = 𝑡 → 𝑡 = 𝑥 𝑦 𝑑𝑥 = − (𝑡 2 ) 𝑑𝑡, entonces:

𝑏 1/𝑏
1 1
∫ 𝑓(𝑥)𝑑𝑥 = ∫ − 𝑓 ( ) 𝑑𝑡
𝑎 1/𝑎 𝑡2 𝑡

Esta ecuación es válida para si uno de los dos extremos es finito y el otro infinito; y 𝑎, 𝑏 > 0. Es
decir 𝑎, 𝑏 deben ser del mismo signo. En el caso que 𝑎, 𝑏 no sean del mismo signo, como 𝑎 = −∞
y 𝑏 > 0, se separa la integral.
𝑏 −𝐴 𝑏
∫ 𝑓(𝑥)𝑑𝑥 = ∫ 𝑓(𝑥)𝑑𝑥 + ∫ 𝑓(𝑥)𝑑𝑥
𝑎 −∞ −𝐴

Donde la integral impropia se resolverá con el cambio de variable anterior, utilizando una fórmula
abierta o semiabierta y la integral definida con alguna fórmula cerrada de Newton Cotes:

Formulas abiertas: Evita evaluar la integral en los límites donde puede haber un punto singulars

Formula semiabiertas: Son usadas cuando uno u otro extremo del intervalo es cerrado, la
preferida:
𝑥𝑛
∫ 𝑓(𝑥)𝑑𝑥 = ℎ[𝑓(𝑥1/2 ) + 𝑓(𝑥3/2 ) + ⋯ + 𝑓(𝑥𝑛−3/2 ) + 𝑓(𝑥𝑛−1/2 )]
𝑥0

Esta fórmula se conoce como regla extendida del punto medio. Se basa en límites de integración
que están ℎ/2 después y antes del último dato respectivamente.

Diferenciación numérica
Se obtendrán aproximaciones de las derivadas por diferencias divididas de alta exactitud
utilizando los términos de la serie de Taylor.
 64

Fórmulas de diferenciación de alta precisión en diferencias divididas

Se parte de la serie de Taylor hacia adelante, con ℎ = (𝑥𝑖+1 − 𝑥𝑖 )

𝑓´´(𝑥𝑖 ) 2
𝑓(𝑥𝑖+1 ) = 𝑓(𝑥𝑖 ) + 𝑓´(𝑥𝑖 )(ℎ) + ℎ +⋯ (1)
2!

De donde se despeja:

𝑓(𝑥𝑖+1 ) − 𝑓(𝑥𝑖 ) 𝑓´´(𝑥𝑖 )

𝑓´(𝑥𝑖 ) = − ℎ + 𝑂(ℎ2 ) (2)
ℎ 2!

Una estimación de la segunda derivada se obtiene usando la serie de Taylor:

𝑓´´(𝑥𝑖 )
𝑓(𝑥𝑖+2 ) = 𝑓(𝑥𝑖 ) + 𝑓´(𝑥𝑖 )(2ℎ) + (2ℎ)2 + ⋯ ; 𝑐𝑜𝑛 2ℎ = (𝑥𝑖+2 − 𝑥𝑖 )
2!

Restándole a esta expresión la ecuación (1) multiplicada por 2 y despejando 𝑓´´(𝑥):

𝒇(𝒙𝒊+𝟐 ) − 𝟐𝒇(𝒙𝒊+𝟏 ) + 𝒇(𝒙𝒊 )

𝑓(𝑥𝑖+2 ) − 2𝑓(𝑥𝑖+1 ) = −𝑓(𝑥𝑖 ) + 𝑓´´(𝑥𝑖 )ℎ2 + ⋯ → 𝒇´´(𝒙𝒊 ) = + 𝑶(𝒉)
𝒉𝟐

Reemplazando esta última expresión en (2):

𝑓(𝑥𝑖+1 ) − 𝑓(𝑥𝑖 ) 𝑓(𝑥𝑖+2 ) − 2𝑓(𝑥𝑖+1 ) + 𝑓(𝑥𝑖 )

𝑓´(𝑥𝑖 ) = −[ + 𝑂(ℎ)]ℎ + 𝑂(ℎ2 )
ℎ 2! ℎ2

*Nota: 𝑂(ℎ) ∗ ℎ = 𝑂(ℎ2 )

Reordenando queda:

−𝒇(𝒙𝒊+𝟐 ) + 𝟒𝒇(𝒙𝒊+𝟏 ) − 𝟑𝒇(𝒙𝒊 )

𝒇´(𝒙𝒊 ) = + 𝑶(𝒉𝟐 )
𝟐𝒉

*Nota: los datos 𝑥𝑖+2 , 𝑥𝑖+1 y 𝑥𝑖 se encuentran igualmente espaciados.

Esta es una aproximación de la primera derivada por diferencias divididas hacia adelante. Las
aproximaciones por diferencias divididas hacia atrás y centrada:
 65

Extrapolación de Richardson

La extrapolación de Richardson usa dos estimaciones de la derivada para obtener una tercera más
exacta.

La extrapolación de Richardson tiene una fórmula similar a la del caso de las integrales, para las
integrales era:

4 1
𝐼 ≅ 𝐼(ℎ2 ) − 𝐼(ℎ1 )
3 3

Para la derivada:

𝟒 𝟏
𝑫≅ 𝑫(𝒉𝟐 ) − 𝑫(𝒉𝟏 )
𝟑 𝟑

Para aproximaciones por diferencias centradas de error 𝑂(ℎ2 ), la aplicación de esta fórmula dará
una nueva estimación de error 𝑂(ℎ4 ).

El caso general de la extrapolación de Richardson, que de dos estimaciones de error 𝑂(ℎ𝑛 )

obtiene una nueva aproximación de 𝑂(ℎ𝑛+2 ) es:

𝑫(𝒉𝟐 ) − 𝑫(𝒉𝟏 )
𝑫 ≅ 𝑫(𝒉𝟐 ) +
𝒉 𝒏
( 𝟏) − 𝟏
𝒉𝟐

Derivadas de datos irregularmente espaciados

Una manera de obtener derivadas con datos irregularmente espaciados consiste en primero
ajustar un polinomio de interpolación de Lagrange de segundo grado a cada conjunto de tres
puntos adyacentes:

(𝑥−𝑥𝑖 )(𝑥−𝑥𝑖+1 ) (𝑥−𝑥 )(𝑥−𝑥 ) (𝑥−𝑥𝑖−1 )(𝑥−𝑥𝑖 )

𝑓(𝑥) = (𝑥 𝑓(𝑥𝑖−1 ) + (𝑥 −𝑥𝑖−1 )(𝑥 −𝑥𝑖+1 ) 𝑓(𝑥𝑖 ) + (𝑥 𝑓(𝑥𝑖+1 )
𝑖−1 −𝑥𝑖 )(𝑥𝑖−1 −𝑥𝑖+1 ) 𝑖 𝑖−1 𝑖 𝑖+1 𝑖+1 −𝑥𝑖−1 )(𝑥𝑖+1 −𝑥𝑖 )

Y luego se deriva analíticamente para obtener:

𝟐𝒙 − 𝒙𝒊+𝟏 − 𝒙𝒊 𝟐𝒙 − 𝒙𝒊+𝟏 − 𝒙𝒊−𝟏

𝒇´(𝒙) = 𝒇(𝒙𝒊−𝟏 ) + 𝒇(𝒙𝒊 )
(𝒙𝒊−𝟏 − 𝒙𝒊 )(𝒙𝒊−𝟏 − 𝒙𝒊+𝟏 ) (𝒙𝒊 − 𝒙𝒊−𝟏 )(𝒙𝒊 − 𝒙𝒊+𝟏 )
𝟐𝒙 − 𝒙𝒊 − 𝒙𝒊−𝟏
+ 𝒇(𝒙𝒊 )
(𝒙𝒊+𝟏 − 𝒙𝒊−𝟏 )(𝒙𝒊+𝟏 − 𝒙𝒊 )

Esta ecuación nos permite estimar la derivada en cualquier punto dentro de un intervalo
determinado por tres puntos y su error es de 𝑶(𝒉𝟐 ).
 66

Unidad 7: Ecuaciones diferenciales ordinarias

Las ecuaciones diferenciales ordinarias expresan la razón de cambio de una variable como una
función de variables y parámetros. Por ejemplo la razón de cambio de la posición con respecto al
tiempo es:

𝑑𝑣 𝑐
=𝑔− 𝑣
𝑑𝑡 𝑚

Donde 𝑔, 𝑐 y 𝑚 son parámetros y 𝑣 es una variable.

Las ecuaciones de orden superior pueden reducirse a un sistema de ecuaciones de menor

orden. Supóngase la ecuación que describe la posición 𝑥 de un sistema masa-resorte:

𝑑2 𝑥 𝑑𝑥
𝑚 2
+𝑐 + 𝑘𝑥 = 0
𝑑𝑡 𝑑𝑡

Se define una nueva variable:

𝒅𝒙
𝒚= (𝟏)
𝒅𝒕

Reemplazando y despejando:

𝑑𝑦 𝒅𝒚 𝒄𝒚 + 𝒌𝒙
𝑚 + 𝑐𝑦 + 𝑘𝑥 = 0 → =− (𝟐)
𝑑𝑡 𝒅𝒕 𝒎

Donde (1) y (2) forman un sistema de ecuaciones diferenciales de primer orden, equivalente a la
ecuación original de segundo orden.

Métodos para resolver EDO sin computadora

Si bien la algunas EDO se puede resolver de forma analítica, hay muchas EDO para las cuales no
están disponibles soluciones exactas.

Método de la Linealización: Como la mayoría de las ecuaciones no lineales no pueden resolverse

de manera exacta, una táctica para resolverlas era linealizarlas. Una ecuación diferencial ordinaria
lineal tiene forma general:

𝑑𝑛 𝑦
𝑎𝑛 (𝑥)𝑦 𝑛 + ⋯ + 𝑎1 (𝑥)𝑦´ + 𝑎0 𝑦 = 𝑓(𝑥), 𝑐𝑜𝑛 𝑦 𝑛 =
𝑑𝑥 𝑛

Ej: Linealización de la ecuación que describe el movimiento de un péndulo oscilante:

𝑑2 𝜃 𝑔
( 2 ) + 𝑠𝑒𝑛(𝜃) = 0
𝑑𝑡 𝑙

Como para ángulos pequeños 𝑠𝑒𝑛(𝜃) ≅ 𝜃 se obtiene:

 67

𝑑2 𝜃 𝑔
( 2 )+ 𝜃 = 0
𝑑𝑡 𝑙

Métodos Runge-Kutta
Los métodos Runge-Kutta se dedican a dar soluciones a las EDO de la forma:

𝑑𝑦
= 𝑓(𝑥, 𝑦)
𝑑𝑥

Un método para resolver esta ecuación es:

𝑦𝑖+1 = 𝑦𝑖 + ∅ℎ

Donde ∅ es una pendiente estimada. Esta ecuación

predice un nuevo valor de "𝒚𝒊+𝟏 " usando la pendiente,
para extrapolar linealmente sobre el tamaño de paso 𝒉.
Aplicando esta ecuación paso a paso se obtienen valores
posteriores, que sirven para trazar una trayectoria de la
solución.

Todos los métodos Runge-Kutta se basan en esta última

expresión y en lo que se diferencian es en la forma en la
que estiman la pendiente.

Método de Euler

Este método toma pendiente al inicio del intervalo como

una aproximación promedio de la pendiente sobre todo el
intervalo.

La primera derivada ofrece una estimación directa de la

pendiente en 𝑥𝑖 :

∅ = 𝑓(𝑥𝑖 , 𝑦𝑖 )

Donde el punto (𝑥𝑖 , 𝑦𝑖 ) es el punto al inicio del intervalo.

Luego la solución de la ecuación diferencial queda:

𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ∗ ℎ
 68

Esta ecuación predice un nuevo valor de "𝑦" usando la pendiente, para extrapolar linealmente
sobre el tamaño de paso ℎ. La estimación obtenida del método de Euler es exacta si la función que
se quiere estimar el lineal, de ahí que se diga que es un método de primer orden.

Error método de Euler: Analizamos primero el error de truncamiento local, por serie de Taylor:

𝑓´(𝑥𝑖 , 𝑦𝑖 ) 2 𝑓 𝑛 (𝑥𝑖 , 𝑦𝑖 ) 𝑛
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )ℎ + ℎ + ⋯+ ℎ + 𝑂(ℎ𝑛+1 )
2! 𝑛!

Comparando con la fórmula de Euler

𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ∗ ℎ

El error de truncamiento verdadero es:

𝑓´(𝑥𝑖 , 𝑦𝑖 ) 2
𝐸𝑡 = ℎ + ⋯ + 𝑂(ℎ𝑛+1 )
2!

Para un tamaño del paso lo suficientemente pequeño los términos de mayor grado pierden peso y
se obtiene:

𝐸𝑎 = 𝑂(ℎ2 )

Se puede demostrar que el error de truncamiento global es 𝑶(𝒉).

El error de redondeo disminuye con el tamaño del paso.

Método de Heun

El método de Heun es una mejora del método de Euler.

Mejora la estimación en un paso, calculando el promedio de la
pendiente al comienzo y al final del intervalo. Para lograr esto
utiliza una ecuación predictora para estimar el valor de la
función al final del intervalo.

𝒚𝟎𝒊+𝟏 = 𝒚𝒊 + 𝒇(𝒙𝒊 , 𝒚𝒊 ) ∗ 𝒉, 𝑐𝑜𝑛 𝑦´𝑖 = 𝑓(𝑥𝑖 , 𝑦𝑖 )

Luego utiliza este valor para calcular la pendiente al final del

intervalo:
0
𝑦´𝑖+1 = 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )

Para poder sacar el promedio:

0
𝑦´𝑖 + 𝑦´𝑖+1 𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦̅´ = =
2 2

Finalmente este promedio se utiliza como pendiente para

 69

extrapolar y calcular 𝑦𝑖+1.

𝒇(𝒙𝒊 ,𝒚𝒊 )+𝒇(𝒙𝒊+𝟏 ,𝒚𝟎𝒊+𝟏 )

𝒚𝒊+𝟏 = 𝒚𝒊 + ∗𝒉
𝟐

Esta ecuación recibe el nombre de ecuación correctora. El

método de Heun es un método predictor-corrector, esto se
refiere a que una estimación anterior se utilizara de forma
repetida para calcular una mejor estimación. Se puede aplicar la
ecuación predictora y correctora de forma iterativa para mejorar
la estimación, el error asociado a cada iteración se puede
calcular como:

𝑗 𝑗−1
𝑦𝑖+1 − 𝑦𝑖+1
|ℰ𝑡 | = | 𝑗
| ∗ 100
𝑦𝑖+1

*Nota: j (iteración actual del corrector), y j-1 (iteración anterior del corrector).

Este es un método de segundo orden, por lo tanto es exacto para funciones cuadráticas.

Error del método de Heun: Hay una estrecha relación entre la regla del trapecio y el método de
Heun, la regla del trapecio se define como:
𝑥𝑖+1
𝑓(𝑥𝑖 ) + 𝑓(𝑥𝑖+1 )
∫ 𝑓(𝑥)𝑑𝑥 ≅ ℎ
𝑥𝑖 2

La ecuación correctora de Heun se puede expresar como:

𝑦𝑖+1 𝑥𝑖+1
∫ 𝑑𝑦 = ∫ 𝑓(𝑥)𝑑𝑥
𝑦𝑖 𝑥𝑖

𝑥𝑖+1
𝑦𝑖+1 = 𝑦𝑖 + ∫ 𝑓(𝑥)𝑑𝑥
𝑥𝑖

En conclusión el error truncamiento del método de Heun está dado por la regla del trapecio:

𝑓´´(ℰ) 3
𝐸𝑡 = − ℎ
12

El método tiene un error local 𝑶(𝒉𝟑 ) y un error global 𝑶(𝒉𝟐 ).

Método de Punto Medio

Esta técnica utiliza el método de Euler para calcular un valor de y en el punto medio del intervalo.
 70

ℎ
𝑦𝑖+1/2 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )
2

Este punto luego es utilizado para calcular la pendiente en el punto medio:

,
𝑦𝑖+1/2 = 𝑓 (𝑥 1, 𝑦 1)
𝑖+ 𝑖+
2 2

Está pendiente la utiliza el método como una aproximación a la pendiente promedio de todo el
intervalo. Luego está pendiente se utiliza para extrapolar linealmente de 𝑥𝑖 hasta 𝑥𝑖+1 :

𝑦𝑖+1 = 𝑦𝑖 + 𝑓 (𝑥 1, 𝑦 1) . ℎ
𝑖+ 𝑖+
2 2

Se puede observar que al no estar el término 𝑦𝑖+1 en cada miembro de la ecuación no se puede
aplicarse de manera iterativa esta ecuación para mejorar la aproximación.

Error del método del punto medio: Al ser la pendiente una aproximación centrada tiene un error
de truncamiento local de 𝑂(ℎ2 ) y en consecuencia los errores local y global del método son 𝑂(ℎ2 )
y 𝑂(ℎ3 ), respectivamente.

Métodos Runge-Kutta de cuarto orden

Hay muchos métodos Runge-Kutta de cuarto orden el que se muestra a continuación es el método
clásico RK de cuarto orden:

𝟏
𝒚𝒊+𝟏 = 𝒚𝒊 + (𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒 )𝒉
𝟔

Donde:

- 𝑘1 = 𝑓(𝑥𝑖 , 𝑦𝑖 )
1 1
- 𝑘2 = 𝑓(𝑥𝑖 + 2 ℎ, 𝑦𝑖 + 2 𝑘1 ℎ)
1 1
- 𝑘3 = 𝑓(𝑥𝑖 + 2 ℎ, 𝑦𝑖 + 2 𝑘2 ℎ)
- 𝑘4 = 𝑓(𝑥𝑖 + ℎ, 𝑦𝑖 + 𝑘3 ℎ)

Este método utiliza múltiples estimaciones de la pendiente para obtener una mejor aproximación
de la pendiente promedio en el intervalo. La anterior ecuación representa una recta cuya
pendiente es el promedio ponderado de las pendientes 𝑘1 , 𝑘2 , 𝑘3 , 𝑘4.

Error del método RK de cuarto orden: El error de truncamiento local del método es 𝑂(ℎ5 ) y el
global 𝑂(ℎ6 ).
 71

Sistemas de ecuaciones

Los sistemas de EDO se presentan como:

𝑑𝑦1 (𝑥)
= 𝑓1 (𝑥, 𝑦1 , 𝑦2 , … , 𝑦𝑛 )
𝑑𝑥
𝑑𝑦2 (𝑥)
= 𝑓2 (𝑥, 𝑦1 , 𝑦2 , … , 𝑦𝑛 )
𝑑𝑥

𝑑𝑦𝑛 (𝑥)
= 𝑓𝑛 (𝑥, 𝑦1 , 𝑦2 , … , 𝑦𝑛 )
𝑑𝑥

La resolución de sistemas de ecuaciones consiste en aplicar la técnica elegida por ecuación antes
de proceder con la siguiente. Para el caso de que la técnica elegida sea el método clásico RK de
cuarto orden este se resuelve de la siguiente manera:

Supongamos que el sistema es de dos ecuaciones y un intervalo (𝑖, 𝑖 + 1)

𝑑𝑦1 (𝑥) 𝑑𝑦2 (𝑥)

= 𝑓1 (𝑥, 𝑦1 , 𝑦2 ) 𝑦 = 𝑓2 (𝑥, 𝑦1 , 𝑦2 )
𝑑𝑥 𝑑𝑥

Primero encontramos las pendientes al inicio del intervalo, 𝑖. (Se conoce 𝑥𝑖 , 𝑦1 𝑖 , 𝑦2 𝑖 )

𝒌𝟏𝟏 = 𝒇𝟏 (𝒙𝒊 , 𝒚𝟏 𝒊 , 𝒚𝟐 𝒊 ) 𝒚 𝒌𝟏𝟐 = 𝒇𝟐 (𝒙𝒊 , 𝒚𝟏 𝒊 , 𝒚𝟐 𝒊 )

Los cuales se utilizan para calcular los primeros valores de 𝑦1 y 𝑦2 en el punto medio del intervalo:

ℎ ℎ
𝑦1 𝑖+1/2 = 𝑦1 𝑖 + 𝑘11 𝑦 𝑦2 𝑖+1/2 = 𝑦2 𝑖 + 𝑘12
2 2

Y así calcular el primer conjunto de pendientes en el punto medio:

ℎ ℎ
𝑘21 = 𝑓1 (𝑥𝑖 + , 𝑦1 1 , 𝑦2 1 ) 𝑦 𝑘22 = 𝑓2 (𝑥𝑖 + , 𝑦1 1 , 𝑦2 1 )
2 𝑖+
2
𝑖+
2 2 𝑖+
2
𝑖+
2

Los cuales se utilizan para obtener una mejor aproximación de los valores de 𝑦1 y 𝑦2 en el punto
medio del intervalo:

ℎ ℎ
𝑦1 𝑖+1/2 = 𝑦1 𝑖 + 𝑘21 𝑦 𝑦2 𝑖+1/2 = 𝑦2 𝑖 + 𝑘22
2 2

Estos se usan para calcular el segundo conjunto de pendientes en el punto medio:

ℎ ℎ
𝑘31 = 𝑓1 (𝑥𝑖 + , 𝑦1 1 , 𝑦2 1 ) 𝑦 𝑘32 = 𝑓2 (𝑥𝑖 + , 𝑦1 1 , 𝑦2 1 )
2 𝑖+
2
𝑖+
2 2 𝑖+
2
𝑖+
2

Con este último conjunto calculamos 𝑦1 y 𝑦2 al final del intervalo, 𝑖 + 1:

 72

𝑦1 𝑖+1 = 𝑦1 𝑖 + 𝑘31 ℎ 𝑦 𝑦1 𝑖+1 = 𝑦1 𝑖 + 𝑘31 ℎ

Que se utiliza para calcular el último conjunto de pendientes al final del intervalo, 𝑖 + 1:

𝑘41 = 𝑓1 (𝑥𝑖 + ℎ, 𝑦1 𝑖+1 , 𝑦2 𝑖+1 ) 𝑦 𝑘42 = 𝑓2 (𝑥𝑖 + ℎ, 𝑦1 𝑖+1 , 𝑦2 𝑖+1 )

Con todo esto finalmente podemos calcular:

𝟏
𝒚𝟏 𝒊+𝟏 = 𝒚𝟏 (𝒙𝒊 + 𝒉) = 𝒚𝟏 (𝒙𝒊 ) + (𝒌𝟏𝟏 + 𝟐𝒌𝟐𝟏 + 𝟐𝒌𝟑𝟏 + 𝒌𝟒𝟏 )𝒉
𝟔
𝟏
𝒚𝟐 𝒊+𝟏 = 𝒚𝟐 (𝒙𝒊 + 𝒉) = 𝒚𝟐 (𝒙𝒊 ) + (𝒌𝟏𝟐 + 𝟐𝒌𝟐𝟐 + 𝟐𝒌𝟑𝟐 + 𝒌𝟒𝟐 )𝒉
𝟔

Caso General

𝑑5 𝑦
= 𝑓(𝑥, 𝑦 𝐼 , 𝑦 𝐼𝐼 , 𝑦 𝐼𝐼𝐼 , 𝑦 𝐼𝑉 )
𝑑𝑥 5

Defino:

𝒛𝟏 ≜ 𝒚𝑰𝑽

Entonces:

𝒅𝒛𝟐
𝒛𝟐 ≜ 𝒚𝑰𝑰𝑰 → 𝒛𝟏 =
𝒅𝒙
𝒅𝒛𝟑
𝒛𝟑 ≜ 𝒚𝑰𝑰 → 𝒛𝟐 =
𝒅𝒙
𝒅𝒛𝟒
𝒛𝟒 ≜ 𝒚𝑰 → 𝒛𝟑 =
𝒅𝒙
𝒅𝒛𝟓
𝒛𝟓 ≜ 𝒚 → 𝒛𝟒 =
𝒅𝒙

Luego:

𝑑𝑧1
= 𝑓1 (𝑥, 𝑧5 , 𝑧4 , 𝑧3 , 𝑧2 , 𝑧1 )
𝑑𝑥
𝑑𝑧2
= 𝑧1 = 𝑓2 (𝑥, 𝑧5 , 𝑧4 , 𝑧3 , 𝑧2 , 𝑧1 )
𝑑𝑥
𝑑𝑧3
= 𝑧2 = 𝑓3 (𝑥, 𝑧5 , 𝑧4 , 𝑧3 , 𝑧2 , 𝑧1 )
𝑑𝑥
𝑑𝑧4
= 𝑧3 = 𝑓4 (𝑥, 𝑧5 , 𝑧4 , 𝑧3 , 𝑧2 , 𝑧1 )
𝑑𝑥
 73

𝑑𝑧5
= 𝑧4 = 𝑓5 (𝑥, 𝑧5 , 𝑧4 , 𝑧3 , 𝑧2 , 𝑧1 )
𝑑𝑥