Calculo Numerico

Universidad Nacional de Santiago del Estero
Facultad de Agronomía y Agroindustrias
Departamento Físico-Matemático
CALCULO NUMERICO
Curso práctico con aplicaciones

a la Ingeniería en Alimentos
Dra. Lucrecia Lucía Chaillou
2008
A mis adorados Abuelos y Padres
con todo mi amor, respeto y admiración
2
Prólogo
Los métodos numéricos constituyen una herramienta muy valiosa para la resolución de
problemas prácticos de Ingeniería, por ello el objetivo de este libro es presentarlos de manera
práctica y sintética.
Los distintos capítulos se diseñaron de acuerdo con las exigencias requeridas para la en-
señanza de cálculo numérico, asignatura de 3º año de la carrera de Ingeniería en Alimentos, des-
de el punto de vista práctico, considerando que los temas desarrollados servirán como base para
estudios más profundos.
Los temas tratados incluyen aspectos teóricos y prácticos sobre modelos y algoritmos;
aproximaciones y errores; solución numérica de ecuaciones; sistemas de ecuaciones lineales;
aproximación polinomial y funcional; simulación; series de Fourier; transformada de Laplace y
ecuaciones diferenciales. Se presentan gráficos aclaratorios, algoritmos de cada método numéri-
co, así como también ejercicios resueltos y propuestos aplicados a la Ingeniería de Alimentos.

Cátedra de Cálculo Numérico
Facultad de Agronomía y Agroindustrias
Universidad Nacional de Santiago del Estero
3
Índice
Pág.
Dedicatoria 2
Prólogo 3
Índice 4-6
Capítulo 1: METODOS NUMERICOS, MODELOS Y ALGORITMOS 7-14
1.1. METODOS NUMÉRICOS 7
1.2. MODELOS MATEMÁTICOS 8
1.2.1. Clasificación de modelos matemáticos 9
1.3. ALGORITMOS 10
1.4. RESOLUCIÓN NUMÉRICA DE UN PROBLEMA REAL 12
EJERCICIOS PROPUESTOS 14
Capítulo 2: APROXIMACIONES Y ERRORES 15-20
2.1. CIFRAS SIGNIFICATIVAS 15
2.2. EXACTITUD Y PRECISIÓN 16
2.3. ERRORES 16-19
2.3.1. Error absoluto y relativo 16
2.3.2. Errores en la resolución numérica 16
2.3.2.1. Error de truncamiento 17
2.3.2.2. Error de redondeo 17
2.3.2.3. Otros tipos de error 19
Capítulo 3: SOLUCION NUMERICA DE ECUACIONES 21-40
3.1. METODO GRAFICO 21
3.2. METODOS NUMERICOS DE CALCULO DE UNA RAIZ 22-35
3.2.1. Métodos cerrados 22-26
3.2.1.1. Método de la Bisección 22
3.2.1.2. Método de la Falsa Posición o Regula Falsi 24
3.2.2. Métodos abiertos 26-35
3.2.2.1. Método de Aproximaciones sucesivas 26
3.2.2.2. Método de Newton-Raphson o de la Tangente 29
3.2.2.3. Método de Newton de segundo orden 31
3.2.2.4. Método de Von Mises 32
3.2.2.5. Método de la secante 34
3.3. RAÍCES DE POLINOMIOS 35-39
3.3.1. Teoremas fundamentales de la Teoría de ecuaciones algebraicas 35
3.3.2. División sintética 36
3.3.3. Regla de los signos de Descartes 37
3.3.4. Raíces racionales 37
3.3.5. Raíces irracionales 37
3.3.5.1. Método de Newton-Raphson 38
Capítulo 4: SISTEMAS DE ECUACIONES LINEALES 41-55
4.1. CONCEPTOS PREVIOS 41-42
4.2. MÉTODOS DE RESOLUCIÓN DE SISTEMAS DE ECUACIONES ALGEBRAI- 42-54
CAS LINELES
4.2.1. Métodos directos 42-48
4.2.1.1. Método de eliminación de Gauss 42
4.2.1.2. Método de Gauss - Jordan 46
4.2.1.3. Partición de matrices 47
4.2.2. Métodos iterativos 48
4.2.2.1. Método de Jacobi 49
4.2.2.2. Método de Gauss–Seidel 52
4
4.2.2.3. Método de Relajación 54
Capítulo 5: APROXIMACION POLINOMIAL Y FUNCIONAL 56-80
5.1. APROXIMACIÓN POLINOMIAL 56-71
5.1.1. Diferencias finitas 56
5.1.2. Diferencias divididas 58
5.1.3. Interpolación con incrementos constantes. Interpolación de Newton 59
5.1.4. Interpolación con incrementos variables. Interpolación de Lagrange 61
5.1.5. Interpolación inversa 63
5.1.6. Derivación numérica 63
5.1.7. Integración numérica 66-71
5.1.7.1. Regla trapecial 67
5.1.7.2. Regla de Simpson 1/3 69
5.1.7.3. Regla de Simpson 3/8 69
5.2. APROXIMACIÓN FUNCIONAL 71-77
5.2.1. Regresión lineal 71
5.2.2. Linealización de relaciones no lineales 73
5.2.3. Regresión polinomial 74
5.2.4. Regresión lineal múltiple 75
Capítulo 6: SIMULACION 81-93
6.1. METODOLOGÍA DE SIMULACIÓN 81
6.1.1. Métodos de Montecarlo 82-89
6.1.1.1. Generación de números aleatorios 82
6.1.1.2. Generación de Números pseudo aleatorios 83
6.1.1.2.1. Método de los cuadrados centrales 83
7.1.1.2.2. Método de los productos centrales 83
6.1.1.2.3. Métodos congruenciales 83
6.1.1.3. Aplicaciones de los métodos de Montecarlo 85
6.1.1.3.1. Paseo aleatorio 85
6.1.1.3.2. Integración por simulación 87
6.1.1.3.2. Línea de espera 88
6.2. Modelos Demográficos y de la Cinética Química 89-93
6.2.1. Modelo demográfico 89
6.2.2. Modelos de la cinética química 90
6.2.2.1. Método diferencial 91
6.2.2.2. Método integral 92
6.2.2.3. Dinámica de sistemas cinetoquímicos 92
Capítulo 7: SERIES DE FOURIER 94-104
7.1. CONSIDERACIONES PREVIAS 94-97
7.1.1. Funciones periódicas 94
7.1.2. Series trigonométricas 95
7.1.3. Funciones seccionalmente continuas 95
7.1.4. Funciones pares e impares 96
7.2. DESARROLLO EN SERIE DE FOURIER 97-101
7.2.1. Cálculo de los coeficientes de Fourier 98
7.2.2. Expresión de la serie de Fourier para funciones de período arbitrario 99
7.2.3. Forma exponencial de la serie de Fourier 100
7.2.4. Consideraciones simplificatorias 100
7.2.5. Espectro de frecuencias 101
7.3. INTEGRALES DE FOURIER 102
Capítulo 8: TRANSFORMADA DE LAPLACE 105-111
8.1. DEFINICIÓN DE TRASFORMADA DE LAPLACE 105
5
8.2. PROPIEDADES IMPORTANTES DE LA TRANSFORMADA DE LAPLACE 106
8.2.1. Transformada de Laplace de operaciones 106
8.3. MÉTODOS PARA CALCULAR TRANSFORMADAS DE LAPLACE 107
8.4. VENTAJAS DEL MÉTODO DE LA TRANSFORMADA DE LAPLACE 107
8.5. TRANSFORMADA DE LAPLACE DE ALGUNAS FUNCIONES IMPORTANTES 108
8.6. TRANSFORMADA INVERSA DE LAPLACE 109
8.7. INTEGRAL DE CONVOLUCION 109
Capítulo 9: ECUACIONES DIFERENCIALES 112-131
9.1. CONCEPTOS PREVIOS 112-113
9.2. ECUACIONES DIFERENCIALES ORDINARIAS 113-125
9.2.1. Solución analítica 114
9.2.1.1. Solución por el método del operador diferencial 114
9.2.2. Solución por métodos numéricos 116
9.2.2.1. Métodos de un paso 116
9.2.2.1.1. Método de la Serie de Taylor 116
9.2.2.1.2. Método de Euler 117
9.2.2.1.3 Métodos de Runge-Kutta 119
9.2.2.1.4 Métodos de Heun 123
9.2.2.2. Métodos de pasos múltiples 125
9.3. ECUACIÓN DIFERENCIAL LINEAL DE SEGUNDO ORDEN 125-129
9.3.1. Métodos de resolución 125
9.3.1.1. Método del operador diferencial 125
9.3.1. 1.1. Solución complementaria 126
9.3.1.1.1.1. Caso sobre-amortiguado 126
9.3.1.1.1.2. Caso crítico 127
9.3.1.1.1.3. Caso oscilatorio amortiguado 128
9.3.1.1.2. Solución particular 128
9.1.1.3. Solución general 129
9.3.1.2. Método de los coeficientes indeterminados 129
9.4. ECUACIONES DIFERENCIALES PARCIALES 131-132
9.4.1. Ecuaciones elípticas 132
9.4.2. Ecuación Parabólica 132
9.4.3. Ecuación Hiperbólica 132
BIBLIOGRAFÍA CONSULTADA 134
6
Cálculo Numérico
Capítulo 1
METODOS NUMERICOS, MODELOS Y ALGORITMOS
La resolución de problemas de Ingeniería está asociada, por lo general, a resultados

numéricos puesto que se requieren respuestas prácticas. Muchos de estos problemas sólo se
pueden resolver de forma aproximada, por ello es importante el estudio de una rama de las Ma-
temáticas denominada Análisis Numérico, esta rama involucra el estudio de Métodos Numéricos.
Su desarrollo estuvo y está notablemente influenciado y determinado por las computadoras digita-
les que permiten realizar los cálculos de manera veloz, confiable y flexible.
1.1. METODOS NUMÉRICOS
Se pueden definir a los métodos numéricos como las técnicas mediante las cuales es
posible formular problemas de manera que puedan resolverse utilizando operaciones aritméti-
cas, ó también como el grupo de conocimientos matemáticos relacionados con el diseño y análi-
sis de algoritmos necesarios para resolver problemas de ciencia e ingeniería. Estos métodos se
caracterizan porque: permiten dar más importancia a la formulación e interpretación de la solución,
los cálculos involucrados están relacionados con cantidades discretas, permiten obtener resulta-
dos aproximados y ayudan a identificar, cuantificar y minimizar los errores.
Existen varios motivos por los cuales deben estudiarse estos métodos:
1. Son herramientas poderosas para la solución de problemas. Permiten manejar sis-

temas de ecuaciones grandes, no linealidades y geometrías complicadas.
2. Su teoría es la base de programas de métodos numéricos.
3. Su conocimiento permite diseñar programas propios.
4. Son un vehículo eficiente para aprender a servirse de las computadoras persona-

les.
5. Son un medio para reforzar la comprensión de las matemáticas.
Estos métodos permiten:
9 Encontrar las raíces de ecuaciones lineales y no lineales
9 Resolver grandes sistemas de ecuaciones algebraicas lineales
9 Encontrar aproximaciones de funciones
7
Cálculo Numérico
9 Realizar interpolaciones para encontrar valores intermedios en tablas de datos
9 Aproximar derivadas de cualquier orden
9 Integrar cualquier función
9 Resolver problemas de valor inicial y de frontera
9 Obtener soluciones numéricas para ecuaciones diferenciales parciales
9 Realizar ajustamiento de curvas a datos
1.2. MODELOS MATEMÁTICOS
El mundo real es naturalmente complejo y en muchas ocasiones los problemas a resolver

resultan difíciles de sintetizar. Para identificar sus aspectos esenciales y expresarlos en términos
precisos se debe realizar un proceso de abstracción. Esa abstracción del problema del mundo real
simplificando su expresión se denomina modelización. La modelización es una de las áreas más
atractivas de la ingeniería y las ciencias aplicadas. De hecho, los ingenieros necesitan construir
modelos para resolver problemas de la vida real.
El objetivo de un modelo consiste en reproducir la realidad de la forma más fiel posible, tra-
tando de entender cómo se comporta el mundo real y obteniendo las respuestas que pueden es-
perarse de determinadas acciones. Su selección es una etapa crucial para obtener una solución
satisfactoria a un problema real. Las estructuras matemáticas asociadas no son arbitrarias, sino
una consecuencia de la realidad misma.
Un modelo mental puede definirse como una representación de la realidad en la que se

consideran los aspectos más relevantes, es decir el modelo representa una parte de la realidad.
Un modelo matemático es una formulación o ecuación que expresa las características

fundamentales de un sistema o proceso físico en términos matemáticos. Los modelos pueden
estar constituidos por simples ecuaciones algebraicas hasta grandes y complicados sistemas de
ecuaciones diferenciales. Sus características son:
1. Describe un sistema o proceso natural en términos matemáticos.
2. Representa una idealización y una simplificación de la realidad.
3. Conduce a resultados predecibles y, en consecuencia, puede utilizarse para

propósitos de predicción.
La generación de un modelo matemático involucra dos etapas fundamentales, la de con-

ceptualización y la de formulación. En la primera se debe caracterizar el contexto del problema
real, definir claramente el propósito y los límites del modelo, identificar y establecer relaciones
8
Cálculo Numérico
entre las variables; en la segunda etapa se deben determinar las ecuaciones asociadas al modelo
y seleccionar y estimar los parámetros del modelo.
El objetivo del modelo es aplicarlo para obtener alguna información del problema o fenó-
meno que se estudia. Frecuentemente sufre modificaciones y a veces es descartado y aunque
contenga errores, puede poner en evidencia componentes esenciales de una realidad compleja.
1.2.1. Clasificación de modelos matemáticos
Los modelos matemáticos pueden clasificarse en función del tratamiento de la incertidum-

bre; del origen de la información; de su campo de aplicación, etc.
a) En función del tratamiento de la incertidumbre
Determinístico: se conoce de manera puntual la forma del resultado ya que no hay incertidumbre.
Además, los datos utilizados para alimentar el modelo son completamente conocidos y determina-
dos.
Estocástico: probabilístico, no se conoce el resultado esperado, sino su probabilidad y existe por

lo tanto incertidumbre.
b) En función del origen de la información utilizada para construirlos
Modelos heurísticos: del griego euriskein, hallar, inventar. Son los que están basados en las
explicaciones sobre las causas o mecanismos naturales que dan lugar al fenómeno estudiado.
Modelos empíricos: del griego empeiricos (experiencia, experimento) Son los que utilizan las
observaciones directas o los resultados de experimentos del fenómeno estudiado.
c) En función de su campo de aplicación
Modelos conceptuales: son los que reproducen mediante fórmulas y algoritmos matemáticos
más o menos complejos los procesos físicos que se producen en la naturaleza.
Modelo matemático de optimización: los modelos matemáticos de optimización son ampliamen-

te utilizados en diversas ramas de la ingeniería para resolver problemas que por su naturaleza son
indeterminados, es decir presentan más de una solución posible.
d) En función del factor tiempo
Modelos estáticos: son independientes del tiempo, consideran situaciones estacionarias.
9
Cálculo Numérico
Modelos dinámicos: son los que describen el comportamiento del sistema en estudio en función
del tiempo.
1.3. ALGORITMOS
Un algoritmo puede definirse como una secuencia lógica de pasos necesarios para la eje-
cución de una tarea específica, tal como la solución de un problema, ó también como una secuen-
cia de instrucciones para alcanzar un resultado deseado en un tiempo finito.
Un buen algoritmo se caracteriza por: terminar luego de una cantidad finita de pasos, ser lo
más general y preciso posible, ser determinístico, no dejar nada al azar y permitir obtener resulta-
dos independientes de quien lo está utilizando.
Para generar un algoritmo se debe seguir una serie de pasos:
1. Determinar el objetivo de la tarea
2. Identificar los datos de entrada y de salida
3. Determinar el proceso involucrado
4. Identificar las variables internas
5. Dividir el proceso en acciones elementales
6. Determinar la secuencia de estas acciones
7. Incorporar estructuras de control
Por lo general, el objetivo del algoritmo será el de implementar un procedimiento numérico

para resolver un problema o para aproximar una solución del problema. Consta de un principio;
de una serie de pasos en los que se deben definir los valores iniciales de las variables del proble-
ma, operar con estos valores hasta llegar a un resultado, proporcionar un resultado y de un final.
Un algoritmo se puede representar mediante un pseudocódigo que especifica los datos

de entrada, la forma de los resultados deseados y los pasos involucrados ó bien mediante un dia-
grama de flujo que es una representación visual o gráfica del algoritmo que emplea una serie de
bloques y flechas. Cada bloque representa una operación particular o un paso en el algoritmo. Las
flechas indican la secuencia en que se implementan las operaciones.
Los símbolos que se utilizan en diagramas de flujo se representan en la Figura 1.1.
10
Cálculo Numérico
Bloque Terminal: indica el inicio o finalización del algoritmo
Bloque de Proceso: representa los cálculos o manipulación de datos
Bloque de Entrada/Salida: indica la entrada o salida de datos e informa-

ción
Bloque de Decisión: representa una comparación o pregunta que deter-

mina alternativas diferentes a seguir
Bloque de Iteración: representa cálculos repetitivos
Conector: indica un truncamiento en el camino del diagrama de flujo

cuando el diagrama es grande y no cabe en una página.
Figura 1.1. Símbolos que se utilizan en diagramas de flujo
¾ Por ejemplo, si se deseara escribir el algoritmo para la solución de un problema simple tal co-
mo, dados dos números x1 y x2, escribir el mayor de ellos, se presentan a continuación el
pseudocódigo y el diagrama de flujo (Figura 1.2) correspondientes:
Algoritmo Mayor (algoritmo que muestra el mayor

Inicio
de dos números)
Las variables de entrada son: x1, x2

Introducir
Paso 1: Introducir x1, x2
x1 y x2
Paso 2: si x1 > x2 entonces escribir x1
Paso 3: si no escribir x2
NO SI
x1 > x2
PARAR
Imprimir Imprimir
x2 x1
FIN
Figura 1.2. Diagrama de flujo del algoritmo Mayor
11
Cálculo Numérico
1.4. RESOLUCIÓN NUMÉRICA DE UN PROBLEMA REAL
Consideremos el problema físico de una encomienda que se deja caer desde un globo ae-
rostático. Se desea determinar la velocidad de caída luego de 12 segundos si la masa del cuerpo
es 70 kg y que el coeficiente de roce es 0,27 kg/m.
En la Figura 1.3 puede observarse un esquema de la situación planteada.
t=0
Fγ y(t)
t=t
Figura 1.3. Representación de las fuerzas que actúan sobre un cuerpo

en caída libre.
El modelo físico que lo rige está dado por la segunda Ley de Newton. Su expresión ma-
temática es
F= M a (1.1)
Este modelo es una idealización y simplificación de la realidad, no incluye los efectos de la

relatividad, donde F corresponde a la fuerza neta que actúa sobre el cuerpo, M es la masa del
objeto y a su aceleración.
Para un cuerpo que cae dentro del perímetro de la Tierra, la fuerza neta está formada por
dos fuerzas opuestas: la atracción hacia abajo debida a la gravedad P y la fuerza hacia arriba de-
bida a la resistencia del aire Fγ, si consideramos el sistema de referencia positivo hacia abajo, esta
última tendrá signo negativo. La resistencia ofrecida por el aire puede expresarse de varias mane-
ras, una aproximación sencilla es suponer que:
Fγ = γ v (1.2)
Además la aceleración puede expresarse como la razón de cambio de la velocidad con

respecto al tiempo (dv/dt), por lo tanto la ecuación (1.1) puede escribirse como
dv
M = Mg − γv (1.3)
dt
Si se divide por M a ambos miembros para normalizarla se llega a la ecuación
dv γ
= g− v (1.4)
dt M
12
Cálculo Numérico
La ecuación (1.4) es un modelo matemático y es una ecuación diferencial puesto que está
escrita en términos de la razón de cambio diferencial (dv/dt). Para resolverla pueden utilizarse dos
métodos: el analítico que aplica las reglas del cálculo diferencial o bien el método numérico.
a) Solución analítica
La solución analítica exacta de la ecuación (1.4) es:
⎛ γ
gM ⎜ − t ⎞⎟
v( t ) = 1− e M (1.5)
γ ⎜⎜ ⎟⎟
⎝ ⎠
Al sustituir los valores de los parámetros en la ecuación (1.5) se obtiene
⎛ 0,27 ⎞
9,8(m / s 2 )70(kg) ⎜ − t
70 ⎟
v( t ) = 1 − e
0,27(kg / m) ⎜⎜ ⎟⎟
⎝ ⎠
b) Solución numérica
Como se mencionó anteriormente, los métodos numéricos permiten reformular el problema

para que se pueda resolver mediante simples operaciones aritméticas. Entonces se aproxima la
razón de cambio de la velocidad con respecto al tiempo por:
dv ∆v v( ti+1) − v( ti )
≈ = (1.6)
dt ∆t ti+1 − ti
Reemplazando en (1.4)
v( ti+1) − v( ti ) γ
= g − v( t i ) (1.7)
ti+1 − ti M
Esta ecuación puede reordenarse para obtener la velocidad en el instante ti+1
⎡ γ ⎤
v( ti+1) = v( ti ) + ⎢g − v( ti )⎥ ( ti+1 − ti ) (1.8)
⎣ M ⎦
De manera que la ecuación diferencial se transforma en una ecuación algebraica, en la

⎡ γ ⎤
que se puede calcular v (ti+1) si se da un valor inicial a v (ti), y donde ⎢g − v( ti )⎥ es la pendiente
⎣ M ⎦
de la recta descripta por esta ecuación, es decir:
Valor de v ) = (Valor de v ) + (Valor de la ) x (Incremento de )

(nuevo anterior pendiente tiempo
Como en el instante inicial la velocidad del cuerpo es 0, se toma éste para calcular la velo-
cidad en t=2 s y así sucesivamente. En la Tabla 1.4 se muestran los valores de velocidad obteni-
dos para la solución analítica y la solución numérica.
13
Cálculo Numérico
Tabla 1.4. Solución analítica y numérica

al problema de un cuerpo que cae
Solución Solución
t (s) analítica numérica
v (m/s)
0 0 0
2 19,5246 19,6000
4 38,8991 39,0488
6 58,1248 58,3476
8 77,2027 77,4975
10 96,1341 96,4996
12 114,9199 115,3552
Puede observarse que por un método numérico la solución se aproxima bastante bien a
solución exacta. Para minimizar las diferencias se puede utilizar un menor intervalo de cálculo, por
ejemplo intervalos de 1 s. Con la ayuda de una computadora digital se puede efectuar un gran
número de cálculos en pocos segundos y modelar con exactitud la velocidad de un cuerpo que
cae sin tener que resolver la ecuación diferencial.
EJERCICIOS PROPUESTOS
1.1. Ejemplifique el proceso de generación de un modelo matemático, a partir de un fenómeno o

problema del mundo real. Detalle los aspectos involucrados en la conceptualización y la formula-
ción del mismo.
1.2. Responda verdadero o falso:
a) Un algoritmo debe ser finito y preciso
b) Los métodos numéricos son aquellos en los que se reformula un problema matemático pa-
ra que pueda resolverse mediante operaciones aritméticas.
c) Un modelo matemático nunca puede ser modificado
1.3. Enuncie las características relevantes del Cálculo Numérico e indique por lo menos cinco pro-
blemas matemáticos que surgen en Ingeniería y pueden resolverse por Métodos Numéricos.
1.4. Calcule analítica y numéricamente la velocidad de caída, a los 20 s, de un cuerpo de 50 kg

que se deja caer desde un aeroplano, considerando que la fuerza de roce es γ v2 (γ=0,27
kg/m).Grafique ambas soluciones.
1.5. Escriba el algoritmo y el diagrama de flujo correspondiente al problema de la multiplicación de

dos números.
14
Cálculo Numérico
Capítulo 2
APROXIMACIONES Y ERRORES
El análisis del error en un resultado numérico es esencial en cualquier cálculo hecho a ma-
no o con una computadora. Los datos de entrada rara vez son exactos puesto que se basan en
ensayos experimentales o bien son estimados y los métodos numéricos introducen errores de
varios tipos, por ello brindan resultados aproximados. En la práctica profesional, los errores son
costosos y en algunos casos letales. Además como los resultados de los métodos numéricos son
aproximaciones, es necesario tener en claro los conceptos de cifras significativas, exactitud y pre-
cisión.
2.1. CIFRAS SIGNIFICATIVAS
La confiabilidad de un valor numérico está dada por sus cifras significativas que se defi-
nen como el número de dígitos, más un dígito estimado que se pueda usar con confianza. Por
ejemplo, si se leen 25 ml en una bureta, que está graduada en 0,1 ml, se puede decir que el nivel
del líquido es mayor que 25,1 y menor que 25,2 ml como puede observarse en la Figura 2.1. Has-
ta puede estimarse con una aproximación de ± 0,05 ml, por lo tanto el volumen vertido es 25,15 ml
que tiene 4 cifras significativas. Los primeros tres dígitos son seguros y el último es una estima-
ción.
25
Figura 2.1. Representación de la sección de una bureta
Un cero puede ser significativo o no, dependiendo de su posición en un número dado. Los
ceros que solamente sitúan la cifra decimal no son significativos, si se escribiera 25,15 ml como
0,02515 l, el número de cifras significativas sigue siendo el mismo. Los ceros al final de un núme-
ro pueden ser significativos o no. Si se dice que un tanque de agua se encuentra a 10,0 m de altu-
ra, significa que la altura se conoce hasta las décimas de metro. Si esa misma altura se da como
15
Cálculo Numérico
1000 cm la expresión es confusa y para mantener el criterio de cifras significativas se utiliza la

notación científica 1,0 x 103.
2.2. EXACTITUD Y PRECISIÓN
La precisión se refiere al número de cifras significativas que representan una cantidad ó a

la extensión en las lecturas repetidas de un instrumento que mide alguna propiedad física.
La exactitud se refiere a la aproximación de un número o de una medida al valor verdade-

ro que se supone representa.
2.3. ERRORES
2.3.1. Error absoluto y relativo
El error se aplica para indicar la inexactitud y la imprecisión de las mediciones.
El error numérico es igual a la diferencia entre el valor verdadero y el aproximado:
E v = valor verdadero – valor aproximado (2.1)
El error relativo fraccional resulta de normalizar el error respecto al valor verdadero:
error
Error relativo fraccional = (2.2)
valor verdadero
Si se expresa en porcentaje:
error verdadero
Ev = x 100 (2.3)
valor verdadero
El error relativo porcentual de aproximación está dado por:
aproximación actual − aproximación previa

Ea = x 100 (2.4)
aproximación actual
Por lo general se toma el valor absoluto del error.
De acuerdo con Scarborough, se tiene la seguridad de que el resultado es correcto en al

menos n cifras significativas si se cumple el siguiente criterio:
∈s = (0,5 x 10 2 − n )% (2.5)
donde ∈s es la tolerancia prefijada.
16
Cálculo Numérico
2.3.2. Errores en la resolución numérica
Las soluciones que resultan de la aplicación de los métodos numéricos son aproximadas
debido a que existen incertidumbres en los datos, puesto que éstos son empíricos; en el modelo
ya que es una idealización y simplificación de la realidad y en la resolución numérica debido a
errores de truncamiento y de redondeo.
2.3.2.1. Error de truncamiento
Estos errores resultan al usar una aproximación en lugar de un procedimiento matemático

exacto y dependen del método numérico empleado.
La serie de Taylor (2.6) puede utilizarse para predecir el valor de la función en x i+1 en
términos de la función y de sus derivadas en la vecindad de un punto xi.
f ' ' ( xi ) f (n)( xi )

f ( xi+1) = f ( xi ) + f ' ( xi )( xi+1 − xi ) + ( xi+1 − xi )2 + L + ( xi+1 − xi )n + Rn (2.6)
2! n!
El término residual considera todos los términos desde n+1 hasta el infinito.
f (n+1) (ξ)
Rn = ( xi+1 − xi )n+1, con ( xi+1 − xi ) = h (paso) (2.7)
(n + 1) !
Si se considera el primer término de la serie, la aproximación es de orden cero; con dos

términos la aproximación es de 1º orden, tres términos corresponden a 2º orden y así sucesiva-
mente. Cuanto mayor sea el número de términos incluidos, menor será el error de truncamiento.
2.3.2.2. Error de redondeo
Estos errores resultan de representar en forma aproximada números exactos, dependen

de la computadora o de quien realice los cálculos. Si se realizan las operaciones algebraicas a
mano se deben tener en cuenta las reglas de redondeo.
Reglas de redondeo
1. Se conservan las cifras significativas y el resto se descarta. El último dígito que se

conserva se aumenta en 1si el primer dígito descartado es mayor que 5. De otra
manera se deja igual. Si el primer dígito descartado es 5 o es 5 seguido de ceros,
entonces el último dígito retenido se incrementa en 1 solo si es impar.
2. En la suma y la resta el redondeo se lleva a cabo de manera que el último dígito re-
tenido en la respuesta corresponda al último dígito más significativo de los números
que están sumando o restando. Un dígito en la columna de las centésimas es más
significativo que en las milésimas.
17
Cálculo Numérico
3. Para la multiplicación y la división el redondeo es tal que la cantidad de cifras del

resultado es igual al número más pequeño de cifras significativas que contienen las
cantidades en operación.
4. En el caso de operaciones combinadas se ejecutan las operaciones entre parénte-

sis y el resultado se redondea antes de proceder con la otra operación.
Si los cálculos se realizan utilizando una computadora se debe considerar que la mayoría
de ellas guardan un número finito de cifras significativas durante un cálculo y resultan críticos en
dos casos:
1. Ciertos métodos requieren cantidades extremadamente grandes para obtener una res-
puesta, además los cálculos dependen entre sí, por lo tanto los cálculos posteriores son
dependientes de los anteriores y por lo tanto el efecto de acumulación en el transcurso de
una gran cantidad de cálculos resulta significativo.
2. Si se realizan operaciones algebraicas con números muy pequeños y muy grandes al

mismo tiempo.
Además la mayoría de las computadoras representan a los números como números de

