CURSO

DE

METODOS

NUMERICOS

Pedro Ochoa Moreno Enero de 2002.

Pedro Ochoa Moreno Enero de 2002.

 CONTENIDO
1. Introducción 1
2 Raíces de ecuaciones 3
Método de Tanteos 4
Series de Taylor 5
Método de Bisección 7
Método de Newton-Raphson 11
Método de la Secante 13
3. Sistemas de Ecuaciones Algebraicas Lineales 15
Métodos Directos 15
Eliminación de Gauss 15
Gauss-Jordán 19
Métodos Iterativos 24
Gauss-Seidel 24
4. Ajuste de Curvas 29
Interpolación 30
Diferencia Dividida para la Interpolación de Polinomios 30
Interpolación Lineal 30
Interpolación Cuadrática 32
Formula General de la Interpolación de Polinomios de
Newton 35
Interpolación de Polinomios de Lagrange 37
Regresión Lineal 40
Criterios para un Mejor Ajuste 41
Cuantificación del Error de una Regresión Lineal 44
5. Integración 47
Regla Trapezoidal 49
Reglas de Simpson 57
6. Ecuaciones Diferenciales Ordinarias 67
Fundamento matemático 68
Método de Euler 71
Método Mejorado del Polígono (Euler Mejorado) 76
Métodos de Runge-Kutta 77

Pedro Ochoa Moreno Enero de 2002.

Bibliografía:

1. Introducción a la Computación para Ingenieros. Steven C.

Chapra y Raymond P. Canale. Mc Graw-Hill, 1986.
2. Métodos Numéricos para Ingenieros. Steven C, Chapra y
Raymond P. Canale. Mc Graw-Hill, 1998. 3ra. Edición.
3. Métodos Numéricos Aplicados a la Ingeniería. Antonio Nieves y
Federico C. Domínguez. CECSA, 1998. 4ta. Edición.

Pedro Ochoa Moreno Enero de 2002.

 Curso de Métodos Numéricos

1. INTRODUCCION

El análisis numérico y sus métodos numéricos son una dialéctica entre el

análisis numérico cualitativo y el análisis numérico cuantitativo; el primero nos
dice por ejemplo que bajo ciertas condiciones algo existe, que es o no único,
etc. Mientras que el segundo complementa al primero, permitiendo calcular
aproximadamente el valor de aquello que existe.

Los métodos numéricos son científicos en el sentido que representan técnicas

sistemáticas para resolver problemas matemáticos. Sin embargo, hay cierto
grado de arte, juicios subjetivos y una concordancia, asociada con su uso
efectivo en la práctica de ingeniería. Para cada uno de los problemas, la
confrontación es con varios métodos numéricos alternativos y con muchos tipos
de computadoras. Por lo tanto, la elegancia y la eficiencia de los diferentes
enfoques de los problemas es muy individualista y se correlaciona con la
habilidad de escoger en forma prudente entre todas las opciones.
Desafortunadamente, como sucede con cualquier proceso intuitivo, los factores
que influyen en esta elección son difíciles de comunicar.

MODELO

FISICO MATEMATICO

ESTATICO DINAMICO ESTATICO DINAMICO

NUMERICO ANALITICO NUMERICO

SIMULACION

ANALOGICA DIGITAL

Pedro Ochoa Moreno Enero, 2002. 1

 Curso de Métodos Numéricos

Razones por las cuales se deben estudiar los métodos numéricos:

1. Los métodos numéricos son herramientas muy poderosas para la

solución de problemas. Son capaces de manejar sistemas de ecuaciones
grandes, no linealidades y geometrías complicadas, comunes en la
práctica de la ingeniería y a menudo, imposibles de resolver
analíticamente. Por lo tanto aumentan la habilidad de quien los estudia
para resolver problemas.
2. Hay muchos problemas que no pueden plantearse al emplear programas
“hechos”. Si esta versado en los métodos numéricos y es un adepto de la
programación de computadoras, entonces tiene la capacidad de diseñar
sus propios programas para resolver los problemas sin tener que
comprar un software costoso.
3. Los métodos numéricos son un vehículo para eficiente para aprender a
servirse de las computadoras. Es bien sabido que un camino efectivo
para aprender programación actualmente es escribir programas de
computadora, porque la mayoría de los métodos numéricos son
diseñados para implementarlos en las computadoras, por lo tanto son
ideales para este propósito.
4. Los métodos numéricos son un medio para reforzar su comprensión de
las matemáticas, ya que una de sus funciones es convertir las
matemáticas superiores a operaciones aritméticas básicas, porque
profundizan en los temas que de otro modo resultarían oscuros. Esta
alternativa aumenta la capacidad de comprensión y entendimiento en la
materia.

Definiciones de error

Error verdadero Et =Valor verdadero − aproximación

Error relativo porcentual Valor verdadero − aproximación

εt = 100%
verdadero Valor verdadero

Error relativo porcentual Aproximación actual − aproximación anterior

εa = 100%
aproximado aproximación actual

Criterio de paro Terminar los cálculos cuando

εa < εs
donde ε s es el error relativo porcentual deseado

Pedro Ochoa Moreno Enero, 2002. 2

 Curso de Métodos Numéricos

2. Raíces de ecuaciones

Los problemas que se verán en esta unidad están relacionados con el valor de
una variable o de un parámetro que satisface una ecuación. Son especialmente
valiosos en proyectos de ingeniería, donde con frecuencia resulta imposible
despejar de manera analítica los parámetros de las ecuaciones de diseño.

Resolver f ( x ) = 0 para x.

Si pretendemos resolver la siguiente ecuación:

f ( x ) = ax 2 + bx + c = 0 1

Recordemos que se nos enseño a usar la formula cuadrática:

− b ± b 2 − 4ac
x= 2
2a

Los valores calculados por la ecuación anterior se les llaman “raíces” de la

ecuación f (x ) . Estas raíces representan los valores de x que harán que la
ecuación f ( x ) = 0.
Como pudimos observar la ecuación cuadrática es útil para resolver la ecuación
f (x ) , pero hay que recordar que existen muchas funciones diferentes que no
se pueden resolver de una manera tan fácil. En estos casos, los métodos
numéricos, nos proporcionaran métodos eficientes para obtener la respuesta.

Pedro Ochoa Moreno Enero, 2002. 3

 Curso de Métodos Numéricos

Método de Tanteos

Este método determina los valores de f (x ) correspondientes a valores

sucesivos de x hasta que se presenta un cambio de signo en f (x ) . El cambio
de signo indica que se ha pasado por una raíz. Entonces se puede obtener una
mayor aproximación al valor de la raíz, volviendo al último valor de x que
precede el cambio de signo ya partir de este valor de x , determinar
nuevamente valores de f (x ) correspondientes a valores sucesivos de x
utilizando ahora un incremento menor que el que se uso inicialmente, hasta
que cambie de nuevo de signo f (x ) . Este procedimiento se repite con
incrementos de x cada vez más pequeños hasta lograr un valor
suficientemente preciso de la raíz. Si se desean raíces adicionales, el
incremento de x puede continuar hasta que se localice la siguiente raíz en
forma aproximada mediante otro cambio de signo de f (x ) , y así
sucesivamente.

Ejemplo:
Obténgase una raíz aproximada de la función: f ( x ) = e − x − x

Solución:

1ra. Aproximación 2da. Aproximación 3ra. Aproximación

x f (x ) x f (x ) x f (x )
0.0 1.0000 0.50 0.1065 0.57 -0.0044
0.2 0.6187 0.52 0.0745 0.565 0.0033
0.4 0.2703 0.54 0.0427 0.567 0.0002
0.6 -0.0512 0.56 0.0112 0.568 -0.0013
0.8 -0.3506 0.58 -0.0201
1.0 -0.6321 0.60 -0.0511

El código utilizado en el entorno MatLab se presenta a continuación:

%Demostración del Método gráfico para la

%función: y=e^a(-x)-x.
%Determine la raíz aproximada de la función
%(Cuando y=0.002, entonces x=0.567)

clf; clear;
x=0:0.001:1;
y(1,:)=exp(-x)-(x)
y(2,:)=0.*(x);
plot(x,y);
axis([0 1 -1 1])
xlabel('x'); ylabel('y=exp(-x)-(x)')
title('Método Gráfico')

Pedro Ochoa Moreno Enero, 2002. 4

 Curso de Métodos Numéricos

La respuesta de este código se presenta en la siguiente gráfica:

Método Gráfico
1

0.8

0.6

0.4

0.2
y=exp(-x)-(x)

-0.2

-0.4

-0.6

-0.8

-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Series de Taylor

Las series de Taylor nos proporcionan una formulación para predecir el valor de
la función en xi +1 en términos de la función y de sus derivadas en una vecindad
al punto x i .
En lugar de presentar en conjunto la serie de Taylor, se obtendrá una mayor
comprensión si se construye término a término. Por ejemplo, el primer término
de la serie es:

f ( x i+1 ) ≅ f ( x i ) 3

Esta igualdad conocida como aproximación de orden cero, indica que el valor
de f en el nuevo punto es el mismo que el valor en el punto anterior. Este
resultado se logra intuitivamente ya que si x i y xi +1 están muy próximos uno
del otro, entonces es igualmente posible que el nuevo valor sea probablemente
similar al anterior.
La ecuación 3 da una estimación perfecta, si la función que se va aproximar es
una constante. Sin embargo, si la función cambia en todo el intervalo, entonces
se requieren términos adicionales de la serie de Taylor para obtener una mejor
aproximación. Por ejemplo, la aproximación a primer orden se obtiene
sumando otro término al anterior para obtener:

Pedro Ochoa Moreno Enero, 2002. 5

 Curso de Métodos Numéricos

f ( xi +1 ) ≅ f ( x i ) + f ' ( xi )( xi +1 − x i ) 4

El termino adicional de primer orden consiste de la pendiente

f ' ( xi ) multiplicada por la distancia entre x i y xi +1 .Por lo tanto, la expresión
ahora representa una línea recta y es capaz de predecir un incremento o
decremento de la función entre x i y xi +1 .
Aunque la ecuación 4 puede predecir un cambio, solo es exacta para una línea
recta o es de dirección lineal. Por lo tanto, se le agrega a la serie un término de
segundo orden para obtener algo más sobre la curvatura de la función si es que
la tiene:
f '' ( x i )
f ( x i+1 ) ≅ f ( x i ) + f ' ( xi )( x i+1 − xi ) + ( xi +1 − x i ) 2 5
2!

De manera similar, se pueden agregar términos adicionales para desarrollar la

expansión completa de la serie de Taylor:

f '' ( x i )
f ( x i+1 ) = f ( xi ) + f ' ( xi )( xi +1 − xi ) + ( xi +1 − x i ) 2
2!
6
f "' ( x i ) f (n ) ( xi )
+ ( x i+1 − xi ) + L +
3
( x i+1 − xi ) n + R
3! n!

Debido a que la ecuación 6 es una serie infinita, el signo igual remplaza el de

aproximación usada en las ecuaciones 3, 4 y 5. Se incluye un término adicional
llamado término residual para considerar todos los términos desde n + 1 hasta
el infinito:

f (n +1) (ξ )
Rn = ( x i+1 ) n+1 7
( n + 1)!

Donde ξ es un valor de x que se encuentra entre x i y xi +1 . Frecuentemente es

conveniente simplificar la serie de Taylor, definiendo un paso h = x i+1 − xi y
expresando la ecuación 6 como se muestra a continuación:

f ' ' ( xi ) 2
f ( x i+1 ) = f ( xi ) + f ' ( xi ) h + h
2!
8
f "' ( x i ) 3 f ( n ) ( xi ) n
+ h +L + h + Rn
3! n!

Entonces el término residual es ahora:

f ( n+1) (ξ ) n+1
Rn = h 9
( n + 1)!

Pedro Ochoa Moreno Enero, 2002. 6

 Curso de Métodos Numéricos

Ejemplo: Usense términos en la serie de Taylor de cero a cuarto orden para

aproximar la función:

f ( x ) = − 0.1x 4 − 0.15 x 3 − 0.5 x 2 − 0.25 x +1.2

Desde x i = 0 con h =1 . Esto es, predecir el valor de la función en x i+1 =1 .

Aproximación grafica de la serie de Taylor para la función

f ( x ) = − 0.1x 4 − 0.15 x 3 − 0.5 x 2 − 0.25 x +1.2 en x =1 .

Método de Bisección

Conocido también como el de corte binario, de partición en dos intervalos

iguales o método Bolzano, es un método de búsqueda incremental en el que el
intervalo se divide siempre en dos. Si la función cambia de signo sobre un
intervalo, se evalúa el valor de la función en el punto medio. 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. El proceso se repite hasta obtener una mejor
aproximación. Esto lo podemos visualizar mediante los pasos del siguiente
algoritmo:

Pedro Ochoa Moreno Enero, 2002. 7

 Curso de Métodos Numéricos

1. Elija los valores iniciales inferior x l y superior x u de forma tal que

la función cambie de signo sobre el intervalo. Esto se puede
verificar asegurándose de que f ( x l ) f ( x u ) < 0 .

2. La primera aproximación a la raíz x r se determina como sigue:

xl + xu
xr =
2

3. Realice las siguientes evaluaciones para determinar en qué

subintervalo cae la raíz.

a) Si f ( x l ) f ( x u ) < 0 , entonces la raíz se encuentra dentro del

subintervalo. Por lo tanto, tome x u = xr y continué con el
paso 2.
b) Si f ( x l ) f ( x u ) > 0 , entonces la raíz se encuentra dentro del
subintervalo superior. Por lo tanto, resuelvase x l = x r , y
continúese con el paso 2.

4. Si f ( x l ) f ( x u ) = 0 , la raíz es igual a x r ; termina el calculo.

En la siguiente figura se muestra un bosquejo gráfico de este método:

Pedro Ochoa Moreno Enero, 2002. 8

 Curso de Métodos Numéricos

Una ilustración de las formas que puede tener una raíz en un intervalo prescrito
por los límites inferior x l y superior x u , se presentan a continuación:

a) Si los x l y x u tienen el mismo

signo, entonces no habrá raíces.

b) Si los intervalos tienen signos

opuestos, por lo cual habrá un
número impar de raíces.

c) Si los intervalos tienen signos

iguales, entonces existen un par
de raíces.

d) Si los intervalos tienen signos

opuestos en los extremos,
entonces habrá un número
impar de raíces dentro del
intervalo.

El siguiente paso es el de desarrollar un criterio objetivo para decidir cuando