punto flotante. La representación de punto flotante de un número está dada por la siguiente ex-
presión:
fl ( x ) = ε 0. a1 a2 a3 L apBb (2.8)
Donde:
ε: es el signo del número, puede ser positivo o negativo
a1a2a3...ap: es la parte fraccionaria significativa
B: es la base, puede ser 2, 10 ó 16
b: es el exponente entero, las computadoras de 12 dígitos tiene un valor de b de ± 999
p: es el número de dígitos significativos (precisión)
Por ejemplo, si se quiere representar el número 24,12 como un número de punto flotante
con B=10 y p=4, este será:
fl(x)= +.2412x102
Si la computadora admite solo p cifras significativas el redondeo se hace al número más

próximo. Dado el número:
x = ε 0. a1 a2 a3 Lapap +1 ap + 2 L10b
El redondeo para este número utilizando el punto flotante fl(x) es:
18
Cálculo Numérico
fl ( x ) = ε 0. a1 a2 a3 Lap10b si ap +1 < 5
fl ( x ) = ε 0. a1 a2 a3 L(ap + 1) 10b si ap +1 ≥ 5
Las cotas del error de redondeo serán:
a) Error absoluto Ea= fl(x)-x con su valor absoluto Ea ≤ 0.5.10b −p

E
b) Error relativo Er= a con su valor absoluto Er ≤ 0.5.10 −p
x
2.3.2.3. Otros tipos de error

Errores por equivocación
Son errores por torpeza o por equivocación, son debidos por lo general a errores humanos.
Las equivocaciones ocurren en cualquier etapa del proceso de modelado matemático y pueden
contribuir con las otras componentes del error.
Errores de formulación
Estos errores degeneran en un modelo matemático incompleto y si se está usando un mo-

delo deficiente, ningún método numérico generará resultados adecuados.
Incertidumbre en los datos
Algunas veces se introducen errores en un análisis debido a la incertidumbre de los datos

físicos sobre los que se basa el modelo. Son errores que muestran inexactitud e imprecisión.
2.1. Suponga que debe cuantificar la cantidad de β-caroteno en lechuga y experimentalmente se

determinó que el valor es 0,042 mg/100g. Si el valor verdadero es 0,048 mg/100g, indique el error
verdadero y el error relativo porcentual.
2.2. Estime el valor de e0.5 utilizando la expansión en serie de Mac Laurin, calculando los errores
relativos porcentuales real y aproximado (considere que el valor real de e0.5 es1, 648721271 des-
pués del agregado de cada término hasta que el valor absoluto del error aproximado sea menor
que el criterio establecido por la fórmula de Scarborough para 3 cifras significativas.
2.3. En la tabla que sigue se muestran las velocidades de formación del compuesto C, mediante
una reacción enzimática, a partir de los reactivos A y B. Se indican las velocidades de formación
con 3, 4, 5 y 6 cifras significativas. Calcule los errores relativos porcentuales para un tiempo t =
12s, considerando que el valor real con 10 cifras significativas es 4984,921508 µg/s.
19
Cálculo Numérico
Velocidad de formación (µg/s)

Tiempo (s)
Cifras significativas
3 4 5 6
4 3200 3200 3200,4 3200,46
8 4489 4489 4489,3 44892,41
12 4980 4984 4984,8 4984,91
2.4. Aplique las reglas de redondeo a:
a) 5,6723 a 3 cifras significativas
b) 10,406 a 4 cifras significativas
c) 7,3500 a 2 cifras significativas
2.5. Evalúe:
2,2 – 1,768 y 0,0642 x 4,80
2.6. Utilice términos en la Serie de Taylor de cero a cuarto orden para aproximar la función f(x)= -
0,2 x4 - 0,35 x3-0,5 x2-0,45 x+1,8 para x= 2 y calcule el error de truncamiento en cada caso. Consi-
dere h=1.
2.7. Represente las siguientes cantidades como números de punto flotante. Considerando base
10 y 4 dígitos significativos.
a) 28,31; b) -0,00144; c) 38000
20
Cálculo Numérico
Capítulo 3
SOLUCION NUMERICA DE ECUACIONES
Uno de los problemas más antiguos y básicos del cálculo numérico es el problema de
búsqueda de la solución de una ecuación, es decir encontrar los valores de la variable x que satis-
facen la ecuación f(x)=0, para una función f dada. Las ecuaciones pueden ser algebraicas (la fun-
ción f es un polinomio), por ejemplo: x2+5x-4=0 o bien trascendentes puesto que están constitui-
das por funciones trascendentes tales como funciones exponenciales, trigonométricas, logarítmi-
cas, etc., por ejemplo: e-x – x; sen x; ln x2 – 1. Solamente en casos muy simples, de ecuaciones
algebraicas, existen fórmulas que permiten resolverlas en términos de sus coeficientes, para el
resto de las ecuaciones se utilizan métodos aproximados que permiten mejorar la solución por
simple repetición del mismo método hasta adquirir el grado de aproximación requerido. Estos
métodos son apropiados para realizarlos utilizando computadoras puesto que comprenden la re-
petición de un proceso, es decir iteración. A continuación se describen métodos numéricos que
permiten calcular las raíces de ecuaciones algebraicas y trascendentes.
3.1. METODO GRAFICO
Este es un método muy simple, consiste en calcular valores de la variable dependiente pa-
ra distintos valores de la variable independiente. A continuación se grafican en un sistema de ejes
coordenados cartesianos y se observa el punto de intersección de la función con el eje de las
abscisas. Este punto proporciona una primera aproximación a la raíz de la ecuación.
¾ Por ejemplo, si se desea determinar, aplicando el método gráfico, los valores aproximados de
las raíces de x2-6 x +1 =0. Para ello se calcula el valor de la función en el intervalo [-2,8] y se
representan los valores en un sistema de ejes cartesianos (Figura 3.1).
y 20
15
10
-5
-10
-2 0 2 4 6 8 x
Figura 3.1. Gráfico del polinomio x2-6 x +1

Se observa en el gráfico que las raíces aproximadas son 0 y 6.
21
Cálculo Numérico
3.2. METODOS NUMERICOS DE CÁLCULO DE UNA RAIZ
Los métodos de cómputo de una raíz real de una ecuación involucra dos pasos, en primer
lugar la determinación del intervalo de búsqueda, es decir el intervalo al que la raíz pertenece,
siempre que la ecuación esté vinculada a un sistema físico y en segundo lugar la selección y apli-
cación de un método numérico apropiado para determinar la raíz con la exactitud adecuada.
Estos métodos se clasifican en dos categorías: métodos abiertos y métodos cerrados. Los
métodos cerrados, tales como el método de la bisección y el de la falsa posición, son aquellos que
usan intervalos, se caracterizan por ser siempre convergentes pero la velocidad de convergencia
es lenta. Los métodos abiertos, método de aproximaciones sucesivas, de Newton-Raphson, de
Newton de 2º orden, de Von Mises, de la secante, requieren información únicamente de un punto,
o de dos pero que no necesariamente encierran a la raíz. La convergencia es más rápida pero
algunas veces divergen.
3.2.1. Métodos cerrados
3.2.1.1. Método de la Bisección
El método de la bisección, conocido también como de corte binario, de partición en dos

intervalos iguales, de búsqueda binaria o de Bolzano se basa en el Teorema del Valor Intermedio
y en el teorema de Bolzano.
Teorema del valor intermedio: Si f∈[a,b] y k es un número cualquiera comprendido entre f(a) y f (b) entonces
existe un c en el intervalo (a,b) tal que f(c)=k.
Teorema de Bolzano: sea f una función contínua en el intervalo [a,b], con f(a)f(b)<0 entonces existe al me-
nos un punto c∈[a,b] tal que f(c)=0
Si se tiene una función f(x) continua en el intervalo [xL,xU], con f(xL) y f(xU) de signos
opuestos, por el teorema anterior, existe un valor x* incluido en el intervalo (xL,xU) tal que f(x*)=0.
El método requiere de dividir el intervalo a la mitad y localizar la mitad que contiene a la raíz. El
proceso se repite y su aproximación mejora a medida que los subintervalos se dividen en interva-
(x L + x U )
los más y más pequeños. La primera aproximación a la raíz, se determina como xM = ,
2
ver Figura 3.2. y
xL (xM,0)
xU x
Figura 3.2. Esquema del método de la Bisección
22
Cálculo Numérico
Para determinar en que subintervalo está la raíz debe considerarse lo siguiente:
Si f (xM) = 0, entonces la raíz es igual a xM.

Si f (xL). f(xM) < 0, la raíz se encuentra en el primer subintervalo (xL, xM)
Si f(xL). f (xM) > 0, la raíz se encuentra en el segundo subintervalo (xM, xU).
Se calcula una nueva aproximación a la raíz en el nuevo subintervalo y se continúa con las ite-
raciones hasta la cota de error fijada de antemano.
Las ventajas y desventajas del método se detallan a continuación:
Ventajas Desventajas
9 Es siempre convergente - Converge muy lentamente
- Si existe más de una raíz en el interva-
lo, el método permite encontrar sólo una
de ellas
A continuación se presenta un algoritmo de este método iterativo.
Algoritmo para el método de Bisección

Permite encontrar una solución a la ecuación f(x)=0, dada la función continua f en el interva-
lo [xL,xU].
Considerando la siguiente notación:
xL: límite inferior del intervalo considerado
xU: límite superior del intervalo considerado
xM: raíz aproximada
E: cota de error o criterio de detención
N: número máximo de iteraciones
Paso 1: Tomar i = 1
Paso 2: Mientras i ≤ N seguir con los pasos 3 a 6
(x L + x U )
Paso 3: Tomar xM= xL+
2
xL + xU
Paso 4: Si < E ó f(xM)=0, SALIDA xM
2
PARAR
Paso 5: Tomar i = i+1
Paso 6: Si f(xL). f(xM) >0, tomar xL=xM,
si no tomar xU=xM
Paso 7: SALIDA (‘Procedimiento completado sin éxito después de N’ iteraciones)
PARAR
23
Cálculo Numérico
¾ Por ejemplo, si se desea determinar, aplicando el método de la bisección, una de las raíces de
la ecuación x3+x2-3x-3=0, considerando que la función cambia de signo en el intervalo (1,2).
La estimación inicial de la raíz se sitúa en el punto medio de este intervalo:
(x L + x U ) 1+ 2
XM = = = 1,5
2 2
Ahora se calcula f(xL). f(xM):
f(1) . f(1,5) = (-4) . (-1,875) = 7,5 > 0 no hay cambio de signo entre a y x1, entonces la raíz se en-
cuentra en el intervalo (1,5, 2).
1,5 + 2
La aproximación a la raíz en la segunda iteración se calcula como xM = = 1,75
2
y f(1,5) . f(1,75) = -1,70312< 0 , por lo tanto la raíz está en (1,5, 1,75), entonces la tercera itera-
ción es:
1,5 + 1,75
xM = = 1,625
2
y así sucesivamente, en la sexta iteración se llega a un valor de x6=1,734 bastante próximo al va-
lor verdadero de la raíz 1,7321
3.2.1.2. Método de la Falsa Posición o Regula Falsi
Este método es similar al de la bisección salvo que la siguiente iteración se toma en la in-
tersección de una recta entre el par de valores x y el eje de las abscisas en lugar de tomar el pun-
to medio. El reemplazo de 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 regla falsa.
Para aplicarlo se eligen los extremos xL y xU del intervalo entre los que se encuentra la raíz,
verificando que se cumpla que f(xL). f(xU) < 0. Si se observa la Figura 3.3, por semejanza de trián-
gulos, puede escribirse la siguiente igualdad:
f(xL ) f(xU )
= (3.1)
xM − xL xM − xU
Y despejando de la expresión (3.1) el valor de xM, que es una aproximación de la raíz, se obtiene
la siguiente fórmula de iteración o recurrencia:
xU − xL
xM = xU − f(xU ) (3.2)
f(xU ) − f(xL )
El valor de xM, calculado con la ecuación (3.2), reemplaza a uno de los dos valores, xL o xU
que produzca un valor de la función que tenga el mismo signo de f(xM). De esta manera los valo-
res xL y xU siempre encierran a la raíz.
Si f(xM)=0 el proceso termina.
Si f(xM) tiene el mismo signo de f(xL), el próximo paso es elegir xL = xM y xU = xU.
Si f(xM) tiene el mismo signo de f(xU) el próximo paso es elegir xL = xL y xU = xM.
24
Cálculo Numérico
El proceso se repite en la misma forma hasta llegar a la cota de error.

En la Figura 3.3 se presenta un esquema del método.
xL (xM,0)
xU x
Figura 3.3. Esquema del método de la Falsa Posición
Un algoritmo para este método iterativo es el que sigue.
Algoritmo para el método de la Falsa Posición

Permite encontrar una solución a la ecuación f(x)=0, dada la función continua f en el interva-
lo [xL,xU], considerando la siguiente notación:
xL: límite inferior del intervalo considerado
xU: límite superior del intervalo considerado
xM: raíz aproximada
f(xL): valor de la función en xL
f(xU): valor de la función en xU
Paso 1: Tomar i = 1
xU − xL
Paso 3: Tomar x M = x U − f ( x U )
f(xU ) − f(xL )
xU − xL
Paso 4: Si f ( x U ) < E ó f(xM)=0, SALIDA xM
f(xU ) − f(xL )
PARAR
Paso 6: Si f(xL). f(xM) >0, tomar xL=xM,
si no tomar xU=xM

PARAR
25
Cálculo Numérico
¾ Por ejemplo, si se desea determinar, aplicando el método de la falsa posición, una de las raí-
ces de la ecuación planteada en el ítem 3.2.1.1, considerando que la función cambia de signo
en el intervalo (1,2).
Se iniciarán los cálculos con los valores iniciales xL = 1 y xU = 2
Primera iteración:
xL = 1, f(xL) = -4
xU = 2 f(xU) = 3
(2 − 1)
xM = 2 − 3 = 1,57142
3 − ( −4)
Segunda iteración:
Como f(xM)= -1,36449 tiene el mismo signo que f(xL), xM se convierte en el límite superior de la
siguiente iteración, xU = 1,57142
xL = 1, f(xL) = -4
xU = 1,57142 f(xU) = -1,36449

(1,57142 − 1)
x M = 1,57142 − 1,36449 = 1,70540
− 1,36449 − ( −4)
Se procede de esta manera hasta que en la quinta iteración el valor de xM es 1,73194 muy
próximo al valor verdadero 1,7321
Las ventajas y desventajas del método son:
9 Es siempre convergente - Si existe más de una raíz en el interva-
lo, el método permite encontrar sólo una
de ellas
9 Converge más rápidamente que el
método de la bisección
3.2.2. Métodos abiertos
3.2.2.1. Método de Aproximaciones sucesivas
El método de aproximaciones sucesivas o iteración de punto fijo es una forma muy útil y
simple de encontrar la raíz de una ecuación de la forma f(x)=0. Para ello se reordena la ecuación
de manera que x sea igual a g(x). Esta transformación se puede llevar a cabo mediante operacio-
nes algebraicas o simplemente agregando x en ambos miembros de la ecuación original. A una
solución de esta ecuación se le llama un punto fijo de la función g. Sin embargo, es muy importan-
te la selección de la función g(x), ya que no siempre converge con cualquier forma elegida de g(x).
26
Cálculo Numérico
En síntesis, sea:
f(x)=0 (3.3)
una ecuación algebraica o trascendente y x = x* una raíz de ella o sea un valor de x tal que la veri-
fique idénticamente, es decir: f(x*)=0. Sumando x a ambos miembros de (3.3) se tiene f(x) +x = x y
llamando g(x)= x + f(x) se tiene que:
x = g(x) (3.4)
El método de aproximaciones sucesivas consiste en sustituir x0, un valor aproximado de la

raíz x* en el segundo miembro de la ecuación (3.4), con lo que se obtiene: x1=g(x0)
Procediendo reiteradamente de esta manera, la i-ésima aproximación o i-ésima iteración

es: xi+1=g(xi) (3.5)
Algoritmo para el método de Iteración de Punto Fijo

Permite encontrar una solución a la ecuación x=g(x), dada una aproximación inicial x0.
x0: aproximación inicial a la raíz
x: aproximación a la raíz
Paso 1: Tomar i = 1
Paso 3: Tomar x= g(x0)
Paso 4: Si x − x 0 < E , SALIDA x
PARAR
Paso 6: Tomar x0=x

PARAR
Si a medida que n crece, xi+1 tiende a x* se dice que el método converge, en caso contrario
diverge. Si el método converge, la diferencia entre dos iteraciones sucesivas será más pequeña a
medida que i aumenta, lo que proporciona un criterio de terminación de aplicación del método. Se
acepta el siguiente teorema sin demostración:
27
Cálculo Numérico
Teorema: El método de aproximaciones sucesivas converge si existe un número fijo m tal

que: f ' ( x ) ≤ m < 1 .
Un planteamiento gráfico diferente es el de separar la ecuación x = g(x) en dos partes, co-

mo y1 = x (recta a 45º) y y2 = g(x), éstas se pueden graficar por separado. Los valores de x co-
rrespondientes a las intersecciones de estas funciones representan las raíces de f(x) = 0. En la
Figura 3.4 se muestra la convergencia (a) y (b) ya que verifican el teorema de la convergencia y la
divergencia (c) y (d) en el método de Aproximaciones sucesivas.
y
(a) (b)
y
x x1 x x* x x0 x2 x* x1 x
(c)
(d)
y
y
x* x1 x2 x1 a x0 x2 x
x
Figura 3.4. Convergencia y divergencia del Método de Iteración de Punto Fijo
¾ Por ejemplo, si se desea determinar, aplicando el método de aproximaciones sucesivas, una

de las raíces de la ecuación x2 - 4x + 2=0, existen muchas formas de cambiar la ecuación a la
forma x=g(x).
x2 + 2
Si se despeja x de la ecuación se tiene: x = , por lo tanto la ecuación de recurrencia
4
( xi )2 + 2
o iteración es xi+1 =
4
28
Cálculo Numérico
En la tabla siguiente se muestran los valores obtenidos, se comienza con una aproxima-
ción x0=1. El valor real de la raíz (0.586) se alcanza luego de cinco iteraciones.
xi xi+1
1 0,75
0,75 0,641
0,641 0,603
0,603 0,591
0,591 0,587
0,587 0,586

9 Es simple - No siempre es convergente, depende de
la forma de la función g(x)
9 Es flexible
3.2.2.2. Método de Newton-Raphson o de la Tangente
En este método si el valor inicial de la raíz es xi, se puede extender una tangente desde el
punto (xi, f(xi)). El punto donde esta tangente corta al eje x representa una aproximación mejorada
de la raíz.
Existen por lo menos tres maneras usuales de introducir el método de Newton – Raphson
puesto que se puede derivar a partir de un método gráfico ó a partir de la de iteración de punto fijo
ó bien utilizando la serie de Taylor. El desarrollo a partir de esta serie es el siguiente:
f ' ' ( ξ)
f ( x i +1 ) = f ( x i ) + f ' ( x i )( x i +1 − x i ) + ( x i +1 − x i ) 2 (3.6)
2
donde ξ se encuentra en alguna parte del intervalo xi y x i+1. Truncando la serie de Taylor después
de la primera derivada, se obtiene:
f(x i+1) ≅ f(xi) + f ’ (xi)(x i+1-xi) (3.7)
donde f ’(xi) es además de la derivada primera, la pendiente de la recta descripta.
En la intersección con el eje x, f(x i+1) debe ser igual a cero, o sea:
0 ≅ f(xi) + f ’ (xi)(x i+1-xi) (3.8)
f ( xi )
Resolviendo para x i+1: x i +1 = x i − (3.9)
f ' ( xi )
La fórmula (3.9) se denomina Fórmula de Newton – Raphson.
29
Cálculo Numérico
Este método definido por el denominador f ’(xi) hace que geométricamente se pase de una
aproximación a la siguiente por la tangente a la curva y = f(x) trazada en el punto correspondiente
a la aproximación presente, esto puede observarse en la Figura 3.5.
y
x* x2 x1 x0 x
Figura 3.5. Método de Newton-Raphson

Algoritmo para el método de Newton-Raphson

Permite encontrar una solución a la ecuación f(x)=0, dada una aproximación inicial x0.
f(x): función en estudio
f’(x): derivada de la función
Paso 1: Tomar i = 1
f(x 0 )
Paso 3: Tomar x = x 0 −
f ' (x 0 )
PARAR
Paso 6: Tomar x0=x

PARAR
¾ Por ejemplo, si se desea determinar, aplicando el método Newton-Raphson, una de las raíces
de la ecuación x2 - 4x + 2=0, se calcula la derivada primera de la función dada.
30
Cálculo Numérico
Como f(xi)= x2 - 4x + 2, su derivada primera es: f’(xi)= 2x - 4. Por lo tanto la fórmula de recu-
rrencia es:
f(xi ) x 2 − 4x + 2
x i +1 = x i − = xi − i
f ' (xi ) 2x i − 4
En la tabla siguiente se muestran los valores obtenidos, se comienza con una aproxima-
ción x0=1. El valor real de la raíz (0.586) se alcanza luego de dos iteraciones.
xi xi+1
1 0,5
0,5 0,583
0,583 0,586

9 Converge más rápido que cual- - No siempre es convergente, depende de
quiera de los métodos analizados la naturaleza de la función
hasta ahora.
- No es conveniente en el caso de raíces
múltiples
- Puede alejarse del área de
interés si la pendiente es cercana a cero
3.2.2.3. Método de Newton de segundo orden
Si en lugar de considerar los dos primeros términos de la serie de Taylor se consideran los
tres primeros términos (3.6), se representa con ∆xi a la diferencia entre x i+1 y xi y se iguala a cero,
se tiene:
( ∆x i ) 2
f ( x i ) + ∆x i f ' ( x i ) + f ' ' (xi ) = 0 (3.10)
2
f(xi )
y sustituyendo ∆xi por − (a partir de la fórmula de Newton-Raphson) queda:
f '(xi )
⎡ 1 f(xi ) ⎤
f ( x i ) + ∆x i ⎢f ' ( x i ) − f ' ' ( x i )⎥ = 0 (3.11)
⎣ 2 f ' (xi ) ⎦
Despejando ∆xi se obtiene:
f(x i )
∆x i = − (3.12)
f(xi )
f ' (xi ) − f ' ' (xi )
2f ' ( x i )
31
Cálculo Numérico
De la ecuación (3.12) se puede despejar el valor de x i+1:
f(xi )
x i +1 = x i − (3.13)
f(xi )
f ' (xi ) − f ' ' (x i )
2f ' ( x i )
Este método considera un mayor número de términos de la serie por lo tanto converge
más rápidamente que el método de Newton-Raphson.
Algoritmo para el método de Newton de segundo orden

f’(x): derivada primera de la función
f’’(x): derivada segunda de la función
Paso 1: Tomar i = 1
f(x 0 )
f(x0 )
f ' (x 0 ) − f ' ' (x 0 )
2f ' ( x 0 )
PARAR
Paso 6: Tomar x0=x

PARAR
3.2.2.4. Método de Von Mises
El método de Newton-Raphson puede ser problemático si se está en puntos alejados de

las raíces y cercanos a puntos donde el valor de f’(xi) sea próximo a cero. Para ello von Mises
sugirió utilizar Newton-Raphson (fórmula 3.9) sustituyendo el denominador f’(xi) por f’(x0), es decir
obtener geométricamente las siguientes aproximaciones por medio de paralelas a la primera tan-
gente. La ecuación de recurrencia es:
32
Cálculo Numérico
f ( xi )
x i +1 = x i − (3.14)
f ' ( x0 )
En la Figura 3.6 se muestra un esquema de la aplicación del Método de von Mises.
x* x3 x2 x1 x0 x
Figura 3.6. Método de von Mises

Algoritmo para el método de Von Mises

f’(x00): valor de la derivada de la función para el valor inicial de aproximación x0
Paso 1: Tomar i = 1
f(x 0 )
f ' ( x 00 )
PARAR
Paso 6: Tomar x0=x

PARAR
3.2.2.5. Método de la secante
Surge como una variación del método de Newton-Raphson, en lugar de tomar la tangente
se toma la secante. De manera que la derivada se aproxima por una diferencia dividida, es decir:
33
Cálculo Numérico
f ( x i −1 ) − f ( x i )
f ' ( xi) ≈ (3.15)
x i −1 − x i
Esto puede sustituirse en la fórmula (3.9), de manera que se llega a:
f ( x i )( x i −1 − x i )
x i +1 = x i − (3.16)
f ( x i −1 ) − f ( x i )
La ecuación (3.16) es la fórmula para el método de la secante. El método requiere de dos valores
iniciales pero como no se requiere que f(x) cambie de signo en el intervalo considerado, no se lo
incluye dentro de los métodos que utilizan intervalos.
Algoritmo para el método de la Secante

Permite encontrar una solución a la ecuación f(x)=0, dadas dos aproximación inicial x0 y x1.
Paso 1: Tomar i = 2
f ( x 1 )( x 0 − x1 )
f ( x 0 ) − f ( x1 )
PARAR
Paso 6: Tomar x0=x1;
x1=x

PARAR
En la Figura 3.7 se muestra un esquema del método.
34
Cálculo Numérico
y
f(xi)
f(xi-1)
x* xi-1 xi x
Figura 3.7. Método de la Secante
3.3. RAÍCES DE POLINOMIOS
Los polinomios son funciones de relevancia en modelos de ciencia e ingeniería, a conti-

nuación se detallan teoremas y reglas necesarias para el cálculo de sus raíces.
3.3.1. Teoremas fundamentales de la Teoría de ecuaciones algebraicas.
El teorema fundamental del Algebra indica:
Teorema Fundamental del Algebra: Toda ecuación algebraica de grado n admite n raíces reales o comple-
jas.
A continuación se enuncian y demuestran algunos teoremas de interés para encontrar las

raíces de polinomios.
Teorema del residuo: el residuo que resulta de dividir el polinomio P(x) entre el binomio (x – a), es igual al
valor del polinomio cuando x = a.
Demostración:
Sea Q(x) el cociente y R el residuo que resulta de dividir P(x) entre x – a, por definición de
división de un polinomio entre un binomio se tiene:
P(x) = (x – a) Q(x) + R
y si: x = a, P(a) = (a-a) Q(a) +R ⇒ P(a) = R con lo que queda demostrado el teorema.
Teorema recíproco: el valor del polinomio P(x) para x = a, es igual al residuo que resulta de dividir P(x) entre
x – a.
35
Cálculo Numérico
Teorema del factor: si x =a es una raíz de la ecuación P(x) = 0, entonces x – a es un factor del polinomio
P(x).
P(x) = (x – a) Q(x) + R
y si: x = a, P(a) = (a-a) Q(a) +R ⇒ P(a) = R
Como P(a) =0 puesto que a es una raíz de P(x) =0, el resto R es igual a cero por lo tanto
puede afirmarse que (x – a) es un factor de P(x).
3.3.2. División sintética
La división de un polinomio P(x) entre (x-a) puede expresarse como:
P(x) = (x – a) Q(x) + R (3.17)
donde Q(x) es el cociente y R el resto o residuo.
Si se considera que
Q ( x ) = A 0 x n − 1 + A1 x n − 2 + A 2 x n − 3 + L + A n − 2 x + A n − 1 es el cociente y R el re-
siduo o resto que resulta de dividir el polinomio: Pn ( x ) = a 0 x n + a1 x n −1 + L + an −1 x + an entre el

binomio (x- a), entonces puede expresarse como:
a 0 x n + a1x n −1 + L + a n −1x + a n = A 0 x n + ( A 1 − aA 0 )x n −1 + ( A 2 − aA 1 )x n − 2 + L +
+ ( A n −1 − aA n − 2 )x + (R − aA n −1 )
Como los polinomios de ambos miembros son iguales los coeficientes de las mismas po-
tencias de x en ambos polinomios deben ser iguales entre sí, luego:
a0=A0
a1=A1 - a A0
...
an-1=A n-1 - a A n-2
an=R - a A n-1
de donde se obtiene:
A0 =a0
A1= a1 + a A0
...
A n-1 = a n-1 + a A n-2
R = an + a A n-1
Estos son los coeficientes del polinomio, cociente y residuo buscados, los cálculos se pue-
den arreglar de la siguiente forma:
36
Cálculo Numérico
a0 a1 ... a n -1 an
a a A0 a A n-2 a A n-1
A0 A1 ... An-1 R
Este es un esquema de la división sintética, se debe ordenar el polinomio P(x) en poten-

cias decrecientes de x, insertando un 0 para todos los términos con coeficientes nulos. Si P(x) es
de grado n, entonces el cociente Q(x) es de grado n-1.
3.3.3. Regla de los signos de Descartes
El número de raíces reales positivas en la ecuación algebraica de coeficiente reales P(x)

=0, es igual al número de cambios de signo en el polinomio P(x) o es menor que este número en
un número par. El número de raíces negativas es igual al número de cambios de signo en el poli-
nomio P(- x) o es menor que este número en un número par.
3.3.4. Raíces racionales
Para determinar las raíces racionales de una ecuación algebraica de coeficientes enteros
o reales si se elimina la parte decimal multiplicándose por un número lo suficientemente grande
pueden establecerse los siguientes pasos:
1. Escribir la ecuación en orden descendente de potencias de x.
2. Separa todas las raíces nulas.
3. Determinar los números máximos de raíces positivas y negativas por la regla de los
signos de Descartes.
4. Establecer las posibles raíces racionales.
5. Probar que una de estas es raíz, aplicando el teorema recíproco del factor.
6. Separar la raíz determinada y estudiar la ecuación reducida obtenida, de manera de

eliminar de la lista original de posibles raíces racionales las que ya no pueden ser.
3.3.5. Raíces irracionales
Si una ecuación algebraica posee raíces irracionales, en primer lugar se deben aplicar los
procedimientos descriptos anteriormente para encontrar y separar las raíces racionales, de forma
que se tenga una ecuación reducida que posea solamente raíces irracionales. Si esta ecuación es
de primer o segundo grado, sus raíces se obtienen por medio de fórmulas, para grados superiores
al segundo se pueden aplicar los métodos detallados anteriormente.
37
Cálculo Numérico
3.3.5.1. Método de Newton-Raphson
Como la función considerada es un polinomio, se puede escribir la fórmula (3.9) como:
P( x i )
x i +1 = x i − (3.18)
P' ( x i )
El polinomio P(x) puede expresarse como la ecuación (3.17), si se toma x = xi, se tiene que:
P(xi)= R (3.19)
El denominador de la ecuación (3.18) puede obtenerse derivando la expresión (3.17) con