debe terminar el método. Como sugerencia inicial podemos decidir la
finalización del cálculo cuando el error se encuentre por debajo de algún nivel
prefijado. Puede decidirse que el método termine cuando se alcance un error
mas bajo, por ejemplo del 0.1%. Esta estrategia es valida siempre y cuando no
se conozca el valor verdadero de la raíz de la función. Por lo tanto, se requiere
estimar el error de manera que no incluya el conocimiento previo de la raíz. Se
puede calcular el error relativo aproximado ε a de la siguiente manera:

x rnuevo − xranterior
εa = 100% 10
x rnuevo

Donde x rnuevo es la raíz de la iteración actual y x ranteriores la raíz de la iteración

anterior. Se usa el valor absoluto, ya que por lo general importa solo la
magnitud de ε a sin considerar su signo. Cuando ε a es menor que un valor

Pedro Ochoa Moreno Enero, 2002. 9

 Curso de Métodos Numéricos

previamente fijado ε s , el calculo termina. Lo anterior se puede expresar como

el criterio de paro:

xrnuevo − xranterior
100% ≤ ε s 11
xrnuevo

Método de falsa posición

Un defecto del método de bisección es que al dividir el intervalo de x l a x u en

mitades iguales, no se toma en cuenta la magnitud de f ( xl ) y f ( x u ) . Por
ejemplo, si f ( xl ) es mucho más cercana a cero que f ( x u ) , es lógico que la raíz
se encuentra más cerca de x l que de x u , como se observa en la siguiente
figura:

Pedro Ochoa Moreno Enero, 2002. 10

 Curso de Métodos Numéricos

Este método alternativo que aprovecha la visualización gráfica es la de unir

f ( xl ) y f ( x u ) con una línea recta. La intersección de esta línea con el eje de las
x representa una mejor estimación de la raíz. El hecho de que se remplace la
curva por una línea recta da una “posición falsa” de la raíz; de aquí el nombre
de método de la falsa posición, o en latín regula falsi. También se le conoce
como el método de interpolación lineal.
Usando triángulos semejantes (como se observan en la figura anterior), la
intersección de la línea recta con el eje de las x puede ser estimado como:

f (x l ) f (x u )
= 12
xr − xl x r − xu

La ecuación 12 puede resolverse de la siguiente manera:

f ( xu )( x l − xu )
xr = xu − 13
f ( xl ) − f ( xu )

Esta es la formula de la falsa posición. El valor de x r , calculado con la ecuación

13 remplaza uno de los valores iniciales x l o x u , que produzca un valor de la
función que tenga el mismo signo de f (x ) . De esta manera, los valores de x l
y x u siempre encierran a la raíz. El proceso se repite hasta que la aproximación
de la raíz sea la adecuada. El algoritmo es idéntico al de la bisección, excepto
que la ecuación 13 se usa en el paso 2. Además, se usan los mismos criterios
de para terminar los cálculos.

Desventajas del método de la falsa posición

Aunque el método de la falsa posición parecería ser siempre la mejor opción de

los que usan intervalos, hay casos donde funciona deficientemente.

Método de Newton-Raphson

Tal vez dentro de las formulas para localizar raíces, la formula de Newton-
Raphson sea la más usada. Si el valor inicial de la raíz es x i , entonces se
puede extender una tangente desde el punto [xi , f ( x ) ] . El punto donde esta
tangente cruza el eje x representa una aproximación mejorada de la raíz.
El método de Newton-Raphson se puede derivar geométricamente (una forma
de hacerlo es mediante el uso de la serie de Taylor).

A continuación se presenta un esquema gráfico del método de Newton-

Raphson:

Pedro Ochoa Moreno Enero, 2002. 11

 Curso de Métodos Numéricos

Como podemos observar en la figura anterior se extrapola una tangente a la

función en el punto x i (esto es, f ' ( xi ) ) hasta el eje x para obtener una
estimación de la raíz en xi +1 . Como en la figura, la primera derivada de x es
equivalente a la pendiente, o sea

f (x i ) − 0
f ' ( xi ) = 14
xi − x i+1

Que se puede reordenar par obtener:

f (x i )
x i+1 = xi − 15
f ' (x i )

La cual es conocida como fórmula de Newton-Raphson.

Como con los otros métodos de localización de raíces, la siguiente ecuación se

puede usar como criterio de paro.

xi +1 − x i
100% ≤ ε s 16
xi+1

Además, el desarrollo del método con base en la serie de Taylor proporciona

un conocimiento teórico relacionado con la velocidad de convergencia
expresado como: Ei+1 = O ( Ei2 ). De esta forma, el error debe ser casi proporcional
al cuadrado del error anterior. En otras palabras, el número de cifras
significativas aproximadamente se duplica en cada iteración.

Pedro Ochoa Moreno Enero, 2002. 12

 Curso de Métodos Numéricos

El algoritmo del método de Newton-Raphson se obtiene al sustituir la ecuación

15 bajo el formato de la formula predictiva x i+1 = g ( xi ) .

Desventajas del método de Newton-Raphson

Aunque el método de Newton-Raphson en general es muy eficiente, hay

situaciones en que se comporta en forma eficiente. Un caso especial – raíces
múltiples.

Método de la Secante

Un problema potencial en la implementación del método de Newton-Raphson

es el de la evaluación de la derivada. Aunque esto no es un inconveniente para
los polinomios y muchas otras funciones, existen algunas funciones cuyas
derivadas pueden ser en extremo difíciles de evaluar. En estos casos, la
derivada se puede aproximar mediante una diferencia dividida finita regresiva.
Esto se muestra en el esquema grafico del método de la secante que se da a
continuación:

El método de la secante usa una diferencia más que una derivada para estimar
la pendiente.

f ( x i−1 ) − f ( xi )
f ' ( xi ) ≅ 17
xi−1 − xi

Sí sustituimos la ecuación 16 en la ecuación 15, podemos obtener la siguiente

ecuación iterativa:

Pedro Ochoa Moreno Enero, 2002. 13

 Curso de Métodos Numéricos

f ( x i )( xi−1 − xi )
x i+1 = xi − 18
f ( x i−1 ) − f ( xi )

La ecuación 17 representa la fórmula para el método de la secante. Observe

que el planteamiento requiere de dos puntos iniciales de x . Sin embargo,
debido a que no se requiere que f (x ) cambie de signo entre estos valores, este
método no es clasificado como aquellos que usan intervalos.

El criterio de paro para el método de la secante es idéntico al usado en el

método de Newton-Raphson y que se da por la ecuación 16.

Diferencia entre los métodos de la secante y de la falsa posición

Observe la similitud entre los métodos de la secante y de la falsa posición. Por

ejemplo, las ecuaciones 17 y 13 son idénticas en todos sus términos. Ambas
usan dos estimaciones iniciales para calcular una aproximación de la pendiente
de la función que se usa para proyectar hacia el eje x una nueva aproximación
a la raíz. Sin embargo, existe una diferencia crítica entre ambos métodos. Tal
diferencia estriba en la forma en que uno de los valores es reemplazado por la
nueva aproximación. Recuerdese que en el método de la falsa posición, la
ultima aproximación de la raíz reemplaza cualquiera de los valores dando una
función con el mismo signo como f ( x r ) . En consecuencia las dos
aproximaciones siempre encierran a la raíz. Por lo tanto, para todos los casos
prácticos, el método siempre converge, ya que la raíz se encuentra dentro del
intervalo. En contraste el método de la secante reemplaza los valores en una
secuencia estricta, con el nuevo valor xi +1 reemplaza a x i y x i reemplaza a
xi −1 . Como resultado de esto, los dos valores pueden caer en un mismo lado de
la raíz. En algunos casos esto puede provocar divergencia.

Pedro Ochoa Moreno Enero, 2002. 14

 Curso de Métodos Numéricos

3. Sistemas de ecuaciones algebraicas lineales

En esencia, estos problemas son similares a los de raíces de ecuaciones en el
sentido de que están relacionados con valores que satisfacen ecuaciones. Sin
embargo, a diferencia de satisfacer una sola ecuación, se busca un conjunto de
valores que satisfaga simultáneamente un conjunto de ecuaciones algebraicas.
Las ecuaciones lineales simultáneas surgen en el contexto de una variedad de
problemas y en todas las disciplinas de la ingeniería. En particular, se originan a
partir de modelos matemáticos de sistemas grandes interrelacionados, como:
estructuras, circuitos eléctricos y redes de flujo; aunque también pueden
encontrarse en otras áreas de los métodos numéricos como el ajuste de curvas
y ecuaciones diferenciales.

Dadas las a y las c resolver

a11 x1 + a12 x 2 = c1
a 21 x1 + a 22 x 2 = c 2

Para las x.

Métodos Directos

Eliminación de Gauss

El procedimiento para llevar a cabo este método consta de dos pasos:

a) Se manejan las ecuaciones para eliminar una incógnita de las

ecuaciones. El resultado de este paso de eliminación es una sola
ecuación con una incógnita.

Pedro Ochoa Moreno Enero, 2002. 15

 Curso de Métodos Numéricos

b) Por consiguiente, esta ecuación se puede resolver directamente y el

resultado se sustituye hacia atrás en las ecuaciones originales

Este comportamiento básico se puede extender a sistemas grandes de

ecuaciones desarrollando un esquema sistemático para eliminar incógnitas y
sustituir hacia atrás. La eliminación Gaussiana es una de las técnicas más
comunes. En esta técnica se utilizan dos fases que son: eliminación hacia
delante y sustitución hacia atrás (inversa).
El procedimiento está diseñado para resolver un conjunto de n ecuaciones:

a11 x1 + a12 x 2 + a13 x3 + L + a1n x n = c1

a 21 x1 + a 22 x 2 + a 23 x3 + L + a 2n x n = c2
19
M M
a n1 x1 + a n2 x 2 + an 3 x 3 + L + a nn xn = c n

El primer paso permite reducir el sistema dado por la ecuación 19, a un

sistema triangular superior. EL paso inicial será el de eliminar la primera
incógnita x1 , desde la segunda hasta la n-ésima ecuación. Este paso se lleva
acabo mediante la normalización de la primera ecuación y posteriormente
multiplicar esta por el coeficiente de la segunda ecuación y el resultado
restárselo a la segunda ecuación para poder eliminar x1 de esta ecuación. El
segundo paso se repite hasta la n-ésima ecuación. A todo este proceso se le
denomina eliminación hacia delante. Posterior a este se lleva a cabo la
sustitución hacia atrás, como se muestra a continuación:

 a11 a12 a13 c1 

a a 22 a 23 c 2 
 21
a 31 a 32 a 33 c3 

Eliminación hacia
adelante

a11 a12 a13 c1 

 '
a 22 '
a 23 c 2' 

 a ''33 c3'' 

x 3= c3'' a33
''

Sustitución hacia
x 2 = (c k' − a 23
'
x3 ) a '22
atrás
x 3 = ( c1 − a12 x 2 − a13 x3 ) a11

Pedro Ochoa Moreno Enero, 2002. 16

 Curso de Métodos Numéricos

Desventajas de los métodos de eliminación

Mientras que hay muchos sistemas de ecuaciones que se pueden resolver con
la eliminación de Gauss simple, existen algunas desventajas que se deben
analizar antes de escribir un programa de computo donde se implemente el
método.

División entre cero

Por ejemplo:

2 x 2 + 3x 3 = 8
4 x1 + 6 x 2 + 7 x 3 = − 3
2 x1 + x 2 + 6 x3 = 5

Para resolver la normalización de la primera fila involucra una división entre

a11 = 0 .

Errores por redondeo

El número de cifras significativas incide en el incremento de errores por

redondeo. El problema de los errores de redondeo puede ser particularmente
importante cuando se trata de resolver un gran número de ecuaciones. Por lo
cual se recomienda usar fracciones para minimizar el problema.

Sistemas mal condicionados

Los sistemas mal condicionados son aquellos en donde pequeños cambios en

los coeficientes generan grandes cambios en la solución. Debido a que los
errores de redondeo pueden inducir pequeños cambios en los coeficientes,
estos cambios artificiales pueden generar grandes errores en su solución para
sistemas mal condicionados.

Ejemplo:

Resuelvase el siguiente sistemas de ecuaciones mediante el método de

eliminación de Gauss.

4 x1 − 2 x 2 − 3 x3 = − 23 1
− 4 x1 + 3x 2 + x3 = 11 2
8x1 − 5 x2 + 4 x3 = 6 3

Representando el sistema en forma matricial, tenemos:

Pedro Ochoa Moreno Enero, 2002. 17

 Curso de Métodos Numéricos

 4 − 2 − 3 − 23
− 4 3 1 11 

 8 − 5 4 6 

1 1 3 23
R1 → R1 →1 − − − → R1
4 2 4 4
R2 → R2 + 4 R1 →4 −2 −3 −23
−4 3 1 11

0 1 −2 −12 → R2

R3 → R3 − 8 R1 → −8 4 6 46
8 −5 4 6

0 −1 10 52 → R3

El nuevo sistema queda de la siguiente manera:

 1 3 23 
 1 − − −
2 4 4
0 1 −2 − 12 
 
0 − 1 10 52 
 

R3 → R3 + R2 →0 1 −2 −12
0 −1 10 52

0 0 8 40 → R3

Entonces el nuevo sistema se reduce al esperado y este queda de la siguiente

manera:

 1 3 23 
1 − 2 −
4
−
4
0 1 −2 − 12 
 
0 0 8 40 
 

Pedro Ochoa Moreno Enero, 2002. 18

 Curso de Métodos Numéricos

Del sistema anterior, mediante la sustitución inversa podemos obtener:

x1 = − 3
x2 = − 2
x3 = 5

Gauss-Jordan

El método de Gauss-Jordan es una variación de la eliminación de Gauss. La

principal diferencia consiste en que cuando una incógnita se elimina en el
método de Gauss-Jordan, esta es eliminada de todas las otras ecuaciones en
lugar de hacerlo solo en las subsecuentes. De esta forma, el paso de
eliminación genera una matriz identidad en vez de una matriz triangular. Por
consiguiente, no es necesario emplear la sustitución hacia atrás para obtener la
solución. El esquema de desarrollo del método de Gauss-Jordan, los asteriscos
indican que el vector de términos independientes se ha modificado varias
veces.

 a11 a12 a13 c1 

a a 22 a 23 c 2 
 21
a 31 a 32 a 33 c3 

1 0 0 c1* 
 *
0 1 0 c 2 
0 0 1 c3* 
 

x1 = c1*
x2 = c *2
x3 = c *3

Pedro Ochoa Moreno Enero, 2002. 19

 Curso de Métodos Numéricos

Ejemplo:

Distribución de recursos (Ingeniería en General)

Un Ingeniero Industrial supervisa la producción de cuatro tipos de

computadoras. Se requieren cuatro clases de recursos: horas-hombre, metales,
plásticos y componentes electrónicos en la producción. En el cuadro 1 se
resumen las cantidades necesarias para cada uno de estos recursos en la
producción de cada tipo de computadora. Si se dispone diariamente de 504
horas-hombre, 1970 Kg. de metal, 970 Kg. de plástico y 601 componentes
electrónicos. ¿Cuántas computadoras de cada tipo se pueden construir por día?

Horas- Componentes
Computadora hombre Metales Kg. Plásticos Kg. unidades
/computadora /computadora /computadora /computadora
1 3 20 10 10
2 4 25 15 8
3 7 40 20 10
4 20 50 22 15
Cuadro 1. Recursos necesarios para producir cuatro tipos de computadoras.

Sean x1 , x 2 , x3 y x 4 las cantidades totales de computadoras producidas

diariamente de cada clase. Se sabe que la cantidad total de horas-hombre
disponibles diariamente es de 504. Por lo tanto, la suma de las distribuciones
de horas-hombre en la producción de cada una de las computadoras debe ser
igual o menor de 504. Por lo tanto, usando los datos del cuadro 1, tenemos:
3x1 + 4 x2 + 7 x3 + 20 x 4 ≤ 504 1

De la misma manera para metales, plásticos y componentes:

20 x1 + 25x 2 + 40 x 3 + 50 x 4 ≤ 1970 2
10 x1 + 15 x 2 + 20 x3 + 22 x 4 ≤ 970 3
10 x1 + 8 x 2 + 10 x 3 + 15 x4 ≤ 601 4

Cada una de las ecuaciones anteriores se debe satisfacer de forma simultánea

de otra manera, se acabaría uno o más recursos necesarios en la producción de
los cuatro tipos de computadoras. Si los recursos disponibles, representados
por el vector de términos independientes de las ecuaciones anteriores, se
reducen todos a cero simultáneamente, entonces se puede remplazar el signo
≤ por el de =. Entonces el sistema de ecuaciones es el siguiente:

Pedro Ochoa Moreno Enero, 2002. 20

 Curso de Métodos Numéricos

3x1 + 4 x2 + 7 x3 + 20 x 4 = 504
20 x1 + 25x 2 + 40 x 3 + 50 x 4 =1970
10 x1 +15 x 2 + 20 x3 + 22 x 4 = 970
10 x1 + 8 x 2 +10 x3 + 15x 4 = 601

Como podemos observar el sistema no es diagonalmente dominante, el método

de Gauss-Seidel puede divergir. Por lo tanto, usaremos el método de Gauss-
Jordan:

 3 4 7 20 504 
20 25 40 50 1970
 
10 15 20 22 970 
 
10 8 10 15 601 

1 1.333 2.333 6.667 168 

0 − 1.667 − 6.667 − 83.33 − 1390 
 
0 1.667 − 3.333 − 44.667 − 710 
 
0 − 5.333 − 13.33 − 51.667 − 1079 

1 0 − 3 − 60 − 944 
0 1 4 50 834 

0 0 − 10 − 128 − 2100 
 
0 0 8 215 3369 

1 0 0 − 21.6 − 314
0 1 0 − 1.2 − 6 

0 0 1 12.8 210 
 
0 0 0 112.6 1689 

1 0 0 0 10  x1 = 10
0 1 0 0 12  x2 = 12

0 0 1 0 18 x3 = 18
 
0 0 0 1 15 x3 = 15

Con base a la información anterior podemos calcular las ganancias totales. Por
ejemplo, supóngase que las ganancias correspondientes a cada computadora

Pedro Ochoa Moreno Enero, 2002. 21

 Curso de Métodos Numéricos

están dadas por p1 , p 2 , p 3 y p 4 . Entonces la ganancia total asociada con un

día de actividad (P) esta dada por

P = p1 x1 + p 2 x 2 + p 3 x 3 + p4 x 4 5

Si las ganancias correspondientes a cada una de las computadoras esta dada

por el cuadro 2.

Ganancia $
Computadora /computadora
1 1000
2 700
3 1100
4 400
Cuadro 2. Ganancia por computadora.

De donde:

P = 1000(10) + 700(12) + 1100(18) + 400(15) = 44,200

P = $ 44,200 Ganancia diaria, para los recursos especificados.

Ejemplo:

Un ingeniero requiere 4,800 m³ de arena, 5,810 m³ de grava fina, 5,690 m³ de

grava gruesa para la construcción de un proyecto. Existen tres bancos donde se
pueden obtener estos materiales. La composición de cada banco se da en el
cuadro 3. ¿Cuántos metros cúbicos se debe tomar de cada banco para cumplir
con las necesidades del ingeniero?

Arena Grava fina Grava gruesa

Banco % % %
1 52 30 18
2 20 50 30
3 25 20 55

Pedro Ochoa Moreno Enero, 2002. 22

 Curso de Métodos Numéricos

Sean x1 , x 2 y x3 la cantidad de material (metros cúbicos) que se debe tomar

de cada banco.

0.52 x1 + 0.20 x 2 + 0.25x 3 ≤ 4800 1

0.30 x1 + 0.50 x 2 + 0.20 x 3 ≤ 5810 2

0.18 x1 + 0.30 x2 + 0.55 x3 ≤ 5690 3

Expresando estas ecuaciones en forma matricial, tenemos:

0.52 0.20 0.25   x1  4800

0.30 0.50 0.20   x  = 5810
  2   
0.18 0.30 0.55   x3  5690

Usando el método de Gauss_Jordan en el calculo de la matriz inversa.

0.52 0.20 0.25 1 0 0

[A] = 0.30 0.50 0.20 0 1 0
0.18 0.30 0.55 0 0 1

1 0.3846 0.4808 1.9230 0 0 

0 0.3846 0.0558 − 0.5769 1 0 
 
0 0.2308 0.4635 − 0.3461 0 1 

1 0 0.4250 2.5 −1 0
0 1 0.1450 − 1.5 2.6 0 
 
0 0 0.4301 0 − 0.6 1 

1 0 0 2.5 − 0.407 − 0.9883

0 1 0 − 1.5 2.8022 − 0.3371
 
0 0 1 0 − 1.395 2.325 

 2.5 − 0.407 − 0.9883

[A] =  − 1.5 2.8022 − 0.3371
−1

 0 − 1.395 2.325 

Pedro Ochoa Moreno Enero, 2002. 23

 Curso de Métodos Numéricos

Multiplicando la matriz inversa por el vector independiente y obtenemos:

x1 = 4800( 2.5) + 5810( −0.407 ) + 5690( −0 .9883) = 4,011.504

x 2 = 4800( −1.5) + 5810( 2.8022) + 5690( −0.3371) = 7,162.683

x 2 = 4800(0) + 5810( −1.395) + 5690( 2.325) = 5,126.830

De donde que las variables son expresadas en m³:

x1 = 4,011.504 m 3
x2 = 7,162.683 m 3
x3 = 5,126.830 m 3

Métodos Iterativos

Gauss-Seidel

Este método es útil en la disminución de errores de redondeo en sistemas, y

esto se debe a que un método de aproximación se puede continuar hasta que
converja dentro de una tolerancia de error previamente especificada. De esta
forma, el redondeo no es un problema, ya que se controla el nivel de error
aceptable. Supóngase que se da un conjunto de n ecuaciones, como las de la
ecuación 19:

a11 x1 + a12 x 2 + a13 x3 + L + a1n x n = c1

a 21 x1 + a 22 x 2 + a 23 x3 + L + a 2n x n = c2
M M
a n1 x1 + a n2 x 2 + an 3 x 3 + L + a nn xn = c n

Si los elementos de la diagonal son diferentes de cero, la primera ecuación se

puede resolver para x1 , la segunda para x 2 , etcétera, lo que lleva a:

c1 − a12 x 2 − a13 x 3 − L − a1n x n

x1 = 20
a11

Pedro Ochoa Moreno Enero, 2002. 24

 Curso de Métodos Numéricos

c 2 − a21 x1 − a31 x3 − L − a 2n x n
x2 = 21
a 22

M M M

c n − a n1 x1 − a n2 x 2 − L − a n, n−1 xn −1
xn = 22
a nn

Ahora, se puede empezar el proceso de solución usando un valor inicial para las
x . La solución trivial puede servir de valor inicial, esto es que todas las x valen
cero. Estos ceros se pueden sustituir en la ecuación 20 que se puede usar para
calcular un nuevo valor de x1 = c1 a . Luego, se sustituye el nuevo valor de x1
11

junto con los valores previos de cero para x3 en la ecuación 21 y calcular el

nuevo valor para x 2 . Este proceso se repite en la ecuación 22 para calcular un
nuevo estimado de x n . Después se regresa a la primera ecuación y se repite
todo el procedimiento hasta que la solución converja lo suficientemente cercana
a los valores reales. La convergencia se puede verificar usando el criterio:

x ij − xij −1
ε a, i = 100% < ε s 23
xij

Para toda i en donde j y j − 1 denotan la iteración actual y anterior.

Mejoramiento en la convergencia usando relajación

La relajación representa una pequeña modificación del método de Gauss-Seidel

y esta diseñada para aumentar la convergencia. Después que cada nuevo valor
de x se calcula usando las ecuaciones 20, 21 y 22, el valor se modifica
mediante un promedio pesado de los resultados de las iteraciones anteriores y
actuales:

x inuevo = λx inuevo + (1 − λ ) xipasado 24

En donde λ es un factor de peso al cual se le asigna un valor entre 0 y 2. Si

λ =1 , (1 − λ ) es igual a cero y el resultado permanece inalterado. Sin embargo,
si a λ se le asigna un valor entre 0 -1, el resultado es un promedio pesado de
los resultados previos y actuales. A este tipo de modificación se le conoce
como sobrerelajación. Por lo general, esta opción se emplea para convertir un
sistema divergente en uno convergente. Si λ se encuentra entre 1 – 2 se
considera otro peso en el valor actual. En este caso, existe una suposición
implícita de que el nuevo valor se mueve en la dirección correcta hacia la
solución real pero con una velocidad muy lenta. Por lo que este tipo de

Pedro Ochoa Moreno Enero, 2002. 25

 Curso de Métodos Numéricos

modificación, al cual se le llama sobrerelajación, esta diseñado para acelerar la

convergencia de un sistema que ya es convergente.

Ejemplo:

Resuelvase el sistema de ecuaciones que se da a continuación, usando el

método de Gauss-Seidel con un criterio de paro de ε s = 10 % .

10 x1 − 3 x2 + 6 x 3 = 24.5 1
x1 + 8x 2 − 2 x 3 = 9 2
− 2 x1 + 4 x2 − 9 x3 = − 50 3

Despejando cada una de las variables, en forma diagonal:

24.5 + 3x 2 − 6 x3
x1 = 1’
10
9 − x1 + 2 x3
x2 = 2’
8
− 50 + 2 x1 − 4 x2
x3 = 3’
−9

A continuación iniciaremos con la primera iteración, suponiendo que x 2 y x3 en

la ecuación 1’ son cero, entonces el valor de x1 es:

x1 = 2.45

Este valor junto con el de x 3 = 0 , pueden sustituirse en la ecuación 2’ y obtener

el valor de x 2 :

x 2 = 0.81875

Los valores de x1 y x 2 se sustituyen en la ecuación 3’, para obtener x3 :

x 3 = 5.375

En la segunda iteración, se repite el proceso, pero utilizando los valores de la

primera iteración, como se observa a continuación:

x1 = − 0.529375
x 2 = + 2.5349
x 3 = 6.7998

Pedro Ochoa Moreno Enero, 2002. 26

 Curso de Métodos Numéricos

Con base a las dos iteraciones anteriores, podemos estimar el error:

− 0.529375 − 2.45
ε a,1 = 100 = 562.8099 % > ε s
− 0.529375

2.5349 − 0.18875
ε a, 2 = 100 = 67.70 % > ε s
2.5349

6.7998 − 5.375
ε a, 3 = 100 = 20.9535 > ε s
6.7998

La tercera iteración, se lleva a cabo dado que ε a, i > ε s .

x1 = − 0.86941
x 2 = + 2.9336
x 3 = 7.05258

Tomando los valores de las iteraciones dos y tres podemos estimar de nuevo
el error:

− 0.86941 − ( −0.529375)
ε a,1 = 100 = 39.11 % > ε s
− 0.86941

2.9336 − 2.5349
ε a, 2 = 100 = 13.59 % > ε s
2.9336

7.05258 − 6.7998
ε a, 3 = 100 = 3.5248 < ε s
7.05258

La cuarta iteración se efectuara, debido a que todavía no se cumple con el

error preestablecido.

x1 = − 0.901468
x 2 = 3.0008285
x 3 = 7.089583

Para estimar el error, tomaremos los valores de la tercera y cuarta iteración, de

donde obtenemos:

− 0.901468 − 0.86941)
ε a,1 = 100 = 3.5562 % < εs
− 0.901468

Pedro Ochoa Moreno Enero, 2002. 27

 Curso de Métodos Numéricos

3.0008285 − 2.9336
ε a, 2 = 100 = 2.240331 % < εs
3.0008285

7.089583 − 7.05258
ε a, 3 = 100 = 0.521934% < εs
7.89583

Como pudimos observar, el error (ε s ) cumple con el criterio de paro por lo

tanto el proceso se iteraciones se detiene. Y la manera de verificar si los
valores de las variables son las adecuadas, es necesario sustituirlas en las
ecuaciones 1,2 y 3.

Pedro Ochoa Moreno Enero, 2002. 28

 Curso de Métodos Numéricos

4. Ajuste de curvas
A menudo se presentara la oportunidad de ajustar curvas a un conjunto de
datos representados por puntos. Las técnicas que se han desarrollado para este
fin pueden dividirse en dos categorías generales: regresión e interpolación. La
primera se desarrolla cuando hay un grado significativo de error asociado a los
datos; con frecuencia los datos experimentales son de esta clase. Para estas
situaciones, la estrategia es encontrar una curva que represente la tendencia
general de los datos sin necesidad de tocar los puntos individuales. En
contraste, la interpolación se maneja cuando el objetivo es determinar valores
intermedios entre datos que estén, relativamente libres de error. Tal es el caso
de la información tabulada. Para estas situaciones, la estrategia es ajustar una
curva directamente mediante los puntos y usar esta curva para predecir valores
intermedios.