respecto a x:
P’(x)=(x-xi)Q’(x)+Q(x) (3.20)
Y haciendo x=xi, se llega a que:
P’(xi)= Q(xi) (3.21)
Y de acuerdo con el teorema recíproco al del residuo Q(x) puede determinarse como el residuo R’
que resulta de dividir Q(x) entre (x-xi) puesto que:
Q(x)=(x-xi) S(x)+R’ (3.22)
Si se toma x=xi, resulta:
P’(xn) = Q(xn) =R’ (3.23)
Sustituyendo en (3.18) se obtiene la expresión de Newton-Raphson para resolver una ecuación

algebraica:
R
x i +1= x i − (3.24)
R'
Realizando consideraciones similares se llega a la fórmula de recurrencia de Newton de

segundo orden para una ecuación algebraica:
1 P ' ( x i ) 1 P' ' ( x i ) R' ' R'

=− + = − (3.25)
∆x i P( x i ) 2 P ' ( x i ) R' R
¾ Por ejemplo, si se desea determinar, aplicando el método Newton-Raphson, una de las raíces
irracionales de la ecuación x3 +2x2-3x-3=0, (se sabe que una de las raíces es 1.4605), se co-
mienza con x=2:
1 2 -3 -3
2 2 8 10
1 4 5 7
2 2 12
1 6 17
38
Cálculo Numérico
Entonces se aplica la fórmula (3.24) para obtener una mejor aproximación de la raíz:
7
x1 = 2 − = 1.5882
17
1 2 -3 -3
1,5882 1,5882 5,6988 4,2862
1 3,5882 2,6988 1,2862
1,5882 1,5882 8,2211
1 5,1764 10,9198
1,2862
x 2 = 1,5882 − = 1,4704
10,9198
Se realiza la división sintética para este valor y así sucesivamente hasta llegar al valor
buscado.
3.1. Calcule la raíz cuadrada negativa de 0.8 utilizando el método de aproximaciones sucesivas.
3.2. Evalúe, aplicando Iteración de Punto Fijo, el factor de fricción f en una tubería por la que circu-
la un fluido con flujo turbulento. Este factor está dado por la ecuación:
1 ⎡e 9,35 ⎤
= 1,14 − 2 log⎢ + ⎥ . Considere que el diámetro D= 0,2 m; el espesor e=0,0035 m y el
f ⎣ D Re f ⎦
número de Reynolds, Re =3x104.
3.3. Resuelva la ecuación: sen 2x = 0, a partir de x0 =1,165, aplicando el método de Newton de 2º

orden.
3.4. Aplique el método de von Mises, y luego compárelo con Newton-Raphson y Newton de 2º
orden, para resolver: 4x3 – 18 x2 + 12 x – 6 = 0
3.5. Considere la pared de un horno, de 0,08 m de espesor, la temperatura del lado interno es 642
K, si las pérdidas de calor desde la superficie externa se producen por convección y radiación,
determine la temperatura del lado externo de la pared (T1). La ecuación que rige esta situación
problemática es:
k
∆x
( )
(T1 − T0 ) + ε σ T14 − Tf 4 + h (T1 − Tf ) = 0
Los datos son:
Conductividad térmica, k =1,33 W/mK; Emisividad, ε = 0.8;Temperatura del lado interno de la pa-
red, T0 = 642 K; Temperatura del lado externo de la pared, T1;Temperatura del aire, Tf =299 K;
39
Cálculo Numérico
Coeficiente convectivo de transferencia calórica, h =18 W/m2K; Constante de Stefan-Boltzmann, σ

=5,67 x10-8 W/m2K4; Espesor de la pared, ∆x=0,08 m.
3.6. Un proyecto de Ingeniería Química requiere que se calcule, exactamente, el volumen molar
(v) de monóxido de carbono a 80 bares de presión y 226 K, de tal forma que se pueda determinar
el tubo adecuado que los contenga. Aplique los métodos de Newton Raphson y von Mises.
Datos: R=0.08314 bar m3/kgmol K; a=1.396 bar m6/(kgmol)2; b= 0,0345 m3/kgmol
3.7. Demuestre que: 1, 2 y -2, son raíces de la ecuación: x3 –x2 - 4 x + 4 =0, utilizando el teorema
del factor.
3.8. Resuelva la ecuación: x3 –7 x2 - 10 x + 16 =0, utilice la regla de los signos de Descartes.
40
Cálculo Numérico
Capítulo 4
SISTEMAS DE ECUACIONES LINEALES
En muchas ocasiones los problemas de matemáticas aplicadas a la ciencia e ingeniería

pueden reducirse a un sistema de ecuaciones algebraicas lineales. Estos sistemas pueden resol-
verse tanto por métodos exactos como por métodos aproximados.
4.1. CONCEPTOS PREVIOS
Una ecuación algebraica lineal es una ecuación en donde en cada término aparece úni-
camente una variable o incógnita elevada a la primera potencia.
Por ejemplo, a11 x1 + a12 x2 + a13 x3 + … + a1n xn = c1, es una ecuación algebraica lineal en
las variables x1, x2, x3, …, xn. Se admite que los coeficientes a11, a12, a13,…, a1n y el término inde-
pendiente c1, son constantes reales.
Un sistema de ecuaciones lineales es un conjunto de ecuaciones que deben resolverse

simultáneamente. Por ejemplo,
a11 x1 + a12 x2 + a13 x3 + … + a1n xn = c1
a21 x1 + a22 x2 + a23 x3 + … + a2n xn = c2
a31 x1 + a32 x2 + a33 x3 + … + a3n xn = c3 (4.1)
an1 x1 + an2 x2 + an3 x3 + … + ann xn = cn
Aplicando la definición de producto de matrices, en este sistema de n ecuaciones algebrai-

cas lineales con n incógnitas, puede escribirse en forma matricial:
⎡a11 a12 a13 K a1n ⎤ ⎡ x1 ⎤ ⎡ c 1 ⎤

⎢a a 22 a 23 a 2n ⎥ ⎢ x ⎥ ⎢c ⎥
⎢ 21 ⎥ ⎢ 2⎥ = ⎢ 2⎥ (4.2)
⎢K K K K⎥ ⎢K ⎥ ⎢K⎥
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎣ a n1 a n2 a n3 a nn ⎦ ⎣ x n ⎦ ⎣c n ⎦
Este sistema de ecuaciones se simboliza como [A]nxn [X]nx1 = [C]nx1, en donde A es la matriz
del sistema, X es el vector incógnita y C es el vector de términos independientes, que en forma
sintética se simboliza como AX=C.
41
Cálculo Numérico
La matriz formada por A, a la que se le agrega el vector de términos independientes como

última columna, se denomina matriz ampliada del sistema que se simboliza como [Aa].
Solución de un sistema de ecuaciones: es un conjunto de valores de las incógnitas que

verifican simultáneamente a todas y cada una de las ecuaciones del sistema.
De acuerdo con su solución, un sistema puede ser compatible, si admite solución; o in-
compatible si no admite solución. Un sistema compatible puede ser determinado, si la solución
es única; o indeterminado, si la solución no es única.
Teorema de Rouchè – Frobenius: Si el rango de la matriz de coeficientes es igual al rango de la matriz am-
pliada (rg (A) = rg (Aa)) entonces A X = C es compatible, y recíprocamente.
El corolario de este teorema es el siguiente:

Un sistema Compatible será determinado (solución única) si el rango de la matriz de coeficientes es
igual al número de incógnitas r(A)=n, y será indeterminado (infinitas soluciones) si el rango de la matriz de
coeficientes es menor que el número de incógnitas r(A) < n
Las soluciones de un sistema compatible de la forma AX=C permanecen invariantes ante
las siguientes operaciones elementales:
Intercambio de dos filas o renglones cualesquiera.
Multiplicación de una fila por un escalar no nulo.
Suma a una fila de una combinación lineal no nula de otro renglón
4.2. MÉTODOS DE RESOLUCIÓN DE SISTEMAS DE ECUACIONES ALGEBRAICAS LINELES
Para la resolución de Sistemas de ecuaciones algebraicas lineales se utilizan dos tipos de

métodos: métodos directos y métodos iterativos.
4.2.1. Métodos directos
Los métodos directos son aquellos que obtiene la solución exacta, salvo errores de redon-
deo en los cálculos, luego de un número finito de operaciones elementales. Pertenecen a este
grupo el método de eliminación de Gauss, el método de Gauss-Jordan, partición de matrices, etc.
A continuación se describen los métodos mencionados.
4.2.1.1. Método de eliminación de Gauss

Consideremos el sistema de ecuaciones algebraicas lineales:
42
Cálculo Numérico
a11 x1 + a12 x2 + a13 x3 + … + a1n xn = c1 (4.3a)
a21 x1 + a22 x2 + a23 x3 + … + a2n xn = c2 (4.3b)
a31 x1 + a32 x2 + a33 x3 + … + a3n xn = c3 (4.3c)
an1 x1 + an2 x2 + an3 x3 + … + ann xn = cn (4.3d)
El procedimiento para resolver este sistema consta de dos pasos:
1. Eliminación hacia adelante de incógnitas.
2. Sustitución hacia atrás
1. Eliminación hacia adelante
En la eliminación hacia adelante, se reduce el conjunto de ecuaciones a un sistema trian-

gular superior. El paso inicial consiste en multiplicar la primera ecuación del sistema (4.3a) por el
a
cociente entre los coeficientes de la primera incógnita de la segunda y primera ecuación, 21 ,
a11
obteniéndose:
⎛ a ⎞ ⎛ a ⎞ c
a 21x 1 + ⎜ a 21 12 ⎟ x 2 + K + ⎜ a 21 1n ⎟ x n = a 21 1 (4.4)
⎝ a11 ⎠ ⎝ a11 ⎠ a11
Como el primer término de la primera ecuación modificada (4.4) es idéntico al primer

término de la segunda ecuación (4.3b), se elimina la primera incógnita de (4.3b) restando la ecua-
ción (4.4) de (4.3b) y se llega a:
⎛ a ⎞ ⎛ a ⎞ c
⎜ a 22 − a 21 12 ⎟ x 2 + K + ⎜ a 2n − a 21 1n ⎟ x n = c 2 − a 21 1 , es decir:
⎝ a11 ⎠ ⎝ a11 ⎠ a11
a’22x2+…+a’2nxn =c’2 (4.5)
el apóstrofe se utiliza para indicar que los coeficientes de las incógnitas han sufrido modificaciones
en sus valores. El proceso se repite hasta que se elimina la primera incógnita de las ecuaciones
restantes dando como resultado el siguiente sistema modificado:
a11 x1 + a12 x2 + a13 x3 + … + a1n xn = c1 (4.6a)
a’22 x2 + a’23 x3 + … + a’2n xn = c’2 (4.6b)
a’n2 x2 + a’n3 x3 + … + a’nn xn = c’n (4.6d)
43
Cálculo Numérico
La ecuación pivotal, es decir la que permanece invariante es la ecuación (4.6a). A conti-

nuación se repite el proceso para eliminar la segunda incógnita (x2) desde la tercera ecuación has-
a' i2
ta la última, restando a cada ecuación la segunda ecuación (4.6b) multiplicada por (recor-
a' 22
dando que i representa al número de fila o renglón, siendo i≥3): Una vez completado este paso se
repite el procedimiento de manera de eliminar las incógnitas iniciales de las ecuaciones subsi-
guientes hasta llegar a la última, transformándose el sistema en un sistema triangular superior:
a11 x1 + a12 x2 + a13 x3 + … + a1n xn = c1 (4.7a)
a’22 x2 + a’23 x3 + … + a’2n xn = c’2 (4.7b)
a’’33 x3 + … + a’’3n xn = c’’3 (4.7c)
a(n-1)nn xn = c(n-1)n (4.7d)
2. Sustitución hacia atrás
La ecuación (4.7c) puede resolverse para xn:
(n −1)
cn
xn = (4.8)
n −1)
a(nn
Este resultado se puede sustituir en la (n-1)-ésima ecuación y resolver ésta para x n-1. El
procedimiento se repite, evaluando las x restantes. Esquemáticamente:
⎡ a 11 a 12 a 13 c1 ⎤
⎢a a 22 a 23 c2⎥
⎢ 21 ⎥
⎣⎢ a 31 a 32 a 33 c 3 ⎦⎥
Eliminación hacia adelante
⇓
⎡ a 11 a 12 a 13 c1 ⎤
⎢ a ' 22 a ' 23 c '2 ⎥
⎢ ⎥
⎣⎢ a ' ' 33 c ' ' 3 ⎥⎦
⇓
c ' '3
x3 =
a ' '3
(c '1 − a ' 23 x 3 )
x2 = Sustitución hacia atrás
a 22
(c 1 − a 12 x 2 − a 13 x 3 )
x1 =
a 11
Una de las desventajas de este método es que durante el proceso en las fases de elimina-
ción y sustitución es posible que ocurra una división entre cero. Por ello se ha desarrollado una
estrategia del pivoteo que evita parcialmente estos problemas.
44
Cálculo Numérico
Si se resuelve un pequeño número de ecuaciones, el error por redondeo es pequeño y

generalmente no afecta sustancialmente la precisión de los resultados, pero si se van a resolver
simultáneamente muchas ecuaciones, el efecto acumulativo del error por redondeo puede introdu-
cir errores relativamente grandes en la solución. Por esta razón el número de ecuaciones simultá-
neas que se puede resolver satisfactoriamente con el método de eliminación de Gauss, utilizando
de 8 a 10 dígitos significativos en las operaciones aritméticas, se limita generalmente a 15 o 20.
A continuación se indica el algoritmo del método.
Algoritmo de Eliminación de Gauss

Considerando el sistema (4.3) y la siguiente notación:
n: número de ecuaciones
aij: elementos de la matriz ampliada Aa (1 ≤ i ≤ n y 1 ≤ j ≤ n + 1)
p: índice del elemento pivote
Fi: fila i
Paso 1: Para i = 1, ..., n-1 seguir los Pasos 2 a 4. (Eliminación hacia adelante)
Paso 2: Sea p el menor entero con i ≤ p ≤ n y api≠0.
Si p no puede encontrarse, SALIDA (‘No existe solución única’)
PARAR
Paso 3: Si p≠i intercambiar la fila p por la fila i
Paso 4: Para j= i+1, ..., n seguir los Pasos 5 a 6
a ij
Paso 5: Hacer m ij =
a ii
Paso 6: Realizar Fj-mijFi e intercambiarla por la fila Fj
Paso 7: Si ann=0 entonces SALIDA (‘No existe solución única’)
PARAR
cn
Paso 8: Hacer x n = (Sustitución hacia atrás)
a nn
⎛ ⎞
1 ⎜ n ⎟
Paso 9: Para i= n-1, ..., 1 tomar xi = ⎜ c i − ∑ a ij x j ⎟
a ii ⎜ j = i +1 ⎟⎟
⎜
⎝ ⎠
Paso 10: SALIDA X (es decir (x1,x2, ...,xn))
PARAR
⎧x1 + x 2 + 2x 3 = 3
⎪⎪
¾ Por ejemplo, se desea resolver el sistema ⎨3 x 1 − x 2 + x 3 = 1 , aplicando el método de eli-
⎪
⎪⎩2x 1 + 3 x 2 − 4 x 3 = 8
minación gaussiana.
45
Cálculo Numérico
a) Eliminación hacia adelante
Como el coeficiente de la primera incógnita es 1, se multiplica la primera ecuación por 3 y

se resta el resultado de la segunda ecuación, luego se multiplica por 2 la primera ecuación y se
resta de la tercera de manera que el sistema queda reducido a:
⎧x 1 + x 2 + 2x 3 = 3
⎪⎪
⎨ − 4 x 2 − 5 x 3 = −8
⎪
⎪⎩ x 2 − 8x 3 = 2
Se procede ahora a eliminar la segunda incógnita de la tercera ecuación, para ello se divi-
de la segunda ecuación por -4 y se multiplica por el coeficiente de la tercera ecuación que en este
5
caso es 1, quedando la segunda como: x 2 + x 3 = 2 y se resta este resultado de la tercera
4
ecuación. El sistema es ahora:
⎧
⎪x 1 + x 2 + 2x 3 = 3
⎪
⎪
⎨ − 4 x 2 − 5 x 3 = −8
⎪
⎪ 37
⎪⎩ − x3 = 0
4
b) Sustitución hacia atrás
Se despeja x3 de la tercera ecuación, en este caso: x3=0, se reemplaza este valor en la se-
gunda ecuación: − 4 x 2 = −8 por lo tanto x2=2 y por último se reemplazan estos valores en la pri-
mera ecuación x 1 + 2 = 3 entonces x1=1.
4.2.1.2. Método de Gauss - Jordan
Es diferente al método de eliminación gaussiana puesto que cuando se elimina una incóg-
nita no sólo se elimina de las ecuaciones siguientes sino de todas las otras ecuaciones. De esta
forma el paso de eliminación genera una matriz identidad en lugar de una matriz triangular. Por
consiguiente no es necesario emplear la sustitución hacia atrás para obtener la solución. A conti-
nuación se esquematiza el método.
46
Cálculo Numérico
⎡ a 11 a 12 a 13 c1 ⎤
⎢a a 22 a 23 c2 ⎥
⎢ 21 ⎥
⎢⎣ a 31 a 32 a 33 c 3 ⎥⎦
⇓
⎡1 0 0 c *1 ⎤
⎢0 1 0 c *2 ⎥
⎢ ⎥
⎢⎣0 0 1 c * 3 ⎥⎦
⇓
x1 = c *1
x2 = c *2
x3 = c *3
¾ Por ejemplo, si se desea resolver el sistema anterior utilizando este método, se escribe el sis-
tema en forma matricial, se trabaja con la matriz ampliada (formada por la matriz de coeficien-
tes a la que se le adiciona una última columna constituida por los términos independientes),
luego se efectúan operaciones elementales en las filas hasta llegar a la matriz identidad que-
dando los valores de las incógnitas en la última columna de la matriz ampliada.
⎡1 1 2 3⎤ −3F1 +F2 ⎡1 1 2 3⎤ ⎡1 1 2 3⎤ −1F2 +F1

⎢3 − 1 1 1⎥ ⎯⎯ ⎯⎯→ 0 − 4 − 5 − 8 ⎯⎯ ⎯⎯→⎢0 1 5 / 4 2⎥ ⎯⎯2⎯⎯
− 2F1 +F3 ⎢ ⎥ −1 / 4F2 3 −F +F3
→
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎣⎢2 3 − 4 8⎦⎥ ⎣⎢0 1 − 8 2 ⎥⎦ ⎣⎢0 1 − 8 2⎦⎥
⎡1 0 3/ 4 1⎤ 4 ⎡1 0 3 / 4 1⎤ −3 / 4F3 +F1 ⎡1 0 0 1⎤
− F3
⎢0 1 5/4 ⎥
2 ⎯⎯ ⎯
37 ⎯→⎢0 1 5 / 4 2⎥ ⎯⎯ ⎯ ⎯ ⎯→⎢0 1 0 2⎥
−5 / 4F3 +F2
⎢ ⎥ ⎢ ⎥ ⎢ ⎥
⎣⎢0 0 − 37 / 4 0⎦⎥ ⎣⎢0 0 1 0⎥⎦ ⎣⎢0 0 1 0⎦⎥
Se llega al mismo resultado que con el método anterior, es decir x1=1, x2=2 y x3=0.
Este no es un método recomendable ya que involucra alrededor de un 50% de cálculos

adicionales, sin que haya más beneficios.
4.2.1.3. Partición de matrices
Este método puede utilizarse para resolver grandes sistemas de ecuaciones lineales, y
consiste en dividir las matrices de la ecuación matricial del sistema en submatrices de orden
menor.
Dado un sistema de ecuaciones, expresado en notación matricial, [A]2x2 [X]2x1 = [C]2x1. Si

A22 representa una submatriz de A con inversa conocida, formada por alguno de los renglones da
A y el mismo número de columnas de A, se puede expresar la ecuación matricial en forma nota-
cional como:
⎡ A 11 A 12 ⎤ ⎡ x1 ⎤ ⎡ c 1 ⎤
⎢A ⎥ ⎢ ⎥= ⎢ ⎥ (4.9)
⎣ 21 A 22 ⎦ ⎣ x 2 ⎦ ⎣c 2 ⎦
47
Cálculo Numérico
En donde fijada la posición de la submatriz A22 en A, quedan obligadas las de las otras
submatrices.
Para resolver el sistema matricial se realizan las operaciones matriciales indicadas en él,
producto e igualdad, teniendo en cuenta que los elementos matriciales son a la vez matrices.
Se tiene:
A11 x1 + A12 x2 = c1 (4.10)
A21 x1 + A22 x2 = c2 (4.11)
Despejando el vector x2 de (4.11), se obtiene: A22 x2 = c2 - A21 x1 y por lo tanto:
x2 = A22-1(c2 - A21 x1) (4.13)
Sustituyendo (4.13) en (4.10) y despejando el vector x1, se llega a:
A11 x1 + A12 A22-1(c2 - A21 x1) = c1
(A11 - A12 A22-1 A21) x1 = c1 - A12 A22-1c2
x1 = (A11 - A12 A22-1 A21)-1(c1 - A12 A22-1c2) (4.14)
En esta última expresión se hace B = A11 - A12 A22-1 A21 y D = A12 A22-1, por lo tanto (4.14)
queda como:
x1 = B-1(c1 - Dc2) (4.15)
Reemplazando la ecuación (4.15) en la ecuación (4.13):
x2 = A22-1 c2 - A22-1 A21 B-1(c1 - Dc2)
x2 = - A22-1 A21 B-1c1 + (A22-1+ A22-1 A21 B-1D) c2 (4.16)
Haciendo en esta última expresión: E = A22-1 A21 y F = A22-1+ E B-1D, se puede escribir
(4.16) como sigue:
x2 = - EB-1 c1 + F c2 (4.17)
En síntesis se llega a:
⎡ x1 ⎤ ⎡ B −1 − B −1D⎤ ⎡ c 1 ⎤
⎢x ⎥ ⎢= ⎥⎢ ⎥ (4.18)
⎣ 2 ⎦ ⎣− EB −1 F ⎦ ⎣c 2 ⎦
4.2.2. Métodos iterativos
Los métodos iterativos o de aproximaciones sucesivas se utilizan cuando los sistemas de

ecuaciones a resolver son de grandes dimensiones o bien son sistemas dispersos (la matriz de
coeficientes posee muchos ceros). Estos métodos construyen una sucesión de aproximaciones a
48
Cálculo Numérico
la solución de las incógnitas hasta obtener una precisión determinada o hasta completar un núme-
ro determinado de iteraciones. Son ejemplos el método de Jacobi, el de Gauss-Seidel, etc.
4.2.2.1. Método de Jacobi
Si se considera un sistema de ecuaciones algebraicas, que puede escribirse en forma ma-

tricial como [A] [X] = [C] y que A = D + R, donde D es una matriz diagonal; es decir, una matriz
cuadrada cuyos elementos sobre la diagonal principal son los únicos diferentes de cero. Entonces
puede escribirse que:
(D + R) X = C
DX = C – RX
X = D-1C – D-1R X (4.19)
Se admite que la diagonal de A no contiene elementos nulos, para que exista la matriz in-
versa D-1. La ecuación (4.19) sugiere el método iterativo:
X k+1= D-1C – D-1R Xk (4.20)
Este es el método iterativo de Jacobi, de iteraciones totales o de desplazamientos simultá-

neos, definido por la ecuación de recurrencia (4.20), significa que del sistema de ecuaciones, se
despeja x1 de la primera ecuación, x2 de la segunda, etc., obteniéndose las siguientes ecuaciones:
c − a12 x 2 − a13 x 3 − K − a1n x n

x1 = 1 (4.21)
a11
c 2 − a 21x 1 − a 23 x 3 − K − a 2n x n
x2 = (4.22)
a 22
c − a 31x1 − a 32 x 2 − K − a 3n x n
x3 = 3 (4.23)
a 33
c n − a n1x1 − a n2 x 2 − K − a n,n −1x n −1

xn = (4.24)
a nn
Para aplicar el método, se considera una primera aproximación al valor de las incógnitas x,
que se denomina X(0) (el supraíndice indica el orden de aproximación). Se sustituye esta primera
aproximación en los segundos miembros de las ecuaciones (4.21) a (4.24), por ejemplo, si se to-
ma la solución trivial, en la ecuación (4.21) se encuentra x1 haciendo x2 hasta xn iguales a cero.
Luego se calcula x2 de la ecuación (4.22) tomando x1, x3, xn iguales a cero y así sucesivamente
hasta llegar a la última ecuación y encontrar xn. Se obtiene de esta manera una nueva aproxima-
49
Cálculo Numérico
ción a los valores de las incógnitas, es la aproximación de orden 1, es decir X(1). El procedimiento
se repite hasta que la solución converja cerca de los valores reales. La convergencia se puede
verificar usando el criterio de error relativo.
Este método es muy poco utilizado debido a que el método de Gauss-Seidel converge más
rápidamente a la solución y además lo hace cuando no se logra que el método de Jacobi conver-
ja.
La condición suficiente para que el método de Jacobi converja es que la matriz de coefi-
cientes sea diagonal dominante, es decir que cada elemento de la diagonal principal es mayor en
valor absoluto que la suma del resto de los elementos de la misma fila en la que se encuentra el
elemento en cuestión.
A continuación se presenta un algoritmo para este método iterativo.
Algoritmo de Jacobi
aij: elementos de la matriz A (i indica el número de fila y j el número de columna en el que se
encuentra el elemento en cuestión)
ci: elementos del vector C
x0i: componentes de la primera aproximación al vector solución (esta primera aproximación
es X0)
xi(k): componentes de la aproximación de orden k al vector solución, k varía de 1 a N (indica
el orden de aproximación o iteración)
Paso 1: Para k = 1
Paso 2: Mientras k ≤ N seguir con los pasos 3 a 6
Paso 3: Para i= 1, ..., n, calcular la aproximación de orden 1 mediante la fórmula:
⎛ ⎞
1 ⎜ n ⎟
xi = ⎜ c i − ∑ a ij x 0 j ⎟
a ii ⎜ j=1 ⎟⎟
⎜
⎝ j≠i ⎠
Paso 4: Si X − X 0 < E , SALIDA X (es decir (x1,x2, ...,xn))
PARAR
Paso 5: Tomar k = k+1
Paso 6: Para i= 1, ..., n tomar x0i=xi
Paso 7: SALIDA (‘Número máximo de iteraciones completado’)
PARAR
50
Cálculo Numérico
⎧2x 1 − x 2 − 0,1x 3 = 4
⎪⎪
¾ Por ejemplo, se desea resolver el sistema ⎨0,1x1 + 7 x 2 − 0,3 = 20 utilizando este método,
⎪
⎪⎩3 x 1 − 2x 2 + 100 x 3 = 450
suponiendo una cota de error del 3 %.
Se despeja x1 de la primera ecuación, x2 de la segunda y x3 de la tercera ecuación.

4 + x 2 + 0,1x 3 20 − 0,1x1 + 0,3 x 3 450 − 3 x 1 + 2 x 2
x1 = ; x2 = ; x3 =
2 7 100
Se supone como primera aproximación la solución trivial, entonces de la primera ecuación

se despeja x1 suponiendo que x2 = x3= 0 por lo tanto x1= 2.
20
Se calcula x2 suponiendo que x1=x3=0, es decir: x 2 = = 2,857 y por último x3 que es:
7
450
x3 = = 4,5 .
100
Se realiza una nueva iteración con los valores: x1= 2; x2= 2.857 y x3= 4.5, es decir:
4 + 2,857 + 0,1( 4,5)

x1 = = 3,653
2
20 − ((0,1) . ( 2)) + 0,3( 4,5)

x2 = = 3,021
7
450 − ((3).( 2)) + 2 (2,857 )

x3 = = 4,497
100
Se calcula el error relativo porcentual de aproximación:
3,653 − 2
E x1 = .100 = 45,25% ;
3,653
3,021 − 2,857
E x2 = .100 = 5,43% ;
3,021
4,497 − 4,5
E x3 = .100 = 0,07%
4,497
Se realiza una nueva iteración, ahora con x1= 3,653; x2= 3,021 y x3= 4,497 y se llega a que
x1=3,735 con Ex1= 2,19 %; x2=2,998 con Ex2= 0,77 %; x3=4,451 con Ex3= 1,03 %.
51
Cálculo Numérico
4.2.2.2. Método de Gauss–Seidel

Es un método iterativo que disminuye el error de redondeo, se denomina también de des-
plazamientos sucesivos o de iteraciones parciales. Si se tiene un conjunto de n ecuaciones (4.1),
que puede escribirse en forma matricial como: [A] [X] = [C] y si los elementos de la diagonal prin-
cipal son diferentes de cero, la primera ecuación se puede resolver para x1, la segunda para x2,
etc., lo que lleva a las ecuaciones (4.21) a (4.24).
Se puede comenzar el proceso de solución utilizando una aproximación inicial X0 a la solu-