Pedro Ochoa Moreno Enero, 2002. 29

 Curso de Métodos Numéricos

Interpolación

Ha menudo hemos tenido la oportunidad de estimar valores intermedios entre

datos precisos. El método más común que se usa para este propósito es la
interpolación del polinomio. Recuerde que la fórmula general para un polinomio
de n-ésimo orden es

f ( x ) = a 0 + a1 x + a 2 x 2 + L + a n xn 24

Para n+1 puntos, hay uno y sólo un polinomio de orden n que pasa a través de
todos los puntos. Por ejemplo, hay sólo una línea recta (es decir, un polinomio
de primer orden) que conecta dos puntos (ver figura 1.a). De manera similar,
únicamente una parábola conecta un conjunto de tres puntos (ver figura 1.b).
Interpolación polinomial consiste en determinar el único polinomio de n-ésimo
orden que ajusta n+1 puntos. Este polinomio entonces proporcionar una
fórmula para calcular valores intermedios. En la figura 1.c se observa una curva
(polinomio de tercer orden) que conecta tres puntos.

Fig. 1. Ejemplos de interpolación polinomial.

Aunque hay uno y solo un polinomio de n-ésimo que ajusta n+1 puntos, existe
una variedad de formatos matemáticos en los cuales este polinomio puede
expresarse.

Diferencia Dividida de Newton para la interpolación de

Polinomios

Interpolación Lineal

El modo más simple de interpolación es conectar dos puntos con una línea
recta. Esta técnica, llamada Interpolación lineal, se ilustra en la figura 2
mediante triángulos semejantes,

f 1 ( x ) − ( x0 ) f ( x1 ) − f ( x 0 )
=
x − x0 x1 − x 0

La cual se puede reordenar para dar

Pedro Ochoa Moreno Enero, 2002. 30

 Curso de Métodos Numéricos

Figura 2. Ilustración gráfica de interpolación lineal.

f ( x1 ) − f ( x 0 )
f 1 ( x ) = f ( x0 ) + ( x − x0 ) 25
x1 − x 0

Que es una fórmula de interpolación lineal. La notación f 1 ( x ) designa que es

una interpolación de polinomios de primer orden. Observe que además de
representar la pendiente de la línea que conecta los puntos, el término
[ f ( x1 ) − f ( x 0 )] /( x1 − x0 ) es una aproximación por diferencia dividida finita de la
primera derivada. En general, cuanta más pequeña sea el intervalo de datos,
mejor será la aproximación. Esto se debe al hecho de que, en tanto el intervalo
disminuya, una función continua se aproximará mejor por una línea recta. Esta
característica se demuestra en el siguiente ejemplo:

Ejemplo:

Estime el logaritmo natural de 2 mediante interpolación lineal. Primero, realice

el calculo por interpolación entre ln 1 = 0 y ln 6 = 1.781759. Después, repita
el procedimiento, pero use un intervalo más pequeño de ln 1 a ln 4 (1.386294).
Observe que el valor real de ln 2 es 0.6931472.

Para resolver este problema utilizaremos la ecuación 25 y una interpolación

lineal para ln(2) de x 0 =1 a x 0 = 6 para dar

1 .791759 − 0
f 1 (2) = 0 + ( 2 −1) = 0. 3583519
6 −1
Que representa un error de ε t = 48.3 % . Con el intervalo más pequeño de, x 0 =1
a x1 = 4 se obtiene

Pedro Ochoa Moreno Enero, 2002. 31

 Curso de Métodos Numéricos

1.386294 − 0
f 1 (2) = 0 + ( 2 −1) = 0.4620981
4 −1

Así, con el intervalo más corto se reduce el error relativo porcentual a

ε t = 33.3 % . Ambas interpolaciones se muestran en la figura3 junto con la
función real.

Figura 3.Dos interpolaciones lineales para estimar ln 2.

Interpolación cuadrática

El error en el ejemplo anterior resulta de nuestra aproximación a una curva con

una línea recta. Por consiguiente, una estrategia para mejorar la estimación es
introducir alguna curvatura en la línea que conecta los puntos. Si tres puntos de
los datos están disponibles, esto puede realizarse con un polinomio de segundo
orden (también conocido como polinomio cuadrático o parabólico). Una forma
en particular conveniente para este propósito es

f 2 ( x) = b0 + b1 ( x − x0 ) + b2 ( x − x 0 )( x − x1 ) 26

Observe que aunque la ecuación 26 parecería diferir del polinomio general

(véase la ecuación 24), las dos ecuaciones son equivalentes. Esto puede
demostrarse al multiplicar los términos de la ecuación 26 para dar

f 2 ( x) = b0 + b1 x − b x 0 + b2 x 2 − b2 x0 x1 − b2 xx0 − b2 x x1 )

O agrupando términos

f 2 ( x) = a0 + a1 x +a 2 x 2

Donde

Pedro Ochoa Moreno Enero, 2002. 32

 Curso de Métodos Numéricos

a 0 = b0 − b1 x0 + b2 x0 x1
a1 = b1 − b2 x 0 − b2 x1
a 2 = b2

Así, las ecuaciones 24 y 26 son formulaciones alternativas equivalentes del

único polinomio de segundo orden que une los tres puntos.
Un procedimiento simple puede usarse para determinar los valores de los
coeficientes. Para b0 , la ecuación 26 con x = x 0 puede ser usada para calcular

b0 = f ( x 0 ) 27

La ecuación 27 puede sustituirse en la ecuación 26, la cual puede evaluarse en

x = x1 para

f ( x1 ) − f ( x 0 )
b1 = 28
x1 − x 0

Por último, las ecuaciones 27 y 28 se pueden sustituir en la 26, la cual puede

evaluarse en x = x 2 y resolver (después de algunas manipulaciones algebraicas)
para

f ( x 2 ) − f ( x1 ) f ( x1 ) − f ( x 0 )
−
x2 − x1 x1− x0
b2 = 29
x 2 − x0

Note que, como fue el caso con la interpolación lineal, b1 todavía representa la
pendiente de la línea que une los puntos x 0 y x1 . Así, los primeros dos términos
de la ecuación 26 son equivalentes a la interpolación lineal de x 0 a x1 , como se
especifico antes en la ecuación 25. El último término, b2 ( x − x0 )( x − x1 ) ,
introduce la curvatura de segundo orden en la formula.
Antes de ilustrar cómo usar la ecuación 26, debemos examinar la forma del
coeficiente b2 . Es muy similar a la aproximación por diferencias divididas finitas
de la segunda derivada. Así, la ecuación 26 comienza a manifestar una
estructura que es muy similar a la serie de expansión de Taylor. Esta
observación será objeto de más exploración en la sección de errores al
interpolar polinomios de Newton. Pero primero, desarrollaremos un ejemplo
que muestra cómo se usa la ecuación 26 para interpolar entre tres puntos.

Ejemplo:

Ajuste los tres puntos usados en el ejemplo anterior, a un polinomio de

segundo orden:

Pedro Ochoa Moreno Enero, 2002. 33

 Curso de Métodos Numéricos

x 0 =1 f (x 0 ) = 0
x1 = 4 f ( x1 ) =1.386294
x2 = 6 f ( x 2 ) =1.791759

Use el polinomio para evaluar ln 2.

Aplicando la ecuación 27 se obtiene

b0 = 0

La ecuación 28 da
1.386294 − 0
b1 = = 0.4620981
4 −1
Y con la ecuación 29 se obtiene

1.791759 −1.386294
− 0.4620981
6−4
b2 = = − 0.0518731
6 −1
Sustituyendo estos valores en la ecuación 26 se obtiene la formula cuadrática

f 2 ( x) = 0 + 0.4620981( x −1) − 0.0518731 ( x − 1)( x − 4)

Que puede evaluarse en x = 2 para

f 2 ( 2) = 0.568444

La cual representa un error relativo de ε t = 18. 4 % . Así, la curvatura introducida

por la fórmula cuadrática (véase figura 4) mejora la interpolación comparada
con el resultado que se obtiene mediante líneas rectas en el ejemplo anterior.

Figura 4. Uso de la interpolación cuadrática para estimar ln 2.

Pedro Ochoa Moreno Enero, 2002. 34

 Curso de Métodos Numéricos

Formula general de la interpolación de polinomios de Newton

El análisis anterior puede ser generalizado para ajustar un polinomio de n-

ésimo orden a n+1 datos. El polinomio de n-ésimo orden es

f n ( x ) = b0 + b1 ( x − x0 ) +L + bn ( x − x0 )( x − x1 ) L ( x − x n−1 ) 30

Como se hizo antes con las interpolaciones lineales y cuadráticas, los puntos de
los datos evaluaban los coeficientes b0 , b1 ,L , bn . Para un polinomio de n-ésimo
orden se requiere n+1 puntos: [x 0 , f ( x0 ) ], [x1 , f ( x1 ) ], L, [x n , f ( x n ) ] . Usamos estos
datos y las siguientes ecuaciones para evaluar los coeficientes:

b0 = f ( x0 ) 31
b1 = f [x1 , x0 ] 32
33
b2 = f [x 2 , x1 , x0 ]
M
bn = f [x n , xn −1 , L, x1 , x 0 ] 34

Donde las ecuaciones de la función puestas entre paréntesis son diferencias

divididas finitas. Por ejemplo, la primera diferencia dividida finita se representa
por lo general como

f ( xi ) − f ( x j )
[ ]
f xi , x j =
xi − x j
35

La segunda diferencia dividida finita, la cual representa la diferencia de las dos

primeras diferencias divididas, se expresa por lo general como

[ ] [
f xi , x j − f x j , xk ]
[ ]
f x i , x j , xk = 36
xi − x k

En forma similar, la n-ésima diferencia dividida finita es

f [x n , x n−1 , L, x1 ]− f [x n−1 , x n−2 , L, x0 ]

f [x n , xn −1 , L, x1 , x 0 ] = 37
x n − x0

Estas diferencias pueden usarse para evaluar los coeficientes en las ecuaciones
31 hasta la 34, las cuales entonces se sustituirán en la ecuación 30 para
obtener el polinomio de interpolación

f n ( x ) = f ( x0 ) + ( x − x 0 ) f [ x1 , x 0 ]+ ( x − x 0 )( x − x1 ) f [x 2 , x1 , x0 ]
38
+ L + ( x − x0 )( x − x1 ) L ( x − xn −1 ) f [x n , xn −1 , L, x0 ]

Pedro Ochoa Moreno Enero, 2002. 35

 Curso de Métodos Numéricos

Que es conocido como polinomio de interpolación por diferencias divididas de