ción que es el vector columna X. La solución trivial puede servir de valor inicial, se supone que x2,
..., xn valen 0. Estos valores se sustituyen en la ecuación (4.21) y de ella se despeja un nuevo va-
c1
lor de x 1 = . Luego se sustituye el nuevo valor de x1 con x3,…, xn iguales a cero en la ecua-
a11
ción (4.22) y se calcula un nuevo valor de x2. Este procedimiento se repite en cada una de las
ecuaciones hasta llegar a la ecuación (4.24) de la que se calcula un nuevo valor de xn. Se regresa
a la primera ecuación y se repite todo el proceso hasta que la solución converja cerca de los valo-
res reales. La convergencia se puede verificar usando el criterio de error relativo. Este método se
diferencia del de Jacobi puesto que una vez que se calcula una aproximación a una incógnita se
utiliza esta aproximación en la misma iteración.
Las condiciones suficientes para que el método de Gauss-Seidel converja es que la matriz
de coeficientes sea diagonal dominante o bien que la matriz de coeficientes sea simétrica y defini-
da positiva. Un algoritmo del método se muestra a continuación.
Algoritmo de Gauss-Seidel
aij: elementos de la matriz A (i indica el número de fila y j el número de columna en el que se
encuentra el elemento en cuestión)
ci: elementos del vector C
x0i: componentes de la primera aproximación al vector solución (esta primera aproximación
es X0)
xi(k): componentes de la aproximación de orden k al vector solución, k varía de 1 a N (indica
el orden de aproximación o iteración)
Paso 1: Para k = 1
Paso 2: Mientras k ≤ N seguir con los Pasos 3 a 6
Paso 3: Para i= 1, ..., n, calcular la aproximación de orden i mediante la fórmula:
1 ⎛ i −1 n ⎞
xi = ⎜ c i − ∑ a ij x − ∑ a ij x ⎟ continua
a ii ⎜ j 0j ⎟
⎝ j =1 j=i+ 1 ⎠
52
Cálculo Numérico
Paso 4: Si X − X 0 < E , SALIDA X (es decir (x1,x2, ...,xn))
PARAR
Paso 5: Tomar k = k+1
Paso 6: Para i= 1, ..., n tomar x0i=xi
Paso 7: SALIDA (‘Número máximo de iteraciones completado’)
PARAR
⎧2x 1 − x 2 − 0,1x 3 = 4
⎪⎪
¾ Por ejemplo, se desea resolver el sistema ⎨0,1x1 + 7 x 2 − 0,3 = 20 utilizando este método,
⎪
⎪⎩3 x 1 − 2x 2 + 100 x 3 = 450
suponiendo una cota de error del 3 %.
Se despeja x1 de la primera ecuación, x2 de la segunda y x3 de la tercera ecuación.

4 + x 2 + 0,1x 3 20 − 0,1x1 + 0,3 x 3 450 − 3 x 1 + 2 x 2
x1 = ; x2 = ; x3 =
2 7 100
Se supone que x2 = x3= 0, por lo tanto x1= 2, con este nuevo valor y con x3= 0 se calcula x2,
es decir:
20 − ((0,1) . ( 2))
x2 = = 2,828
7
con este valor y x1= 2 se calcula x3,
450 − ((3).(2)) + 2 (2,828 )

x3 = = 4,497
100
Con estos nuevos valores de x2 y x3 se determina un nuevo valor de x1
4 + 2,828 + 0.1( 4,497 )

x1 = = 3,639
2
Con este valor y x3= 4,497 se calcula x2:
20 − ((0,1) . (3,639 )) + 0,3( 4,497 )

x2 = = 2,998
7
Y con este valor y x1=3,639 se encuentra el nuevo valor de x3
53
Cálculo Numérico
450 − ((3).(3,639 )) + 2 (2,998 )

x3 = = 4,451
100
Se calcula el error relativo porcentual de aproximación:
3,639 − 2
E x1 = .100 = 45,04% ;
3,639
2,998 − 2,828
E x2 = .100 = 5,67% ;
2,998
4,451 − 4,497
E x3 = .100 = 1,03%
4,451
En la tercera iteración se llega a que x1=3,722 con Ex1= 2,23 %; x2=2,995 con Ex2= 0,1 %;
x3=4,448 con Ex3= 0,07 %.
4.2.2.3. Método de Relajación

El método de relajación fue ideado por el ingeniero británico Richard Southwell, converge
más rápidamente que el de Gauss-Seidel. Consiste en tomar nueva aproximación a la solución
como una combinación lineal de la solución de la etapa anterior, es decir aplicando la fórmula
(4.24) correspondiente al método de Gauss-Seidel pero con la incorporación de un factor de rela-
jación denominado w:
1 ⎛ i −1 n ⎞
x (ik ) = (1 − w )x (ik −1) + w ⎜ c i − ∑ a ij x − ∑ a ij x ⎟
⎜ j 0j ⎟ (4.25)
a ii ⎝ j =1 j=i+ 1 ⎠
En esta fórmula el supraíndice k indica que es la k-ésima iteración y a w se le asigna un
valor entre 0 y 2, este valor se determina por prueba y error y el método se denominará:
Método de sub-relajación, si 0<w<1, es aplicable en los casos en que el método de Gauss-
Seidel diverge
Método de sobre-relajación, si 1<w<2, se aplica para acelerar la convergencia del método de
Gauss-Seidel
Método de Gauss-Seidel si w=1
4. 1. Utilice el método de eliminación gaussiana y el de Gauss-Jordan para resolver el siguiente

sistema de ecuaciones:
54
Cálculo Numérico
3 x1 – 0,1 x2 –2 x3 = 14
0,1 x1 + 8 x2 – 2,3 x3 = - 21
0,4 x1 – 0,8 x2 + 10 x3 = 75
4.2. En una fábrica de embutidos, se elaboran tres tipos diferentes de estos productos que res-
ponden a las formulaciones A, B y C. En el cuadro que sigue se presentan las proporciones nece-
sarias de las materias primas principales de cada formulación.
Si se dispone diariamente de 140 tn de hígado, 130 tn de carne porcina y 180 tn

de tocino, indique la producción diaria de cada tipo de producto en kg/d.
Tipo de Hígado Carne Tocino

embutido
A 3 1 2
B 1 3 2
C 3 2 4
4.3. Calcule el valor de las incógnitas x1 y x2 del sistema de 5 ecuaciones con 5 incógnitas, apli-
cando el método de partición de matrices.
x1 + x2 + x3 + x4 + x5 = 18
x1 + 2 x2 + 2 x3 +3 x4 + 4 x5 = 27
2 x1 + x2 - 2 x3 + 2 x4 - 3 x5 = - 32
3 x1 + 2 x2 + 3 x3 + 4 x4 + x5 = 6
- x1 + x2 - 4 x3 + 4 x4 + 2 x5 = - 24
4.4. Utilice el método de Jacobi para resolver el sistema indicado a continuación, calculando para
cada iteración el error relativo porcentual. Utilice como criterio de parada un error del 6 %, compa-
re con la solución obtenida anteriormente.
16 x1 +2 x2 – 2 x3 = 16
4 x1 + 24 x2 +18 x3 = 24
2 x1 – 14 x2 + 34 x3 = -16
4.5. Resuelva el sistema anterior utilizando el método de Gauss-Seidel.

4.6. Resuelva el sistema anterior por el método de sobre-relajación, utilizando w=1.4
55
Cálculo Numérico
Capítulo 5
APROXIMACION POLINOMIAL Y FUNCIONAL
En la práctica es frecuente tratar funciones que no son del tipo de las elementales, además
de funciones definidas de manera tabular o gráfica, de las que se desconoce su expresión analíti-
ca y de las que se necesita conocer valores de la variable que no están tabulados.
Existen casos de funciones expresadas en forma tabular en los que se requiere una alta
aproximación y para ello existen métodos numéricos que por lo general utilizan funciones raciona-
les enteras (polinomios), de manera que la curva descripta por los mismos toque todos los puntos
definidos. Si no se requiere gran aproximación se deriva una curva simple que represente el com-
portamiento general de los datos.
5.1. APROXIMACIÓN POLINOMIAL
Si se tiene una función definida en forma tabular de la que se desconoce su expresión

analítica puede afirmarse que para n+1 datos o puntos existe uno y solo un polinomio de n-ésimo
grado que pasa a través de todos los puntos y existe una variedad de fórmulas matemáticas que
permiten expresar a este polinomio.
5.1.1. Diferencias finitas
El cálculo de las diferencias finitas permite encontrar el grado del polinomio por el cual
puede describirse una función tabular.
Dada la función y=f(x) definida en forma tabular como la que se presenta en la Tabla 5.1, y
suponiendo que los valores de la variable independiente xn, están igualmente espaciados entre sí,
es decir que el incremento o paso es igual a un valor constante denominado h.
Tabla 5.1. Función tabular de paso constante
x y
x0 y0
x1= x0+h y1
x2= x0+2h y2
x3= x0+3h y3
… …
xn= x0+nh yn
Se denominan primeras diferencias hacia adelante y se representan con ∆yi a las dife-
rencias entre dos valores consecutivos de y, es decir:
56
Cálculo Numérico
a0 = y1 – y0 (5.1)
a1 = y2 – y1 (5.2)
a2 = y3 – y2 (5.3)
…
an-1 = yn – yn-1 (5.4)
Las diferencias de las primeras diferencias se llaman segundas diferencias hacia adelan-
te, ∆2yi:
b0 = a1 – a0 = y2 – 2y1 + y0 (5.5)
b1 = a2 – a1 = y3 – 2y2 + y1 (5.6)
b2 = a3 – a2 = y4 – 2y3 + y2 (5.7)
…
bn-2 = a n-1 – an-2 = yn – 2yn-1 + yn-2 (5.8)
Las diferencias de las segundas diferencias se llaman terceras diferencias hacia adelan-
3
te, ∆ yi:
c0 = b1 – b0 = y3 – 3y2 + 3y1 – y0 (5.9)

c1 = b2 – b1 = y4 – 3y3 + 3y2 – y1 (5.10)
c2 = b3 – b2 = y5 – 3y4 + 3y3 – y2 (5.11)
…
cn-3 = b n-2 – bn-3= yn – 3yn-1 +3 yn-2 – yn-3 (5.12)
Siguiendo este proceso se definen las cuartas, quintas, etc., diferencias hacia adelante.
Todas las diferencias pueden arreglarse en una tabla de diferencias (ver Tabla 5.2), en donde
cada diferencia se indica entre los dos elementos que la producen, ver la tabla siguiente:
Tabla 5.2. Tabla de diferencias
xi yi ∆yi ∆2yi ∆3yi …

x0 y0
a0 = y1 – y0
x1= x0+h y1 b 0 = a1 – a 0
a1 = y2 – y1 c 0 = b1 – b 0
x2= x0+2h y2 b 1 = a2 – a 1
a2 = y3 – y2 c 1 = b2 – b 1
x3= x0+3h y3 b 2 = a3 – a 2
… … … c n - 3 = b n - 2 – bn - 3
… … … bn - 2 = an-1 – an - 2
an-1 = yn – yn -1
xn= x0+nh yn
57
Cálculo Numérico
Si una de estas diferencias se vuelve constante (o aproximadamente constante), puede

decirse que los valores tabulados pueden describirse por un polinomio de grado igual al orden de
la diferencia constante (o aproximadamente constante).
5.1.2. Diferencias divididas
Si se considera la función y=f(x) definida en forma tabular, parecida a la presentada en la

Tabla 5.1, pero sin que los valores de la variable independiente tengan paso constante, puede
escribirse un polinomio de grado n-ésimo que pase por todos los (n+1) puntos definidos de la fun-
ción tal como:
Pn ( x ) = a 0 + ( x − x 0 )a1 + ( x − x 0 )( x − x1 )a 2 + L + ( x − x 0 )( x − x 1 )L( x − x n −1 )a n (5.13)
Los coeficientes ai pueden determinarse fácilmente si se utilizan las diferencias divididas

de los valores tabulados.
La diferencia dividida de orden cero se define como:
f [xr]= f(xr) (5.14)
esta diferencia se puede denotar también como yr ó bien fr.
La diferencia de primer orden o diferencia de orden uno es igual a:
f(x s ) − f(xr )
f [x r , x s ] = (5.15)
x s − xr
Las diferencias de orden superior se definen en términos de diferencias de orden inferior,

por ejemplo una diferencia de orden n se define como:
f [ x 1,K, x n ] − f [ x 0 , x1,K, x n −1 ]
f [ x 0 , x 1,K, x n ] = (5.16)
xn − x0
Una vez definidas las diferencias, los coeficientes se determinan como sigue:
a0=f(x0) (5.17)
f ( x1 ) − f ( x 0 )
a1=f[x0,x1] = (5.18)
x1 − x 0
f [ x1, x 2 ] − f [ x 0 , x1 ] [ f ( x 2 ) − f ( x1 )] − [ f ( x1 ) − f ( x 0 )]
a2=f[x0, x1, x2] = = (5.19)
x2 − x0 x2 − x0
f [ x 1,K, x n ] − f [ x 0 , x1,K, x n −1 ]
an= f [ x 0 , x 1,K, x n ] = (5.20)
xn − x0
58
Cálculo Numérico
Por lo tanto si se reemplaza en el polinomio dado por (5.13) se llega a que:
Pn ( x ) = f ( x 0 ) + ( x − x 0 )f [ x 0 , x1 ] + ( x − x 0 )( x − x1 )f [ x 0 , x1, x 2 ] + L +
(5.21)
+ ( x − x 0 )( x − x1 )L ( x − x n −1 )f [ x 0 ,L, x n ]
Este polinomio se conoce como polinomio de interpolación con diferencias divididas de

Newton.
5.1.3. Interpolación con incrementos constantes. Interpolación de Newton
Si se desea encontrar un valor incluido entre dos valores consecutivos de una función ta-
bular puede utilizarse la interpolación de Newton. Para ello se utilizan las diferencias finitas defini-
das en el ítem 5.1.1.
Dada la función y = f(x), definida en la Tabla 5.1, para encontrar un valor de x incluido entre
dos valores consecutivos de la tabla mencionada, x k < x< x k+1, se supone que la función f(x) se
aproxima a un polinomio Pn(x) de grado n, que pasa por todos los puntos que definen a la función
(puesto que la diferencia de orden n es aproximadamente constante). Recordando la definición de
diferencias hacia adelante y que las diferencias de orden superior se definen en función de las
diferencias de orden inferior pueden calcularse los valores de la variable dependiente y en función
de estas diferencias como se indica a continuación:
y1 = y0 + a0 (5.22)
y2 = y1 + a1= (y0 + a0) + (a0 + b0) = y0 + 2a0 + b0 (5.23)
y3 = y2 + a2= y0 + 2a0 + b0 + (a1 + b1)= (y0 + 2a0 + b0) + (a0 + b0+ b0+c0)= y0 + 3a0 + 3b0+c0 (5.24)
En estas expresiones puede verse que aparecen las primeras de las distintas diferencias
de órdenes sucesivos a partir de y0, afectadas por los coeficientes del desarrollo del binomio de
Newton. Suponiendo que esto es verdadero para cualquier valor de y, puede establecerse que:
k(k − 1) k(k - 1)(k - 2) k(k − 1)(k − 2)(k − 3)

y k = y 0 + ka 0 + b0 + c0 + d0 + ... (5.25)
2! 3! 4!
Y como: a0=∆y0, b0=∆2y0, c0=∆3y0, d0=∆4y0…, puede escribirse:
k(k − 1) 2 k(k - 1)(k - 2) 3 k(k − 1)(k − 2)(k − 3) 4

y k = y 0 + k∆y 0 + ∆ y0 + ∆ y0 + ∆ y 0 + .. . (5.26)
2! 3! 4!
Esta fórmula es verdadera para todo valor entero positivo de k, se denomina fórmula de interpo-
lación de Newton y es aplicable para cualquier valor de xk correspondiente o no a la tabla. En
esta fórmula, yk es un valor aproximado (interpolado) de la función obtenida para x = xk; y0 es el
valor inicial de y, el cual se considera inmediato al valor que se trata de interpolar; ∆y0, ∆2y0, ∆3y0,
59
Cálculo Numérico
∆4y0, …, son las diferencias hacia adelante de órdenes sucesivas correspondientes a y0; y k se
determina como sigue:
x −x 0
xk = x0 + kh ⇒ k = k (5.27)
h
Si se considera el polinomio interpolante de Newton en función de las diferencias divididas,

el error para un polinomio de grado n es:
f [ x n + 1, x n , x n −1,K, x 0 ]( x − x 0 )( x − x1 )( x − x 2 )L ( x − x n ) (5.28)
El algoritmo utilizando estas diferencias puede escribirse como:
Algoritmo de Interpolación de Newton

n+1: número de datos o puntos definidos de la función tabular
k: subíndice que indica el número de orden de los datos, varía de 2 hasta n
xk: valores de la variable independiente
yk: valores de la variable dependiente
xs: valor de la variable para el que se desea conocer el valor de f(xs)=ys
D: diferencias divididas
Paso 1: Tomar j= 2, ..., n
Paso 2: Para k= j, ..., n, calcular las diferencias divididas con la fórmula:
Dk, j −1 − Dk −1, j −1
D k, j =
xk − xk − j
Paso 3: Para j=1, ..., n hacer C=D(1,j)

Paso 4: Hacer Pa=1
Paso 5: Para j=1, ..., n, tomar Pa= Pa*(xs-xj)
Paso 6: Tomar ys=0
n
Paso 7: Para j= 1, ..., n , calcular ys como ys= ys + ∑ C* Pa
j =1
Paso 9: SALIDA ys
PARAR
¾ Por ejemplo si se desea encontrar el valor de la variable independiente para x=6,2 de la fun-
ción definida por la tabla
x 0 4 8 12 16 20
y 4 16 124 424 1012 1984
60
Cálculo Numérico
Se calculan las diferencias hacia adelante:
x y ∆y ∆2y ∆3y
0 4
4 16 12
8 124 108 96
12 424 300 192 96
16 1012 588 288 96
20 1984 972 384 96
Como puede observarse las diferencias de orden 3 ó terceras diferencias se mantienen
constantes, por lo tanto la función puede describirse por un polinomio de tercer grado. Se aplica la
fórmula (5.26) para encontrar el valor deseado. Siendo:
x0=4, y0=16; xk=6,2; h=4; ∆y=108 ; ∆2y=192 ; ∆3y=96; k=0,55
0,55(0,55 − 1) 0,55(0,55 - 1)(0,55 - 2)

y k = 16 + 0,55.108 + 192 + 96 = 62,1
2! 3!
O bien puede resolverse el problema encontrando el polinomio de interpolación de New-

ton, para ello se calculan las diferencias divididas hasta el orden quinto, es decir:
x y f[x0,x1] f[x0,x1,x2] f[x0,x1,x2, x3] f[x0,x1,x2, x3,x4] f[x0,x1,x2, x3,x4]

0 4
4 16 3
8 124 27 3
12 424 75 4 0,083
16 1012 147 4,5 0,031 -3,25 x10-3
20 1984 243 4,8 0,015 -8 x10-4 -2,025 x10-4
El polinomio interpolador es un polinomio de quinto grado:
P(x) = 4 + 3(x-0) + 3(x-0)(x-4) + 0,083(x-0)(x-4)(x-8) - 3,25 x10-3(x-0)(x-4)(x-8)(x-12) +

- 2,025 x10-4(x-0)(x-4)(x-8)(x-12)(x-16)
Por lo tanto si x=6,2, se reemplaza en el polinomio de Newton y el valor de la variable dependiente
es: 61,1.
5.1.4. Interpolación con incrementos variables. Interpolación de Lagrange

Este método se utiliza para funciones tabulares en las cuales los valores de x no son equi-
distantes. Para realizar la interpolación, se busca un polinomio que pase por todos los puntos. Si
se tienen n puntos el polinomio debe ser de grado n-1, o sea:
y = a 0 x n −1 + a1x n − 2 + L + a n − 2 x + a n −1 (5.29)
Este polinomio puede escribirse en la forma:
y = A1(x-x2) (x-x3)… (x-xn) + A2(x-x1) (x-x3)… (x-xn) + . . . + An(x-x1) (x-x2)… (x-xn-1) (5.28)
61
Cálculo Numérico
el grado del polinomio es n-1. Los coeficientes A0, A1, A2,… An se determinan de manera que la
gráfica del polinomio pase por todos los puntos especificados.
Si x = xn, se tiene y = yn, entonces reemplazando en la fórmula (5.28) se llega a que:
yn = An(xn – x1) (xn – x2) (xn – x3) . . . (xn – xn-1) (5.30)
y despejando el coeficiente An:
yn
An = (5.31)
( x n − x1 )( x n − x 2 )( x n − x 3 )L( x n − x n −1 )
Si se sustituyen los coeficientes dados por la ecuación (5.31) en la ecuación (5.29) se tiene
la fórmula de interpolación de Lagrange:
( x − x 2 )( x − x 3 )L ( x − x n ) ( x − x 1 )( x − x 3 )L ( x − x n )
y= y1 + y2 +
( x 1 − x 2 )( x 1 − x 3 )L ( x 1 − x n ) ( x 2 − x1 )( x 2 − x 3 ) L ( x 2 − x n )
(5.32)
( x − x1 )( x − x 2 )L ( x − x n −1 )
+L+ yn
( x n − x1 )( x n − x 2 )L ( x n − x n −1 )
Esta fórmula suele expresarse también como:
n n
Pn −1( x ) = ∑ L j ( x )f ( x i ) = ∑ L j ( x )y i (5.33)
j =1 j =1
n x − xk
con: L j ( x ) = ∏
k =1 x j − x k
k≠j
La fórmula de error para el polinomio de Lagrange de grado n-1es similar a la fórmula de

error para el polinomio de Taylor salvo que contiene un producto de n términos:
f (n) ( ξ( x ))
( x − x 1 )( x − x 2 )L ( x − x n −1 ) (5.34)
n!
El algoritmo para la interpolación de Lagrange es:
Algoritmo de Interpolación de Lagrange

n: número de datos o puntos definidos de la función tabular
x: valores de la variable independiente
y: valores de la variable dependiente
xs: valor de la variable para el que se desea conocer el valor de f(xs)=ys
L: es el coeficiente polinómico de Lagrange
Paso 1: Tomar yi = 0
Paso 2: Para i= 1, ..., n+1, tomar L =1
62
Cálculo Numérico
continua
Paso 3: Para j = 1, ..., n+1, si i≠j calcular el coeficiente polinómico mediante la

xs − x j
fórmula: L j = L
xi − x j
n +1
Paso 4: Calcular ys mediante la fórmula ys=yi+ ∑ Ljyi
j =1
Paso 5: SALIDA ys
PARAR
¾ Por ejemplo si se deseara encontrar el valor de la función tabular dada por la tabla para x=4
aplicando la interpolación de Lagrange:
x 2 8 12
y 0,69315 2,07944 2,48491
Se utiliza interpolación de orden 2 puesto que se tiene 3 puntos definidos:
( x − x 2 )( x − x 3 ) ( x − x 1 )( x − x 3 ) ( x − x 1 )( x − x 2 )
P2 ( x ) = y1 + y2 + y 3 y reemplazan-
( x 1 − x 2 )( x 1 − x 3 ) ( x 2 − x 1 )( x 2 − x 3 ) ( x 3 − x 1 )( x 3 − x 2 )
do por los valores correspondientes:
( 4 − 8)( 4 − 12) ( 4 − 2)( 4 − 12) ( 4 − 2)( 4 − 8)

P2 ( 4) = 0,69315 + 2,07944 + 2,48491 = 1,25899
(2 − 8)(2 − 12) (8 − 2)(8 − 12) (12 − 2)(12 − 8)
5.1.5. Interpolación inversa
El problema de interpolación inversa consiste en determinar el valor de la variable inde-

pendiente x conocido el valor de la función f(x). Se resuelve utilizando la fórmula de interpolación
de Lagrange y formando una tabla con los valores de la variable dependiente como valores de x y
los de la independiente como los de y.
5.1.6. Derivación numérica
Dada la función y = f(x) se trata de calcular el valor de sus derivadas en algunos de los
puntos x = x0, x1, x2, …, xn. Si se acepta aproximar la función con la fórmula de interpolación de
Newton (5.26):
k(k − 1) 2 k(k - 1)(k - 2) 3 k(k − 1)(k − 2)(k − 3) 4

y k = y 0 + k∆y 0 + ∆ y0 + ∆ y0 + ∆ y 0 + ...
2! 3! 4!
Derivando ambos miembros de la ecuación (5.26) con respecto a x y teniendo en cuenta que el
segundo miembro es una función compuesta de x, se tiene:
63
Cálculo Numérico
d d k2 − k 2 k 3 − 3k 2 + 2k 3 dk
f(x) = [y 0 + k∆y 0 + ∆ y0 + ∆ y0 + L ] (5.35)
dx dk 2! 3! dx
dk 1
y de (5.27), = , entonces la ecuación (5.35) queda reducida a:
dx h
d 1 2k − 1 2 3k 2 − 6k + 2 3
f ( x ) = [∆y 0 + ∆ y0 + ∆ y0 + L ] (5.36)
dx h 2 6
Derivando esta ecuación con respecto a x, se obtiene la segunda derivada.
d2 1 d 2k − 1 2 3k 2 − 6k + 2 3 dk
f ( x) = [∆y 0 + ∆ y0 + ∆ y0 + L ] por lo que:
dx h dk 2 6 dx
d2
2
f ( x) =
1
2
[∆2 y 0 + (k − 1)∆3 y 0 + L ] (5.37)
dx h
Si se deriva de nuevo con respecto a x, se obtiene la tercera derivada:
d3
3
f ( x) =
1
3
[∆3 y 0 + L] (5.38)
dx h
Considerando que si x=x0 entonces k=0, se puede escribir una fórmula muy simple para
estimar la primera derivada de la función en ese punto xi como se indica a continuación:
d 1⎡ 1 1 ∆n y 0 ⎤
f ( x ) = f ' ( x i ) = ⎢∆y 0 − ∆2 y 0 + ∆3 y 0 + L ( −1)n −1 ⎥ (5.39)
dx h ⎣⎢ 2 3 n ⎦⎥
Si se utiliza el primer término de la fórmula (5.39) esto significa que la interpolación es line-
al y la fórmula, considerando el error (último término) es:
1 1
f ' (xi ) = [∆y 0 ] − h∆2 f ( 2) (ξ) (5.40)
h 2
Con dos términos se tiene:
1⎡ 1 1
f ' (xi ) = ∆y 0 − ∆2 y 0 ⎤⎥ + h 2 f (3 ) (ξ) (5.41)
h ⎢⎣ 2 ⎦ 3
La fórmula de la ecuación (5.39) se denomina aproximación de diferencias hacia ade-

lante debido a que todas las diferencias se realizan con valores de la función que están delante
de yi. Si se utilizan diferencias entre valores de la función que están antes de yi, se tienen la
aproximación de diferencias hacia atrás y si se utilizan intervalos de diferencias en los que el
valor considerado se encuentra en el centro del intervalo se tiene la aproximación de diferencias
centrales. En la Tabla 5.1 se presentan fórmulas para calcular derivadas de primer a tercer orden
64
Cálculo Numérico
utilizando uno a cuatro términos de la ecuación (5.39) con diferencias hacia adelante, hacia atrás y
centrales. Para facilitar la escritura se hace f(xi) = fi.
Tabla 5.3. Fórmulas para calcular derivadas mediante aproximaciones de diferencia

Primera derivada
Aproximaciones de diferencia hacia adelante
f −f 1
fi' = i +1 i + E, E ≈ − h fi' '
h 2
− f i + 2 + 4 fi + 1 − 3 f i 1
fi' = + E, E ≈ − h 2 fi' ' '
2h 3
2fi + 3 − 9fi + 2 + 18fi +1 − 11fi 1
fi' = + E, E ≈ − h 3 fi' ' ' '
6h 4
Aproximaciones de diferencia hacia atrás
f −f 1 ''
fi' = i i −1 + E, E≈ h fi
h 2
3fi − 4fi −1 + fi − 2 1 2 '''
fi' = + E, E≈ h fi
2h 3
11fi − 18 fi −1 + 9fi − 2 − 2fi − 3 1 3 ''''
fi' = + E, E≈ h fi
6h 4
Aproximaciones de diferencias centrales
f −f h2 '''
fi' = i +1 i −1 + E, E≈− fi
2h 6
− fi + 2 + 8 f i + 1 − 8 f i − 1 + f i − 2 1 4 (v )
fi' = + E, E≈ h fi
12h 30
Segunda derivada
f − 2fi +1 + fi
fi' ' = i + 2 + E, E ≈ −h fi' ' '
h2
− fi + 3 + 4fi + 2 − 5fi +1 + 2fi 11 2 ' ' ' '
fi' ' = + E, E≈ h fi
h2 12
f − 2 fi − 1 + fi − 2
fi' ' = i + E, E ≈ h fi' ' '
h2
− fi + 4fi −1 − 5fi − 2 + 2fi − 3 11 2 ' ' ' '
fi' ' = + E, E≈ h fi
2 12
h
f − 2 fi + fi − 1 h2 ''''
fi' ' = i +1 + E, E≈ fi
h2 12
65
Cálculo Numérico
− fi + 2 + 16fi + 1 − 30fi + 16 fi −1 − fi − 2 1 4 ( vi)