Newton. Debe observarse que no es necesario que los datos utilizados en la
ecuación 38 sean igualmente espaciados o que los valores de la abscisa deban
estar en orden ascendente, como se ilustra en el siguiente ejemplo. También,
observe cómo las ecuaciones 35 a 37 son recursivas (es decir, las diferencias
de órdenes mayores se calculan al tomar diferencias de orden menor,
obsérvese el cuadro que se da a continuación.

i xi f ( xi ) Primero Segundo Tercero

0 x0 f (x 0 ) f ( x1 , x0 ) f ( x 2 , x1 , x 0 ) f ( x 3 , x2 , x1 , x 0 )
1 x1 f ( x1 ) f ( x 2 , x1 ) f ( x 3 , x 2 , x1 )
2 x2 x
f (x 2 ) f ( x 3 , x2 )
i

3 x3 f (x 3 )

Cuadro 1. Ilustración gráfica de la naturaleza de las diferencias divididas.

Ejemplo:

En el ejemplo anterior, los datos en x 0 =1 , x1 = 4 y x 2 = 6 se utilizaron para

estimar ln 2 con una parábola. Ahora, agregando un cuarto punto
( x 3 = 5 ; f ( x3 ) =1.609438 ), calcule el ln 2 con una interpolación del polinomio de
Newton de tercer orden.

El polinomio de tercer orden utilizando la ecuación 30 con n=3, es

f 3 ( x ) = b0 + b1 ( x − x 0 ) + b2 ( x − x 0 )( x − x1 ) + b3 ( x − x0 )( x − x1 )( x − x 2 )

Las primeras diferencias divididas son

1.386294 − 0
f [x1 , x0 ] = = 0.4620981
4− 0

1.791759 −1.386294
f [ x2 , x1 ] = = 0.2027326
6− 4

1.609438 −1.386294
f [x 3 , x 2 ] = = 0.0.1823216
5−6
Las segundas diferencias divididas son

0.2027326 − 0.4620981
f [x 2 , x1 , x 0 ] = = − 0.05187311
6 −1

Pedro Ochoa Moreno Enero, 2002. 36

 Curso de Métodos Numéricos

0.1823216 − 0.2027326
f [x 3 , x 2 , x1 ] = = − 0.02041100
5− 4

La tercera diferencia dividida es

− 0.02041100 − (−0.05187311)
f [x 3 , x 2 , x1 , x0 ] = = 0.007865529
5 −1

Los resultados para f [x1 , x0 ] , f [x 2 , x1 , x0 ] y f [ x3 , x 2 , x1 , x0 ] representan los

coeficientes b1 , b2 y b3 de la ecuación 30. Junto con b0 = f ( x0 ) = 0.0 , la
ecuación 30 es

f 3 ( x ) = 0 + 0.4620981 ( x − 1) − 0.05187311( x − 1)( x − 4) + 0.007865529 ( x − 1)( x − 4)( x − 6)

Que puede usarse para evaluar f 3 ( 2) = 0.6287686, el cual representa un error

relativo de ε t = 9.3 %. El polinomio cúbico completo se muestra en la figura 5.

Figura 5. Uso de la interpolación cúbica para estimar ln 2.

Interpolación de Polinomios de Lagrange

La interpolación de polinomios de Lagrange es simplemente una reformulación

del polinomio de Newton que evita el cálculo por diferencias divididas. Se puede
expresar de manera concisa como

Pedro Ochoa Moreno Enero, 2002. 37

 Curso de Métodos Numéricos

n
f n (x ) = ∑ L (x ) f (x )
i i 39
i= 0

Donde

n x − xj
Li ( x ) = ∏ 40
j= 0 xi − x j
j ≠1

Donde ∏ designa el “producto de”. Por ejemplo, la versión lineal (n=1) es

x − x1 x − x0
f1 (x ) = f (x 0 ) + f ( x1 ) 41
x 0 − x1 x1 − x 0

Y la versión de segundo orden es

( x − x1 )( x − x2 ) ( x − x0 )( x − x 2 )
f 2 ( x) = f (x 0 ) + f ( x1 )
( x 0 − x1 )( x 0 − x2 ) ( x1 − x0 )( x1 − x2 )
42
( x − x0 )( x − x1 )
+ f (x 2 )
( x 2 − x0 )( x 2 − x1 )

La ecuación 39 se puede derivar de manera directa a partir del polinomio de

Newton. Sin embargo, el racional resaltado de la formulación de Lagrange se
puede captar de manera directa al darse cuenta que cada término Li (x ) será 1
en x = x i y 0 en todos los puntos de la muestra como se muestra en la figura 6.

Figura 6. Una ilustración visual del racional detrás del polinomio de Lagrange.

Pedro Ochoa Moreno Enero, 2002. 38

 Curso de Métodos Numéricos

De esta manera, cada producto Li ( x ) f ( xi ) toma el valor de f ( x i ) en el

puntote muestra x i . En consecuencia, la sumatoria de todos los productos
designados para la ecuación 39 es el único polinomio de n-ésimo orden que
pasa de manera exacta a través de todos los n+1 puntos.

Ejemplo:

Use una interpolación del polinomio de Lagrange de primer y segundo orden

para evaluar ln 2 con base en los datos dados en el ejemplo anterior.

x0 = 1 f (x 0 ) = 0
x1 = 4 f ( x1 ) = 1.386294
x2 = 6 f ( x 2 ) = 1.791760

El polinomio de primer orden se puede usar para obtener la estimación en

x = 2,

2−4 2 −1
f 1 (2) = ( 0) + (1.386294) = 0.4620981
1− 4 4 −1

De manera similar, el polinomio de segundo orden se desarrolla como

( 2 − 4)( 2 − 6) ( 2 − 1)( 2 − 6)
f 2 ( 2) = (0) + (1.386294)
(1 − 4)(1 − 6) ( 4 − 1)( 4 − 6
( 2 − 1)( 2 − 4)
+ (1.791769) = 0.5658444
( 6 − 1)(6 − 4)
Como se esperaba, ambos resultados concuerdan con los que se obtuvieron
antes al usar la interpolación para el polinomio de Newton.

Pedro Ochoa Moreno Enero, 2002. 39

 Curso de Métodos Numéricos

Regresión Lineal

Su primera situación en el ajuste de curvas podría haber sido determinar

valores intermedios a partir de datos tabulados (por ejemplo, de tablas de
interés para ingeniería económica, o a partir de tablas de vapor para
termodinámica).
El ejemplo más simple de una aproximación por mínimos cuadrados es
mediante el ajuste de un conjunto de pares de observaciones:
( x1 , y1 ), ( x2 , y 2 ),L , ( x n , y n ) , a una línea recta. La expresión matemática para
esta última es
y = a0 + a1 x + e 43

Donde a 0 y a1 son coeficientes que representan el intercepto y la pendiente,

respectivamente, y e es el error, o residuo, entre el modelo y las
observaciones, las cuales se pueden representar al reordenar la ecuación 43
como
e = y − a0 − a1 x

Así, el error o residuo es la discrepancia entre el valor real de y y el valor

aproximado, a 0 + a1 x , predicho por la ecuación lineal.
En la figura 7 podemos observar datos que exhiben un error significativo y sus
diferentes tipos de ajuste: 7.b, ajuste polinomial oscilando más allá del rango
de datos; 7.c, ajuste por mínimos cuadrados.

Figura 7. Ajuste de datos.

Pedro Ochoa Moreno Enero, 2002. 40

 Curso de Métodos Numéricos

Criterios para un “mejor” ajuste

Una estrategia para ajustar a la ¿”mejor”? línea a través de los datos podría ser
minimizar la suma de los errores residuales para todos los datos disponibles,
como en

n n

∑ ei = ∑ ( yi − a 0 − a1 xi ) 44
i =1 i =1

Donde n = número total de puntos. Sin embargo, éste es un criterio

inadecuado, como lo ilustra la 8, la cual muestra el ajuste de una línea recta de
dos puntos.

Figura 8. Ejemplos de algunos criterios para “el mejor ajuste” inadecuado para
regresión: a) minimiza la suma de los residuos, b) minimiza la suma de los
valores absolutos de los residuos y c) minimiza el error máximo de cualquier
punto individual.

Obviamente, el mejor ajuste es la línea que conecta los puntos. Sin embargo,
cualquier línea recta que pasa a través del punto medio que conecta la línea
(excepto una línea perfecta vertical) resulta en un valor mínimo de la ecuación
44 igual a cero debido a los errores que se cancelan.

Pedro Ochoa Moreno Enero, 2002. 41

 Curso de Métodos Numéricos

Por lo tanto, otro criterio lógico podría ser minimizar la suma de los valores
absolutos de las discrepancias, como en

n n

∑e i = ∑ ( yi − a0 − a1 xi )
i =1 i =1

La figura 8.b demuestra por qué este criterio es también inadecuado. Para los
cuatro puntos expuestos, cualquier línea recta que esté dentro de las líneas
punteadas minimizará el valor absoluto de la suma. Así, este criterio tampoco
da un único mejor ajuste. Una tercera estrategia para ajustar a la mejor línea
recta es el criterio de minmax. En esta técnica, la línea se elige de manera que
minimice la máxima distancia que tenga un punto individual desde la línea.
Como se ilustra en la figura 8.c, tal estrategia no es adecuada para regresión,
ya que tiene una excesiva influencia en puntos fuera del conjunto; es decir, un
solo punto con un gran error. Debe hacerse la observación que el principio
minmax es en algunas ocasiones muy adecuado para ajustar una función
simple a una función complicada.
Una estrategia que supera los defectos de los procedimientos mencionados es
minimizar la suma de los cuadrados de los residuos ente la y medida y la
y calculada con el modelo lineal

n n n
S r = ∑ ei = ∑ ( yi ,medida − y i, modelo) = ∑ ( y i − a0 − a1 xi )
2 2 2
45
i =1 i =1 i =1

Este criterio tiene varias ventajas, entre ellas el hecho de que se obtiene una
línea única para cierto conjunto de datos. Antes de analizar esas propiedades,
presentaremos una técnica para determinar los valores de a 0 y a1 que
minimizan la ecuación 45.

Ajuste por mínimos cuadrados de una línea recta

Para determinar los valores de a 0 y a1 , la ecuación 45 es diferenciada con

respecto a cada coeficiente:

∂S r
= −2∑ ( y i − a 0 − a1 x i )
∂a 0

∂S r
= −2∑ ( yi − a0 − a1 xi ) xi
∂a1

Observe que hemos simplificado los símbolos de la sumatoria; a menos que se

indique otra cosa, todas las sumatorias son de i = 1 hasta n . Al fijar esas
derivadas igual a cero, resultara en un mínimo S r . Si se hace esto, las
ecuaciones se pueden expresar como

0 = ∑ yi − ∑ a0 − ∑ a1 xi

Pedro Ochoa Moreno Enero, 2002. 42

 Curso de Métodos Numéricos

0 = ∑ yi x i − ∑ a 0 x i − ∑ a1 xi2

Ahora, si hacemos que ∑a 0 = na 0 , podemos expresar las ecuaciones como un

conjunto de dos ecuaciones lineales con dos incógnitas ( a 0 y a1 ):

na 0 + (∑ xi )a1 = ∑ y i 46

(∑ x )a + (∑ x )a = ∑ x y
i 0 i
2
1 i i 47

Estas son llamadas ecuaciones normales, y pueden ser resueltas en forma

simultánea

n∑ x i y i − ∑ x i ∑ yi
a1 =
n∑ x i2 − (∑ xi )
48

Este resultado, entonces, se puede usar en conjunto con la ecuación 46 y

resolver para

a 0 = y − a1 x 49

Donde x y y son las medias de y y x , respectivamente.

Ejemplo: (Regresión Lineal)

Ajuste a una línea recta los valores de x y y en las dos primeras columnas del
Cuadro 2.

xi yi ( y i − y) 2 ( y i − a 0 − a1 x i ) 2
1 0.5 8.5765 0.1687
2 2.5 0.8622 0.5625
3 2.0 2.0408 0.3473
4 4.0 0.3265 0.3265
5 3.5 0.0051 0.5896
6 6.0 6.6122 0.7972
7 5.5 54.2908 0.1993
Suma 24.0 22.7143 2.9911

Cuadro 2. Cálculos para un análisis de error de ajuste lineal.

Se calculan las siguientes cantidades:

Pedro Ochoa Moreno Enero, 2002. 43

 Curso de Métodos Numéricos

n =7 ∑x yi = 119.5
i ∑x2
i = 140
28
∑x i = 28 x=
7
=4
24
∑y i = 28 y=
7
= 3.428571

Mediante las ecuaciones 48 y 49

7 (119.5) − 28( 24)

a1 = = 0.8392857
7 (140) − ( 28) 2
a 0 = 3. 428571 − 0,8392857 (4) = 0.07142857

Por lo tanto, el ajuste por mínimos cuadrados es

y = 0.07142857 + 0.8392857 x

La línea junto con los datos, son mostrados en la figura 8.c.

Cuantificación del error de una regresión lineal

Cualquier otra línea que la calculada en el ejemplo anterior resulta en una gran
suma de cuadrados de los residuos. Así, la línea es única y en términos de
nuestro criterio elegido en una línea “mejor” a través de los puntos. Un número
adicional de propiedades de este ajuste se puede dilucidar al examinar más de
cerca la forma en que se calcularon los residuos. Recuerde que la suma de los
cuadrados se define como

n n
S r = ∑ ei = ∑ ( yi − a 0 − a1 xi )
2 2
50
i =1 i =1

Obsérvese que el cuadrado del residuo representa el cuadrado de la

discrepancia entre los datos y una sola estimación de la medida de tendencia
central (la media). En la ecuación 50, el cuadrado de los residuos representa el
cuadrado de la distancia vertical entre los datos y otra medida de tendencia
central: la línea recta, observe esto en la figura 9.

Pedro Ochoa Moreno Enero, 2002. 44

 Curso de Métodos Numéricos

Figura 9. El residuo en la regresión lineal representa la distancia vertical entre

un dato y una línea recta.

Pedro Ochoa Moreno Enero, 2002. 45

 Curso de Métodos Numéricos

Pedro Ochoa Moreno Enero, 2002. 46

 Curso de Métodos Numéricos

5. Integración
Tal como se representa, una interpretación física de la integración numérica es
la determinación del área bajo la curva. La integración tiene muchas
aplicaciones en la práctica de la ingeniería, empezando por la determinación de
centroides de objetos con formas extravagantes, hasta el cálculo de cantidades
totales basadas en conjuntos de medidas discretas. Adicionalmente, las
formulas de integración numérica juegan un papel importante en la solución de
ecuaciones diferenciales.

I = ∫ f ( x ) dx ≅ ∫ f n ( x ) dx
b b
51
a a

Donde f n (x ) es un polinomio de la forma

f n ( x) = a0 + a1 x + L + an −1 x n−1 + a n x n

Donde n es el orden del polinomio. Por ejemplo, en la figura 1, se usa un

polinomio de primer orden (una línea recta) como una aproximación. En la
figura 1.b, se emplea una parábola para el mismo propósito.

Figura 1. La aproximación de una integral como el área bajo a) una sola línea
recta y b) una sola parábola.

La integral se puede también aproximar mediante una serie de polinomios

aplicada por pedazos a la función o datos sobre segmentos de longitud
constante. Por ejemplo, en la figura 2, se usan tres segmentos de línea recta
para aproximar la integral.
Pueden utilizarse polinomios de orden superior para los mismos propósitos. Con
este antecedente, reconocemos que el “método de barras” en la figura 3 se
emplea una serie de polinomios de orden cero (es decir, constantes) para
aproximar la integral.

Pedro Ochoa Moreno Enero, 2002. 47

 Curso de Métodos Numéricos

Figura 2. La aproximación de una integral como el área bajo tres segmentos de

línea recta.

Figura 3. Uso de rectángulos, o barras, para aproximar la integral.

Se dispone de formas cerradas y abiertas de fórmulas de Newton-Cotes. Las

formas cerradas son aquellas donde los datos al inicio y al final de los límites de
integración son conocidos, véase figura 4.a. Las formas abiertas tienen límites
de integración que se extienden más allá del rango de los datos, como se ve en
la figura 4.b. Las formas abiertas de Newton-Cotes no se usan por lo general
para integración definida. Sin embargo, se utilizan para evaluar integrales
impropias y la solución de ecuaciones diferenciales ordinarias.

Pedro Ochoa Moreno Enero, 2002. 48

 Curso de Métodos Numéricos

Figura 4. La diferencia entre fórmulas de integración a) cerrada y b) abierta.

Regla Trapezoidal

La regla trapezoidal es la primera de las fórmulas de integración cerrada

Newton-Cotes. Corresponde al caso donde el polinomio en la ecuación 51 es de
primer orden:

I = ∫ f ( x ) dx ≅ ∫ f n ( x ) dx
b b

a a

Pero si recordamos de la unidad anterior (interpolación) que una línea recta se

representa como:

f (b) − f ( a)
f 1 ( x ) = f ( a) + ( x − a) 52
b−a

El área bajo esta línea recta es un estimado de la integral de f (x ) entre los

límites a y b :

b
 f (b) − f ( a) 
I = ∫  f (a ) + ( x − a) dx
a 
b−a 

El resultado de la integración es

f ( a) + f ( b)
I = (b − a ) 53
2

La cual se denomina regla trapezoidal.

Geométricamente, la regla trapezoidal es equivalente a aproximar el área

trapezoide bajo la línea recta que conecta f (a ) y f (b) en la figura 5.

Pedro Ochoa Moreno Enero, 2002. 49

 Curso de Métodos Numéricos

Figura 5. Ilustración gráfica de la regla trapezoidal.

Recuerde de geometría que la fórmula para calcular el área de un trapezoide es

la altura por el promedio de las bases, véase la figura 6.a. En nuestro caso, el
concepto es el mismo, pero el trapezoide está sobre su lado, véase la figura
6.b.

Figura 6. a) La fórmula para calcular el área de un trapezoide: altura por el

promedio de las bases. b) Para la regla trapezoidal, el concepto es el mismo
pero ahora el trapezoide está sobre su lado.

Por lo tanto, la integral se representa como

I ≅ ( ancho)( altura promedio ) 54

I ≅ (b − a) (altura promedio ) 55

Pedro Ochoa Moreno Enero, 2002. 50

 Curso de Métodos Numéricos

Donde, para la regla trapezoidal, la altura promedio es el promedio de los

valores de la función en los puntos extremos, o [ f (a ) + f (b) ] 2.
Todas las fórmulas de Newton-Cotes pueden expresarse en el formato general
de la ecuación 55. De hecho, solo la difieren con respecto a la formulación de la
altura promedio.

Error de la regla trapezoidal

Cuando empleamos la integral bajo un segmento de línea recta para aproximar

la integral bajo la curva, obviamente podemos incurrir en un error que puede
ser sustancial, véase la figura 7.

Figura 7. Ilustración gráfica del uso de una sola aplicación de la regla

trapezoidal para aproximar la integral de
f ( x ) = 0.2 + 25 x − 200 x + 675 x − 900 x + 400 x de x = 0 a 0.8.
2 3 4 5

Una estimación para el error de truncamiento local de una sola aplicación de la

regla trapezoidal es

1 ''
Et = − f (ξ )( b − a ) 3 56
12

Donde ξ está en algún lugar en el intervalo de a y b . La ecuación 56 indica

que si la función sujeta a integración es lineal, la regla trapezoidal será exacta.
De otra manera, para funciones con derivadas de segundo orden y superior (es
decir, con curvatura), puede ocurrir algún error.

Pedro Ochoa Moreno Enero, 2002. 51

 Curso de Métodos Numéricos

Ejemplo: Aplicación simple de la regla trapezoidal

Use la ecuación 53 para integrar numéricamente

f ( x ) = 0.2 + 25 x − 200 x 2 + 675 x 3 − 900 x 4 + 400 x 5

Desde a = 0 a b = 0.8. Recuerde que el valor exacto de la integral se puede

determinar de forma analítica y es 1.640533.

Los valores de la función son entonces

f (0) = 0.2
f (0.8) = 0.232

Pueden sustituirse en la ecuación 53 para dar

(0.2 + 0.232)
I ≅ ( 0.8) = 0.1728
2

El cual representa un error de

Et = 1,640533 − 0.1728 = 1.467733

Que corresponde a error relativo porcentual de ε t = 89.5% . La razón para este

error tan grande es evidente de la gráfica de la figura 7. Observe que el área
bajo la línea recta no toma una porción significativa de la integral que esta por
arriba de la línea. En situaciones actuales, podríamos no conocer previamente
el valor real. Por lo tanto, se requiere una aproximación del error estimado.
Para obtener dicha estimación, la segunda derivada de la función sobre el
intervalo podría calcularse al diferenciar dos veces la función original para dar

f ' ( x) = −400 + 4050 x − 10800 x 2 + 8000 x 3

El valor promedio de la segunda derivada se puede calcular mediante la

ecuación de la media para datos continuos o discretos, y gráficamente esta se
puede observar en la figura 8.

Media =
∫
a
f ( x )dx
b−a

Pedro Ochoa Moreno Enero, 2002. 52

 Curso de Métodos Numéricos

∫ (−400 + 4050 x − 10800 x + 8000 x 3 )dx

2

f '' ( x) = a = − 60
0 .8 − 0

Que podría sustituirse en la ecuación 56 para obtener

1
Ea = − ( −60)(0.8) 3 = 2.56
12

La cual es del mismo orden de magnitud y signo que el error real. Sin embargo,
existe una discrepancia, ya que, de hecho, para un intervalo de este tamaño el
promedio de la segunda derivada no es necesariamente una aproximación
exacta de f '' (ξ ). Así, indicamos que el error es aproximado mediante la
notación E a , más exacto que usar Et .

Figura 8. Una ilustración de la media para datos continuos.

Aplicación múltiple de la regla trapezoidal

Una forma de mejorar la exactitud de la regla trapezoidal es dividir el intervalo

de integración desde a a b en un número de segmentos y aplicar el método a
cada uno de ellos, véase la figura 9. El área de segmentos individuales se
puede entonces agregar para dar la integral para todo el intervalo. Las
ecuaciones resultantes son llamadas fórmulas de integración, de múltiple
aplicación o compuestas.
La figura 10 muestra el formato general y la nomenclatura que usaremos para
caracterizar integrales de múltiple aplicación. Hay n+1 puntos base igualmente
espaciados ( x0 , x1 , x 2 ,L , xn ) . En consecuencia, hay n segmentos de igual
anchura:
b−a
h= 57
n
Si a a b son designados como x 0 y x n , respectivamente, la integral total se
representará como

Pedro Ochoa Moreno Enero, 2002. 53

 Curso de Métodos Numéricos

Figura 9. Ilustración de la regla trapezoidal de aplicación múltiple: a) Dos

segmentos, b) tres segmentos, c) cuatro segmentos y d) cinco segmentos.

I = ∫ f ( x) dx + ∫ f ( x ) dx + L + ∫
x1 x2 xn
f ( x) dx
x0 x1 x n− 1

Al sustituir la regla trapezoidal para cada integral se obtiene

Figura 10. El formato general y nomenclatura para integrales de aplicación

múltiple.

Pedro Ochoa Moreno Enero, 2002. 54

 Curso de Métodos Numéricos

f ( x 0 ) − f ( x1 ) f ( x1 ) − f ( x2 ) f ( x n−1 ) − f ( x n )
I =h +h +L+h 58
2 2 2

0 mediante la agrupación de términos,

h n

I = 
2
f ( x0 ) + 2 ∑ f ( xi ) + f ( xn )  59
i =1 

O, usando la ecuación 57 para expresar la ecuación 59 en el formato general de

la ecuación 55

n
f ( x0 ) +2∑ f ( xi ) + f ( x n )
i =1
I = (b − a ) 60
2n

Ancho Altura promedio

Como la sumatoria de los coeficientes de f (x ) en el numerador dividido entre

2n es igual a 1, la altura promedio representa un promedio ponderado de los
valores de la función. De acuerdo con la ecuación 60, los puntos interiores dan
dos veces el peso de los dos puntos extremos f ( x 0 ) y f ( x n ) .
Un error para la regla trapezoidal de múltiple aplicación se puede obtener al
sumar los errores individuales de cada segmento para dar

(b − a ) 3 n
Et = −
12n 3
∑f
i =1
''
(ξ i ) 61

Donde f '' (ξ i ) es la segunda derivada en un punto ξ i localizado en el

segmento i . Este resultado se puede simplificar al estimar la media o valor
promedio de la segunda derivada para todo intervalo como

∑f
i =1
''
(ξ i )
f ≅
''
62
n

Por lo tanto, ∑f ''

(ξ i ) ≅ n f '' y la ecuación 61 se puede reescribir como

(b − a ) 3
Ea = − f ''
63
12n 2

Así, si el número de segmentos se duplica, el error de truncamiento disminuirá

a un cuarto. Observe que la ecuación 63 es un error aproximado debido a la
naturaleza aproximada de la ecuación 62.

Pedro Ochoa Moreno Enero, 2002. 55

 Curso de Métodos Numéricos

Ejemplo: regla trapezoidal de múltiple aplicación

Use dos segmentos de la regla trapezoidal para estimar la integral de

f ( x ) = 0.2 + 25 x − 200 x 2 + 675 x 3 − 900 x 4 + 400 x 5

Desde a = 0 a b = 0.8. Recuerde que el valor exacto de la integral se puede

determinar de forma analítica y es 1.640533.

Entonces, n = 2(h = 0.4) :

f (0) = 0.2 f (0.4) = 2.456 f (0.8) = 0.232

0.2 + 2( 2.456) + 0.232

I = (0.8) =1.0688
4

Et = 1.640533 − 1.0688 = 0.57173 ε t = 34.9%

( 0.8) 3
Ea = − (−60) = 0.64
12( 2) 2

Donde -60 es el promedio de la segunda derivada, determinada en el ejemplo

anterior.
Los resultados del ejemplo anterior, junto con valores de n = 3 hasta n = 10 se
presentan en la tabla1. Observe como el error disminuye en tanto que el
número de segmentos aumenta.

n h I ε t (%)
2 0.4 1.0688 34.9
3 0.2667 1.3695 16.5
4 0.2 1.4848 9.5
5 0.16 1.5399 6.1
6 0.1333 1.5703 4.3
7 0.1143 1.5887 3.2
8 0.1 1.6008 2.4
9 0.0889 1.6091 1.9
10 0.08 1.6150 1.6

Pedro Ochoa Moreno Enero, 2002. 56

 Curso de Métodos Numéricos

Reglas de Simpson

Además de aplicar la regla trapezoidal con segmentación fina, otra forma de

obtener una estimación más exacta de una integral es con el uso de polinomios
de orden superior para conectar los puntos. Por ejemplo, si hay un punto extra
a la mitad del camino entre f (a ) y f (b) , los tres puntos se pueden conectar
con una parábola, véase figura 11.a. Si hay dos puntos igualmente espaciados
entre f (a ) y f (b) , los cuatro puntos se pueden conectar con un polinomio de
tercer orden, como se observa en la figura 11.b. Las fórmulas que resultan al
tomar las integrales bajo esos polinomios son conocidas como reglas de
Simpson.

Figura 11. a) Ilustración gráfica de la regla de Simpson 1/3: consiste en tomar

el área de una parábola que conecta tres puntos. b) Ilustración gráfica de la
regla de Simpson 3/8: consiste en tomar el área bajo una ecuación cúbica que
conecta cuatro puntos.

Regla de Simpson 1/3

La regla de Simpson 1/3 resulta cuando una interpolación polinomial de

segundo orden es sustituida en la ecuación 51:

I = ∫ f ( x ) dx ≅ ∫ f 2 ( x) dx
b b

a a

Si a y b se designan como x 0 , x 2 y f 2 ( x) es representada por un polinomio

de Lagrange de segundo orden, entonces la integral se transforma en

x2  ( x − x )( x − x ) ( x − x0 )( x − x 2 ) ( x − x0 )( x − x1 ) 
I = ∫x  1 2
f (x 0 ) + f ( x1 ) + f ( x 2) dx
0
 ( x 0 − x1 )( x0 − x 2 ) ( x1 − x0 )( x1 − x 2 ) ( x 2 − x0 )( x 2 − x1 ) 

Después de la integración y manejo algebraico, resulta la siguiente fórmula:

Pedro Ochoa Moreno Enero, 2002. 57

 Curso de Métodos Numéricos

I ≅
h
[ f ( x 0 ) + 4 f ( x1 ) + f ( x 2 )] 64
3

Donde, para este caso, h = ( b − a) / 2. Esta ecuación es conocida como regla de

Simpson 1/3. Es la segunda fórmula de integración cerrada de Newton-Cotes.
La especificación “1/3” surge del hecho de que h está dividida entre 3 en la
ecuación 64. La regla de Simpson 1/3 se expresa también con el uso del
formato de la ecuación 55.

f ( x 0 ) +4 f ( x1 ) + f ( x 2 )
I = (b − a ) 65
6

Ancho Altura promedio

Donde a = x0 , b = x 2 y x = punto a la mitad del camino entre a y b , que está

dado por (b + a) / 2. Observe que, de acuerdo a la ecuación 65, el punto medio
está ponderado por dos tercios y los dos puntos extremos por un sexto.
Se puede mostrar que una aplicación de un solo segmento de la regla de
Simpson 1/3 tiene un error de truncamiento, como se observa a continuación:

f ( x0 ) +4 f ( x1 ) + f ( x 2 ) 1 ( 4)
I = (b − a ) − − f (ξ ) h 5 66
6 90

Regla de Simpson 1/3 Error de truncamiento

De donde,

1 5 ( 4)
Et = − h f (ξ )
90

o, como h = ( b − a) / 2,

(b − a) 5 ( 4)
Et = − f (ξ ) 67
2880

Donde ξ está en algún lugar en el intervalo desde a y b . Así, la regla de

Simpson 1/3 es más exacta que la regla trapezoidal. Sin embargo, en
comparación con la ecuación 67 indica que es más exacta de lo esperado. En
lugar de ser proporcional a la tercera derivada, el error es proporcional a la
cuarta derivada. Esto es porque, si analizamos el termino del coeficiente de
tercer orden este tendera a cero durante la integración polinomial. En
consecuencia, la regla de Simpson 1/3 tiene una precisión de tercer orden aun

Pedro Ochoa Moreno Enero, 2002. 58

 Curso de Métodos Numéricos

cuando se base solo en tres puntos. En otras palabras, da resultados exactos

para polinomios cúbicos aun cuando se derive de una parábola.

Ejemplo: Aplicación simple de la regla de Simpson 1/3

Use la ecuación 65 para integrar

f ( x ) = 0.2 + 25 x − 200 x 2 + 675 x 3 − 900 x 4 + 400 x 5

Desde a = 0 a b = 0.8. Recuerde que la integral exacta es 1.640533.

f (0) = 0.2 f (0.4) = 2.456 f (0.8) = 0.232

Por lo tanto la ecuación 65 se puede usar para calcular

0.2 + 4( 2.456) + 0.232

I = (0.8 − 0) =1.367467
6

Que representa un error exacto de

Et =1.640533 − 1.367467 = 0.2730667 ε t =16.6 %

El cual es aproximadamente 5 veces más preciso que una aplicación de la regla

trapezoidal.
El error estimado se calcula usando la ecuación 67 y este es

(0.8) 5
Ea = − ( −2400) = 0.2730667
2880

Donde -2400 es el promedio de la cuarta derivada para el intervalo. Como

ocurrió en el ejemplo anterior a este, el error es aproximado ( Ea ) y es debido a
que el promedio de la cuarta derivada no es una estimación exacta de f ( 4) (ξ ).
Sin embargo, como en este caso tiene que ver con polinomios de quinto orden,
el resultado concuerda.

La regla de Simpson 1/3 de aplicación múltiple

Así, como con la regla trapezoidal, la regla de Simpson se puede mejorar al

dividir el intervalo de integración en un número de segmentos de igual anchura
como se observa en la figura 12.

b−a
h= 68
n

La integral total se puede representar como

Pedro Ochoa Moreno Enero, 2002. 59

 Curso de Métodos Numéricos

I = ∫ f ( x) dx + ∫ f ( x) dx + L + ∫
x2 x4 xn
f ( x) dx
x0 x2 x n −2

Figura 12. Representación gráfica de la regla de Simpson 1/3 de aplicación

múltiple. Nótese que el método se puede emplear sólo si el número de
segmentos es par.

Al sustituir la regla de Simpson 1/3 para la integral individual se obtiene