fi' ' = + E, E≈ h fi
2 90
12h
Tercera derivada
f − 3fi + 2 + 3fi +1 − fi 3
fi' ' ' = i + 3 + E, E ≈ − h 2 fi' ' ' '
h3 2
f − 3fi −1 + 3fi − 2 − 2fi − 3 3 2 ''''
fi' ' ' = i + E, E≈ h fi
h3 2
f − 2fi + 1 + 2fi − 1 − 2 fi − 2 1
fi' ' = i + 2 + E, E ≈ − h 2 fi( v )
2h 3 4
¾ Por ejemplo si se deseara encontrar el valor de la segunda derivada para x=4 para la función
tabular:
x 2 3 4 5
y 3,010 4,771 6,021 6,990
Se utiliza la fórmula dada en la tabla 5.1:

fi+1 − 2fi + fi−1 6,990 − 2 ∗ 6,021 + 4,771
fi'' = = = −0,281
h2 1
5.1.7. Integración numérica
Dada la función f(x) cuya aproximación está dada por la ecuación (5.26) que pasa por los
n+1 puntos pivotes para los que x = x0, x1, x2,…, xn, todos ellos igualmente espaciados, entonces
se podrá tener una aproximación a la integral de la función:
xn xn
⎡ k(k − 1) 2 k(k − 1)(k − 2) 3 ⎤
∫ f ( x )dx = ∫ ⎢ y 0 + k∆y 0 + ∆ y0 + ∆ y 0 + L⎥dx (5.42)
x0 x0 ⎣ 2! 3! ⎦
Haciendo un cambio de variable en el que x es igual a x0+kh, y diferenciando, se tiene que

dx =hdk. Además si x= x0 entonces k=0 y si x = xn = x0 +nh, entonces k=n. Sustituyendo estos va-
lores en la ecuación (5.42) se llega a que:
xn n⎡ k2 − k 2 k 3 − 3k 2 + 2k 3 ⎤
∫ f ( x )dx = ∫ ⎢ y 0 + k∆y 0 + ∆ y0 + ∆ y 0 + L⎥hdk =
x0 0 ⎢⎣ 2 6 ⎥⎦
0 (5.43)
⎡ 2 ⎛ k3 k2 ⎞ 2 ⎛ 4 3 2⎞ ⎤
k ⎟∆ y 0 + ⎜ k − k + k ⎟∆3 y 0 + L⎥ n
= h ⎢ky 0 + ∆y 0 + ⎜ −
2 ⎜ ⎟ ⎜ ⎟
⎣⎢ ⎝ 6 4 ⎠ ⎝ 24 6 6 ⎠ ⎦⎥
Por lo tanto, reemplazando los límites se llega a que:
66
Cálculo Numérico
xn ⎡ n2 ⎛ n3 n2 ⎞ 2 ⎛ 4 3 2⎞ ⎤
⎟ ∆ y 0 + ⎜ n − n + n ⎟∆3 y 0 + L⎥ (5.47)
∫ f ( x )dx = h ⎢ny 0 + ∆y 0 + ⎜ −
⎢⎣ 2 ⎜ 6 4 ⎟⎠ ⎜ 24 6 6 ⎟⎠ ⎥⎦
x0 ⎝ ⎝
A continuación se presentan algunas reglas de aproximación a la integral de una función:

regla del trapecio, regla de Simpson1/3 y regla de Simpson de 3/8.
5.1.7.1. Regla trapecial
Si la interpolación se limita al primer orden, y la integral solo se calcula entre los dos prime-
ros valores de x, es decir entre x0 y x1 (n=1), aplicando (5.47), se obtiene:
x1
12 x1
h
∫ f ( x )dx = h( y 0 + ∆y 0 ) y como ∆y 0 = y 1 − y 0 , ∫ f ( x )dx = ( y 0 + y1 ) (5.48)
x0 2 x 2
0
Geométricamente, el primer miembro de la expresión anterior corresponde al área A1 bajo

la curva y= f(x) entre las rectas x =x0 y x = x1 y el eje x. El segundo miembro, que es el valor
aproximado de la integral representa al área del trapecio formado por las tres rectas antes men-
cionadas y la que une los puntos pivotes usados en la interpolación lineal. De manera semejante
a como se obtuvo la expresión (5.48) para el área elemental A1, se pueden obtener expresiones
xn
h
para A2, A3, …, An (ver Figura 5.1). Por ejemplo, para An se tiene: ∫ f ( x )dx = ( y n −1 + y n ) .
x n −1 2
Sumando miembro a miembro las expresiones individuales se obtiene que:
xn
h
∫ f ( x )dx = (y 0 + y n + 2( y1 + y 2 + y 3 + L + y n −1) (5.48)
x n −1 2
El segundo miembro de esta expresión, que proporciona un valor aproximado a la integral

del primer miembro, recibe el nombre de fórmula o regla trapecial, la cual se expresa en la for-
h
ma: A 1 = [y 0 + y n + 2 ∑ resto de las ordenadas] (5.49)
2
2
A1 An
A3 A4
A2
x0 x1 x
Figura 5.1 Aproximación al área bajo la curva con la fórmula trapecial
67
Cálculo Numérico
Para mejorar la exactitud de la regla del trapecio, se puede dividir el intervalo de integra-
ción en h segmentos y aplicar el método a cada uno de ellos. El ancho de segmento se calcula
x − x0
como h = n y la regla se denomina trapecial de segmentos múltiples.
n
Algoritmo para Integración por la regla trapecial de segmentos múltiples

m: número de segmentos
h: ancho de segmentos
x0: límite inferior de integración
xn: límite superior de integración
yi: valor de la variable dependiente
SU: valor de y para la subrutina de cálculo
A: área del segmento
IN: valor de la integral
Paso 1: Ingresar n
Paso 2: Hacer m=n-1
Paso 3: Ingresar x0,xn
Paso 4: Calcular h como h = x n − x 0

m
Paso 5: Ingresar y1
Paso 6: Hacer SU= y1
Paso 7: Para i=2 hasta m hacer SU=SU+2*yi
Paso 8: Calcular el área mediante la fórmula A=(SU+yn)/2*m
Paso 9: Calcular IN mediante la fórmula IN = ( x n − x 0 ) * A
Paso 10: la SALIDA es IN

PARAR
5.1.7.2. Regla de Simpson 1/3
Si la interpolación es limitada de segundo orden y la integral solo se calcula entre los tres
primeros valores de x (n=2) en (5.47) se obtiene:
x2 ⎡ 22 ⎛ 23 2 2 ⎞ 2 ⎛ 4 3 2⎞ ⎤
⎟∆ y 0 + ⎜ 2 − 2 + 2 ⎟∆3 y 0 + L⎥ (5.50)
∫ f ( x )dx = h ⎢2y 0 + ∆y 0 + ⎜ −
2 ⎜ ⎟ ⎜ ⎟
x0 ⎣⎢ ⎝ 6 4 ⎠ ⎝ 24 6 6 ⎠ ⎦⎥
y como ∆y0 = y1 - y0 y ∆2y0 = y2 – 2 y1 + y0 se tiene:
x2
⎡ 1 ⎤ h
∫ f ( x )dx = h ⎢2y 0 + 2( y1 − y 0 ) + ( y 2 − 2y 1 + y 0 )⎥ = ( y 0 + 4 y 1 + y 2 ) , en forma general:
x0 ⎣ 3 ⎦ 3
68
Cálculo Numérico
xn
h
∫ f ( x )dx = 3 ( y n − 2 + 4y n −1 + y n ) (5.51)
x 0−2
Sumando miembro a miembro se obtiene:
xn
h
∫ f ( x )dx = [y 0 + y n + 2( y 2 + y 4 + y 6 + L + y n − 2 ) + 4( y1 + y 3 + y 5 + L + y n −1 )] (5.52)
x0 3
El segundo miembro de la ecuación se denomina fórmula de Simpson del 1/3 y se expre-

sa como:
h
A1 = [y 0 + y n + 2 ∑ ordenadas de orden par + 4 ∑ ordenadas de orden impar ] (5.53)
3
3
En la Figura 5.2 se muestra un esquema de aplicación de esta fórmula:
A1
x1 x2 x
Figura 5.2 Aproximación al área bajo la curva con la fórmula

de Simpson de 1/3
5.1.7.3. Regla de Simpson 3/8
Si la interpolación es de tercer orden y se toma n =3, se obtiene:
x3 ⎡ 32 ⎛ 33 32 ⎞ 2 ⎛ 4 3 2⎞ ⎤
⎟∆ y 0 + ⎜ 3 − 3 + 3 ⎟∆3 y 0 + L⎥
∫ f ( x )dx = h ⎢3 y 0 + ∆y 0 + ⎜ − (5.54)
2 ⎜ ⎟ ⎜ ⎟
x0 ⎢⎣ ⎝ 6 4 ⎠ ⎝ 24 6 6 ⎠ ⎥⎦
x3
⎡ ⎤
∫ f ( x )dx = h ⎢⎣3y 0 + 2 ( y1 − y 0 ) + 4 ( y 2 − 2y1 + y 0 ) + 8 (y 3 − 3y 2 + 3y1 − y 0 )⎥⎦ =
9 9 3
x0
3
=
h ( y 0 + 3y 1 + 3y 2 + y 3 )
8
en forma general:
xn
3
∫ f ( x )dx = h( y n − 3 + 3 y n − 2 + 3 y n −1 + y n ) (5.55)
x0 −3 8
Sumando miembro a miembro:
69
Cálculo Numérico
xn
3 ⎡ y + y + 2( y + y + y + L + y − 3 ) + ⎤
∫ f ( x )dx = h ⎣⎢+03( y n+ y + 3y + 6y + 9y + y +nL + +
(5.56)
x0 8 1 2 4 5 7 8 y n−2 y n −1 ⎥⎦
)
El segundo miembro de la ecuación se denomina fórmula de Simpson del 3/8 y se expresa co-
mo:
3
A3 = h[y 0 + y n + 2 ∑ ordenadas de orden múltiplo de 3 + 3 ∑ resto de ordenadas ] (5.57)
8
8
En la Figura 5.3 se muestra un esquema de aplicación de esta fórmula:
x0 x1 x2 x3 x
Figura 5.3 Aproximación al área bajo la curva con la fórmula

de Simpson de 3/8
¾ Por ejemplo si se deseara encontrar el valor aproximado de la integral de la función tabular:
x 2 4 6 8 10 12 14 16
y 8 4 2 1 2 6 9 10
Utilizando la regla trapecial y las de Simpson se obtiene:
2[8 + 10 + 2( 4 + 2 + 1 + 2 + 6 + 9)]
A1/ 2 = = 66
2
2[8 + 10 + 2(2 + 2 + 9) + 4( 4 + 1 + 6)]

A1/ 3 = = 58,66
3
2.3[8 + 10 + 2(1 + 9) + 3( 4 + +2 + 2 + 6)]

A3 / 8 = = 60
8
70
Cálculo Numérico
5.2. APROXIMACIÓN FUNCIONAL
En este tipo de aproximación se trata de encontrar la ecuación de una curva que, aunque
no pase por todos los puntos, tenga pocas variaciones, es decir sea suave y pase lo más cerca
posible de todos ellos, para ello es necesario aplicar el criterio de mínimos cuadrados. Antes de
aplicar este criterio, debe escogerse la forma de la curva que se va a ajustar al conjunto de puntos
dado y su ecuación puede obtenerse desde un conocimiento previo del problema, es decir por su
interpretación física o en forma arbitraria observando que ecuación conocida describe aproxima-
damente a esta curva.
Dentro de la aproximación funcional por mínimos cuadrados se distinguen tres tipos de re-
gresión: lineal, polinomial, lineal múltiple.
5.2.1. Regresión lineal
El ejemplo más simple de aproximación por mínimos cuadrados es el ajuste de un conjun-

to de datos a una línea recta.
La expresión matemática de una recta es:
y = a0 + a1x + E (5.58)
en donde a0 y a1 son coeficientes que representan la intersección con el eje de las ordenadas y la
pendiente, respectivamente y E es el error o residuo entre el modelo y las observaciones. Reorde-
nando, se puede calcular el error como
E = y - a0 - a1x (5.59)
es decir, es la diferencia entre el valor real de y y el valor aproximado, a0 + a1x que predice la
ecuación lineal.
Una forma de obtener un mejor ajuste es minimizar la suma de cuadrados de los residuos,
Sr, de la siguiente manera:
n n
Sr = ∑ Ei2 = ∑ ( yi − a0 − a1x i )2 (5.50)
i =1 i =1
Para encontrar los valores de a0 y a1 que minimicen la ecuación (5.50) se debe derivar esta
ecuación con respecto a los coeficientes indicados, es decir:
∂S r n
= −2 ∑ ( y i − a 0 − a1x i ) (5.51)
∂a 0 i =1
∂S r n
= −2 ∑ ( y i − a 0 − a1x i )x i ] (5.52)
∂a1 i =1
71
Cálculo Numérico
Para generar un mínimo, se igualan estas derivadas a cero y se expresan como un con-
junto de dos ecuaciones lineales con dos incógnitas a0 y a1.
na0 + Σ xia1= Σ yi (5.53)
Σxia0 + Σ xi2a1= Σ xi yi
Si se resuelve este sistema se obtiene:
n∑ x i y i − ∑ x i ∑ y i
a1 = (5.54)
n∑ x i 2 − ( ∑ x i ) 2
a0 = y − a1x (5.55)
donde: y y x son las medias aritméticas de y y x respectivamente.
El error estándar de aproximación Sy/x, que indica el error para los valores predichos de y
correspondientes a los valores particulares de x y permite cuantificar la dispersión alrededor de la
línea de regresión, se calcula mediante la siguiente ecuación:
Sr
Sy / x = (5.56)
n−2
Para cuantificar la eficiencia del ajuste, que es particularmente útil en la comparación de

varias regresiones, se utilizan el coeficiente de determinación r2 y el de correlación r, que es la raíz
cuadrada del coeficiente anterior.
El coeficiente de determinación se calcula como sigue:
S t − Sr
r2 = (5.57)
St
donde: S t = ∑ ( yi − y ) 2 es la cantidad de dispersión en la variable dependiente que existe antes de

la regresión.
La diferencia entre las dos cantidades cuantifica la mejora en la reducción del error debido
al modelo de la línea recta.
Para un ajuste perfecto Sr = 0 y r2 =1, así la línea recta explica un 100 % de la variabilidad.
El algoritmo para regresión lineal es:
72
Cálculo Numérico
Algoritmo para Regresión lineal

xi: valor de la variable independiente (i varía de 1 a n)
yi: valor de la variable dependiente
Sx: suma de los valores de xi
Sy: suma de los valores de yi
X^2: suma de cuadrados de xi
Y^2: suma de cuadrados de yi
XY: suma del producto
xm: valor medio de x
ym: valor medio de y
a1: pendiente
a0: término independiente
Paso 1: Para i= 1, ..., n, seguir los pasos 2 a 4
Paso 2: calcular la suma de los valores de x como: Sx= Sx+xi
Paso 3: Para i= 1, ..., n, calcular la suma de los valores de y como: Sy= Sy+yi
Paso 4: Para i= 1, ..., n, calcular la suma de cuadrados de x como: X^2= xi+ xi xi
Paso 5: Para i= 1, ..., n, calcular la suma de cuadrados de y como: Y^2= yi+ yi yi
Paso 6: Hacer xm= Sx/n
Paso 7: Hacer ym= Sy/n
Paso 8: Calcular a1 mediante la fórmula: ym =(n*XY- Sx Sy) /(n* X^2- Sx Sy)
Paso 9: Calcular a0 mediante la fórmula: a0= ym- a1* xm
Paso 10: la SALIDA es a0 y a1
PARAR
5.2.2. Linealización de relaciones no lineales
En el caso de tener relaciones no lineales se pueden hacer transformaciones que expresen

los datos de manera que sean compatibles con la regresión lineal. A continuación se presentan
algunos ejemplos.
Modelo exponencial
y = d1e b1x (5.58)
en donde d1 y b1 son constantes. Para linealizar este modelo se aplican logaritmos naturales, es
decir:
ln y = ln d1 + b1 ln e (5.59)
73
Cálculo Numérico
En una gráfica semilogarítmica de ln y vs. x se genera una línea recta con pendiente b1 y
ordenada al origen: ln d1.
Modelo de ecuación elevada a una potencia
y = d 2 x b2 (5.60)
En donde d2 y b2 son coeficientes, puede linealizarse mediante logaritmos en base 10, o sea:
log y = b2 log x + log d2 (5.61)
de forma que en una gráfica logarítmica de log y vs log x se genera una línea recta con pendiente
b2 y ordenada al origen log d2.
Modelo de crecimiento a saturación
x
y = d3 (5.62)
b3 + x
Los coeficientes d3 y b3 son constantes y puede linealizarse si se invierte la ecuación (5.62), es

decir:
1 b3 1 1
= + (5.63)
y d3 x d3
Una gráfica de 1/y contra 1/x será lineal, con pendiente b3/d3 y ordenada al origen 1/d3.
5.2.3. Regresión polinomial
El procedimiento de mínimos cuadrados se puede extender fácilmente y ajustar los datos a

un polinomio de m- ésimo grado.
y = a 0 + a1x + a 2 x 2 + L + a m x m (5.64)
En este caso la suma de cuadrados de los residuos es:
n
S r = ∑ ( y i − a 0 − a1x i − a 2 x i 2 − L − a m x im ) 2 (5.65)
i =1
Siguiendo el procedimiento anterior, se deriva la ecuación con respecto a cada uno de los
coeficientes del polinomio, para obtener:
∂S r n
= −2 ∑ ( y i − a 0 − a1x i − a 2 x i 2 − L − a m x im ) (5.66)
∂a 0 i =1
∂S r n
= −2 ∑ x i ( y i − a 0 − a1x i − a 2 x i 2 − L − a m x im ) (5.67)
∂a1 i =1
74
Cálculo Numérico
...
∂Sr n
= −2 ∑ x im ( y i − a 0 − a1x i − a 2 x i 2 − L − am x im ) (5.68)
∂a m i =1
Si estas ecuaciones se igualan a cero y se reordenan se obtiene un conjunto de ecuacio-

nes normales:
na0 + a1Σ xi + a2 Σ xi2 + . . . +am Σ xim = Σ yi
a0Σ xi + a1Σ xi2 + a2 Σ xi3 + . . . +am Σ xim+1 = Σ xi yi
a0Σ xi2 + a1Σ xi3 + a2 Σ xi4 + . . . +am Σ xim+2 = Σ xi2 yi (5.69)

. . . . .
. . . . .
. . . . .
m m+1 m+2 2m
a0Σ xi + a1Σ xi + a2 Σ xi + . . . +am Σ xi = Σ xim yi
El error de regresión polinomial se calcula con:
Sr
Sy / x = (5.70)
n − (m + 1)
Y el coeficiente de determinación con
S v − Sr
r2 = (5.71)
Sv
5.2.4. Regresión lineal múltiple
Se utiliza si y es una función de dos o más variables. Por ejemplo si y es función lineal de
x1 y x2 de la forma:
y = a0 + a1x1+ a2 x2 (5.72)
La suma de los residuos es:
n
S r = ∑ ( y i − a 0 − a1x1,i − a 2 x 2,i ) 2 (5.73)
i =1
Si se deriva con respecto a cada uno de los coeficientes se obtiene:
∂S r n
= −2 ∑ ( y i − a 0 − a1x1,i − a 2 x 2,i ) (5.74)
∂a 0 i =1
∂Sr n
= −2 ∑ x1,i ( y i − a 0 − a1x1,i − a 2 x 2,i ) (5.75)
∂a1 i =1
75
Cálculo Numérico
∂S r n
= −2 ∑ x 2,i ( y i − a 0 − a1x1,i − a 2 x 2,i ) (5.76)
∂a 2 i =1
Los coeficientes que generan la suma mínima de los cuadrados de los residuos se obtie-
nen igualando a cero cada una de las derivadas parciales y expresando las ecuaciones como un
conjunto de ecuaciones lineales simultáneas, de la forma:
na0 + a1Σ x1,i + a2 Σ x2,i = Σ yi
a0Σ x1,i + a1Σ x1,i2 + a2 Σ x1,ix 2,i = Σ x1,i yi (5.77)
a0Σ x2,i2 + a1Σ x1,ix 2,i + a2 Σ x2,i2 = Σ x2,i yi
¾ Por ejemplo se desean ajustar los datos experimentales de la tabla que se presenta a conti-
nuación a un modelo exponencial.
x 1 2 2,5 4 6
y 0,4 0,7 0,8 1,0 1,2
La ecuación que describe estos datos es la (5.59):
y = d1e b1x
Aplicando logaritmos neperianos a ambos miembros de la ecuación se llega a que:
ln y = ln d1 + b1 ln e
Por lo tanto se puede llamar Y a ln y; a0 al ln d1 y a1 a b1, de manera que la ecuación con
esa transformación es una ecuación lineal:
Y=a0+a1x
En la tabla que sigue se presenta la transformación de los datos
x y Y=ln y x.Y x2
1 0,4 -0,916 -0,916 1
2 0,7 -0,365 -0,712 4
2,5 0,8 -0,223 -0,557 6,25
4 1 0 0 16
6 1,2 0,182 1,092 36
Σ 15,5 -1,313 -1,093 63,25
El valor medio de x es 3,1 y el de Y es -0,262
76
Cálculo Numérico
El cálculo de la pendiente en base a la fórmula (5.54) es:

5( −1,0935 ) − 15,5( −1,313 )
a1 = = 0,195
5(63,25) − (15,5) 2
a0= -0,262 – (0,195. 3,1)= -0,867
Como a0 = ln d1, entonces d1=e-0,867=0,420
Como a1 = b1, entonces b1es 0,195
La ecuación ajustada es:
y = 0,420 e 0,195 x
EJERCICOS PROPUESTOS
5.1. Obtenga, aplicando el método de diferencias finitas, el grado del polinomio que puede repre-
sentar al conjunto de valores de la siguiente tabla:
x 0 6 12 18 24 30
y 2 8 62 212 506 992
5.2. Utilice el método de interpolación que corresponda para:
a) conocer la producción de arroz parboilizado durante el año 95 si los datos de pro-

ducción en toneladas se presentan en la tabla que sigue:
Año 90 92 94 96 98
Producción (tn) 1229 1302 1100 1200 1220
b) determinar la conductividad térmica de las frutillas a -25 ºC si su conductividad en

W/m K a diferentes temperaturas es:
T (ºC) -40 -34 -30 -24 -20 -15 -10 -5

k(W/m K) 1.489 1.450 1.413 1.375 1.388 1.299 1.255 1.183
5.3. Para la función definida en la tabla de la distribución de velocidades de leche en una tubería,
calcule:
a) la primera derivada en x = 2, utilizando las fórmulas de interpolación limitadas al primero,

segundo y tercer orden.
b) La primera y segunda derivadas para x = 4, aplicando la interpolación de segundo orden.
X 2 3 4 5 6 7
y 6,010 8,144 12,022 14,990 16,781 20,451
5.4. Si se tiene una curva de crecimiento de Clostridium sporogenes en la que se representa la

velocidad de crecimiento en función del tiempo, descripta por los datos que se presentan en la
tabla, indique el valor de la pendiente de la curva en un tiempo t = 300 min.
77
Cálculo Numérico
t(min) 0 100 200 300 400

v
( N/min) x 10
6 0 30 45 50 50
5.5. Las dos primeras columnas de la tabla que se presenta a continuación indican las relaciones
tiempo temperatura en el centro de una lata de habas en salsa de tomates que se está esterili-
zando en un autoclave discontinuo. La tercera columna presenta el valor letal de las temperaturas
superiores a 100 ºC. Como el área bajo la curva de letalidad vs. tiempo corresponde al valor este-
rilizante de este tratamiento, determine si el mismo es suficiente considerando que el valor apro-
piado es 6 minutos.
t (min) T (ºC) L t (min) T (ºC) L
-7
0 54,44 2,208 x 10 14 116,77 0,378
-6
1 62,78 1,506 x 10 15 117,22, 0,419
-6
2 68,54 5,675 x 10 16 118,05 0,507
-5
3 75,55 2,851 x 10 17 118,61 0,577
-4
4 81,66 1,164 x 10 18 119,05 0,638
-4
5 87,22 4,188 x 10 19 119,5 0,708
-4
6 92,77 1,503 x 10 20 119,61 0,726
-3
7 97,77 4,753 x 10 21 117,77 0,475
8 102,22 0,013 22 108,33 0,054
-3
9 105,55 0,028 23 93,33 1,710 x 10
-5
10 108,33 0,054 24 80,00 7,943 x 10
-6
11 110,55 0,090 25 68,88 6,138 x 10
12 112,77 0,150
13 114,44 0,221
5.6. Determine la ecuación de la recta de regresión para los datos de la tabla que sigue que re-
presentan la cantidad de β-caroteno (mg) extraído de granos de maíz (kg). Grafique la función
tabular y el ajuste. Calcule la media, desviación estándar y el coeficiente de variación.
Maíz (kg) 0,1 0,4 0,5 0,7 0,8 0,9
β-caroteno (mg) 0,61 0,92 0,99 1,52 1,47 2,03
5.7. Ajuste a un polinomio de segundo orden los datos que siguen:
x 0 1 2 3 4 5
y 4,1 14,7 26,6 52,2 80,9 122,1
78
Cálculo Numérico
Grafique los datos y el polinomio de ajuste, calcule el error de aproximación y el coeficiente de

correlación
5.8. Considere que p representa a la población de levaduras del género Saccharomyces que se
utilizan en la elaboración de cerveza y que f es la concentración de la fuente de carbono que ser-
virá de sustrato a estos microorganismos. Las medidas de k (tasa de crecimiento específico en
-1
días contra f (mg/L) se muestran en la tabla. Determine los valores de kmáx (tasa máxima de cre-
cimiento para valores de f) y K (constante de semi-saturación), empleando la linealización del mo-
delo de promedio de crecimiento a saturación.
f
La ecuación que rige el modelo es: k = k máx
K +f
x 7 9 15 25 40 75 100 150
y 0,28 0,36 0,56 0,74 0,86 0,97 0,99 1,14
5.9. La relación entre la concentración de una sustancia orgánica en función del tiempo está ex-
presada por la siguiente tabla de datos experimentales:
Tiempo (min) 0,1 0,4 0,5 0,7 0,7 0,9
Concentración (mgL-1) 0,61 0,92 0,99 1,52 1,47 2,03
Establezca la relación empírica C=f (t) que mejor se ajuste a los datos y grafique los datos expe-
rimentales y el ajuste.
5.10. La relación entre el contenido en CO2 de un determinado vapor (W en ppm) y la conductivi-

dad específica del condensado (C en mhos cm-1) está dada por la siguiente tabla:
C (mhos cm-1) 0,49 1,56 2,32 3,11 3,74 4,44 5,36
W (ppm) 0,78 3,09 5,25 8,44 11,88 15,79 20,18
Determine la ecuación polinomial W=f (C) que se ajuste a los datos y grafique.
79
Cálculo Numérico
Capítulo 6
SIMULACION
El término simulación se refiere a la construcción de una representación simplificada de un

proceso o un sistema físico con el fin de facilitar su análisis, se caracteriza por el hecho de no in-
cluir todas las propiedades del sistema real; es decir, que el propósito de la simulación es mostrar
el efecto de ciertos factores particulares que se están investigando.
Se utiliza simulación si no existe una formulación matemática o métodos de resolución

analíticos (o son demasiado complejos); se desea experimentar con el sistema antes de su cons-
trucción e interesa controlar las condiciones del sistema y resulta imposible experimentar sobre el
mismo. Sus objetivos son observar, predecir, modificar y optimizar un sistema. De esta manera es
posible estudiar el comportamiento del sistema, imaginar situaciones para mejorarlo, con costos
bajos o poco peligrosos.
6.1. METODOLOGÍA DE SIMULACIÓN
Debido al carácter tan complejo de un sistema real, el estudio que se hace para obtener un
modelo o simulación se considera dividido en cuatro etapas principales, que son las siguientes:
1. Observación del sistema: representa la fase experimental del estudio en la que se hace
un análisis preliminar del sistema, con el fin de determinar la interacción con otros sistemas, sus
restricciones, las variables que interactúan dentro del mismo y sus interrelaciones, las medidas de
efectividad que se van a utilizar para definirlo y estudiarlo y los resultados que se esperan obtener
del estudio que fundamentarán el diseño del modelo.
2. Establecimiento del sistema matemático: de las conclusiones inferidas en la etapa ante-

rior se fijan los postulados sobre los que se fundamenta el modelo matemático y en general puede
consistir en un sistema de ecuaciones diferenciales o de diferencias. En la formulación del modelo
es necesario definir todas las variables que forman parte de él, sus relaciones lógicas y los dia-
gramas de flujo que describan en forma completa el modelo.
3. Solución del modelo matemático: esta proporciona información que predice el compor-
tamiento del sistema real, el cual sirve para conocer el grado de validez del modelo. Si no se pue-
den resolver las ecuaciones que forman el modelo matemático, ya sea por la cantidad de cálculos
matemáticos que requieren, por falta de capacidad de la computadora o por no poder obtenerse
80
Cálculo Numérico
exactamente la solución, se puede hacer un análisis cualitativo del modelo matemático con el fin
de resolverlo; el planteamiento y solución de estos modelos cualitativos es lo que constituye los
métodos de Montecarlo.
4. Verificación y ajuste del modelo: Los resultados obtenidos en la etapa anterior se com-
paran con los resultados del comportamiento que sigue el sistema real y con los cuales se esta-
bleció el modelo de la etapa 1.
6.1.1. Métodos de Montecarlo
Los métodos de Montecarlo abarcan una colección de técnicas que permiten obtener solu-
ciones de problemas matemáticos o físicos por medio de pruebas aleatorias repetidas. En la
práctica, las pruebas aleatorias se sustituyen por resultados de ciertos cálculos realizados con
números aleatorios.
El problema crucial de la aplicación de estos métodos es hallar los valores de una variable
aleatoria (discreta o continua) con una distribución de probabilidad dada por la función p(x) a partir
de una transformación de los valores de una variable aleatoria uniformemente distribuida en el
intervalo [0,1), proporcionada por la computadora o por una rutina incorporada al programa (esto
se denomina generación de números aleatorios).
6.1.1.1. Generación de números aleatorios
En los experimentos de simulación, es necesario aplicar un mecanismo que proporcione

números aleatorios. La forma de obtenerlos es variada, entre los primeros mecanismos se tenían
las ruletas, dados, cartas, etc., actualmente pueden generarse mediante técnicas que utilizan ta-
blas, computadoras analógicas y digitales. En la Tabla 6.1 se presenta un cuadro comparativo de
estas técnicas.
Tabla 6.1. Técnicas de generación de números aleatorios
CARACTERISTICAS
TECNICA DE VENTAJAS DESVENTAJAS
OPERACION
Es lenta y puede quedar limitada si
se trata de una muestra grande.
Los números han sido generados Se puede repetir la misma
TABLAS Los números con los que se trabaja
por algún otro medio. secuencia de números.
son siempre los mismos para cual-
quier problema.
Los números son verdade-
COMPUTADORAS Generan los números por algún No se puede repetir la misma se-
ramente aleatorios con alta
ANALOGICAS proceso físico. cuencia.
frecuencia.
La generación interna por
medio de alguna relación
En el caso de tablas la generación
Los números se generan por de recurrencia se puede
se vuelve lenta y ocupa un gran
medio de tablas con algún adita- considerar como una
espacio en la memoria y en el caso
COMPUTADORAS mento analógico que establezca operación normal de la
de aditamentos analógicos no se
DIGITALES algún proceso físico o bien por computadora y se pueden
puede reproducir una misma se-
generación interna utilizando una lograr grandes muestras
cuencia de números.
relación de recurrencia. de números en poco tiem-
po y es posible reproducir
el proceso.
81
Cálculo Numérico
Cabe aclarar que las relaciones de recurrencia producen secuencias en las que es posible
identificar cualquier número de antemano, estableciendo una secuencia determinística y perdien-
do su carácter aleatorio. Estas secuencias se consideran aleatorias si las sucesiones de números
generados tienen una distribución uniforme y son estadísticamente independientes y se denomi-
nan números pseudoaleatorios o números aleatorios producidos por métodos secuenciales.
6.1.1.2. Generación de Números pseudoaleatorios
Para generarlos pueden utilizarse los métodos de cuadrados centrales, de los productos
centrales y métodos congruenciales.
6.1.1.2.1. Método de los cuadrados centrales
Se parte de un número de n cifras que se eleva al cuadrado y del resultado obtenido, que
por lo general es un número de 2n cifras, se toman las n cifras de la parte central, que a su vez se
vuelven a elevar al cuadrado y así sucesivamente. La desventaja es que tiende a degenerar rápi-
damente hacia valores constantes dependiendo del número inicial que se elija.
¾ Por ejemplo se desean generar números entre 0 y 999 y se elige como número inicial a 204,
que genera la secuencia:
2042=41616→161
1612=25921→592
5922=350464→504
6.1.1.2.2. Método de los productos centrales
Se efectúa el producto entre dos números seguidos de la sucesión de números pseudoa-

leatorios y del resultado tomar los números del medio para obtener el número siguiente de la su-
cesión. También tiende a degenerar a valores constantes.
¾ Por ejemplo se desean generar 3 números aleatorios y se eligen como números iniciales a
1425 y 1915, que generan la secuencia:
(1425)(1915)=2728875→7288
(7288)(1915)=13956520→9565
(9565)(7288)=69709720→7097
6.1.1.2.3. Métodos congruenciales
Dentro de los métodos secuenciales, los más importantes son los de tipo congruencial y
responden a una expresión general de la forma:
Ui+1 = (αUi + k) mód m (6.1)
82
Cálculo Numérico
donde:
Ui+1: número aleatorio generado en el paso i+1

Ui: número aleatorio generado en el paso i
α: multiplicador
k: incremento
m: módulo
Significa que (αUi + k) se divide entre m y el residuo que se obtenga es el valor de Ui+1.
Estos números quedan limitados ya que U0 y k deben ser mayores o iguales a 0, pero menores
que m. El multiplicador α debe ser mayor que cero y menor que m.
Pueden distinguirse los siguientes métodos congruenciales:
a) Método congruencial mixto
Se lo denomina así si k es una constante diferente de cero. Para asegurar que el período
de una secuencia sea mayor a la cantidad de números aleatorios requeridos es importante la se-
lección del módulo, incremento, multiplicador y del valor inicial.
Se puede obtener una muestra con un período igual al módulo m si se cumplen las si-
guientes condiciones:
1. k y m son números primos relativos
2. α = pt +1, siendo p un factor primo de m y t cualquier entero no negativo.
3. En el caso de que 4 sea un factor de m, se aplica la condición:
α = 4t +1, donde t es cualquier entero no negativo.
Como las computadoras emplean números binarios o decimales, se utilizan criterios dife-
rentes para la selección de m, por ejemplo para computadoras decimales conviene que m valga
10d, siendo d el número de dígitos decimales que se tiene en una palabra. U0 puede ser cualquier
número decimal positivo, k un número impar positivo no divisible entre 5, y α= 20t +1, donde t es
cualquier entero no negativo o bien α =10s +1, con s>1.
b) Método congruencial multiplicativo
Es el más empleado, se obtiene cuando k =0, o sea:
Ui+1 = (αUi) mód m (6.2)
c) Método congruencial aditivo
En este caso α =1 y se tiene la fórmula:
Ui+1 = (Ui + Ui+s) mód m (6.3)
83
Cálculo Numérico
y si s es igual a 1, U0 = y U1 =1, se obtiene la secuencia de Fibonacci.
6.1.1.3. Aplicaciones de los métodos de Montecarlo
Como aplicaciones de estos métodos se pueden mencionar: paseo aleatorio, integración

por simulación y línea de espera.
6.1.1.3.1. Paseo aleatorio
Para explicar este método supongamos que se desea saber cual es la distancia D que re-
corre un niño en una plaza después de haber dado n pasos, cada paso tiene la misma longitud P,
variando exclusivamente la dirección hacia la que se mueve dentro de un ángulo α tal que 0º ≤ α ≤
360º, como se muestra en la Figura 6.1.
y1
y3
D
y2
P
y1
αi
x1 x2 x3
x
Figura 6.1. Paseo aleatorio
La distancia esperada se define como:
Desp= P n (6.4)
Para realizar la simulación de la distancia recorrida por el niño se genera un número alea-
torio Ui, para obtener el ángulo αi en radianes se multiplica al número por 2π.
Las componentes de cada paso dado se calculan mediante las expresiones:
xi = P cos αi (6.5)
yi = P sen αI (6.6)
La distancia después de n pasos, con la ecuación:
2 2
⎛n ⎞ ⎛n ⎞
D esp = ⎜⎜ ∑ X i ⎟⎟ + ⎜⎜ ∑ Yi ⎟⎟ (6.7)
⎝ i =1 ⎠ ⎝ i =1 ⎠
Un algoritmo para el paseo aleatorio es:
84
Cálculo Numérico
Algoritmo para el Paseo aleatorio

Para determinar la distancia esperada luego de un número N de simulaciones, para
número de pasos y longitud de paso determinadas
n: número de pasos
P: longitud del paso
N: número máximo de simulaciones
x y y: componentes del paso
Sx: sumatoria de las componentes x; Sy: sumatoria de las componentes y
α: ángulo simulado en radianes
U: número aleatorio de 0 a 1
D(j): distancia acumulada hasta el paseo j-ésimo
DE: distancia esperada que se calcula como P n

DM: distancia media obtenida por simulación
E: diferencia entre la DE y la distancia obtenida por simulación
Paso 1: Tomar j=1
Paso 2: Hacer D(j)=DE
Paso 3: Tomar j=2
Paso 4: Tomar Sx=0;
Sy=0
Paso 5: Si j≤N seguir los pasos 6 a 11
Paso 6: Tomar i=1
Paso 7: Si i≤n seguir los pasos 8 a 11
Paso 8: Generar U=rand(1)
Paso 9: Hacer α=2*π*U;
x= P *cos α
y = P *sen α
Paso 10: Hacer Sx=Sx+x;
Sy=Sy+y
Paso 11: Tomar i=i+1
Paso 12: Si I>N ,calcular la distancia acumulada mediante la fórmula
D( j) = D( j − 1) + Sx ^2 + Sy ^2
Paso 13: Tomar j=j+1
Paso 14: Repetir los pasos 7a 11
Paso 15: Si j>N entonces hacer DM=D(j)/j;
DE − DM
E= .100
DE
Paso 16: la SALIDA es DM, E
PARAR
85
Cálculo Numérico
6.1.1.3.2. Integración por simulación
Una de las aplicaciones más frecuentes que han tenido los métodos de Montecarlo ha sido
el cálculo de áreas o volúmenes.
Si se desea calcular el área bajo la curva de una función, en el intervalo (a,b), como la que
se muestra en la Figura 6.2., se deben seguir ciertos pasos.
a b
x
Figura 6.2. Gráfico de una función cuya área

se desea calcular
Estos pasos son:
1. Definir un área conocida que comprende el área que se desea integrar, en este caso un
rectángulo.
2. Considerar esta área como un blanco con dos zonas, sobre la que se lanzan dardos, que
se distribuyen sobre toda el área. Cada vez que se lanza el dardo, se observa si quedó en
el área bajo la curva, contabilizando como un acierto en el caso de que así ocurra y no in-
dicando nada en caso contrario.
3. Repetir este proceso N veces
4. Se aproxima el número de aciertos, n, a :
n ≅ pa N (6.8)
donde pa es la probabilidad de acertar.
5. A partir de la consideración de que la distribución de dardos es uniforme, se puede pensar

que n es un número asociado al área buscada y N el valor asociado al área del rectángulo,
por lo cual se tiene:
Area ≅ pa (área rectangular), esta probabilidad es pa=n/N de la ecuación (6.8), por lo tanto
Área ≅ n/N (área rectangular) (6.9)
El procedimiento descripto puede hacerse por generación de números aleatorios.
86
Cálculo Numérico
6.1.1.3.2. Línea de espera
Se puede utilizar el método de Montecarlo para calcular el tiempo medio de espera en una
situación dada.
El problema se estudia para distintas frecuencias de llegada y se simula el transcurso del tiempo
por medio del generador de números aleatorios, el cual proporciona un número cada minuto y de
acuerdo al número se tendrá o no el arribo de una muestra para hacer uso del servicio.
El tiempo de servicio, que por lo general no depende de cuando ocurre el período sino más
bien de la longitud del intervalo, es bastante común en las situaciones cotidianas de líneas de es-
pera. Para calcular el tiempo se aplican los siguientes pasos:
1. Cálculo de la probabilidad de que el sistema esté ocupado Pw
λ
Pw = =ρ (6.10)
µ
dónde λ es la tasa promedio de llegadas y µ es la tasa promedio de servicio.
2. Cálculo de la probabilidad de que el sistema no esté trabajando, o esté vacío:
P0 = 1 – Pw = 1 - ρ (6.11)
3. Cálculo de la probabilidad de que haya n muestras en el sistema Pn:
Pn = P0ρn (6.12)
donde n es cualquier entero no negativo.
4. Cálculo del número promedio de muestras que se encuentran en el sistema, L, ya sea

esperando o siendo analizadas:
ρ
L= (6.13)
1− ρ
5. Cálculo del número promedio de muestras que esperan ser analizadas, Lq.
ρ2
Lq = L − ρ = (6.14)
1− ρ
6. Cálculo del tiempo de espera, en el se distinguen:
Tiempo promedio de una unidad en el sistema,
L
W= (6.15)
λ
y el Tiempo de espera de una muestra para ser analizada:
87
Cálculo Numérico
Lq
Wq = (6.16)
λ
¾ Por ejemplo, si se quiere calcular el tiempo de espera en un laboratorio de análisis de agua al

que llegan 20 muestras por hora y el tiempo promedio de análisis es 30 muestras por hora.
λ 20 2
Se calcula la probabilidad de que el sistema esté ocupado Pw = = = , la probabilidad de
µ 30 3
que no lo esté es 1/3.
El número de muestras que se encuentran en el sistema siendo analizadas o en espera es:

2/3 2 4
L= = 2 , por lo tanto el número que espera ser analizado es L q = 2 − = .
1/ 3 3 3
2 1
El tiempo promedio de una unidad en el sistema es W = = h = 6 min y el tiempo pro-
20 10
4/3 1
medio de espera es Wq = = h = 15 min
20 15
6.2. Modelos Demográficos y de la Cinética Química
Como ejemplos de modelos de la dinámica de sistemas pueden mencionarse los modelos

demográficos y los de la cinética química.
6.2.1. Modelo demográfico
La relación entre la población en un instante de tiempo P(T) y un instante anterior P(T-1),

donde H es el tiempo transcurrido entre ambos, y N y D son los nacimientos y las defunciones
ocurridos en ese tiempo, pueden modelarse por la ecuación:
P(T) = P(T-1) + H*(N - D) (6.17)
Esta ecuación representa un modelo matemático simple de la población y permite determi-

nar la población en el instante T.
Además, por lo general, los nacimientos y las defunciones se expresan en función de las
tasas de natalidad, TN, y de defunción TD, entonces se tienen las siguientes expresiones:
N = P(T-1)*TN (6.18)
D= P(T-1)*TD (6.19)
Se tienen tres ecuaciones, que se pueden interpretar como tres instrucciones de progra-
mación que pueden integrarse en un programa de una computadora. Para generarlo deben reali-
zarse las siguientes tareas:
88
Cálculo Numérico
1. Encabezamiento: se identifica el nombre del programa y su versión, año, autores, finali-

dad y otros datos de interés para el usuario.
2. Inicialización: se definen todos los elementos que participan en el programa, por ejemplo,
variables, parámetros, constantes, etc. Es conveniente especificar el nombre o el rol de
cada variable.
3. Dimensionamiento de las variables: se especifica la dimensiones de las variables que

requieren exigen la reserva de memoria en la computadora y se almacenan en vectores
para su posterior impresión en tablas.
4. Ingreso de datos: el programa debe ser versátil, es decir debe permitir la modificación de
los datos de entrada.
5. Módulo de cálculo: en el se incluye una estructura iterativa mediante la cual se ordena a

la computadora que realice la tarea incluida entre las instrucciones FOR T =… y NEXT T,
hasta alcanzar una determinada condición.
6.2.2. Modelos de la cinética química
La cinética química tiene por objetivo el estudio de las velocidades de las reacciones quí-
micas en sistemas homogéneos y heterogéneos, así como los factores que determinan esas velo-
cidades. El estudio de estos últimos permite formular conclusiones sobre el mecanismo que rige
en esas reacciones.
En el estudio de estos modelos se pueden considerar dos aspectos: el fenomenológico y el

teórico.
Con respecto al primero se estudia la forma en que las velocidades dependen de la con-
centración.
La velocidad de una reacción, puede definirse como la cantidad de sustancia transformada

en la unidad de tiempo. El estudio se realiza sobre una reacción conocida es decir una en la que
ya están establecidas las reacciones estequiométricas entre reactantes y productos,
aA + bB + … mM + nN + … (6.20)
en la que el segundo miembro puede ser reemplazado por P (productos):
aA + bB + … P (6.21)
La teoría demuestra que, en general, la velocidad de reacción en un instante t es propor-
cional a los productos de las concentraciones molares actuales (masas activas) de los reactivos,
elevadas a una potencia igual a su respectivo coeficiente estequiométrico. La cantidad de sustan-
cia transformada puede expresarse en número de unidades de reacción convertidos en productos.
Si esta cantidad es x en el tiempo t, la velocidad es:
89
Cálculo Numérico
dx
v= = k[A ]a [B]b [C]c L (6.22)
dt
esta expresión constituye la ley básica de la cinética química.
La velocidad de una reacción puede expresarse en función de la variación de la concentra-

ción molar de uno de los reactivos (que se consume) o de uno de los productos (que se forma). La
sustancia se selecciona según convenga por su significación, por su facilidad de ser valorada, etc.
La velocidad de consumo de A, puede expresarse como:
d[A ]
= −k[A ]a [B]b [C]c L (6.23)
dt
La constante k es la constante de velocidad o velocidad específica, se determina isotérmi-

camente. El orden total de la reacción, N, se define como el número total de moléculas o átomos
de las sustancias reaccionantes cuyas concentraciones, al variar en función del tiempo, participan
en la expresión de la velocidad de la reacción.
La experiencia prueba que no siempre las reacciones ocurren de acuerdo a lo expresado

en las expresiones anteriores y la ley fundamental de la cinética se expresa en función de los
órdenes de reacción νi con respecto a la concentración que ponderan, es decir:
d[A ]
= −k[A ]ν a [B]νb [C]ν c L (6.24)
dt
El objetivo de los modelos es determinar esos exponentes (coeficientes estequiométricos),

que se llaman x, y , z, ….Pueden emplearse dos métodos: el método diferencial y el integral.
6.2.2.1. Método diferencial
Se determinan los valores de la derivada de la concentración con respecto al tiempo.
Una forma de obtenerlos es a través de las concentraciones con el tiempo y representando

[…] frente a t, ya que la pendiente es la derivada buscada. Estas se comparan con la ecuación
cinética.
La determinación de los órdenes de reacción se efectúa determinando el valor inicial de la

velocidad, para lo que se emplean tiempos de reacción cortos de manera que no haya cambios
apreciables en la concentración y se mide la velocidad. Por ejemplo si se mantienen constantes
todas las concentraciones y al duplicar la de A se observa que la velocidad [A]’ se cuadruplica,
entonces x = 2 y así se determinan todos los exponentes. A continuación se determina k emple-
ando valores medidos, en un experimento de [A], [B], [C], ….
90
Cálculo Numérico
6.2.2.2. Método integral
La ecuación diferencial se transforma por integración en otra que exprese la concentración

en función del tiempo. Para ello se mide la velocidad de reacción durante un largo período de
tiempo y luego se ajustan los datos a la ecuación integrada.
Por ejemplo, se supone que la cinética es de primer orden,
d [ A]
= −[ A] (6.25)
dt
la ecuación integrada es
[A] = [A]0 e –kt (6.26)
por lo tanto una representación de log

[A] frente al tiempo t es una línea recta de pendiente –k,
[A]0
por lo que a partir de las medidas experimentales de [A]con el tiempo se puede obtener el valor de
la constante de velocidad k a partir de la pendiente de la gráfica.
El método integral tiene mucho de artesanal, ya que se formula una hipótesis sobre la
ecuación de velocidad y los órdenes de reacción en reactantes y productos y se valida por con-
traste con los resultados experimentales. Si no existe buena concordancia se formula otra hipóte-
sis y se repiten los cálculos hasta llegar a un resultado satisfactorio.
Este método es el más usado, su principal inconveniente es que se supone que los órde-
nes de reacción son números enteros y esto no siempre es así, además se aplica a una sola curva
de cinética en determinadas condiciones.
El método diferencial es teóricamente más correcto, su inconveniente deriva de la determi-

nación de las pendientes.
6.2.2.3. Dinámica de sistemas cinetoquímicos
Considerando una reacción simple:
A+B C+D (6.27)
Si una molécula de A entra en colisión con una de B, pueden reaccionar y formar C y D.

Suponiendo coeficientes estequiométricos iguales a 1, la velocidad de reacción se puede expresar
como:
d[ ]
= −k[A ][B] (6.28)
dt
Si la variación de fuese lineal, para un tiempo T transcurrido entre T – 1 y T, con H = T –

(T-1), se obtiene:
91
Cálculo Numérico
[L](T ) − [L](T −1)

= −K[A ](T −1) [B](T −1) (6.29)
H
Entonces, el comportamiento temporal de la concentración del reactante A se expresa co-

mo:
[A ](T ) = [A ](T −1) − K[A ](T −1) [B](T −1) H (6.30)
Reemplazando A en el primer miembro por B, C y D, se obtiene las restantes expresiones

y se obtiene la velocidad por cualquiera de los métodos estudiados.
6.1. Encuentre 7 números aleatorios generados por el método de cuadrados centrales, entre 0 y
999, considerando como número inicial 612.
6.2. Calcule los cinco primeros números aleatorios aplicando el método de los productos centrales
si U0 =1422 y U1= 1963
6.3. A una línea de espera de un laboratorio de análisis de miel llegan 30 muestras por hora y el
tiempo promedio de análisis es 10 muestras por hora, realice un análisis de esta línea de espera.
6.4. Considere la sección de empaque de una gran fábrica de galletitas a donde llega el producto
listo para ser envasado. Si llegan a una tasa de 90 grupos a ser envasados por hora y existen 10
líneas de envasado si la tasa de servicio es 12 por hora calcule el tiempo de espera promedio y la
utilización de la instalación.
6.5. Considere una reacción irreversible de primer orden y que la velocidad de reacción se ha se-
guido midiendo la absorbancia de la disolución en distintos tiempos, con los siguientes resultados:
Tiempo 0 18 57 130 240 337 398

(min)
Absorbancia 1,39 1,26 1,03 0,706 0,398 0,251 0,18
Si [A]0= 1,39 y [B]0=0, encuentre el valor de k y los valores de [A]t por aplicación del mode-
lo y por integración directa para t = 0 hasta70 minutos con H = 10 min.
92
Cálculo Numérico
Capítulo 7
SERIES DE FOURIER
En numerosos problemas de ingeniería aparecen funciones periódicas discontinuas y fun-
ciones aperiódicas que pueden representarse por funciones senoidales, esto conduce a la Teoría
de las Series e Integrales de Fourier.
La idea básica de las series de Fourier es que toda función periódica de período T puede
ser expresada como una suma trigonométrica de senos y cosenos del mismo período T.
Neugebauer en 1952 descubrió que los Babilonios utilizaron una forma primitiva de las se-
ries de Fourier para predecir eventos celestiales. La historia moderna de las series de Fourier co-
menzó en 1747con D'Alembert y su tratado sobre las oscilaciones de las cuerdas de un violín.
Euler en 1748 propuso que la solución podía ser expresada en una serie senoidal. Las mismas
ideas fueron luego expuestas por D. Bernoulli (1753) y Lagrange (1759). La formula para calcular
los coeficientes apareció por primera vez en un artículo escrito por Euler en 1777. La contribución
de Fourier comenzó en 1807 con sus estudios del problema del flujo del calor.
7.1. CONSIDERACIONES PREVIAS
7.1.1. Funciones periódicas
Una función f(x) es periódica con período T si está definida para cualquier valor de x real y
si existe un número positivo T tal que:
f(x+T)= f(x) (7.1)
En la Figura 7.1 se muestra un ejemplo de función periódica. Además de la definición pue-

de deducirse que si el período es T entonces cualquier múltiplo entero de T también lo es.
f(x)
T
2T
Figura 7.1. Función periódica con período T
93
Cálculo Numérico
El valor más pequeño de T para el cual se cumple la ecuación (7.1) se denomina período
fundamental de f(x).
Puede demostrarse que si f(x) y g(x) son funciones de período T, la función: h(x) = f(x) +
g(x) también tiene período T.
7.1.2. Series trigonométricas
Una serie trigonométrica puede representarse como:
a0+a1cos x+b1 sen x+a2 cos 2x + b2 sen 2x+ ... (7.2)
donde an y bn son coeficientes de la serie y cada término tiene período 2π. Si la serie converge, su
suma es una función de período 2π.
7.1.3. Funciones seccionalmente continuas
Una función de valor real f es seccionalmente continua en un intervalo [a, b] si:
a) f es definida y continua en todos, excepto en un número finito de puntos de [a, b] y
b) los siguientes límites existen en cada punto de discontinuidad x0 en [a, b]:
f ( x 0 + ) = lim f ( x 0 + h) (7.3)
h→0
f ( x 0 − ) = lim f ( x 0 − h) (7.4)
h→0
Las discontinuidades aceptadas se denominan discontinuidades por salto, de magnitud:
Salto = f(x0+) – f(x0-)
Propiedades
a) Si f(x) ∈ C [a, b] entonces existe la integral de la función en este intervalo y es indepen-
diente de los valores que toma la función f en sus puntos de discontinuidad. En particular si f(x) y
g(x) son idénticas en [a, b] excepto en sus puntos de discontinuidad, entonces puede escribirse la
siguiente igualdad:
b b
∫ f ( x ) dx = ∫ g( x ) dx (7.5)
a a
94
Cálculo Numérico
b) Si f(x) y g(x) ∈ C [a, b] también lo es su producto interno, entonces la integral del producto de
dos funciones seccionalmente continuas siempre existe.
Definido el producto interno en C [a, b] puede demostrarse que el conjunto de funciones
seccionalmente continuas en el intervalo mencionado es un espacio euclidiano denotado por C

[a, b] y contiene a C [a, b] como subespacio.
7.1.4. Funciones pares e impares
Las funciones pares e impares se caracterizan por su simetría con respecto al eje de las
ordenadas y al origen de coordenadas, respectivamente.
Analíticamente las funciones pares son aquellas para las cuales se verifica que:
f(-x)= f(x) (7.6)
mientras que las funciones impares satisfacen:
f(-x)= - f(x) (7.7)
En la Figura 7.2 (a) y (b) se esquematizan estas funciones.
y y (b)
(a)
x x
Figura 7.2. Funciones con simetría: par (a) e impar (b)
Son ejemplos de funciones pares x2, cos nx, x2n, x y de funciones impares x, x3, sen nx,
x2n+1.
95
Cálculo Numérico
Estas funciones tienen propiedades que son muy útiles para el desarrollo de funciones en
serie de Fourier, a continuación se enuncian las más interesantes:
9 La suma o diferencia y el producto o cociente de dos funciones pares da como resultado

una función par.
9 La suma o diferencia de dos funciones impares es impar.
9 El producto o cociente entre dos funciones impares es par.
9 La suma o diferencia de una función impar y una función par no es par ni impar.
9 El producto o cociente entre una función impar y una función par es impar.
7.2. DESARROLLO EN SERIE DE FOURIER
Sea f(x) una función de la variable real x, que en el intervalo (θ, θ +2π) satisface las condi-
ciones de Dirichlet, es decir:
1) f(x) tiene un número finito de discontinuidades,
2) f(x) tiene un número finito de máximos y mínimos
3) f(x) es acotada
entonces, f(x) pertenece a la clase de funciones seccionalmente continuas o continuas por partes
o por tramos, f(x) ∈ C (θ, θ +2π) y admite como modelo matemático la Serie trigonométrica de
Fourier:
b ∞ ∞
f ( x ) = 0 + ∑ a n sen nx + ∑ b n cos nx , con n = 1, 2, 3, … (7.8)
2 n =1 n =1
Siendo an y bn los coeficientes de Fourier que se calculan con las fórmulas siguientes:
1 π
an = ∫ f ( x ) sen nx dx (7.9)
π −π
1 π
bn = ∫ f ( x ) cos nx dx (7.10)
π −π
b0 1 π
= ∫ f ( x ) dx (7.11)
2 2π − π
96
Cálculo Numérico
El intervalo de integración podría extenderse de 0 a 2π o desde x0 a (x0+2π) para todo x0,

sin afectar el valor de los coeficientes an y bn puesto que la función f(x) se repite idénticamente al
cabo del periodo 2π independientemente del límite inferior del intervalo. La elección del intervalo
fundamental [-π,π] es a los fines de facilitar las discusiones de simetría que se presentan más ade-
lante.
7.2.1. Cálculo de los coeficientes de Fourier
Para demostrar la validez de las ecuaciones (7.9) y (7.11) se utilizan las siguientes propie-
dades de las funciones senoidales:
π ⎧⎪0 para n ≠ m
a) ∫ sen nx sen mx dx = ⎨
−π ⎪⎩π para n = m
π ⎧⎪0 para n ≠ m
b) ∫ cos nx cos mx dx = ⎨
−π ⎪⎩π para n = m
π
c) ∫ sen nx cos mx dx = 0
−π
Para determinar an, se multiplican ambos miembros de la serie (7.8) por sen mx y se inte-
gra en un período es decir:
π π ⎡b0 ∞ ∞ ⎤
∫ f ( x ) sen mx dx = ∫ sen mx ⎢ + ∑ a n sen nx + ∑ b n cos nx ⎥ dx (7.12)
−π −π ⎣ 2 n =1 n =1 ⎦
El resultado de esta integral es:
π
∫ a n sen nx sen mx dx = a n π (7.13)
−π
Por lo tanto,
π 1 π
∫ f ( x ) sen mx dx = a n π ∴ a n = ∫ f ( x ) sen nx dx (7.14)
−π π −π
Un cálculo análogo multiplicando por cos mx e integrando en un período completo permite

obtener bn.
Significado de b0/2
97
Cálculo Numérico
Su expresión, dada por la fórmula (7.11), se denomina coeficiente constante, es el valor

medio de la función en el intervalo (-π,π). Se puede evaluar si se calcula el área encerrada por f(x)
en un período y dividirlo entre este período.
7.2.2. Expresión de la serie de Fourier para funciones de período arbitrario
En la práctica las funciones periódicas rara vez tiene un período igual a 2π, por lo general
tiene cualquier período T, por lo general expresado en unidades de tiempo. Para solucionar esto
se utiliza un cambio de escala.
Suponiendo una función f(t) de periodo T, puede introducirse una nueva variable x tal que
f(t) como función de x tenga período 2 π. Si se hace que:
T
t= x (7.15)
2π
puede deducirse que x es:
2π
x= t (7.16)
T
T
Por lo tanto si x=±π corresponde a t = ± . Esto significa que f(x) tiene período 2π. De
2
donde si f tiene una serie de Fourier, será de la forma:
T ⎞ b0 ∞ ∞
f ( t ) = f ⎛⎜ x⎟ = + ∑ a n sen nx + ∑ b n cos nx (7.17)
⎝ 2π ⎠ 2 n =1 n =1
con coeficientes de la forma:
1 π T 1 π T b0 1 π T
an = ∫ f ( x ) sen nx dx ; b n = ∫ f ( x ) cos nx dx ; = ∫ f ( x ) dx
π − π 2π π − π 2π 2 2 π − π 2π
Si se diferencia la ecuación (7.16) se tiene que:
2π
dx = dt (7.18)
T
T T
Y el intervalo de integración sobre el eje x corresponde al intervalo − ≤ t ≤ , por lo tanto consi-
2 2
2π
derando que ω = y reemplazando en las expresiones de la serie y de los coeficientes se tiene
T
que:
b ∞ ∞
f ( t ) = 0 + ∑ a n sen nω t + ∑ b n cos n ω t con n = 1, 2, 3, … (7.19)
2 n =1 n =1
98
Cálculo Numérico
T
2 2
an = ∫ f ( t ) sen n ω t dt (7.20)
T T
−
2
T
2 2
bn = ∫ f ( t ) cos n ω t dt (7.21)
T T
−
2
T
1 2
b0 = ∫ f ( t ) dt (7.22)
T T
−
2
7.2.3. Forma exponencial de la serie de Fourier
Si se consideran las relaciones de Euler: e±ix = cos x ± i sen x la serie puede expresarse
como sigue:
∞
f ( t ) = ∑ β n e in ω t (7.23)
n = −∞
T
1 2
− inω t
Con β n = ∫ f (t) e dt , n= 0, ±1, ±2, ±3, … (7.24)
T T
−
2
7.2.4. Consideraciones simplificatorias
Tomando en cuenta la geometría de las funciones periódicas pueden realizarse algunas

simplificaciones que ayudan al cálculo de la expresión en serie de Fourier.
9 Si f(x) es una función par en el intervalo considerado (-π, π), su desarrollo en serie de Fou-
rier es:
b ∞ 1 π 2π
f ( x ) = 0 + ∑ b n cos nx y además b n = ∫ f ( x ) cos nx dx = ∫ f ( x ) cos nx dx
2 n =1 π −π π0
9 Si f(x) es una función impar en el intervalo considerado (-π, π), su desarrollo en serie de
Fourier es:
∞ 1 π 2π
f ( x ) = ∑ a n sen nx y además a n = ∫ f ( x ) sen nx dx = ∫ f ( x ) sen nx dx
n =1 π −π π0
99
Cálculo Numérico
El desarrollo en el período (-π, π) es igual al desarrollo en el período (0, π) multiplicado por

dos.
7.2.5. Espectro de frecuencias
Si la función a desarrollar es función del tiempo y si T es el período, se puede expresar de

acuerdo a la expresión (7.19), ésta puede escribirse en forma más compacta como:
b ∞
f ( t ) = 0 + ∑ C n cos(nωt − ϕn ) (7.25)
2 n =1
con Cn = a n 2 + b n 2 (7.26)
a
y con ϕn = tg −1 n (7.27)
bn
La función f(t) queda desarrollada en una suma de componentes armónicas en las que Cn
y ϕn son la amplitud y el ángulo de fase respectivamente, como se observa en la Figura 7.3.
Cn
ϕn
bn
an
nωt
Figura 7.3. Diagrama fasorial de la Serie de Fourier
La ecuación (7.25) muestra claramente que una función periódica que satisface las condi-
ciones de Dirichlet puede descomponerse en un valor medio y componentes senoidales armóni-
camente relacionadas cuyas frecuencias son múltiplos enteros de la frecuencia fundamental 1/T.
La amplitud y fase de las componentes están dadas por las expresiones (7.26) y (7.27), respecti-
vamente.
Si se representa gráficamente el valor absoluto de Cn en función de nω, se obtiene el es-

pectro discreto de frecuencias, si se desea representar también ϕn, se deben ampliar las deno-
minaciones a espectro discreto de amplitud (Figura 7.4.(a)) y el espectro discreto de fase (Fi-
gura 7.4 (b)), respectivamente, en el que cada componente de fase o de amplitud está espaciada
en ω y consiste en líneas discretas.
100
Cálculo Numérico
Cn (a) ϕn (b)
ω 2ω 3ω 4ω nω ω 2ω 3ω 4ω
Figura 7.4. Espectro discreto de amplitud (a) y de fases (b)
7.3. INTEGRALES DE FOURIER
Las integrales de Fourier pueden interpretarse como el resultado de una generalización del
método de las Series de Fourier para incluir funciones no periódicas. La forma más frecuentemen-
te adoptada es del par unilateral de integrales de Fourier, que se representa como sigue:
1 ∞ iωt −1
f (t) = ∫ g(ω) e dω = F (g(ω)) (7.28)
2π − ∞
∞
g(ω) = ∫ f ( t ) e −iωt dt = F( f ( t )) (7.29)
0
Siendo la expresión (7.28) la Integral inversa de Fourier y (7.29) la integral de Fourier.
¾ Por ejemplo si se desea encontrar el desarrollo en serie de Fourier de la función:

-k para -π < x < 0
f(x)=
k para 0 < x < π
Considerando que f(x + 2 π) = f(x)
En primer lugar se realiza un gráfico de la función:
f(x)
-π
π x
-k
101
Cálculo Numérico
Por inspección gráfica se puede afirmar que es una función impar puesto que f(x)=-f(-x), y
tomando en cuenta las consideraciones simplificatorias el desarrollo en Serie de Fourier está dado
∞
por: f ( x ) = ∑ a n sen nx y el coeficiente an se calcula como:
n =1
π
1 π 2π 2π 2⎡ k ⎤ 2k
an = ∫ f ( x ) sen nx dx = ∫ f ( x ) sen nx dx = ∫ k sen nx dx = ⎢− cos nx ⎥ = (− cos nπ + cos 0 ) =
π −π π0 π0 π⎣ n ⎦ 0 nπ
2k
= (1 − cos nπ)
nπ
Puede deducirse que si n es impar (1-cos nπ) es igual a 2 y si n es par (1-cos nπ) es igual a cero,
por lo tanto la serie toma valores solamente para n impar, esto puede expresarse tomando (2n-1)
como coeficiente del argumento. En síntesis, la serie en forma compacta es:
∞ 4k
f ( x) = ∑ sen(2n − 1)x
n =1 ( 2n − 1)π
Una representación de la función y de la 1º, 3º, 5 y 7º armónica se muestra en la figura que

sigue.
16/ π (sin(x)+1/3 sin(3 x)+1/5 sin(5 x)+1/7 sin(7 x))

6
-2
-4 f(x)
1º arm
3º arm
5º arm
7º arm
-6
-3 -2 -1 0 1 2 3
x
7.1. Demuestre que las formas trigonométricas dadas a continuación son expresiones de la serie
de Fourier:
∞ ∞
a) f (t ) = c 0 + ∑
n =1
c n cos( nω t − ϕ n ) b) f (t ) = c 0 + ∑c
n =1
n sen( nω t + φ n )
7.2. A partir de la forma básica de la Serie de Fourier deduzca su expresión exponencial.
102
Cálculo Numérico
7.3. Encuentre el desarrollo en serie de Fourier de la función definida por:

0, en –T/6 < t < -T/12
f(t)= A, en –T/12 < t < T/12
0, en T/12 < t < T/6
7.4. Desarrolle en Serie de Fourier la función:

-8 para -π < x < 0
f(x)=
8 para 0 < x < π
Considere que f(x + 2 π) = f(x)
7.5. Desarrolle f(x) = 3x2, en serie de Fourier si el período es 2 π.
7.6. Sea f(x) = x periódica en el intervalo (-4, 4), encuentre los coeficientes y escriba su desarrollo
en serie de Fourier.
103
Cálculo Numérico
Capítulo 8
TRANSFORMADA DE LAPLACE
El método de la Transformada de Laplace constituye una valiosa herramienta para la reso-

lución de numerosos problemas que se presentan en ingeniería, tales como la resolución de
ecuaciones diferenciales, evaluación de integrales, solución de ecuaciones integrales, etc.
En el caso particular de la resolución de ecuaciones diferenciales, este método tiene nota-

bles ventajas sobre los otros métodos de resolución ya que se pueden transformar las ecuaciones
diferenciales en algebraicas, cualesquiera sean las condiciones iniciales dadas, estas se incorpo-
ran al problema algebraico y el uso de tablas de transformadas de Laplace facilita notablemente la
resolución.
8.1. DEFINICIÓN DE TRASFORMADA DE LAPLACE
Sea f(t) una función de t definida para t ≥ 0, la Transformada de Laplace de f(t), que se
denota por L [f(t)] ó por F(s), está definida por la ecuación:
∞
L [f(t)] = F(s) = ∫ e − st f ( t )dt (8.1)
0
donde s es un parámetro complejo, pero se restringe a valores reales. El símbolo L se denomina

operador de la Transformada de Laplace. La integral impropia de la ecuación (8.1) puede definirse
como:
m
− st
lím ∫ e f ( t )dt (8.2)
m→∞ 0
Si este límite existe, la integral converge y puede asegurarse la existencia de la Transfor-

mada de Laplace. A continuación se indica el teorema para la existencia de la transformada.
Teorema para la existencia de la transformada de Laplace: Si f(t) es seccionalmente continua en el intervalo
0 ≤ t ≤ A, para cualquier A positiva y además f (t ) ≤ Ke at donde t ≥M y K, a y M son constantes reales y K

y M necesariamente positivas, entonces existe la Transformada de Laplace de f(t) para s > a.
Algunas características de la función f(t) son las siguientes:
104
Cálculo Numérico
9 la variable real no necesariamente representa al tiempo, puede utilizarse el símbolo ade-

cuado para el sistema en consideración,
9 f(t) debe ser una función unívoca para que tenga una única F(s). Sin embargo, se admite
que la función tenga un número finito de discontinuidades,
9 si existe discontinuidad para f(t) en t=0, entonces el límite inferior de la transformada se
considera como una aproximación por el lado positivo, es decir f(0+).
8.2. PROPIEDADES IMPORTANTES DE LA TRANSFORMADA DE LAPLACE
Se detallan, en los apartados que siguen, propiedades de la transformada de Laplace que

facilitan notablemente su cálculo.
a) Linealidad
Teorema: Si c1 y c2 son constantes y f1(t) y f2(t) son funciones cuyas transformadas de Laplace son, respec-
tivamente, F1(s) y F2(s), entonces:
L[c1 f1(t) + c2 f2(t)] = c1 L[f1(t)] + c2 L[f2(t)] = c1 F1(s) + c2 F2(s) (8.3)
El símbolo L es el operador transformada de Laplace y debido a la propiedad de linealidad
expresada en este teorema se dice que L es un operador lineal.
b) Traslación o desplazamiento de frecuencia
Teorema: La transformada de Laplace de e −α t veces una función es igual a la transformada de Laplace de

la función con s reemplazado por (s + α).
−α t
L[ e f(t)] = F(s +α) (8.4)
c) Desplazamiento temporal
Teorema: Si L [f(t)] = F(s), entonces: L [f(t – t0) u (t – t0)] = e − t0s F(s) (8.5)
d) Cambio de escala
1 ⎛s⎞
Teorema: Si L [f(t)] = F(s), entonces: L[f(at)] = F⎜ ⎟ (8.6)
a ⎝a⎠
8.2.1. Transformada de Laplace de operaciones
Al aplicar la transformada de Laplace a la solución de ecuaciones diferenciales se deben

transformar no solamente funciones sino también operaciones, por ello se enuncian a continua-
ción los teoremas pertinentes
105
Cálculo Numérico
Teorema de la derivación: Si una función f(t) y su primera derivada son ambas Laplace transformables y si
d f (t )
L [f(t)] = F(s), entonces: L = L[D f(t)] = s F(s) – f(0+) (8.7)
dt
Esta propiedad se extiende para las derivadas de orden superior.
⎡t ⎤ F (s )
Teorema de la integración: Si L [f(t)] = F(s), entonces: L ⎢∫ f (t )dt ⎥ = (8.8)
⎣⎢ 0 ⎦⎥ s
Teorema de la multiplicación por t: Si L [f(t)] = F(s), entonces: L [t f(t)] = - d F (s ) (8.9)

ds
∞
⎡1 ⎤
Teorema de la división por t: Si L
⎣t ⎦ s
∫
[f(t)] = F(s), entonces: L ⎢ f (t )⎥ = F (s ) ds (8.10)
8.3. MÉTODOS PARA CALCULAR TRANSFORMADAS DE LAPLACE
Existen varios métodos para encontrar las transformadas de Laplace:
9 Método directo: utilizando la definición dada en la expresión (8.1).
9 Método de las series: si f(t) tiene un desarrollo en serie dado por:
∞
9 F(t) = a0 + a1t + a2t2 + ... = ∑a t
n =0
n
n
, su transformada de Laplace puede calcularse como la
suma de la transformada de cada uno de los sumandos de la serie, es decir:
∞
a 0 a1 2! a 2 n! a n
9 L [f(t)] =
s s
+ 2 + 3 + L=
s n =0 s
∑
n +1
9 Método de las ecuaciones diferenciales: consiste en encontrar una ecuación diferencial

que sea satisfecha por f(t) y aplicar luego los teoremas anteriores.
9 Mediante el uso de tablas.
8.4. VENTAJAS DEL MÉTODO DE LA TRANSFORMADA DE LAPLACE
9 Simplifica funciones: pues la transformada convierte funciones que ocurren frecuentemen-

te, tales como funciones exponenciales o funciones trascendentes y sus combinaciones en
simples funciones algebraicas, fáciles de operar.
106
Cálculo Numérico
9 Simplifica operaciones: la transformada de Laplace convierte la derivación en multiplica-

ción y la integración en división.
9 Determina constantes automáticamente: al resolver la ecuación diferencial de un sistema

dado por el método de la transformada de Laplace, la solución surge completa y no es ne-
cesario operar adicionalmente para determinar las constantes.
9 Permite generalizar respuestas de sistemas: si se determina la respuesta de un sistema a

las funciones escalón o impulso como excitación, es posible determinar fácilmente la res-
puesta del mismo sistema a cualquier otro tipo de función excitación.
8.5. TRANSFORMADA DE LAPLACE DE ALGUNAS FUNCIONES IMPORTANTES
En la Tabla 8.1 se muestran las transformadas de Laplace para algunas funciones impor-
tantes.
Tabla 8.1. Transformadas de Laplace de algunas funciones utilizadas con frecuencia
Nombre de la Función Expresión Transformada de Laplace

1
Escalón unitaria u(t)
s
1
Exponencial de t e at
s −a
ω
Seno sen ωt
s +ω2
2
s
Coseno cos ωt
s2 +ω 2
b
Seno hiperbólico senh bt
s2 − b2
s
Coseno hiperbólico cosh bt
s2 − b2
n!
Potencias positiva de t tn
sn +1
¾ Por ejemplo, si se desea calcular la transformada de Laplace de la siguiente función:
f(t)= 8 e 5 t + 12 t 3 − 6 sen 4 t + 4 cos 2 t
Utilizando la propiedad de linealidad se puede escribir la transformada como sigue:
L [ 4 e 5 t + 6 t 3 − 3 sen 4 t + 2 cos 2 t ] = L [ 4 e 5 t ] + L [ 6 t 3 ] - L [ 3 sen 4 t ]+L [2 cos 2 t ]
Utilizando la Tabla 8.1, el cálculo se simplifica notablemente y se llega que:
107
Cálculo Numérico
8 72 24 4s
L[f(t)]= + − +
s−5 s 4 2 2
s + 16 s + 4
8.6. TRANSFORMADA INVERSA DE LAPLACE
Si se necesita encontrar una función f(t) cuya trasformada de Laplace se conoce, es decir
F(s), lo que se está buscando es en realidad una función que se denomina transformada inversa
de Laplace que se simboliza como L-1 [f(t)], donde L-1 se denomina operador de la transformada
inversa de Laplace.
Para evaluar estas transformadas inversas se pueden utilizar varios métodos, tales como
el uso de tablas, de teoremas sobre transformadas inversas de Laplace, el método de las convo-
luciones.
¾ Por ejemplo si se desea determinar la función que tiene como transformada de Laplace a:
8 9
F(s) = − , utilizando la columna de la derecha de la Tabla 8.1, se puede deducir que:
(s − 1) (s − 3 )
L-1 F(s)=8 et - 9 e3t
8.7. INTEGRAL DE CONVOLUCION
Si F(s) y G(s) son las transformadas de Laplace de las funciones f(t) y g(t), respectivamen-
te, entonces:
L-1 [F(s)±G(s)]= f(t) ± g(t) (8.11)
Es interesante saber si existe alguna expresión para la transformada inversa de Laplace

del producto F(s) G(s) en términos de las funciones originales, para ello se utiliza el siguiente teo-
rema:
Teorema de convolución:
t
Si L-1[F(s)]= f(t) y L-1[G(s)]= g(t), entonces L-1 [F(s).G(s)]= ∫ f (u) g( t − u) du = f * g .
0
Donde f *g se denomina la convolución de f y g.
En forma equivalente se tiene que:
108
Cálculo Numérico
⎧t ⎫
L-1[f*g]= L ⎨ ∫ f (u)g( t − u)du⎬ = F(s)G(s)
⎩0 ⎭
s
¾ Por ejemplo si se desea encontrar la transformada inversa de Laplace de: , se
(s 2 + 1) 2
⎧⎪ s ⎫⎪ ⎧⎪ 1 ⎫⎪
tiene que: L-1 ⎨ ⎬ = cos t y que L-1 ⎨ ⎬ = sent . Por el teorema de convolución se llega
⎪⎩ (s 2 + 1) ⎪⎭ ⎪⎩ (s 2 + 1) ⎪⎭
a que:
L-1
⎧⎪ s 1 ⎫⎪ t
L-1 ⎨ ⋅ ⎬ = cos t * sent = ∫ cos u sen( t − u) du =
⎪⎩ (s 2 + 1) (s 2 + 1) ⎪⎭ 0
t
= ∫ cos u [sent cos u − cos t senu]du =
0
t t 1 1 1 1
= sent ∫ cos 2 u du − cos t ∫ senu cos u du = sent ⎡⎢ t + sent cos t ⎤⎥ − cos t ⎡⎢ sen 2 t ⎤⎥ = t sent
0 0 ⎣2 2 ⎦ ⎣2 ⎦ 2
8. 1. Demuestre, suponiendo que s > 0 y s > a, que:
1 a s
a) L [ eat ]= ; b)L [sen at] = ; c)L [cos at] =
s −a 2
s +a 2
s + a2
2
8.2. Encuentre la transformada de Laplace de las siguientes funciones:
senω t
a) t2 e3t; b) e-2t sen 4t; c)
t
8.3. Encuentre la transformada inversa de Laplace de:
1 1 4s + 12
a) ; ;
s − a s 2 + 9 s 2 + 8s + 16
8.4. Encuentre la Transformada Inversa de Laplace de las ecuaciones diferenciales:

a) y’’+ y = 16 cos t; y(0)=0, y’(0)=0
b) y’’ + 2 y’+ 5 y = 0; y(0)=3, y’(0)= -7
109
Cálculo Numérico
8.5. Considere dos tanques de salmuera conectados como se indica en la Figura 2. El tanque 1
contiene x kg de sal en 100 l de salmuera y el tanque 2 y kg da sal en 200 l de salmuera. Esta se
conserva uniforme por agitación y se bombea de un tanque a otro a la tasa indicada. Además al
tanque 1 fluye agua pura a una tasa de 7,08 l/s y del tanque 2 salen 7,08 l/s por lo que el volumen
de salmuera en ambos tanques permanece constante. Encuentre la concentración de sal en cada
tanque en función del tiempo, considerando que la concentración inicial es 0,4 kg/l en el tanque 1
y 0,3 kg/l en el tanque 2.
110
Cálculo Numérico
Capítulo 9
ECUACIONES DIFERENCIALES
La importancia de las ecuaciones diferenciales en matemáticas aplicadas se debe al hecho

que la mayor parte de las leyes científicas se expresan en términos de la rapidez de variación de
una variable con respecto otra. Además proporcionan una herramienta esencial para modelar mu-
chos problemas en Ingeniería, Física, Economía y Biología, puesto que estos, por lo general, re-
quieren la determinación de una función que satisface a una ecuación diferencial.
9.1. CONCEPTOS PREVIOS
Una ecuación diferencial es una ecuación que involucra derivadas de una función desco-
nocida de una o más variables.
Una ecuación diferencial ordinaria es aquella en la que la función desconocida depende

de una variable (las derivadas son derivadas ordinarias). Un ejemplo se muestra en la ecuación
(9.1).
dy
= 4 x − y o y' = 4 x − 4 o bien D y = 4 x − 4 (9.1)
dx
Una ecuación diferencial parcial es aquella en la que la función desconocida depende de

más de una variable (las derivadas son derivadas parciales), tal como se muestra en la ecuación
(9.2).
∂2v ∂2v
+ =v (9.2)
∂ t2 ∂ x2
El orden de una ecuación diferencial es el orden de la derivada más alta que aparece en la
ecuación.
El grado de una ecuación diferencial es el grado algebraico de la derivada de mayor orden

que aparece en la ecuación.
111
Cálculo Numérico
Ecuación diferencial lineal: es la ecuación en la que no aparecen potencias de la variable

dependiente y ni de sus derivadas, ni productos de la variable dependiente por sus derivadas o
productos entre derivadas, puede escribirse como:
a 0 ( x ) y (n) + a1( x ) y (n −1) + ... + a n −1( x ) y'+a n ( x ) y = f ( x ) (9.3)
donde f(x) y los coeficientes a0(x), a1(x),…, an(x) son funciones de x y a0(x) no es idéntica a 0.
Solución de una ecuación diferencial: es cualquier relación funcional que no incluya de-
rivadas o integrales de funciones desconocidas y que implique a la propia ecuación diferencial, en
el sentido que la verifique idénticamente por sustitución directa.
Ecuación o condiciones homogéneas: se dan si son satisfechas por una función parti-
cular y(x) y también por cy(x), donde c es una constante arbitraria. Por ejemplo, una condición
homogénea puede obtenerse del requerimiento de que una función o una de sus derivadas (o
alguna combinación lineal de la función y/o ciertas de sus derivadas) sea nula.
Problemas de valor inicial y de frontera
Dada una ecuación diferencial ordinaria de orden n y cualquier grado, se establece que en
su solución general deben aparecer n constantes arbitrarias, entonces la solución general es
una función g de la forma: g(x, y, c1, c2, …, cn) =0 y los valores de estas constantes se determinan
sujetando la solución general a n condiciones independientes.
Dependiendo de cómo se establezcan estas condiciones, se distinguen dos tipos de pro-

blemas: los de valores iniciales y los de valores de frontera.
9 Un problema de valor inicial es un problema que busca determinar una solución a una
ecuación diferencial sujeta a condiciones sobre la función desconocida y sus derivadas
especificadas en un valor de la variable independiente. Estas condiciones se denominan
condiciones iniciales.
9 Un problema de valor de frontera es un problema que busca determinar una solución a

una ecuación diferencial sujeta a condiciones sobre la función desconocida especificadas
en dos ó más valores de la variable independiente. Estas condiciones se denominan
condiciones de frontera.
9.2. ECUACIONES DIFERENCIALES ORDINARIAS LINEALES DE PRIMER ORDEN
Una ecuación diferencial ordinaria puede simbolizarse como
dy
= f ( x, y ) (9.4)
dx
112
Cálculo Numérico
Dentro del tema de las ecuaciones diferenciales ordinarias, se destaca, dentro de las ma-
temáticas aplicadas, el problema del valor inicial.
El conjunto constituido por una ecuación diferencial ordinaria lineal y la condición inicial
constituyen el problema del valor inicial, que puede expresarse como sigue:
⎧⎪ d y ⎫
a1 + a o ( x ) y( x ) = f ( x )⎪
⎨ dx ⎬ (9.5)
⎪⎩y(0) = y 0 ⎪⎭
donde a1 y a2 son coeficientes, y es la variable dependiente, x la variable independiente, f(x) es

una función. Desde el punto de vista sistémico, se puede decir que f(x) es la excitación o entrada
del sistema y y es la respuesta o salida del sistema.
9.2.1. Solución analítica
Si se normaliza la ecuación, es decir se dividen ambos miembros por el coeficiente de la

derivada de mayor orden, se tiene:
⎧⎪ d y ⎫
+ b o ( x ) y( x ) = h( x )⎪
⎨d x ⎬ (9.6)
⎪⎩y(0) = y 0 ⎪⎭
a f (x)
en la que b 0 = 0 y h( x ) =
a1 a1
Recordando lo estudiado en análisis matemático puede afirmarse que la solución general

del problema de valor inicial es:
y( x ) = K e − ∫ b0 ( x ) dx + e − ∫ b0 ( x ) dx ∫ h( x ) e ∫ b0 ( x ) dx dx (9.7)
La solución general está constituida por el primer sumando que corresponde a la solución
transitoria o función complementaria yc(x), que es independiente de la función excitación y el
segundo sumando que corresponde a la solución permanente o forzada yp(x) que depende de
la función excitación.
Por lo general se trabajan con coeficientes constantes, por lo tanto la solución general de
una ecuación diferencial lineal de 1º orden puede expresarse como sigue:
b0 x
y( x ) = K e − b 0 x + e − ∫ h( x ) e
b0 x
dx (9.8)
La constante K se evalúa con la condición inicial y(0) = y0.
9.2.1.1. Solución por el método del operador diferencial
Si se considera ahora una ecuación diferencial en la que la variable independiente es el

tiempo, se tiene:
113
Cálculo Numérico
du
+ α u( t ) = h( t ) (9.9)
dt
u( t 0 ) = u 0
La solución general aplicando la fórmula (9.7) es:
u( t ) = K e −α t + e −α t ∫ h( t ) e α t dt (9.10)
Como se observa, esta solución es la suma de la suma de la solución complementaria y de la par-

ticular:
u(t) = uc(t) + up(t) (9.11)
Si la función excitación h(t) es nula, el problema descripto por (9.9) se puede escribir co-
mo:
du c
+ α uc (t) = 0
dt (9.12)
uc (t 0 ) = u0
y su solución es:
-α t
uc(t) = K e
En forma abreviada puede expresarse que:

-αt
D uc(t) + α uc(t) = 0 ⇒ uc(t) = K e
La inferencia que puede extraerse de la expresión anterior es:

-α t
I1: (D + α) uc(t)= K e (9.13)
Si existe función excitación la ecuación (9.9) puede escribirse en forma abreviada como:
D u(t) +α u(t)= h(t) (9.14)
y su solución particular es:
up ( t ) = e −α t ∫ h( t ) e α t d t (9.15)
por lo tanto puede inferirse que:
1
I2: up ( t ) = h( t ) ⇒ up ( t ) = e −α t ∫ h( t ) e α t d t (9.16)
D+α
1
Donde (D+α) es el operador derivada polinómico y es el operador derivada polinómico in-
D+α
verso.
114
Cálculo Numérico
9.2.2. Solución por métodos numéricos
Dentro de los métodos numéricos de solución aproximada a ecuaciones diferenciales ordi-

narias, se distinguen los métodos de un paso y los de pasos múltiples.
9.2.2.1. Métodos de un paso
Con los métodos de Paso Único, la solución numérica a una ecuación diferencial, se en-
cuentra utilizando la información de la derivada en el punto xi para extrapolar a partir de un valor
anterior yi a un valor actual yi+1 en una distancia determinada, denotada con h. Este procedimiento
se aplica paso a paso para encontrar el valor de y y trazar la trayectoria de la solución. En cada
paso de extrapolación se comete error. Puede simbolizarse como:
y i +1 = y i + φ h (9.17)
Donde yi+1 es el nuevo valor, yi es el valor anterior, φ es la función de incremento o pendiente y h

el paso.
9.2.2.1.1. Método de la Serie de Taylor
Este método no es estrictamente un método numérico pero por lo general se utiliza junto
con los esquemas numéricos, es de aplicabilidad general y sirve como punto de partida para otros
métodos.
Dada una ecuación diferencial como la explicitada en la ecuación (9.4), la relación entre x y
y se obtiene al encontrar los coeficientes de la Serie de Taylor. Desarrollando la serie para un va-
lor xi se tiene:
f ' ' (xi ) f (n ) ( x i )

f ( x i+1 ) = f ( x i ) + f ' ( x i )( x i+1 − x i ) + ( x i+1 − x i ) 2 + L + ( x i+1 − x i )n (9.18)
2! n!
Si se trunca la serie en el tercer término, el método se denomina de 2º orden y reempla-

zando f(x) por y(x) y (xi+1-xi) por h, se obtiene:
y ' ' (xi ) 2

y( x i + 1 ) = y i + y ' ( x i ) h + h (9.19)
2!
El valor de y’(x) se conoce a partir de la ecuación diferencial y y’’(x) se obtiene por diferen-
ciación de la ecuación original. Incrementando el orden del método se obtiene mayor precisión
pero las expresiones se tornan complicadas, por lo general se utiliza muy poco en la práctica.
¾ Por ejemplo si se desea encontrar la solución del siguiente problema de valor inicial para x=1:
⎧⎪ dy
=x+y
⎨ dx
⎪⎩y(0) = 1
115
Cálculo Numérico
Para aplicar la fórmula (9.19) las derivadas son:
y’(xi)=xi+y(xi) y derivando esta expresión con respecto a x, se tiene
y’’(xi)= 1+ y’(xi)=1+ xi+y(xi)
Por lo tanto la expresión de recurrencia a aplicar reemplazando estos valores en (9.18) es:
h2
y( x i + 1 ) = y i + ( x i + y i ) h + (1 + x i + y i )
2!
Se comienza con x=0 para el que y(0)=1 y suponiendo un paso h=0,1 en la tabla que sigue se
explicitan los valores obtenidos. Se incluyen también los valores que se obtiene utilizando la solu-
ción analítica que es: y(x)=2ex-x-1.
Método de la Método Error

Serie de analítico verdadero
Taylor
xi yi yi
0 1,000 1 0
0,1 1,110 1,110 0
0,2 1,240 1,243 0,003
0,3 1,390 1,399 0,09
0,4 1,560 1,584 0,024
0,5 1,750 1,797 0,047
0,6 1,960 2,044 0,084
0,7 2,190 2,327 0,137
0,8 2,440 2,651 0,211
0,9 2,710 3,019 0,309
1 3,000 3,436 0,436
9.2.2.1.2. Método de Euler
Si se trunca la serie de Taylor a partir de su segundo término, se tiene la siguiente expre-

sión:
y( x i + 1 ) = y i + y ' ( x i ) h (9.20)
Por lo general y por razones de simplicidad se denota a y’(xi) como f(xi,yi). A la expresión (9.20) se
la conoce como método de Euler o método de Euler-Cauchy o de pendiente puntual. Se caracteri-
za porque se predice un nuevo valor de y usando la pendiente, es decir la primera derivada en el
valor anterior de x para extrapolar linealmente sobre un tamaño de paso constante h. En la figura
9.1 se muestra un esquema del método.
116
Cálculo Numérico
y Valor predicho
Error
Valor real
xi xi+1
x
Figura 9.1. Esquema del método de Euler

Puede decirse que el método puede sintetizarse como:
Valor nuevo = valor anterior + pendiente * paso
La solución de una ecuación diferencial ordinaria por este método incluye dos tipos de
error: error de truncamiento y de redondeo. El error de truncamiento se compone de dos partes, el
error de truncamiento local que resulta de aplicar el método en un paso y el error de propagación
que resulta de las aproximaciones producidas en los pasos anteriores, la suma de ambos da el
error de truncamiento global. A partir de la serie de Taylor puede deducirse que el error de trun-
camiento local es:
f ' (xi, yi ) 2
Et = h + L + Οh(n + 1) (9.20)
2
Para un paso lo suficientemente pequeño los errores decrecen a medida que el orden crece y por
lo general se calcula el error de truncamiento local aproximado:
f ' (xi, yi ) 2
Eta = h = Ο(h 2 ) (9.21)
2
El algoritmo correspondiente es:
Algoritmo para aproximar la solución de un problema de valor inicial por el

método de Euler
y' = f ( x, y ), x 0 ≤ x ≤ x n , y( x 0 ) = y 0 : problema de valor inicial
f(x,y):función que especifica la ecuación diferencial
X: valor de la variable independiente
x0: valor inicial de la variable independiente
xn: valor final de la variable independiente
Y: valor de la variable dependiente
y0: valor inicial de la variable dependiente
h: tamaño de paso
continua
117
Cálculo Numérico
N: número entero que indica el número de segmentos igualmente espaciados (es

igual al número de datos menos 1)
Paso 1: Ingresar x0,xn, N, y0
(x 0 − xn )
Paso 2: Hacer h =
N
Paso 3: Hacer X= x0, Y= y0
Paso 3: Para i=1, 2, ..., N seguir los pasos 4 y 5
Paso 4: Tomar Y=Y+h*f(x,y)
X=X+h
Paso 5: la SALIDA es (X,Y)
PARAR
¾ Por ejemplo si se desea encontrar la solución del problema de valor inicial planteado en el ítem
anterior la expresión que surge de aplicar la ecuación (9.19) es:
y( x i + 1 ) = x i + ( x i + y i ) h
Si se desea calcular y(x0)= y0 + (x0+ y0)h=1+(0+1)0,1=1,1 y así sucesivamente.
Método de Método Error

Euler analítico verdadero
xi yi yi
0 1,000 1 0
0,1 1,100 1,110 0,010
0,2 1,220 1,243 0,023
0,3 1,362 1,399 0,037
0,4 1,528 1,584 0,056
0,5 1,721 1,797 0,076
0,6 1,943 2,044 0,101
0,7 2,197 2,327 0,130
0,8 2,487 2,651 0,164
0,9 2,816 3,019 0,203
1 3,187 3,436 0,249
9.2.2.1.3 Métodos de Runge-Kutta
Se trata de una familia de métodos y si bien se estima el valor de la función utilizando la

expresión general dada por la ecuación (9.17), en lugar de calcular derivadas de orden superior,
se evalúa la función en un mayor número de puntos, tratando de igualar la precisión del método
de la serie de Taylor.
La fórmula general está dada por la expresión:
y i + 1 = y i + φ( x i , y i , h) h (9.22)
Donde φ(xi,yi,h) es la función de incremento y puede interpretarse como el promedio de la pen-

diente sobre el intervalo.
118
Cálculo Numérico
Esta función de incremento puede escribirse como:
φ = α1b1+ α2b2 + ... + αnbn (9.23)
los valores de α son constantes y los de b son en realidad ecuaciones recurrentes de la forma:
b1=f(xi,yi) (9.23a)
b2=f(xi + l1h, yi + m11b1h) (9.23b)
b3=f(xi + l2h, yi + m21b1h + m22b2h ) (9.23c)
bn=f(xi + ln-1h, yi + mn-1,1b1h + mn-1,2b2h + ...+ mn-1,n-1bn-1h) (9.23c)
Una vez que se establece n los valores de α, l y m, se evalúan igualando la ecuación

(9.22) a la serie de Taylor (9.18). El método de Runge-Kutta de 1º orden (n=1) es el método de
Euler, los métodos de Runge-Kutta de 2º orden utilizan una función de incremento con dos térmi-
nos y debido a que se desprecian los términos con h3 y de orden superior durante la derivación, el
error local de truncamiento es Ο(h3) y el error global es Ο(h2). Si n=3 el método es de 3º orden, si
n=4 es de 4º orden.
9 Método de Runge-Kutta de 2º orden
En este caso las ecuaciones (9.22), (9.23a y 9.23b) pueden expresarse como:
y i + 1 = y i + (α1b1 + α 2b 2 ) h (9.24)
b1=f(xi,yi) (9.24a)
b2=f(xi + l1h, yi + m11b1h) (9.24b)
Para obtener el valor de los coeficientes se debe igualar la ecuación (9.24) a la serie de
Taylor de 2º orden que es:
h2
y( x i + 1 ) = y i + f ( x i , y i ) h + f ' ( x i , y i ) (9.25)
2!
donde f’(xi,yi) debe determinarse utilizando la regla de la cadena, es decir:
∂f ∂f dy
f ' (xi, yi ) = + (9.26)
∂x ∂y dx
y sustituyendo (9.26) en (925) se llega a que:
⎛ ∂f ∂f dy ⎞ h 2
y( x i + 1 ) = y i + f ( x i , y i ) h + ⎜ + ⎟ (9.27)
⎝ ∂x ∂y dx ⎠ 2 !
119
Cálculo Numérico
Para tratar de igualar la expresión (9.24b) a la serie se utiliza la expresión (9.27) pero apli-
cada a una función de dos variables, es decir:
∂f ∂f
f ( x i + l 1h, y i + m11b1h) = f ( x i , y i ) + l 1h + m11b1h + Ο(h 2 ) (9.28)
∂x ∂y
Por lo tanto si sustituimos este resultado en (9.24) se obtiene:
∂f ∂f
y i + 1 = y i + α1h f ( x i , y i ) + α 2 h f ( x i , y i ) + α 2 l 1h 2 + α 2 m11h 2 f ( x i , y i ) + Ο(h 3 ) (9.29)
∂x ∂y
Y reordenando términos en la ecuación (9.29) se llega a que:
⎡ ∂f ∂f ⎤
y i + 1 = y i + [α 1 f ( x i , y i ) + α 2 f ( x i , y i )]h + ⎢α 2 l 1 + α 2m11f ( x i , y i ) ⎥h 2 + Ο(h 3 ) (9.30)
⎣ ∂x ∂y ⎦
Comparando (9.30) con (9.27) se determina que para que las dos expresiones sean equivalentes,
se debe cumplir que:
α1 + α 2 = 1
1
α 2 l 1= (9.31)
2
1
α 2 m 11=
2
Volviendo a la expresión (9.24), para calcular los coeficientes con las expresiones (9.31) se
debe suponer el valor de una de ellas para calcular el resto. Suponiendo que se especifique el
valor de α2, las ecuaciones se resuelven para:
α1 = 1 − α 2
(9.32)
1
l 1= m 11=
2α 2
Si se considera que α2 es igual a ½, las ecuaciones (9.32) pueden resolverse y se obtiene

α1=1/2 y l1=m11=1, sustituyendo estos valores en la ecuación (9.24) se tiene:
1 1
y i + 1 = y i + ( b1 + b 2 ) h (9.33)
2 2
donde:
b1=f(xi,yi) (9.33a)
b2=f(xi + h, yi + b1h) (9.33b)
Un algoritmo para este método es el que sigue.
120
Cálculo Numérico
Algoritmo para aproximar la solución de un problema de valor inicial por el

método de Runge-Kutta de 2º orden suponiendo α2 = ½,
y' = f ( x, y ), x 0 ≤ x ≤ x n , y( x 0 ) = y 0 : problema de valor inicial
f(x,y):función que especifica la ecuación diferencial
X: valor de la variable independiente
x0: valor inicial de la variable independiente
xn: valor final de la variable independiente
Y: valor de la variable dependiente
y0: valor inicial de la variable dependiente
h: tamaño de paso
N: número entero que indica el número de segmentos igualmente espaciados (es
igual al número de datos menos 1)
Paso 1: Ingresar x0,xn, N, y0
(x 0 − xn )
Paso 2: Hacer h =
N
Paso 3: Hacer X= x0; Y= y0
Paso 3: Para i=1, 2, ..., N seguir los pasos 4 y 5
Paso 4: Tomar b1= f(x,y);
b2=f(X + h, Y+ b1h)
h
Paso 5: Tomar Y = Y + (b1 + b 2 ) Y=Y+h*f(x,y)
2
X=X+h
Paso 6: la SALIDA es (X,Y)
PARAR
Si se considera que α2 es igual a 2/3 las ecuaciones (9.32) pueden resolverse y se obtiene
α1=1/3 y l1=m11=3/4, sustituyendo estos valores en la ecuación (9.24) se tiene:
1 2
y i + 1 = y i + ( b1 + b 2 ) h (9.34)
3 3
donde:
b1=f(xi,yi) (9.34a)
3 3
b2= f(xi + h, yi + hb1 ) (9.34b)
4 4
121
Cálculo Numérico
Una versión común es la dada por:
1
y i + 1 = y i + ⎡⎢ (b1 + 4b 2 + b 3 )⎤⎥h (9.35)
⎣6 ⎦
b1=f(xi,yi) (9.35a)
1 1
b2= f(xi + h, y i + hb1 ) (9.35b)
2 2
b3= f(xi + h, y i − hb1 + 2hb 2 ) (9.35c)
Una versión común es la dada por:
1
y i + 1 = y i + ⎡⎢ (b1 + 2b 2 + 2b 3 + b 4 )⎤⎥h (9.36)
⎣6 ⎦
b1=f(xi,yi) (9.36a)
1 1
b2= f(xi + h, y i + hb1 ) (9.36b)
2 2
1 1
b3= f(xi + h, y i + hb 2 ) (9.36c)
2 2
b4= f(xi + h, y i + hb 3 ) (9.36d)
9.2.2.1.4 Método de Heun
Este es un método que permite mejorar la aproximación de la pendiente e implica el cálcu-

lo de dos derivadas, se promedian las mismas y se obtiene una aproximación mejorada de la pen-
diente en el intervalo completo. Puede decirse que el método de Heun utiliza un esquema predic-
tor-corrector.
El valor de la pendiente, es decir
y’i=f(xi,yi) (9.37)
se toma para extrapolar linealmente a yi+1, obteniéndose la ecuación predictora:
y pi +1 = y i + f ( x i , y i ) h (9.38)
este valor de y se utiliza para calcular una pendiente aproximada al final del intervalo:
y' i + 1 = f ( x i + 1, y pi + 1 ) (9.39)
Se promedian las pendientes obtenidas en las ecuaciones (9.37 y 9.39) y se tiene:
122
Cálculo Numérico
p
y '+ y ' f ( xi, yi ) + f ( xi+1, yi+1)
y' = i i+1 = (9.40)
2 2
Esta pendiente se utiliza en la ecuación correctora que permite calcular el valor de la varia-
ble dependiente de la función:
f ( xi, yi ) + f ( xi+1, ypi+1)

yi+1 = yi + h (9.41)
2
El algoritmo correspondiente es el especificado para el método de Runge-Kutta de 2º or-

den con α2 = ½.
¾ Por ejemplo si se desea encontrar la solución del problema de valor inicial, ya resuelto por el
método de la serie de Taylor y de Euler para x=1:
⎧⎪ dy
=x+y
⎨ dx
⎪⎩y(0) = 1
La pendiente en (x0,y0) es 1
Se toma este resultado para extrapolar linealmente a y1, considerando h=0,1se obtiene la ecua-
ción predictora:
y 1p = y 0 + 1h = 1 + 0,1 = 1,1
este valor de y se utiliza para calcular una pendiente aproximada al final del intervalo:
y' i + 1 = f ( x i + 1, y pi + 1 ) =0,1+1,1=1,2
Se promedian las pendientes obtenidas en las ecuaciones y se tiene:
1 + 1,2
y' = = 1,1
2
Esta pendiente se utiliza en la ecuación correctora que permite calcular el valor de la varia-
ble dependiente de la función:
y 1 = y 0 + y' h = 1 + 1,1 * 0,1 = 1,11 y así se va repitiendo el cálculo hasta llegar al valor deseado de y,
a continuación se muestra una tabla con los valores obtenidos aplicando este método.
123
Cálculo Numérico
Método de Método Error

Heun analítico verdadero
xi yi yi
0 1,000 1 0
0,1 1,110 1,110 0
0,2 1,242 1,243 0,001
0,3 1,398 1,399 0,001
0,4 1,582 1,584 0,002
0,5 1,795 1,797 0,002
0,6 2,041 2,044 0,003
0,7 2,323 2,327 0,004
0,8 2,646 2,651 0,005
0,9 3,012 3,019 0,007
1 3,428 3,436 0,008
9.2.2.2. Métodos de pasos múltiples
Los métodos de pasos múltiples se basan en que una vez que los cálculos han comenza-
do, los datos obtenidos en puntos anteriores sirven de guía. Se caracterizan porque no se inician
por sí mismos, requieren de valores previos que deben ser obtenidos por otros métodos, los cam-
bios de paso son complicados.
9.3. ECUACIÓN DIFERENCIAL LINEAL DE SEGUNDO ORDEN
Una ecuación diferencial de 2º orden con coeficientes constantes puede expresarse

como sigue:
d2 y dy
a2 + a1 + a 0 y( t ) = f ( t ) (9.42)
2 dt
dt
Normalizando, es decir dividiéndola por el coeficiente de la derivada de mayor orden se

tiene:
d2 y dy
+ b1 + b 0 y( t ) = h( t ) (9.43)
dt2 dt
9.3.1. Métodos de resolución
9.3.1.1. Método del operador diferencial
Si se utiliza el operador diferencial se puede escribir esta ecuación de la siguiente manera:
D2 y(t) + b1 D y(t) + b0 y(t) = h(t) (9.44)
124
Cálculo Numérico
y factoreando:
(D2 + b1 D + b0) y(t) = h(t) (9.45)
Si se reemplaza por λ a este operador polinómico, y considerando la función excitación h(t)

igual a 0, se llega a la (ecuación característica):
λ2 + b1λ + b0 =0 (9.46)
que se puede interpretar como una ecuación cuadrática y si sus raíces son: λ1 y λ2, el paréntesis
de la ecuación (9.45) puede expresarse en función de las mismas:
(D - λ1) (D - λ2) y(t) = h(t) (9.47)
Las raíces de la ecuación característica se calculan mediante la ecuación (9.48):
2
b ⎛b ⎞
λ 1,2 = − 1 ± ⎜ 1 ⎟ − b 0 (9.48)
2 ⎝ 2 ⎠
que puede expresarse también como λ1,2= α ± β
Se debe recordar que la solución general es la suma de la solución complementaria (yc(t))

y de la particular (yp(t)):
y(t) = yc(t) + yp(t) (9.50)
9.3.1.1.1. Solución complementaria
Considerando la ecuación (9.47), si la función excitación es nula, se obtiene la solución

complementaria:
(D - λ1) (D - λ2) yc(t) = 0 (9.51)
y según el tipo de raíces λ1 y λ2 se tienen tres casos:
1. Caso sobre-amortiguado, raíces reales y distintas
2. Caso oscilatorio amortiguado, raíces complejas y distintas
3. Caso crítico, raíces reales e iguales.
9.3.1.1.1.1. Caso sobre-amortiguado
Las raíces λ1 y λ2 son reales y distintas, la solución es:
y c ( t ) = K 1e λ 1 t + K 2 e λ 2 t (9.52)
125
Cálculo Numérico
donde K1 y K2 son constantes de integración que pueden determinarse a partir de condiciones

iniciales. En la Figura 9.2 se muestra un esquema del caso sobreamortiguado considerando las
constantes K1 y K2 positivas.
K 1e λ1t + K 2 e λ 2 t
K 2 e λ 2t
K 1e λ1t
Figura 9.2. Esquema de un caso sobreamortiguado
9.3.1.1.1.2. Caso crítico
Las raíces λ1 y λ2 son reales e iguales, es decir: λ1=λ2, la solución es:
y c ( t ) = (K 1 + K 2 t ) e λ1 t (9.53)
En la Figura 9.3 se muestra un esquema del caso crítico considerando las constantes K1 y
K2 positivas.
K 1e λ1t + K 2 te λ1t
K 1e λ1t
K 2 te λ1t
Figura 9.3. Esquema de un caso crítico
126
Cálculo Numérico
9.3.1.1.1.3. Caso oscilatorio amortiguado
Las raíces λ1 y λ2 son complejas y distintas, es decir: λ1,2= α ± i β, la solución puede ex-
presarse de cuatro formas diferentes:
(
y c ( t ) = e −α t K1 eiβ t + K 2 e −iβ t ) (9.54)
y c ( t ) = e −α t ( A cos β t + B sen β t ) (9.55)
y c ( t ) = C e −α t cos(β t − ϕ) (9.56)
y c ( t ) = C e −α t sen(β t + φ) (9.57)
En la figura 9.4 se muestra un esquema de este caso.
K 1e λ1t
K 1e λ1t cos β t
− K 1e λ1t
Figura 9.4. Esquema de un caso oscilatorio amortiguado
9.3.1.1.2. Solución particular
Esta solución se expresa como:
2
yp ( t ) = eλ1 t ∫ eλ 2 t ∫ h( t ) e −( λ1 + λ 2 ) t (d t ) (9.58)
9.3.1.1.3. Solución general
Como se expresó anteriormente, la solución general es la suma de la solución complemen-

taria y de la particular, esta suma puede esquematizarse como:
127
Cálculo Numérico
⎧k1eλ1t + k 2eλ 2 t ⎫
⎪ ⎪ 2
y( t ) = ⎨(k1 + k 2 t )eλ1t λ t λ t −( λ + λ ) t
⎬ + e 1 ∫ e 2 ∫ h( t ) e 1 2 (d t ) (9.59)
⎪Ce− αt cos(β t − ϕ)⎪
⎩ ⎭
9 Por ejemplo si se desea calcular la solución del siguiente problema de valor inicial:
y’’ – 4 y’ + 4y =0 con y(0) = 3; y’(0) = 1.
Como es una ecuación homogénea, no tiene función excitación, se necesita calcular sola-
mente la solución complementaria. Aplicando la ecuación (9.48) se obtienen las raíces λ1 y λ2 que
son reales e iguales, es decir: λ1=λ2=2, por lo tanto el caso que se presenta es el caso crítico y la
solución es:
y c ( t ) = (K 1 + K 2 t ) e 21 t
Para calcular las constantes K1 y K2, se aplican las condiciones iniciales:
y(0) = (K 1 + K 2 0 ) e 2. 0 = 3 ⇒ K 1 = 3
Para aplicar la segunda condición inicial se debe derivar la solución, es decir:
(
y' ( t ) = 2K 1e 21 t + K 2 e 21 t + 2 K 2 te 21 t )
Y particularizando para y’(0) = 1, se tiene:
( )
y' (0) = 2 ⋅ 3e 0 + K 2 e 0 + 0 = 1 ⇒ K 2 = 1 − 6 = −5
La solución complementaria que también es la solución general se escribe como:
y c ( t ) = (3 − 5t )e 2t
9.3.1.2. Método de los coeficientes indeterminados
Este método se emplea en ecuaciones diferenciales en las cuales la función excitación

contiene un polinomio, términos de la forma sen r x, cos r x, erx donde r es constante, o combina-
ciones de sumas y productos de estos.
Consiste en hacer una hipótesis inteligente de la forma de la solución particular y sustituir

esta función, que en general involucra uno o más coeficientes desconocidos, en la ecuación dife-
rencial. Para que este método sea aplicable con éxito, se deben determinar los coeficientes des-
conocidos, de tal manera que la función satisfaga a la ecuación diferencial. Por ello se dice que se
basa en la habilidad para suponer la forma general de la solución particular.
Si la ecuación es de la forma:
128
Cálculo Numérico
d2 y dy
a +b + c y = h( x ) (9.60)
2 dx
dx
La función excitación puede tener alguna de las formas que se presentan a continuación:
Pn(x) = a0xn + a1x n-1 + … + an
eαx Pn(x)
h(x) = eαx Pn(x) sen βx (9.61)
eαx Pn(x) cos βx
9
Si h(x) = a0xn + a1x n-1 + … + an, la solución particular se supone como:
yp(x) = A0 xn + A1 x n-1 + An (9.62)
9
Si h(x) = eαx Pn(x) la solución particular se supone como:
yp(x) = eαx (A0 xn + A1 x n-1 + An) (9.63)
cos βt
αx n n-1
9 Si h(x) = e (a0x + a1x + … + an) sen βt , la solución es:
yp(x) = eαx (A0 xn + A1 x n-1 + An) cos βt + eαx (A0 xn + A1 x n-1 + An) sen βt (9.64)
Por ejemplo si se desea calcular la solución particular de la ecuación diferencial:

y’’ + 4y’ + 9y = x2 + 3x
Se supone que la solución particular es:
yp(x) = A0 x2 + A1 x + A2
Se deriva esta solución y se llega a:
y’p(x) = 2A0 x + A1
y se deriva de nuevo:
y’’p(x) = 2A0
Se reemplazan estos valores en la ecuación diferencial original y se tiene:
2A0 +4(2A0 x + A1)+9(A0 x2 + A1 x + A2)= x2 + 3x
Distribuyendo y agrupando términos se obtiene:
9A0 x2 + (8A0+ 9A1) x +2A0 + 4 A1 + A2= x2 + 3x
129
Cálculo Numérico
Por igualdad de polinomios se llega a que:

A0=1/9; A1=19/81; A2= - 94/729
La solución particular es:
1 2 19 94
yp(x) = x + x- A2
9 81 729
9.4. ECUACIONES DIFERENCIALES PARCIALES
Una forma de representar una ecuación diferencial parcial de segundo orden con dos va-
riables independientes es:
∂ 2u ∂ 2u ∂ 2u ∂u ∂u
A +B +C +E +F + G u( x, y ) = H( x, y ) (9.65)
∂ x2 ∂ x ∂y ∂y 2 ∂x ∂y
Estas ecuaciones se clasifican de acuerdo a la relación de los coeficientes A, B, y C, en

los siguientes tipos:
9 Elíptica si B2 – 4AC < 0
9 Parabólica si B2 – 4AC = 0
9 Hiperbólica si B2 – 4AC > 0
A continuación se presentan los ejemplos típicos de cada uno de estos tipos:
9.4.1. Ecuaciones elípticas
Ecuación de Laplace: Ecuación de Poisson:
∂ 2u ∂ 2u ∂ 2u ∂ 2u
+ =0 + = F(x, y )
∂ x2 ∂ y2 ∂ x2 ∂ y2
u(a, y ) = f 1 ( y )
u(b, y ) = f 2 ( y )
u( x, c ) = g 1( x )
u( x, e ) = g 2 ( x )
9.4.2. Ecuación Parabólica
Ecuación de calor:
∂ 2u 1 ∂u
=
∂x 2
α 2 ∂t
u(0, t ) = T1
u(L, t ) = T2
u(x,0 ) = f ( x )
130
Cálculo Numérico
9.4.3. Ecuación Hiperbólica
Ecuación de onda:
∂ 2u 1 ∂ 2u
=
∂x2 α ∂ t2
u(0, t ) = 0
u(L, t ) = 0
u( x,0 ) = f ( x )
∂u
(x,0 ) = g(x )
∂t
Para solucionar estas ecuaciones se utiliza el Método de separación de variables, por

ejemplo si se tienen una ecuación diferencial parcial, de segundo orden, con dos variables inde-
pendientes, u(x,t) la misma puede expresarse como el producto de dos funciones, una que de-
pende solo de x y otra que depende solo de t.
1.1. En una población aislada p(t), la rapidez de crecimiento de la población dp/dt es igual a una
dp
función de la población; así: = f ( p)
dt
La expresión más simple para f(p) es f(p) = ε p(t), donde ε es una constante positiva. Determine: la
población en función del tiempo si p(0) = p0
1.2. Se determinó experimentalmente que un cierto material radioactivo decae con una rapidez
que es proporcional a la cantidad presente. Un bloque de este material tiene originalmente una
masa de 100 g y se observa que al cabo de 20 años es de 75 g.
a) Determine una expresión para la masa en función del tiempo.
b) Determine el tiempo de vida media del material.
c) Represente gráficamente masa en función del tiempo.
1.3. Un cuerpo de masa M se lanza verticalmente con una velocidad v0. Se asume que el aire pre-
senta una fuerza de roce proporcional a la velocidad con un coeficiente de roce γ.
Determine la expresión que permite evaluar la altura máxima y el tiempo que tarda en alcanzarla.
Represente gráficamente.
1.4. De acuerdo con la Ley de enfriamiento de Newton, la velocidad a la que se enfría una sustan-
cia al aire libre es proporcional a la diferencia entre la temperatura de la sustancia y la del aire. Si
la temperatura del aire es 28 ºC y la de la sustancia se enfría desde 100 ºC a 75ºC en 15 minutos,
indique ¿cuánto tardará dicha sustancia en alcanzar 50ºC? Resuelva analíticamente y por los
métodos de la serie de Taylor, Euler, Heun y Runge-Kutta de 4º orden.
131
Cálculo Numérico
1.5. Un trozo de carne de 2,3 kg, originalmente a 10ºC, se lleva a un horno a 300 ºC a las 20
horas, después de 75 minutos se determinó mediante termocupla la temperatura de la carne que
fue de 168 ºC ¿a qué hora se alcanzarán 212ºC , temperatura necesaria para que la carne este
cocida término medio? Resuelva analíticamente y por los métodos de la serie de Taylor, Euler,
Heun y Runge-Kutta de 4º orden.
1.6. Utilizando el método de los coeficientes indeterminados encuentre la solución particular para:
a) y’’ + 4 y = 9 e2x
b) y’’ + 4 y’ +48 y = x2 +2x
c) y’’ + 12 y’ +16 y = 24 sen 8x
1.7. Se encontró experimentalmente que un peso de 8 N estira un resorte 8 cm. Si el peso se lleva
5 cm por debajo de la posición de equilibrio y se suelta:
a) Establezca la ecuación diferencial y las condiciones asociadas que describan el movimien-
to;
b) Encuentre la posición del peso como una función del tiempo;
c) Determine la posición, velocidad, y aceleración del peso 0,6 s después de haberlo soltado.
132
Cálculo Numérico
BIBLIOGRAFIA CONSULTADA
AYRES, F Jr. 1991. Ecuaciones diferenciales. Editorial Mc Graw – Hill/ Interamericana de

México, S.A. México.
BOYCE, W. E.; DI PRIMA, R. C. 1974. Ecuaciones diferenciales y problemas con valores

en la frontera. Editorial LIMUSA. México.
BURDEN, R. L.; FAIRES, J. D. 1985. Análisis Numérico. Grupo editorial Interamericana de

México, S.A. de C. V. México.
CHAPRA, S. C.; CANALE, R. P. 1994. Métodos Numéricos para Ingenieros. Editorial Mc

Graw – Hill/ Interamericana de México, S.A. México.
EDWARDS, C. H. Jr., PENNEY, D. E. 1993. Ecuaciones diferenciales aplicadas. Editorial

Prentice – Hall Hispanoamericana S.A. México. 3º edición. México.
KREYSZIG., E. 1991. Matemáticas avanzadas para Ingeniería. Vol I y II. Editorial LIMUSA,
S. A. de C. V. México.
LASBAINES A. 2003. Introducción al Cálculo Numérico. Cuaderno de Laboratorio de Teor-

ía de Sistemas Nº 6. UNSE. Santiago del Estero. Argentina
LUTHE, R.; OLIVERA, A.; SCHUTZ, F. 1995. Métodos Numéricos. Editorial LIMUSA, S. A.
de C. V. México.
PITA, C. 1989. Ecuaciones diferenciales: una introducción con aplicaciones. Editorial LI-
MUSA, S. A. de C. V. México.
SPIEGEL, M. R. 1997. Ecuaciones diferenciales aplicadas. Editorial Prentice – Hall Hispa-

noamericana S.A. México.
SPIEGEL, M. R. 1997. Transformadas de Laplace. Editorial Mc Graw – Hill/ Interamericana

de México, S.A. México.
WILLIAMS, J. 1975. Transformadas de Laplace. Vol. 10. Editorial LIMUSA. México.
133

Calculo Numerico

Cargado por

Copyright:

Formatos disponibles

Calculo Numerico

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Calculo Numerico

Cargado por

Copyright:

Formatos disponibles

Universidad Nacional de Santiago del Estero

Facultad de Agronomía y Agroindustrias

Curso práctico con aplicaciones

Dra. Lucrecia Lucía Chaillou

Dra. Lucrecia Lucía Chaillou

METODOS NUMERICOS, MODELOS Y ALGORITMOS

La resolución de problemas de Ingeniería está asociada, por lo general, a resultados

1.1. METODOS NUMÉRICOS

1. Son herramientas poderosas para la solución de problemas. Permiten manejar sis-

2. Su teoría es la base de programas de métodos numéricos.

3. Su conocimiento permite diseñar programas propios.

4. Son un vehículo eficiente para aprender a servirse de las computadoras persona-

5. Son un medio para reforzar la comprensión de las matemáticas.

Estos métodos permiten:

9 Encontrar las raíces de ecuaciones lineales y no lineales

9 Resolver grandes sistemas de ecuaciones algebraicas lineales

9 Encontrar aproximaciones de funciones

9 Realizar interpolaciones para encontrar valores intermedios en tablas de datos

9 Aproximar derivadas de cualquier orden

9 Integrar cualquier función

9 Resolver problemas de valor inicial y de frontera

9 Obtener soluciones numéricas para ecuaciones diferenciales parciales

9 Realizar ajustamiento de curvas a datos

1.2. MODELOS MATEMÁTICOS

El mundo real es naturalmente complejo y en muchas ocasiones los problemas a resolver

Un modelo mental puede definirse como una representación de la realidad en la que se

Un modelo matemático es una formulación o ecuación que expresa las características

1. Describe un sistema o proceso natural en términos matemáticos.

2. Representa una idealización y una simplificación de la realidad.

3. Conduce a resultados predecibles y, en consecuencia, puede utilizarse para

La generación de un modelo matemático involucra dos etapas fundamentales, la de con-

1.2.1. Clasificación de modelos matemáticos

Los modelos matemáticos pueden clasificarse en función del tratamiento de la incertidum-

a) En función del tratamiento de la incertidumbre

Estocástico: probabilístico, no se conoce el resultado esperado, sino su probabilidad y existe por

b) En función del origen de la información utilizada para construirlos

c) En función de su campo de aplicación

Modelo matemático de optimización: los modelos matemáticos de optimización son ampliamen-

d) En función del factor tiempo

Modelos estáticos: son independientes del tiempo, consideran situaciones estacionarias.

Para generar un algoritmo se debe seguir una serie de pasos:

1. Determinar el objetivo de la tarea

2. Identificar los datos de entrada y de salida

3. Determinar el proceso involucrado

4. Identificar las variables internas

5. Dividir el proceso en acciones elementales

6. Determinar la secuencia de estas acciones

7. Incorporar estructuras de control

Por lo general, el objetivo del algoritmo será el de implementar un procedimiento numérico

Un algoritmo se puede representar mediante un pseudocódigo que especifica los datos

Los símbolos que se utilizan en diagramas de flujo se representan en la Figura 1.1.

Bloque Terminal: indica el inicio o finalización del algoritmo

Bloque de Proceso: representa los cálculos o manipulación de datos

Bloque de Entrada/Salida: indica la entrada o salida de datos e informa-

Bloque de Decisión: representa una comparación o pregunta que deter-

Bloque de Iteración: representa cálculos repetitivos

Conector: indica un truncamiento en el camino del diagrama de flujo

Figura 1.1. Símbolos que se utilizan en diagramas de flujo

Algoritmo Mayor (algoritmo que muestra el mayor