f ( x0 ) − 4 f ( x1 ) + f ( x 2 ) f ( x 2 ) + 4 f ( x3 ) + f ( x 4 )
I ≅ 2h + 2h
6 6
f ( x n−2 ) + 4 f ( x n−1 ) + f ( x n )
+ L + 2h
6

o combinando términos y usando la ecuación 68,

n −1 n −2
f ( x 0 ) +4 ∑ f ( xi ) + 2
i =1, 3, 5
∑ f (x ) + f (x
j =2 , 4 , 6
j n )
I ≅= (b − a ) 69
3n

Ancho Altura promedio

Observe que, como se ilustra en la figura 12, se debe utilizar un número par de
segmentos para implementar el método. Además los coeficientes “4” y “2” en la
ecuación 65 podrían parecer peculiares a primera vista. Sin embargo, siguen en
forma natural la regla de Simpson 1/3. Los puntos nones representan el
término medio para cada aplicación y, por lo tanto, el peso de 4 en la ecuación
65. Los puntos pares son comunes en las aplicaciones adyacentes y por lo
tanto, se cuentan dos veces.
Un error estimado para la aplicación de la regla de Simpson se obtiene en la
misma forma que para la regla trapezoidal, sumando los errores individuales de
los segmentos y sacando el promedio de la derivada para dar

Pedro Ochoa Moreno Enero, 2002. 60

 Curso de Métodos Numéricos

(b − a ) 5 ( 4)
Ea = − f 66
180n 4

donde f ( 4) es el promedio de la cuarta derivada para el intervalo.

Ejemplo: Versión de la regla de Simpson 1/3 de aplicación múltiple

Use la ecuación 65 con n = 4 para estimar la integral de

f ( x ) = 0.2 + 25 x − 200 x 2 + 675 x 3 − 900 x 4 + 400 x 5

desde a = 0 a b = 0.8. Recuerde que la integral exacta es 1.640533.

n = 4 (= 0.2) :

f (0) = 0.2 f (0.2) = 1.288

f (0.4) = 2.456 f (0.6) = 3.464
f (0.8) = 0.232

A partir de la ecuación 65,

0.2 + 4(1.288 + 3.464) + 2(2.456) + 0.232

I = (0.8 − 0) =1.623467
12

Et =1.640533 − 1.623467 = 0.017067 ε t =1.04 %

El error estimado con base a la ecuación 66, y este es

(0.8) 5
Ea = − ( −2400) = 0.017067
180( 4) 4

Regla de Simpson 3/8

En una manera similar a la derivación de la regla de trapezoidal y de Simpson

1/3, un polinomio de Lagrange de tercer orden se puede ajustar a cuatro
puntos e integrarse:

I = ∫ f ( x ) dx ≅ ∫ f 3 ( x ) dx
b b

a a

Para obtener

I ≅
3h
[ f ( x 0 ) + 3 f ( x1 ) + 3 f ( x2 ) + f ( x3 ) ]
8

Pedro Ochoa Moreno Enero, 2002. 61

 Curso de Métodos Numéricos

Donde h = ( b − a) / 3. Esta ecuación se llama regla de Simpson 3/8 debido a que

h se multiplica por 3/8. Esta es la tercera formula de integración cerrada de
Newton-Cotes. La regla de 3/8 se puede expresar también en la forma de la
ecuación 55:

f ( x 0 ) +3 f ( x1 ) + 3 f ( x 2 ) + f ( x 3 )
I ≅= (b − a ) 70
8

Ancho Altura promedio

Así, los dos puntos interiores tienen pesos de tres octavos, mientras los puntos
extremos pesan un octavo. La regla de Simpson 3/8 tiene un error de

3 5 5
Et = − h f (ξ )
80

o, ya que h = ( b − a) / 3,

(b − a ) 5 5
Et = − f (ξ ) 71
6480

Como el denominador de la ecuación 71 es mayor que el de la ecuación 66, la

regla 3/8 es algo más exacta que la de 1/3.

Ejemplo: Regla de Simpson 3/8

a) Use la regla de Simpson 3/8 para integrar

f ( x ) = 0.2 + 25 x − 200 x 2 + 675 x 3 − 900 x 4 + 400 x 5

Desde a = 0 a b = 0.8.

b) Usela en conjunto con la regla de Simpson 1/3 a fin de integrar la misma

Función para cinco segmentos.

a) una sola aplicación de la regla de Simpson 3/8 requiere de cuatro puntos

igualmente espaciados:

f (0) = 0.2 f (0.2667 ) = 1.432724

f (0.53334) = 3.487177 f (0.8) = 0.232

Usando la ecuación 70

Pedro Ochoa Moreno Enero, 2002. 62

 Curso de Métodos Numéricos

0.2 + 3(1.432724 + 3.487177) + 0.232

I ≅ ( 0 .8 − 0 ) =1.519170
8

Et =1.640533 − 1.519170 = 0.01213630 ε t = 7.4 %

(0.8) 5
Ea = − ( −2400) = 0.1213630
6480

b) Los datos necesarios para la aplicación de cinco segmentos ( h = 0.16 ) es

f (0) = 0.2 f (0.16) = 1.296919

f (0.32) = 1.7433933.487177 f (0.48) = 0.3.186015
f (0.64) = 3.181929 f (0.80) = 0.232

La integral para los dos primeros segmentos se obtiene usando la regla

de Simpson 1/3:

0.2 + 4(1.296919) + 1.743393

I ≅ ( 0.32) = 0.3803237
6

Para los últimos tres segmentos, la regla de 3/8 se puede usar para
obtener:

1.743393 + 3( 3.186015 + 3.181929) + 0.232

I ≅ ( 0.48) =1.264754
8

La integral total se calcula al sumar los dos resultados:

I = 0.3803237 + 1.264754 = 1.65077

Et =1.640533 − 1.645077 = − 0.0045383 ε t = − 0.28 %

Pedro Ochoa Moreno Enero, 2002. 63

 Curso de Métodos Numéricos

Algunos problemas de aplicación de la integral:

Figura 13. Una sección transversal de una corriente.

Figura 15. Un campo límitado por dos caminos y una cresta.

Pedro Ochoa Moreno Enero, 2002. 64

 Curso de Métodos Numéricos

Figura 16. Agua que ejerce presión sobre la cara de una presa corriente arriba:
a) vista lateral donde se muestra que la fuerza aumenta linealmente con la
profundidad; b) vista frontal en la cual se muestra el ancho de la presa en
metros.

Figura 17. a) Una barra bajo carga axial y b) la curva resultante esfuerzo-
deformación unitaria donde el esfuerzo en kilolibras por pulgada cuadrada
(10³lb/pulg²) y la deformación unitaria es adimensional.

Pedro Ochoa Moreno Enero, 2002. 65

 Curso de Métodos Numéricos

Pedro Ochoa Moreno Enero, 2002. 66

 Curso de Métodos Numéricos

6. Ecuaciones diferenciales ordinarias

Cuando una función incluye una variable dependiente y una independiente, la

podemos llamar Ecuación Diferencial Ordinaria (EDO). Esta se clasifica de
acuerdo a su orden o grado. Por ejemplo,

dx c
1. =g − v Ecuación de 1er. orden
dt m

d2x dx
2. m 2 + b + kx = 0 Ecuación de 2do. Orden
dt dt

Las leyes fundamentales de la física, la mecánica, la electricidad y la

termodinámica se basan en general en observaciones empíricas que explican la
variación de las propiedades físicas y estados de los sistemas. En lugar de
describir el estado de los sistemas físicos directamente, las leyes se expresan
en cambios del tiempo y el espacio. En el siguiente cuadro se muestran varios
ejemplos:

Expresión Variables y
Ley matemática parámetros
Segunda ley dv F Velocidad (v),
=
de Newton dt m Fuerza (F), masa( m)
Ley del calor ∂T Conductividad térmica (k)
Flujo de calor = k
de Fourier ∂x y temperatura (T)
Ley de difusión ∂c Coeficiente de difusión (D)
Flujo de masa = − D
Fick ∂x y concentración (c)
di Inductancia(L) y
Caida de voltaje= L
Ley de Faraday dt corriente(i)

Conservación de dc Volumen (V) y

Acumulación = V
la masa dt concentración (c)

Cuadro. Ejemplos de las leyes fundamentales escritas en términos del promedio

de cambio de las variables (t = tiempo y x = posición).

Pedro Ochoa Moreno Enero, 2002. 67

 Curso de Métodos Numéricos

Figura 1. La secuencia de eventos en la aplicación de EDO para resolver

problemas de ingeniería. El ejemplo mostrado es la velocidad de un paracaidista
en caída.

Fundamentos matemáticos

Una solución de ecuación diferencial ordinaria es una función específica de la

variable independiente y de sus parámetros que satisfacen la ecuación
diferencial original. Para ilustrársete concepto, se tiene la siguiente función:

y = − 0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1 72

La cual es un polinomio de cuarto orden. Ahora, si se deriva la ecuación 72, se

obtiene la EDO:

dy
= − 2 x 3 + 12 x 2 − 20x + 8.5 73
dx

Esta ecuación representa o describe el comportamiento del polinomio, pero de

manera diferente al de la ecuación 72. En vez de representar explícitamente los
valores de y para cada uno de los valores de x , la ecuación 73 proporciona la
relación del cambio de y respecto a x (esto es, la pendiente) para cada valor
de x . Esto se puede visualizar en la figura 2.
Nótese que los valores cero de la derivada corresponden a un punto donde la
función originales plan, esto es tiene una pendiente cero. También, los valores
absolutos máximos alcanzados de las derivadas son los extremos del intervalo
en donde las pendientes de una función son más grandes. Aunque, como ya se
ha demostrado, se puede determinar una ecuación diferencial dando a la
función original, el objetivo aquí es determinar la función original dada la
ecuación diferencial. La función original representa entonces la solución. En
este caso, se puede determinar la solución en forma analítica, integrando la
ecuación 73,

Pedro Ochoa Moreno Enero, 2002. 68

 Curso de Métodos Numéricos

Figura 2. a) y contra x y b) dy dx contra x de la función

y = − 0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1 .

[ ]
y = ∫ − 2 x 3 + 12 x 2 − 20 x + 8.5 dx

Aplicando la regla de integración,

u n =1
∫ u du = +c ; n ≠1
n

n +1

de donde, la solución es:

y = − 0.5 x 4 + 4 x 3 − 10 x 2 + 8.5 x + c 74

que es idéntica a la función original con una notable excepción. En el acto de

derivación e integración, se pierde el valor de la constante (1) en la ecuación
original y se gana el valor c. Esta es conocida como la constante de integración.
El hecho de que aparezca una constante indica que la solución no es única. De
hecho, esta es solo una de un número infinito de funciones posibles que
satisfacen a la ecuación diferencial. Por ejemplo, en la figura 3, se muestran
seis posibles funciones que satisfacen la ecuación 73.
Por lo tanto, para especificar la solución completamente, una ecuación
diferencial se acompaña de condiciones auxiliares. Para EDO’s de primer orden,
a un tipo de condición auxiliar se le llama valor inicial y es necesaria para
determinar la constante y obtener una solución única. Por ejemplo, la ecuación

Pedro Ochoa Moreno Enero, 2002. 69

 Curso de Métodos Numéricos

Figura 3. Seis posibles soluciones de la integral − 2 x 3 + 12 x 2 − 20 x + 8.5 . Cada

una conforma un valor diferente de la constante c.

73 puede ir acompañada de una condición inicial en que x = 0, y = 1 . Estos

valores pueden sustituirse en la ecuación 74:

1 = −0.5(0) 4 + 4( 0) 3 − 10( 0) 2 + 8.5( 0) + c 75

para determinar c = 1. Por lo tanto, la solución única que satisface a la ecuación

diferencial y a la condición inicial especificada se obtiene sustituyendo c = 1. en
la ecuación 74 para obtener:

y = − 0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1 76

de esta manera, se ha considerado en la ecuación 74 la condición inicial y al

hacerlo, se ha obtenido una solución única para la EDO y se ha completado el
ciclo hasta la función inicial.
Las condiciones iniciales por lo común tienen interpretaciones muy tangibles en
ecuaciones diferenciales de problemas físicos. Cuando se trata con ecuaciones
diferenciales de n-ésimo orden se requieren de n-condiciones para obtener una
solución única. Si todas las condiciones se especifican en el mismo en el mismo
valor de la variable independiente (por ejemplo, en x o t = 0 ), entonces al
problema se le conoce como problema de valor inicial. Esto contrasta con los
problemas de valor en la frontera en donde las especificaciones de las
condiciones ocurren en valores diferentes de la variable independiente.

Esta unidad se orienta a la resolución de ecuaciones diferenciales de la forma

dy
= f ( x, y)
dx

Para resolver esta ecuación podemos utilizar el método de un paso que se da

en forma general a continuación:

Pedro Ochoa Moreno Enero, 2002. 70

 Curso de Métodos Numéricos

Nuevo valor = valor anterior + pendiente x tamaño de paso

o en términos matemáticos

y i+1 = y i + φh 77

De acuerdo con esta ecuación, la pendiente estimada φ se usa para extrapolar

desde un valor anterior yi a un nuevo valor y i+1 en una distancia h , véase la
figura 4. Esta fórmula se puede aplicar paso a paso para calcular el valor en el
futuro y, por lo tanto trazar la trayectoria de la solución.

Figura 4. Ilustración gráfica del método de un paso.

Todos los métodos de un paso se pueden expresar en esta forma general, que
sólo va a diferir en la manera en la cual se estima la pendiente. Entonces la
pendiente al inicio del intervalo es tomada como una aproximación de la
pendiente promedio sobre todo el intervalo. Este procedimiento, llamado
método de Euler, se analiza a continuación. Después continuaremos con otros
métodos de un paso que cumplen estimaciones de pendientes en forma alterna,
y cuyas resultantes serán predicciones más exactas. Todas estas técnicas se
conocen por lo general como métodos de Runge-Kutta.

Método de Euler

La primera derivada proporciona una aproximación directa a la pendiente en x i ,

como se muestra en la figura 5:

dy
φ = f ( xi , yi ) ⇒ f ( x, y)
dx

donde f ( xi , yi ) es la ecuación diferencial evaluada en xi y y i . Esta

aproximación se sustituye en la ecuación 77:

Pedro Ochoa Moreno Enero, 2002. 71

 Curso de Métodos Numéricos

Figura 5. Método de Euler.

y i+1 = y i + f ( x i , y i ) h 78

A esta fórmula se le conoce como el método de Euler (o método de Euler-

Cauchy o de pendiente puntual). Se predice un nuevo valor de y usando la
pendiente (igual a la primera derivada en el valor original de x ) para extrapolar
linealmente sobre el tamaño del paso h (como se observa en la figura 5).

Ejemplo:

Utilicese, el método de Euler para integra numéricamente la ecuación:

dy
= − 2 x 3 + 12 x 2 − 20x + 8.5
dx

de x = 0 hasta x = 4 con un tamaño de paso de 0.5. La condición inicial en

x = 0 es y = 1.

y = − 0.5x 4 + 4 x 3 − 10 x 2 + 8.5x + 1

Se puede usar la ecuación 78, para implementar el método de Euler.

y i+1 = y i + f ( x i , y i ) h

y (0.5) = y (0) + f ( 0,1)(0.5)

de donde y (0 ) = 1. y la aproximación a la pendiente en x = 0 es:

f (0,1) = − 2(0) 3 + 12(0) 2 − 20(0) + 8.5 = 8.5

Pedro Ochoa Moreno Enero, 2002. 72

 Curso de Métodos Numéricos

Por lo tanto,

y ( 0.5) = 5 .25 ; y ( 0) = 1

Valor aproximado

Pero la solución verdadera en x = 0.5 es:

y (0.5) = − 0.5( 0.5) 4 + 4( 0.5) 3 − 10(0.5) 2 + 8.5(0.5) + 1 = 3.21875

Valor verdadero

Por lo tanto el error es:

Et = 3.21875 − 5.25 = − 2.03125

Expresado como relativo porcentual Et = − 63.1%

En el segundo paso:

y (1.0) = y( 0.5) + f (0.5, 5.25)(0.5) = 5.875

La solución en x = 1.0 es 3.0 y por lo tanto, el error relativo porcentual es -

95.8%. El cálculo se repite y los resultados se compilan en la tabla 1 que se
muestra a continuación:

x y verdadero y Euler Error relativo porcentual

Global Local
0.0 1.00000 1.00000
0.5 3.21875 5.25000 -63.1 -63.1
1.0 3.00000 5.87500 -95.8 -28.0
1.5 2.21875 5.12500 131.0 -1.41
2.0 2.00000 4.50000 -125.0 20.5
2.5 2.71875 4.75000 -74.7 17.3
3.0 4.00000 5.87500 46.9 4.0
3.5 4.71785 7.12500 -51.0 -11.3
4.0 3.00000 7.00000 -133.3 -53.0

Tabla 1. Comparaciones de los valores verdadero y aproximado de la integral

y ' = − 2 x 3 + 12 x 2 − 20 x + 8.5 , con la condición inicial de y = 1 en x = 0 .

En la tabla 1, los valores aproximados se calcularon mediante el método de

Euler con un tamaño de paso de 0.5. El error local se refiere al error incurrido
sobre un solo paso. Se calcula con una expansión de la serie de Taylor. El error
global es la discrepancia total debida a los pasos anteriores así como a los

Pedro Ochoa Moreno Enero, 2002. 73

 Curso de Métodos Numéricos

actuales. Observe que aunque el cálculo captura la tendencia general de la

solución verdadera, el error es considerable. Una alternativa para reducir este
error es la usar un tamaño de paso más pequeño, esto se vera en la siguiente
sección.
En la figura 6, tenemos el resultado gráfico del problema anterior.

Figura 6. Comparación de la solución verdadera con una solución numérica

usando el método de Euler.

Análisis del error en el método de Euler

La solución numérica de una ecuación diferencial ordinaria incluye dos tipos de

error:

1. Errores de truncamiento causados por la naturaleza de los métodos

empleados en la aproximación de y .
2. Errores de redondeo causados por el número limitado de dígitos o
cifras que puede retener la computadora.
Los errores de truncamiento se componen de dos partes. L a primera es un
error de truncamiento local que resulta al aplicar el método en cuestión en un
paso. El segundo es un error de programación que resulta de las
aproximaciones producidas durante los pasos anteriores. La suma de los dos
errores es el error de truncamiento global.
El conocimiento de la magnitud y propiedades del error de truncamiento se
puede obtener derivando el método de Euler directamente de la expansión de
Taylor. Con el fin de hacer esto, recuerdese que la ecuación diferencial que se
esta integrando será de la forma general:

y ' = f ( x, y) 79

donde y ' = dy dx , x y y son las variables independientes y dependientes

respectivamente. Si la solución, esto es, la función que describe el
comportamiento de y tiene derivadas continuas, esta se puede representar

Pedro Ochoa Moreno Enero, 2002. 74

 Curso de Métodos Numéricos

mediante una expansión de la serie de Taylor alrededor del punto de inicial

( xi , y i ), como en:

yi'' 2 y ( n)
y i+1 = yi + y 1i h + h +L + i h n + Rn 80
2! n!

donde h = xi+1 − x i y Rn es el término residual definido como

y n+1 n +1
Rn = h 81
( n + 1)!

Sustituyendo la ecuación 79 en la ecuación 80, obtenemos

(n −1 )
f ' (x i , y i ) 2 f ( xi , yi ) n
y i+1 = y i + f ( x i , y i ) h + h +L + h + 0 (h n+1 ) 82
2! n!

Comparando las ecuaciones 78 y 82, puede verse que el método de Euler

corresponde a la serie de Taylor truncada hasta el término f ( xi , y i )h . Por
ejemplo, el error de truncamiento en el método de Euler es atribuible a los
términos restantes de la expansión que no se incluyen en la ecuación 78.
Restando la ecuación 79 de la 82 se obtiene

f ' (x i , y i ) 2
Et = h + L + 0( h n+1 ) 83
2!

donde Et es el error de truncamiento local. Para un h lo suficientemente

pequeña, los errores en los términos de la ecuación 83 decrecen por lo común a
medida que el orden crece y el resultado a menudo se representa como:

f ' (x i , y i ) 2
Ea = h 84
2!

E a = 0( h 2 ) 85

donde Ea es el error de truncamiento local aproximado.

Ejemplo:

Utilicese la ecuación 83 para aproximar el error del primer paso del ejemplo
anterior. Usese también esa ecuación para determinar el error ocasionado por
cada uno de los términos de orden superior de la expansión de Taylor.

Pedro Ochoa Moreno Enero, 2002. 75

 Curso de Métodos Numéricos

f ' ( xi , y i ) 2 f '' ( x i , y i ) 3 f ''' ( x i , y i ) 4

Et = h + h + h a
2! 3! 4!

Si f ( x i , y i ) = − 2 x 3 + 12 x 2 − 20 x + 8.5

f ' ( xi , yi ) = − 6 x 2 + 24 x − 20 b
f '' ( xi , y i ) = − 12 x + 24 x c
f ''' ( x i , y i ) = − 12 d

Por lo anterior, podemos decir que se pueden omitir los términos adicionales,
dado que estos son cero. En este caso las ecuaciones a,…,d definen
completamente el error de truncamiento de una aplicación del método de Euler.
Entonces el error debido al truncamiento del segundo término se puede calcular
como:

− 6(0) 2 + 24( 0) − 20
Et , 2 = ( 0.5) 2 = − 2.5
2

− 12(0) + 24
Et ,3 = ( 0.5) 2 = 0.5
6
− 12
Et , 4 = ( 0.5) 2 = − 0.03125
24

Estos tres valores se pueden sumar para obtener el error total de truncamiento:

Et = Et , 2 + Et, 3 + Et, 4 = − 2.5 + 0.5 − 0.03125 = −2.03125

Observe, como Et , 2 > Et , 3 > Et , 4 , el cual soporta la aproximación por la ecuación

84.

Método mejorado del polígono (Euler modificado)

En la figura 7 se ilustra la modificación simple del método de Euler. Este

método usa el de Euler para predecir un valor de y en el punto medio del
intervalo:

h
y i+1/ 2 = y i + f ( x i , y i ) 86
2

Entonces este valor predecido se usa en la aproximación de la pendiente en el

punto medio:

y 'i +1 / 2 = f ( xi +1/ 2 , yi +1/ 2 ) 87

Pedro Ochoa Moreno Enero, 2002. 76

 Curso de Métodos Numéricos

Lo cual, se supone, representa una aproximación valida de la pendiente

promedio en el intervalo completo. Esta pendiente se usa para extrapolar
linealmente de x i a xi +1 usando el método de Euler:

y i+1 = y i + f ( x i+1/ 2 , y i+1/ 2 ) h 88

Nótese que debido a que y i+1 no esta en ambos lados, la correctora (ecuación
88) no se puede aplicar iterativamente para mejorar la solución.

Figura 7. Esquema gráfico del polígono mejorado: a) ecuación 86, y b) ecuación

88.
El método del polígono mejorado es superior al método de Euler ya que este
utiliza una aproximación de la pendiente en el punto medio del intervalo de
predicción. Los errores local y global del método del polígono mejorado son
0( h 3 ) y 0( h 2 ) respectivamente.

Métodos de Runge-Kutta

Los métodos de Runge-Kutta tienen la exactitud del esquema de la serie de

Taylor sin necesitar del cálculo de derivadas superiores. Existen muchas
variaciones pero todas ellas se pueden ajustar a la forma general de la ecuación
89:

y i+1 = y i + φ ( x i , y i , h) h 89

Pedro Ochoa Moreno Enero, 2002. 77

 Curso de Métodos Numéricos

donde φ ( xi , yi , h) se le llama función de incremento y puede interpretarse como

el promedio de la pendiente sobre el intervalo. La función de incremento se
puede escribir en la forma general como:

φ = a1k1 + a2 k 2 + L + a n k n 90

en donde las a son constantes y las k son:

k 1 = f ( xi , y1 )
k 2 = f ( x i + p1h, yi + q11k 1h)
k 3 = f ( xi + p2 h, y i + q 21k1h + q22 k 2 h)
M
k n = f ( xi + p n− h, yi + qn −1,1k1h + qn −1, 2 k 2 h + L + q n−1, n−1k n−1h )

Obsérvese que las k son relaciones recurrentes. Esto es k1 aparece en la

ecuación k 2 , que aparece en la ecuación de k 3 , etc. Esta recurrencia hace de
los métodos de RK eficientes para su cálculo en computadora.
Nótese, que el método de RK de primer orden con n = 1 , es de hecho. El
método de Euler.

Método de Runge-Kutta de segundo orden

La versión de segundo orden de la ecuación 89 es:

y i+1 = yi + (a1k1 + a 2 k 2 ) h 91

donde

k 1 = f ( xi , y1 )
k 2 = f ( x i + p1h, yi + q11k 1h)

los valores de a1 , a 2 , a1 y q11 se evalúan igualando la ecuación 91 a la

expansión de la serie de Taylor hasta el segundo término. Haciendo esto, se
obtienen tres ecuaciones para evaluar las cuatro incógnitas constantes. Estas
tres ecuaciones son:

a1 + a2 = 1 92
a 2 p1 = 1 / 2 93
a 2 q11 = 1 / 2 94

debido a que se tienen tres ecuaciones con cuatro incógnitas se debe suponer
el valor de una de las incógnitas para determinar las otras tres. Supóngase que

Pedro Ochoa Moreno Enero, 2002. 78

 Curso de Métodos Numéricos

se especifica el valor de a 2 . Entonces las ecuaciones 92 y 94 se resuelven

simultáneamente para:

a1 = 1 − a 2 95
1
p1 = q11 = 96
2 a2

ya que se puede escoger una cantidad infinita de valores de a 2 , existe un

número infinito de métodos de RK de segundo orden. Cada versión llevaría
exactamente al mismo resultado si la solución de la EDO es cuadrática, lineal o
constante.

Método de Runge-Kutta de tercer orden

Se puede llevar a cabo una derivación análoga a la del método de segundo

orden, para n = 3 . El resultado de esta derivación es de seis ecuaciones con
ocho incógnitas. Por lo tanto, se deben de dar a priori los valores de dos
incógnitas para determinar los parámetros restantes. Una versión común que
resulta es:

1 
y i+1 = y i +  ( k1 + 4k 2 + k 3 )  h 97
6 

donde

k 1 = f ( xi , y1 )
1 1
k 2 = f (x i + h, y i + k1 h)
2 2
k 3 = f ( xi + h, y i − k1h + 2k 2 h )
Obsérvese que si la derivada es una función solo de x , este método de tercer
orden se reduce a la regla de Simpson de 1/3.
Los métodos de RK de tercer orden tienen errores globales de 0( h 4 ) y 0( h 3 ) ,
respectivamente y llevan a resultados exactos cuando la solución es de orden
cúbico.

Ejemplo: Método de RK de tercer orden

Utilicese la ecuación 97 para integrar

a) Una EDO que es exclusivamente una función de x

dy
= − 2 x 3 + 12 x 2 − 20x + 8.5
dx

Pedro Ochoa Moreno Enero, 2002. 79

 Curso de Métodos Numéricos

con y (0) = 1 y con un tamaño de paso igual a 0.5.

Se pueden usar las ecuaciones que resultan de la ecuación 97 para calcular:

k 1 = −2(0) 3 + 12( 0) 2 − 20( 0) + 8.5 = 8.5

k 2 = −2(0.25) 3 + 12(0.25) 2 − 20(0.25) + 8.5 = 4.21875
k 3 = −2( 0.5) 3 + 12( 0.5) 2 − 20( 0.5) + 8.5 =1.25

que se pueden sustituir en la ecuación 97 para obtener:

1 
y (0.5) = 1 +  [8.5 + 4( 4.21875) + 1.25]( 0.5) = 3.21875
6 

b) Una EDO que es una función de x y y.

dy
= 4e 0.8 x − 0.5 y
dx

con y (0) = 2 desde x = 0 a 1 con un tamaño de paso de 1.

Se pueden usar las ecuaciones resultantes de la ecuación 97, para calcular:

k 1 = 4e 0.8( 0) − 0.5( 2) = 3
[ ]
k 2 = 4e 0.8( 0.5) − 0.5 2 + 0.5(1) 3 = 4.21729879
k 3 = 4e 0.8(1.0)
− 0.5[2 − 1(3) + 2(1)( 4.21729879) ] = 5.184864924

que se puede sustituir en la ecuación 97 y obtener:

1 
y (1.0) = 2 +  [3 + 4( 4.21729879) + 5.184864924](1) = 6.17567681
6 

que representa un Et = 0.31%, (el valor verdadero es 6.194).

Pedro Ochoa Moreno Enero, 2002. 80