ApuntesMN2016 2

Descargar como pdf o txt
Descargar como pdf o txt
Está en la página 1de 84

Universidad Autónoma de Baja California

Facultad de Ingenierı́a, Arquitectura y Diseño


Bioingenierı́a

Apuntes de Métodos Numéricos


Dra. Dora-Luz Flores
Ensenada, Baja California, México 2016
Métodos Numéricos

Dra. Dora-Luz Flores

10 de enero de 2017
Índice general

1. Conceptos básicos 2
1.1. Uso de los métodos numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2. Errores numéricos y propagación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1. Errores inherentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2. Errores de truncamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.3. Errores de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.4. Cifras significativas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3. Exactitud y precisión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.1. Convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.2. Recursividad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.3. Series y sucesiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.4. Criterio de convergencia y divergencia . . . . . . . . . . . . . . . . . . . . . 6
1.4. Modelos matemáticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1. Algoritmos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2. Estabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.3. Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.4. Serie binomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.5. Serie de McLaurin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2. Solución numérica de ecuaciones de una sola variable 9


2.1. Método gráfico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1. Ejemplo de aproximación gráfica . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2. Método de bisecciones sucesivas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1. Ejemplo del método de bisección . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3. Método de interpolación lineal (Regla Falsa) . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1. Ejemplo del método de interpolación lineal . . . . . . . . . . . . . . . . . . . 13
2.4. Método de Newton-Raphson de primer orden . . . . . . . . . . . . . . . . . . . . . . 13
2.5. Método de Newton-Raphson modificado para raı́ces múltiples . . . . . . . . . . . . 14
2.5.1. Ejemplo del método de Newton-Raphson para raı́ces múltiples . . . . . . . . 15
2.6. Método de Von Mises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.6.1. Ejemplo de von Mises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6.2. Ejercicio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.7. Método de Birge-Vieta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

i
Dra. Dora-Luz Flores Métodos Numéricos

2.7.1. Método de Horner (división sintética) . . . . . . . . . . . . . . . . . . . . . . 19


2.7.2. Método de Birge-Vieta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.7.3. Ejemplo de Birge-Vieta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.7.4. Ejercicio de Birge-Vieta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3. Solución numérica de sistemas de ecuaciones 24


3.1. Método de eliminación de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.1.1. Algoritmo de eliminación de Gauss . . . . . . . . . . . . . . . . . . . . . . . 29
3.1.2. Código en Matlab del método de Gauss . . . . . . . . . . . . . . . . . . . . . 29
3.2. Método de eliminación de Gauss-Jordan . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.1. Algoritmo de eliminación de Gauss-Jordan . . . . . . . . . . . . . . . . . . . 31
3.2.2. Diagrama de flujo del método de Gauss-Jordan . . . . . . . . . . . . . . . . 31
3.3. Inversa de matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.3.1. Procedimiento para inversa de una matriz . . . . . . . . . . . . . . . . . . . 33
3.3.2. Ejemplo de la inversa de una matriz . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.3. Ejercicios propuestos para inversión de matrices . . . . . . . . . . . . . . . . 34
3.4. Método iterativo de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.4.1. Algoritmo del método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4.2. Ejemplo del método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5. Método de eliminación de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5.1. Algoritmo de eliminación de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . 37
3.5.2. Ejemplo de eliminación de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . 39
3.5.3. Ejercicios propuestos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4. Aproximación polinomial y funcional 42


4.1. Método de interpolación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.1. Diferencias finitas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.1.2. Diferencias divididas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2. Método de interpolación de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.2.1. Ejemplo del polinomio de Newton . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.2. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3. Método de interpolación de Lagrange de primer orden . . . . . . . . . . . . . . . . . 46
4.3.1. Ejemplo del método de interpolación de Lagrange . . . . . . . . . . . . . . . 47
4.4. Métodos de interpolación mediante polinomios de grado n . . . . . . . . . . . . . . 48
4.5. Método de mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.5.1. Regresión lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.5.2. Ejemplo de regresión lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.5.3. Linealización de regresiones . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.5.4. Regresión polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.5.5. Ejemplo de regresión polinomial . . . . . . . . . . . . . . . . . . . . . . . . . 54

ii
Dra. Dora-Luz Flores Métodos Numéricos

5. Integración numérica 55
5.1. Método analı́tico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.2. Método de la Regla del Trapecio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.2.1. Ejemplo del método de la regla del trapecio . . . . . . . . . . . . . . . . . . 57
5.2.2. Aplicación múltiple de la regla trapezoidal . . . . . . . . . . . . . . . . . . . 58
5.2.3. Ejemplo de la aplicación múltiple de la regla del trapecio . . . . . . . . . . . 60
5.3. Método Simpson 1/3 y 3/8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3.1. Regla de Simpson 1/3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.3.2. Ejemplo del método Simpson 1/3 . . . . . . . . . . . . . . . . . . . . . . . . 62
5.3.3. Regla de Simpson 3/8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.3.4. Ejemplo del método Simpson 3/8 . . . . . . . . . . . . . . . . . . . . . . . . 63
5.4. Método de diferenciación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.4.1. Diferenciación hacia adelante, hacia atrás y centrada . . . . . . . . . . . . . 64
5.4.2. Ejemplo de diferenciación numérica . . . . . . . . . . . . . . . . . . . . . . . 66

6. Solución numérica de ecuaciones diferenciales 67


6.1. Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.1.1. Ejemplo del método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.1.2. Análisis del error para el método de Euler . . . . . . . . . . . . . . . . . . . 70
6.1.3. Ejemplo del cálculo del error del método de Euler . . . . . . . . . . . . . . . 72
6.2. Método de Euler mejorado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.2.1. Método de Heun . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.2.2. Ejemplo del método de Heun . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3. Métodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.1. Métodos de Runge-Kutta de segundo orden . . . . . . . . . . . . . . . . . . 77
6.3.2. Métodos de Runge-Kutta de tercer orden . . . . . . . . . . . . . . . . . . . . 78
6.3.3. Métodos de Runge-Kutta de cuarto orden . . . . . . . . . . . . . . . . . . . 78
6.3.4. Ejemplo del método de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . 79

1
Unidad 1

Conceptos básicos

1.1. Uso de los métodos numéricos


En el campo de las ciencias e ingenierı́a, existen infinidad de fenómenos que requieren repre-
sentarse mediante modelos matemáticos. Desafortunadamente, la gran mayorı́a de estos modelos
no tiene una solución exacta ó no es fácil encontrarla. Es estos casos es en donde los métodos
numéricos proporcionan una solución aproximada al problema original. Un método numérico es
aquel que obtiene números que se aproximan a los que se obtendrı́an aplicando la solución analı́tica
de un problema.
Los métodos numéricos son herramientas extremadamente poderosas para la solución de pro-
blemas. Son capaces de manejar sistemas de ecuaciones grandes, no linealidades geométricas com-
plicadas que son comunes en la practica de la ingenierı́a y que, a menudo, son imposibles de resolver
analı́ticamente.

1.2. Errores numéricos y propagación


Deberá ser suficiente, en la práctica de la ingenierı́a y de las ciencias, contar con una solución
aproximada a un problema por las siguientes razones:
Los modelos matemáticos son aproximados esto es, simplificaciones al problema real. No se
toman en cuenta todos los factores que afectan a un fenómeno. Por ejemplo, en el caso del tiro
parabólico, se suele despreciar la resistencia del aire, sin embargo, esta puede ser importante.
Los modelos matemáticos requieren de parámetros, los cuales la mayorı́a de las veces provie-
nen de mediciones experimentales y estas, solo tienen una precisión limitada, que depende
del instrumento de medición. Por ejemplo la constante de los gases ideales. También pueden
provenir de cálculos y estos tienen una precisión limitada que depende tanto del método
como del instrumento de cálculo que se utilicen. Por ejemplo π.
Los modelos matemáticos resultantes son imposibles de resolver por métodos analı́ticos y se
debe de aproximar la solución numéricamente. Por ejemplo una ecuación de quinto grado. Por lo
anterior, se debe aceptar que siempre se tendrán errores presentes, estos pueden clasificarse en:

2
Dra. Dora-Luz Flores Métodos Numéricos

Errores inherentes.
Errores de truncamiento.
Errores de redondeo.

1.2.1. Errores inherentes


Los errores inherentes son aquellos que tienen los datos de entrada de un problema, y son
debidos principalmente a que se obtienen experimentalmente, debiéndose tanto al instrumento de
medición, como a las condiciones de realización del experimento. Por ejemplo, sı́ el experimento
es a temperatura constante y no se logra esto mas que en forma aproximada. También pueden
deberse a que se obtengan
√ de cálculos previos. Por ejemplo el valor calculado es el de un número
irracional como π o 2.

1.2.2. Errores de truncamiento


Los errores de truncamiento se originan por el hecho de aproximar la solución analı́tica de un
problema, por medio de un método numérico. Por ejemplo al evaluar la función exponencial por
medio de la serie de Taylor, se tiene que calcular el valor de la siguiente serie infinita:

x2 x3 xN X xN
ex = 1 + x + + + ... + =
2! 3! N ! N =0 N !
Ante la imposibilidad de tomar todos los términos de la serie, se requiere truncar después de
cierto número de términos. Esto nos introduce ciertamente un error, que es el error de truncamiento.
Este es independiente de la manera de realizar los cálculos. Solo depende del método numérico
empleado.

1.2.3. Errores de redondeo


Los errores de redondeo, se originan al realizar los cálculos que todo método numérico o analı́tico
requieren y son debidos a la imposibilidad de tomar todas las cifras que resultan de operaciones
aritméticas como los productos y los cocientes, teniendo que retener en cada operación el número
de cifras que permita el instrumento de cálculo que se este utilizando. Por ejemplo al calcular el
valor de , tenemos que conformarnos solo con la mayor cantidad de cifras 3, que maneje nuestro
instrumento de calculo.
Los errores anteriores también suelen denominarse como las fuentes de error. La magnitud
del error generada por alguna o todas las fuentes de error mencionadas anteriormente, se puede
cuantificar con ayuda de los siguientes parámetros:
Error.
Error relativo.
Error porcentual.

3
Dra. Dora-Luz Flores Métodos Numéricos

Error
El error se define como la diferencia entre el valor real (Vr ) y una aproximación a este valor Va :

e = Vr − Va

Error relativo
El error relativo se define como el cociente del error entre el valor real Vr (si Vr 6= 0):
e Vr − Va
er = =
Vr Vr
En ciertos métodos numéricos se utilizan esquemas iterativos para calcular resultados. En tales
esquemas, se hace una aproximación en base a la aproximación anterior. Este proceso se repite
varias veces, o de forma iterativa, para calcular sucesivamente más y mejores aproximaciones.
En tales casos, el error a menudo se calcula como la diferencia entre aproximación previa y la
actual por lo tanto, el error relativo porcentual o error porcentual esta dado por el error porcentual.

Error porcentual
El error porcentual es simplemente el error relativo expresado en por ciento ( %).

Vr − Va
ep = ∗ 100 %
Vr
En 1966 Scarberough demostró que si el siguiente criterio se cumple puede tenerse la seguridad
de que el resultado es correcto en al menos n cifras significativas.

Es = 0.5 × 102−n

1.2.4. Cifras significativas


El concepto de cifras significativas se ha desarrollado para designar formalmente la confiabilidad
de un valor numérico. El número de cifras significativas es el número de dı́gitos que se puede usar
con plena confianza. Por ejemplo se puede calcular un número irracional con varias cifras, pero de
ellas no todas, sobre todo las últimas pueden tomarse con plena confianza de que son correctas.
Por otro lado, los ceros no siempre son cifras significativas ya que pueden usarse solo para ubicar al
punto decimal. Por ejemplo los siguientes números tienen todos 4 cifras significativas: 0.00001985,
0.0001985, 0.001985, 1985, 19.85. Para asegurar que un cero nos represente una cifra significativa,
es común emplear la notación cientı́fica. Por ejemplo los siguientes números tienen 3, 4 y 5 cifras
significativas: 4.53 × 10−5 , 4.530 × 10−5 y 4.5300 × 10−5 . También se suele poner explı́citamente
los ceros. Los siguientes números tienen 5 cifras significativas: 19850, 0.019850, 19.850.

4
Dra. Dora-Luz Flores Métodos Numéricos

1.3. Exactitud y precisión


Los errores asociados con los cálculos y mediciones se pueden caracterizar observando su pre-
cisión y exactitud. La mayorı́a de la gente piensa que estos términos son sinónimos, pero no es ası́.
La precisión se refiere al número de cifras significativas que representan una cantidad. La exactitud
se refiere al grado de aproximación que se tiene de un número o de una medida al valor verdadero
que se supone representa, es decir, que tan cerca estamos del valor buscado. Por ejemplo, sı́ leemos
la velocidad del velocı́metro de un auto, esta tiene una precisión de 3 cifras significativas y una
exactitud de ± 5 Kph.

1.3.1. Convergencia
Velocidad de convergencia (rapidez o razón de convergencia): Es el número de iteraciones que
requiere un cálculo o algoritmo para converger o aproximarse a un valor. Es decir, la convergencia
se refiere al hecho de que los métodos numéricos obtienen n términos de una sucesión de valores.
Comenzamos con un valor inicial que sea una aproximación de la solución de un problema x0 .
Aplicando un método numérico se obtiene otra aproximación x1 . Se repite el procedimiento para
obtener x2 y ası́ sucesivamente, es decir, se generar la sucesión x0 , x1 , · · · , xn (todos los términos
son aproximaciones a la solución del problema). Sı́ la sucesión obtenida al cabo de n iteraciones
tiende a un lı́mite se dice que el método es convergente o divergente en caso contrario.

1.3.2. Recursividad
Fórmula recursiva: Relaciona términos sucesivos de una sucesión particular de números, fun-
ciones o polinomios, para proporcionar medios para calcular cantidades sucesivas en términos de
las anteriores.

1.3.3. Series y sucesiones


Serie infinita: 2, 4, 6, 8, · · ·
Serie finita: 2, 4, 6, 8, 10.

Sucesión aritmética
a + (a + d) + (a + 2d) + · · · + (a + (N − 1)d) = n(a + l) donde:
l = (N − 1)d y representa el último término de la sucesión.

Sucesión geométrica
a(1 − xN )
a + ax + ax2 + ax3 + · · · + axN −1 =
1−x
Se dice que una sucesión es creciente si:
an−1 ≤ an ∀n
y es decreciente si:

5
Dra. Dora-Luz Flores Métodos Numéricos

an−1 ≥ an ∀n

1.3.4. Criterio de convergencia y divergencia



X
Sea una serie infinita dada y sea {Sn } la sucesión de sumas parciales que definen esta serie
n=1
infinita. Entonces si el lı́m Sn , existe y es igual a S entonces se dice que la serie converge y que S
n→∞
es la suma infinita dada. Si coexiste al lı́m Sn , entonces se dice que la serie diverge o no converge
n→∞
y S no tiene valor. Una serie infinita es convergente si y solo si, la secuencia correspondiente es
convergente.

1.4. Modelos matemáticos


1.4.1. Algoritmos
Un algoritmo es una secuencia de pasos lógicos necesarios para llevar a cabo una tarea especı́fica,
generalmente los algoritmos se describen mediante un pseudocódigo.
Ejemplo. Algoritmo hecho en pseudocódigo del promedio de n números.

1. Pedir datos

2. Contar datos: n = números de datos.

3. Sumar los datos: suma = suma + dato(i)

4. Dividir suma entre n: prom = suma/n

5. Imprimir el prom.

Los algoritmos pueden ser estables e inestables.

1.4.2. Estabilidad
Algoritmos estables: Son aquellos en los que los cambios pequeños en los datos de entrada
generan cambios pequeños al final o a la salida. Algoritmos inestables: Son aquellos en los que los
cambios pequeños en la entrada producen grandes cambios en la salida.
Por ejemplo si en es un error en alguna etapa de un proceso y k es una constante independiente
de n el número de etapa, entonces si el error después de n operaciones se puede representar por
f (n) = knε , se dice que el crecimiento del error es lineal. Sı́ en cambio el error se representa por
f (n) = k n ε para k > 1, el crecimiento del error se dice que es exponencial.
El crecimiento del error lineal es por lo general inevitable, y cuando k y n son pequeños, los
resultados son aceptables. El crecimiento del error exponencial debe ser evitado, ya que el término
k n será grande, aun para valores relativamente pequeños de n. Por lo tanto si el crecimiento del
error es lineal el método es estable y sı́ es exponencial es inestable.

6
Dra. Dora-Luz Flores Métodos Numéricos

1.4.3. Serie de Taylor


La serie de Taylor permite predecir o calcular el valor de una función en un punto en términos
del valor de la función y sus derivadas en otro punto. Esto quiere decir que cualquier función suave
puede ser aproximada mediante un polinomio.

f 00 (a)(x − a)2 f 000 (a)(x − a)3 f n (a)(x − a)n


f (x) = f (a) + f 0 (a)(x − a) + + + ··· + + Rn
2! 3! n!
donde Rn es el término residual

f n−1 (ξ)hn+1
(n + 1)!
Algunas series tı́picas de Taylor son las siguientes:

x2 x3 x4 xn
1. ln(1 + x) = x − + − + ··· + para −1 ≤ x ≤ 1
2 3 4 n
 2  3  n
x−1 1 x−1 1 x−1 1 x−1
2. ln(x) = + + + ··· + para x ≥ 1/2
x 2 x 3 x n x

x3 x5 x7
3. sin(x) = x − + − + · · · para −∞ < x < ∞
3! 5! 7!

x2 x4 x6
4. cos(x) = 1 − + − + · · · para −∞ < x < ∞
2! 4! 6!

1 2 17 7
5. tan(x) = x + x3 + X 5 + x + · · · para x < π/2
3 15 315

1 x3 1 × 3 x5 1 × 3 × 5 x7
6. sin−1 (x) = x + + + + · · · para x < 1
2 3 2×4 5 2×4×6 7
3 5 7
 
π −1 π 1 x 1 × 3 x 1 × 3 × 5 x
7. cos−1 (x) = − sin (x) = − x + + + + · · · para x < 1
2 2 2 3 2×4 5 2×4×6 7

x3 x5 x7
8. sinh(x) = x + + + + · · · para −∞ < x < ∞
3! 5! 7!

x2 x4 x6
9. cosh(x) = 1 + + + + · · · para −∞ < x < ∞
2! 4! 6!

7
Dra. Dora-Luz Flores Métodos Numéricos

1.4.4. Serie binomial


n(n − 1)an−2 2 n(n − 1)(n − 2)an−3 3
(a + x)n = an + nan−1 x + x + x + ···
2! 3!

1.4.5. Serie de McLaurin


En matemáticas a menudo se pueden representar funciones mediante una serie infinita por
ejemplo la función exponencial se puede utilizar usando

x 2 x3 xn
ex = 1 + x + + + ··· +
2! 3! n!
que es conocida como expansión de serie de McLaurin, que es una modificación de la serie de
Taylor para cuando xi = 0.

8
Unidad 2

Solución numérica de ecuaciones de una


sola variable

Las soluciones de una ecuación f (x) = 0, se llaman ceros o raı́ces de f (x). En algunos casos
las raı́ces pueden ser obtenidas con métodos directos, por ejemplo, para una ecuación cuadrática
se utiliza la fórmula general. Aunque existen ecuaciones que no se pueden resolver directamente,
por ejemplo una función tan simple como f (x) = e−x − x. Para estos casos, la única alternativa es
una técnica de solución numérica.

2.1. Método gráfico


Un método simple para obtener una aproximación a la raı́z de la ecuación f (x) = 0, consiste
en graficar la función y observar en donde cruza al eje x. Este punto representa el valor de x para
el cual f (x) = 0 y proporciona una aproximación de la raı́z de la función f (x).

2.1.1. Ejemplo de aproximación gráfica


Usando la aproximación gráfica para obtener el coeficiente de razonamiento c, necesario para
que un paracaidista de masa = 68.1kg tenga una velocidad de 40m/s después de una caı́da libre
de 10 segundos, donde la función que representa este hecho esta dada por:
−c
 
gm  t
v(t) = 1−em  (2.1)
c

En este caso se reescribe la función de tal manera que sea igualada a cero, lo cual queda:
−c
 
gm  t
f (c) = 1−em −v (2.2)
c

Se obtiene la tabla 2.1.1, y con esos valores se genera la gráfica 2.1.

9
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 2.1: Valores generados de la ecuación 2.2


c f (c)
2 45.00718
4 34.19047
6 25.20892
8 17.71226
10 11.42152
12 6.11394
14 1.61112
16 -2.23026
18 -5.52565
20 -8.36838
22 -10.83416

Figura 2.1: Aproximación gráfica

2.2. Método de bisecciones sucesivas


Este método consiste en encerrar una raı́z entre un intervalo en el cual la función debe cruzar al
eje horizontal (eje x), e ir dividiendo el intervalo a la mitad hasta encontrar la mejor aproximación,
como se puede observar en la figura 2.2.
El algoritmo para determinar la raı́z por el método de bisección se describe a continuación:

1. Elegir lı́mites superior xu e inferior xi .


xi + xu
2. Obtener la aproximación a la raı́z xr = .
2

10
Dra. Dora-Luz Flores Métodos Numéricos

35

30

25

20
xi + xu
15 xr =
2
10

0
4 6 8 10 12 14 16 18 20 22
-5
xi xr xu
xi xr xu
-10
xi xu
-15

Figura 2.2: Método de bisección

3. Si f (xi )f (xr ) < 0, entonces xi = xi y xu = xr . Si no, si f (xu )f (xr ) < 0, entonces xi = xr y


xu = xu .

4. Si f (xu )f (xr ) = 0, la raı́z es igual a xr ; termina el cálculo.

2.2.1. Ejemplo del método de bisección


Resolviendo el ejemplo mostrado en el método de aproximación gráfica con la ecuación 2.2,
ahora utilizando el método de bisección con xi = 12 y xu = 16 se tienen las siguientes iteraciones
y luego se genera la tabla 2.2:
Primera iteración.
xi + xu 12 + 16
xr = = = 14
2 2
−ct (−12)(10)
   
gm  (9.81)(68.1) 
f (xi ) = f (12) = 1−e m −v = 1 − e 68.1  − 40 = 6.1139431
c 12

2.3. Método de interpolación lineal (Regla Falsa)


Un defecto del método de bisección es que al dividir el intervalo de xi a xu en mitades, no se
considera la magnitud de f (xi ) y de f (xu ). Por ejemplo, si f (xi ) esta más cercano a 0 que f (xu ),

11
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 2.2: Valores generados de la ecuación 2.2


i xi xu xr f (xi ) f (xu ) f (xr ) ker k
1 12.0000000 16.0000000 14.0000000 6.1139431 -2.2302607 1.6111164 -
2 14.0000000 16.0000000 15.0000000 1.6111164 -2.2302607 -0.3844581 6.6666667
3 14.0000000 15.0000000 14.5000000 1.6111164 -0.3844581 0.5936984 3.4482759
4 14.5000000 15.0000000 14.7500000 0.5936984 -0.3844581 0.0998300 1.6949153
5 14.7500000 15.0000000 14.8750000 0.0998300 -0.3844581 -0.1434972 0.8403361
6 14.7500000 14.8750000 14.8125000 0.0998300 -0.1434972 -0.0221312 0.4219409
7 14.7500000 14.8125000 14.7812500 0.0998300 -0.0221312 0.0387748 0.2114165
8 14.7812500 14.8125000 14.7968750 0.0387748 -0.0221312 0.0083032 0.1055966
9 14.7968750 14.8125000 14.8046875 0.0083032 -0.0221312 -0.0069187 0.0527704

es lógico pensar que la raı́z se encuentra más cerca de xi que de xu . El método de falsa posición
aprovecha la visualización gráfica de unir f (xi ) y f (xu ) con una recta, donde la intersección de
esta recta con el eje x representa una mejor estimación a la raı́z, como se muestra en la gráfica 2.3.

Figura 2.3: Método de la falsa posición

De la gráfica se observan triángulos semejantes:

f (xi ) f (xu )
= (2.3)
x r − xi xr − xu
de la ecuación anterior se despeja xr y se obtiene:

f (xu )(xi − xu )
xr = xu − (2.4)
f (xi ) − f (xu )

12
Dra. Dora-Luz Flores Métodos Numéricos

El algoritmo es igual al del método de bisección, lo único que cambia es la ecuación para xr .

2.3.1. Ejemplo del método de interpolación lineal


Tomando el mismo ejemplo mostrado en los métodos de aproximación gráfica y bisecciones
sucesivas, con la ecuación 2.2, ahora utilizando el método de interpolación lineal con xi = 12 y
xu = 16 se tienen las siguientes iteraciones y luego se genera la tabla 2.3:
Primera iteración.
xi + xu 12 + 16
xr = = = 14
2 2
−ct (−12)(10)
   
gm  (9.81)(68.1) 
f (xi ) = f (12) = 1−e m −v = 1 − e 68.1  − 40 = 6.1139431
c 12

Tabla 2.3: Valores generados de la ecuación 2.2


i xi xu xr f (xi ) f (xu ) f (xr ) ERP
1 12.0000000 16.0000000 14.0000000 6.1139431 -2.2302607 1.6111164 -
2 14.0000000 16.0000000 15.0000000 1.6111164 -2.2302607 -0.3844581 6.6666667
3 14.0000000 15.0000000 14.5000000 1.6111164 -0.3844581 0.5936984 3.4482759
4 14.5000000 15.0000000 14.7500000 0.5936984 -0.3844581 0.0998300 1.6949153
5 14.7500000 15.0000000 14.8750000 0.0998300 -0.3844581 -0.1434972 0.8403361
6 14.7500000 14.8750000 14.8125000 0.0998300 -0.1434972 -0.0221312 0.4219409
7 14.7500000 14.8125000 14.7812500 0.0998300 -0.0221312 0.0387748 0.2114165
8 14.7812500 14.8125000 14.7968750 0.0387748 -0.0221312 0.0083032 0.1055966
9 14.7968750 14.8125000 14.8046875 0.0083032 -0.0221312 -0.0069187 0.0527704

2.4. Método de Newton-Raphson de primer orden


Los métodos de Bisección y falsa posición son llamados métodos por intervalos en los cuales
los valores iniciales deben encerrar a la raı́z deseada. Los métodos siguientes son llamados métodos
de intervalo abierto dado que las condiciones iniciales no necesariamente tienen que contener a la
raı́z. Si el valor inicial de la raı́z es xi , entonces se puede trazar una tangente del punto f (xi ); el
punto donde esta tangente cruza al eje x representa una aproximación de la raı́z.
Si la pendiente en un punto dado se llama primera derivada de la función f 0 (x) y la pendiente
y2 − y1
de una recta es m = se puede tener que:
x2 − x1
y2 − y1 f (xi )
f 0 (x) = m = = (2.5)
x 2 − x1 xi − xi+1
Despejando xi+1 se tiene la ecuación de Newton-Raphson

13
Dra. Dora-Luz Flores Métodos Numéricos

Figura 2.4: Método de Newton-Raphson

f (xi )
xi+1 = xi − (2.6)
f 0 (xi )

2.5. Método de Newton-Raphson modificado para raı́ces


múltiples
Una raı́z múltiple corresponde a un punto donde una función es tangencial al eje horizontal x,
por ejemplo si f (x) es formada por una multiplicación de binomios iguales se encontrará una raı́z
repetida, como se puede observar en la figura 2.5.
La ecuación anterior tiene una raı́z doble porque un valor de x hace que dos términos de la
ecuación sean iguales a cero en x = 1, se observa que en la curva toca al eje x pero no lo cruza.
Para la ecuación f (x) = x4 − 6x3 + 12x2 − 10x + 3 se observa en la figura 2.6 que la función es
tangente al eje x. En general, la multiplicación impar de raı́ces cruza al eje horizontal x, mientras
que la multiplicidad par no lo hace. Las raı́ces múltiples ofrecen ciertas dificultades a los métodos
anteriormente expuestos.
Para el método de Newton-Raphson modificado, es necesario la obtención de la primera y
segunda derivada de la función y la obtención de la aproximación de la raı́z, esta dada por:

f (xi )f 0 (xi )
xi+1 = xi − (2.7)
[f 0 (xi )]2 − f (xi )f 00 (xi )
La manera en que se realizan las iteraciones es de la misma forma que el del método de Newton-
Raphson.

14
Dra. Dora-Luz Flores Métodos Numéricos

f (x)

Raíz Simple
Raíz Doble
1 2 3 x

Figura 2.5: Raı́ces múltiples, f (x) = x3 − 5x2 + 7x − 3 = (x − 3)(x − 1)2

f (x)

Raíz Simple
Raíz Triple
1 2 3 x

Figura 2.6: Raı́ces múltiples, f (x) = x4 − 6x3 + 12x2 − 10x + 3 = (x − 3)(x − 1)3

2.5.1. Ejemplo del método de Newton-Raphson para raı́ces múltiples


Utilizar el método de Newton-Raphson modificado para evaluar la raı́z múltiple de la ecuación
f (x) = x3 − 5x2 + 7x − 3 con un valor inicial de x0 = 0.
Obteniendo la primera y segunda derivada de la función se tiene que:

f 0 (x) = 3x2 − 10x + 7


f 00 (x) = 6x − 10
Primera iteración

15
Dra. Dora-Luz Flores Métodos Numéricos

x0 = 0

f (x0 ) = x3 − 5x2 + 7x − 3 = −3

f 0 (x0 ) = 3x2 − 10x + 7 = 7

f 00 (x0 ) = 6x − 10 = −10

(−3)(7)
x1 = 0 − = 1.1052
[7]2 − (−3)(−10)

Segunda iteración

x1 = 1.1052

f (x1 ) = x3 − 5x2 + 7x − 3 = −0.0209

f 0 (x1 ) = 3x2 − 10x + 7 = −0.3878

f 00 (x1 ) = 6x − 10 = −3.3684

(−0.0209)(−0.3878)
x2 = 1.1052 − = 1.0030
[−0.3878]2 − (−0.0209)(−3.3684)

1.0030 − 1.1052
er = ∗ 100 = 10.1867
1.0030

Tercera iteración

x2 = 1.0030

f (x2 ) = x3 − 5x2 + 7x − 3 = −0.000019

f 0 (x2 ) = 3x2 − 10x + 7 = −0.01229

f 00 (x2 ) = 6x − 10 = −3.9815

(−0.000019)(−0.01229)
x3 = 1.0030 − = 1.0000024
[−0.01229]2 − (−0.000019)(−3.9815)

1.0000024 − 1.0030817
er = ∗ 100 = 0.3079275
1.0000024

16
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 2.4: Valores generados de la ecuación f (x) = x3 − 5x2 + 7x − 3 usando el método de Newton-
Raphson para raı́ces múltiples
i xi f (xi ) f 0 (xi ) f 00 (xi ) xi+1 |er |
0 0.0000000 -3.0000000 7.0000000 -10.0000000 1.1052632 -
1 1.1052632 -0.0209943 -0.3878116 -3.3684211 1.0030817 10.1867572
2 1.0030817 -0.0000190 -0.0122982 -3.9815100 1.0000024 0.3079275
3 1.0000024 0.0000000 -0.0000095 -3.9999857 1.0000000 0.0002381
4 1.0000000 0.0000000 0.0000000 -4.0000000 1.0000000 0.0000000

2.6. 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 0 (xi ) sea próximo a cero. Para ello Von Mises sugirió
f (xi )
utilizar Newton-Raphson, xi+1 = xi − 0 , sustituyendo el denominador f 0 (xi ) por f 0 (x0 ), es
f (xi )
decir obtener geométricamente las siguientes aproximaciones por medio de paralelas a la primera
tangente, como se muestra en la figura 2.7.

Figura 2.7: Método de Von Mises

Por lo anterior, la ecuación de Von Mises es:

f (xi )
xi+1 = xi − (2.8)
f 0 (x0 )

17
Dra. Dora-Luz Flores Métodos Numéricos

2.6.1. Ejemplo de von Mises


Utilizar el método de von Mises para encontrar la raı́z de e−x − ln x empleando un valor inicial
de x0 = 1.
1
Se tiene que f (x) = e−x − ln(x) y su derivada evaluada en x0 es f 0 (x0 ) = f 0 (1) = −e−1 − =
x
−1.36787944.

Primera iteración

i=0

x0 = 1

f (x0 ) = e−1 − ln(1) = 0.36787944

f (x0 ) 0.36787944
x1 = x0 − 0
=1− = 1.26894142
f (x0 ) −1.36787944

Segunda iteración

i=1

x1 = 1.26894142

f (x1 ) = e−1.26894142 − ln(1.26894142) = 0.04294604

f (x1 ) 0.04294604
x2 = x1 − 0 = 1.26894142 − = 1.30033749
f (x0 ) −1.36787944

1.30033749 − 1.26894142
|er | = ∗ 100 = 2.41445529
1.30033749

Ası́ sucesivamente hasta completar la información que se muestra en la tabla 2.5.

Tabla 2.5: Iteraciones para evaluar la función f (x) = e−x − ln(x) con el método de von Mises
i xi f (xi ) xi+1 |er |
0 1 0.36787944 1.26894142 -
1 1.26894142 0.042946035 1.30033749 2.4144554
2 1.30033749 0.00981599 1.307513555 0.54883309

18
Dra. Dora-Luz Flores Métodos Numéricos

2.6.2. Ejercicio
Utilizar el método de von Mises para calcular la raı́z de la función f (x) = cos x − x con un
valor inicial de x0 = 1. Realizar al menos 6 iteraciones.

2.7. Método de Birge-Vieta


Un caso especial de importancia práctica es encontrar las raı́ces de la ecuación f (x) = 0 cuando
f (x) es un polinomio en x. En esta sección se verá el método Birge-Vieta que encuentra todas las
raı́ces reales de un polinomio.

2.7.1. Método de Horner (división sintética)


Suponer dos polinomios P (x) y Q(x) de la forma

P (x) = a1 xn + a2 xn−1 + . . . + an x + an+1 (2.9)

Q(x) = b1 xn + b2 xn−1 + . . . + bn x + bn+1 (2.10)


donde a1 6= 0. Si la relación entre P (x) y Q(x) esta dada por

P (x) = (x − x0 )Q(x) + bn+1 , (2.11)


se tiene que b1 = a1 , bn+1 = P (x0 ), y

bk = ak + bk−1 x0 , (2.12)
para k = 2, 3, . . . , n + 1.

Lo anterior puede realizarse mediante una tabla de la siguiente manera

x0 a1 a2 a3 ... an an+1
b1 x 0 b2 x 0 b2 x 0 ... bn−1 x0 bn x 0
b 1 = a1 b 2 = a2 + b 1 x 0 b 3 = a3 + b 2 x 0 ... bn = an + bn−1 x0 bn+1 = an+1 + bn x0

El polinomio P (x),

P (x) = a1 xn + a2 xn−1 + . . . + an x + an+1


puede ser representado por el vector de sus coeficientes,

a = [a1 a2 . . . an an+1 ]

de la misma manera Q(x) puede ser representado por el vector b(1 : n)

b = [b1 b2 . . . bn ].

19
Dra. Dora-Luz Flores Métodos Numéricos

Dado que

P (x) = (x − x0 )Q(x) + bn+1 ,


P 0 (x) = Q(x) + (x − x0 )Q0 (x).
Por lo tanto

P 0 (x0 ) = Q(x0 ),
es decir, que P (x0 ) puede evaluarse obteniendo el residuo de la división de Q(x) por (x − x0 )
y evaluando Q(x0 ).

2.7.2. Método de Birge-Vieta


Un polinomio de la forma,

P (x) = a1 xn + a2 xn−1 + . . . + an−1 x + an (2.13)


puede ser factorizado en la forma

P (x) = (x − p1 )(x − p2 ) . . . (x − pn ) (2.14)


donde pi es un cero (o raı́z) del polinomio porque P (pi ) = 0.

El método Birge-Vieta aplica Newton-Raphson para encontrar una raı́z del polinomio P (x).
Dado un punto xk , evalúa P (xk ) y P 0 (xk ) mediante división sintética. Cuando encuentra una raı́z
pi , elimina el factor (x − pi ) mediante división sintética y continúa trabajando sobre el polinomio
resultante. El proceso se repite hasta encontrar todas las raı́ces del polinomio.

2.7.3. Ejemplo de Birge-Vieta


Sea P (x) = x3 − 2x2 − 5x + 6. Valor inicial x = −(−5)/6 = 0.8333.

Primera iteración

0.8333 1 −2 −5 6
0.8333 −0.9722 −4.9769
1 −1.1667 −5.9722 1.0231

0.8333 1 −1.1667 −5.9722


0.8333 −0.2778
1 −0.3333 −6.2500

1.0231
x = 0.8333 − = 0.9970
−6.2500

20
Dra. Dora-Luz Flores Métodos Numéricos

0.9970 1 −2 −5 6
0.9970 −1.0000 −5.9822
1 −1.0030 −6.0000 0.0178

0.9970 1 −1.0030 −6.0000


0.9970 −0.0059
1 −0.0059 −6.0059

0.0178
x = 0.9970 − =1
−6.0059

x = 1 es la primera raı́z.

Se encuentra el siguiente polinomio eliminando el factor (x − 1) mediante división sintética.

1 1 −2 −5 6
1 −1 −6
1 −1 −6 0

El polinomio Q(x) resultante es x2 − x − 6.

Segunda iteración

Continuando con el polinomio x2 − x − 6 y con un valor inicial x = −(−1)/(−6) = −0.1667

−0.1667 1 −1 −6
−0.1667 0.1944
1 −1.1667 −5.8056

−0.1667 1 −1.1667
−0.1667
1 −1.3333

−5.8056
x = −0.1667 − = −4.5208
−1.3333

21
Dra. Dora-Luz Flores Métodos Numéricos

−4.5208 1 −1 −6
−4.5208 24.9588
1 −5.5208 18.9588

−4.5208 1 −5.5208
−4.5208
1 −10.0417

18.9588
x = −4.5208 − = −2.6328
−10.0417

−2.6328 1 −1 −6
−2.6328 9.5646
1 −3.6328 3.5646

−2.6328 1 −3.6328
−2.6328
1 −6.2656

3.5646
x = −2.6328 − = −2.0639
−6.2656

−2.0639 1 −1 −6
−2.0639 6.3237
1 −3.0639 0.3237

−2.0639 1 −3.0639
−2.0639
1 −5.1278

0.3237
x = −2.0639 − = −2.0008
−5.1278

−2.0008 1 −1 −6
−2.0008 6.0040
1 −3.0008 0.0040

−2.0008 1 −3.0008
−2.0008
1 −5.0016

22
Dra. Dora-Luz Flores Métodos Numéricos

0.0040
x = −2.0008 − = −2
−5.0016

x = −2 es la segunda raı́z.

Se encuentra el siguiente polinomio eliminando el factor (x + 2) mediante división sintética.

−2 1 −1 −6
−2 6
1 −3 0

El polinomio Q(x) resultante es x − 3. La tercera raı́z es x = 3.

2.7.4. Ejercicio de Birge-Vieta


Encontrar las raı́ces del polinomio x4 + 4x3 − 6x2 − 4x + 5.
Las raı́ces son: 1, -1, 1, -5.

23
Unidad 3

Solución numérica de sistemas de


ecuaciones

Conceptos básicos de matrices


Una matriz consiste de un arreglo rectangular de elementos representado por un solo sı́mbolo.
 
a11 a12 · · · a1m
 a21 a22 · · · a2m 
A =  .. (3.1)
 
.. . . .. 
 . . . . 
an1 an2 · · · anm
A es de n renglones por m columnas. Su dimensión es de n × m donde An×m . Si n = 1, se le
conoce como vector renglón;
 
b = b1 b 2 b 3 · · · bm
si m = 1, se le conoce como vector columna.
 
c1
 c2 
 
c =  c3 
 
 .. 
.
cn
Si m = n se le llama matriz cuadrada.
 
a11 a12 a13 a14
a21 a22 a23 a24 
A4×4 =
a31
 (3.2)
a32 a33 a34 
a41 a42 a43 a44
A los elementos a11 , a22 , a33 , a44 y en general a aii se les conoce como diagonal principal.

24
Dra. Dora-Luz Flores Métodos Numéricos

Tipos de matrices cuadradas


Matriz simétrica
Es una matriz donde aij = aji ∀i, j.
 
5 1 2
A = 1 3 7
2 7 8

Matriz diagonal
Es una matriz donde todos los elementos fuera de la diagonal principal son cero.
 
a11 0 0
A =  0 a22 0
0 0 a33

Matriz identidad
Es aquella matriz cuyos elementos de la diagonal principal son 1 y los demás 0.
 
1 0 0
A= 0 1 0
0 0 1

Matriz triangular superior


Es aquella matriz en la que todos los elementos por debajo de la diagonal principal son 0.
 
a11 a12 a13
A =  0 a22 a23 
0 0 a33

La matriz triangular inferior


Es aquella matriz en la que todos los elementos por encima de la diagonal principal son 0.
 
a11 0 0
A = a21 a22 0 
a31 a32 a33

Matriz banda, tribanda o tridiagonal


La matriz banda tiene todos los elementos igual a 0 excepto en una banda centrada sobre la
diagonal principal.

25
Dra. Dora-Luz Flores Métodos Numéricos

 
a11 a12 0 0
a21 a22 a23 0 
A= 
 0 a32 a33 a34 
0 0 a43 a44

Operaciones con matrices


Dos matrices de m×n son iguales, si y solo si, cada elemento se encuentra en la segunda matriz.
A = B si aij = bij ∀i, j.

Suma de matrices

C =A+B
(
i = 1, 2, 3, · · · , n
cij = aij + bij para
j = 1, 2, 3, · · · , m

Resta de matrices

C =A−B
(
i = 1, 2, 3, · · · , n
cij = aij − bij para
j = 1, 2, 3, · · · , m
La suma y la resta de matrices son conmutativas y asociativas.

Multiplicación por un escalar

C = gA

donde g es el escalar.

Multiplicación de dos matrices


La multiplicación de dos matrices se puede realizar si y solo si la primera matriz tiene el mismo
número de columnas que el número de renglones de la segunda matriz. El producto de dos matrices
se presenta como:

C = AB
Xn
cij = aik bkj
k=1

26
Dra. Dora-Luz Flores Métodos Numéricos

Inversa de una matriz


Aunque la multiplicación de matrices es posible, la división de matrices no esta definida. Sin
embargo, si la matriz es cuadrada y no singular (determinante diferente de cero) existe una matriz
A−1 llamada inversa de la matriz para la cual

AA−1 = I.

La multiplicación de una matriz por la inversa es análoga a la división, en el sentido de que un


número dividido por si mismo es igual a 1.

Matriz transpuesta
Si los renglones y las columnas de una matriz A se intercambian, entonces la matriz resultante
de n × m se conoce como la transpuesta de A y se denota por AT .

aij = aji
La traza de una matriz es la suma de los elementos de su diagonal principal.
n
X
tr[A] = aii
i=1

Determinante de una matriz


Para matrices de 2 × 2, la determinante es:

det(A) = a11 a22 − a12 a21

Ejemplo de sistemas de ecuaciones lineales


Muchos problemas pueden ser descritos mediante un sistemas de ecuaciones lineales. Por ejem-
plo, considere el circuito eléctrico mostrado en la figura 3.1.

Figura 3.1: Circuito eléctrico que puede ser descrito por un sistema de ecuaciones lineales

Las ecuaciones de malla que describen a este circuito son las siguientes:

27
Dra. Dora-Luz Flores Métodos Numéricos

15i1 − 5i2 = 20
−5i1 + 15i2 − 5i3 = 0
− 5i2 + 20i3 = 0
A partir de las ecuaciones de malla se pueden obtener todas las corrientes, voltajes y potencial
de los elementos del circuito.
Utilizando la notación de matrices se puede obtener una matriz aumentada con este ejemplo.
Se define R, i y v:
 
15 −5 0
R = −5 15 −5
0 −5 20
 
i1
i = i2 
i3
 
20
v=0
0
Se puede expresar el juego de ecuaciones como:
    
15 −5 0 i1 20
Ri = v = −5 15 −5  i2  =  0 
0 −5 20 i3 0
Y que a su vez puede representarse por la matriz aumentada:
 
15 −5 0 20
Rb = −5 15 −5 0 
0 −5 20 0

3.1. Método de eliminación de Gauss


Este método consiste en expresar el sistema como una matriz aumentada de la forma:

Ax = b (3.3)
La idea del método es llevar el sistema a la forma triangular superior y de allı́ despejar una
variable a la vez partiendo de la última. Él último paso se conoce como sustitución en reversa. Para
lograr llevar el sistema a la forma triangular superior, se emplean las operaciones elementales de
matrices como son el intercambio de renglones, división entre un escalar a cada renglón, ası́ como
suma y resta entre renglones.

28
Dra. Dora-Luz Flores Métodos Numéricos

3.1.1. Algoritmo de eliminación de Gauss


n : número de ecuaciones
aij : elementos de la matriz aumentada A
p : ı́ndice del elemento pivote
fi : fila i

1. Para i = 1, · · · , n − 1 seguir los pasos 2 a 4.

2. Sea p el menor entero con i ≤ p ≤ n y api 6= 0


Si p no puede encontrarse, SALIDA(’No existe solución única’)
PARAR

3. Si p 6= i intercambiar la fila p por la fila i

4. Para j = i + 1, · · · , n seguir los pasos 5 a 6


aij
5. Hacer mij =
aii
6. Realizar fj − mij fi e intercambiar por la fila fj

7. Si ann = 0 entonces SALIDA(’No existe solución única’)


PARAR
cn
8. Hacer xn =
ann
1  
ci − nj=i+1 aij xj
P
9. Para i = n − 1, · · · , 1 tomar xi =
aii
10. SALIDA x
PARAR

3.1.2. Código en Matlab del método de Gauss


1 function r e s = EGaussiana (A, v a r a r g i n )
2
3 % Implementación d e l método de e l i m i n a c i ó n g a u s s i a n a para r e s o l v e r s i s t e m a s
4 % de e c u a c i o n e s l i n e a l e s y e n c o n t r a r l a i n v e r s a de una m a t r i z .
5 %
6 % Inv = EGaussiana (A)
7 %
8 % Regresa Inv , l a i n v e r s a de A.
9 %
10 % x = EGaussiana (A, b )
11 %

29
Dra. Dora-Luz Flores Métodos Numéricos

12 % Regresa x , l a s o l u c i ó n d e l s i s t e m a Ax=b .
13
14 MAXPIVOTE = 1 ;
15
16 % Se c r e a l a m a t r i z aumentada
17 [ n ,m] = s i z e (A ) ;
18 i f length ( v a r a r g i n )>=1
19 b = varargin {1};
20 A = [A b ] ;
21 e l s e i f n==m
22 A = [A eye ( n ) ] ;
23 end
24 [ n ,m] = s i z e (A ) ;
25
26 i f n>m
27 error ( ’ n>m en l a m a t r i z aumentada ’ ) ;
28 end
29
30 for i =1:n
31 i f MAXPIVOTE
32 % Encontrar r e n g l o n de maximo p i v o t e
33 k = find ( abs (A( : , i ))==max( abs (A( i : n , i ) ) ) , 1 , ’ l a s t ’ ) ;
34
35 i f k˜= i
36 % Intercambiar renglones
37 r P i v o t e = A( k , : ) ;
38 A( k , : ) = A( i , : ) ;
39 A( i , : ) = r P i v o t e ;
40 end
41 end
42
43 % S a l i r s i determinante es cero
44 i f abs (A( i , i ) ) == 0
45 error ( ’ Determinante e s i g u a l a c e r o ’ )
46 end
47
48 % Hacer c e r o s en l a columna i d e b a j o de l a f i l a i
49 for j=i +1:n
50 A( j , : ) = A( j , : ) − A( i , : ) ∗ A( j , i ) /A( i , i ) ;
51 end
52 end
53
54 for i=n : −1:1

30
Dra. Dora-Luz Flores Métodos Numéricos

55 % Hacer uno e l e l e m e n t o i , i
56 A( i , : ) = A( i , : ) / A( i , i ) ;
57
58 % Hacer c e r o s en l a columna i a r r i b a de l a f i l a i
59 for j=i −1: −1:1
60 A( j , : ) = A( j , : ) − A( i , : ) ∗ A( j , i ) ;
61 end
62 end
63
64 r e s = A( : , n+1:m) ;

3.2. Método de eliminación de Gauss-Jordan


Jordan propuso una modificación al método de Gauss. En vez de llevar el sistema a la forma
triangular superior y de allı́ usar la sustitución en reversa, él pensó que serı́a más fácil continuar el
procedimiento de eliminación de elementos, es decir, él propuso eliminar los elementos tanto arriba
como abajo del pivote hasta llegar a la matriz identidad. De esta manera la solución del sistema
se puede leer directamente de la última columna de la matriz aumentada.

3.2.1. Algoritmo de eliminación de Gauss-Jordan


1. Determinar la primer columna (a la izquierda) no cero.

2. Si el primer elemento de la columna es cero, intercambiarlo por un renglón que no tenga cero.
Multiplicando apropiadamente el renglón igual a 1. Este primer 1 será llamado pivote.

3. Obtener ceros arriba y abajo del pivote sumando múltiplos adecuados a los renglones debajo
de renglón pivote en la matriz completa.

4. Cubrir la columna y el renglón de trabajo y repetir el proceso comenzando en el paso 1 con


la columna siguiente.

Es importante observar que en el método de Gauss-Jordan:

De forma general, la matriz se va escalonando y reduciendo a la vez.

En el paso 2, si el elemento no es cero no se realiza intercambio.

En el paso 3, los elementos que se hacen cero no solo son los inferiores al pivote (Eliminación
Gaussiana) sino también los superiores.

3.2.2. Diagrama de flujo del método de Gauss-Jordan


En la figura 3.2 se muestra el diagrama de flujo del método de eliminación de Gauss-Jordan.

31
Dra. Dora-Luz Flores Métodos Numéricos

Figura 3.2: Diagrama de flujo del método de Gauss-Jordan

3.3. Inversa de matrices


Este método es más teórico. Consiste en expresar el sistema como una ecuación matricial de la
forma Ax = b y despejar el vector columna x. Dado que no esta definida la división de matrices,
se usa la matriz inversa A−1 . Multiplicando por la matriz inversa ambos lados se tiene
A−1 Ax = A−1 b
de donde
Ix = A−1 b
y finalmente
x = A−1 b
El problema se reduce a hallar la matriz inversa para multiplicarla por el vector columna b y
ası́ hallar x.

32
Dra. Dora-Luz Flores Métodos Numéricos

3.3.1. Procedimiento para inversa de una matriz


Para hallar la matriz inversa se puede utilizar el siguiente procedimiento.

1. Se coloca la matriz A junto a una matriz identidad I del mismo tamaño, es decir,
 
a11 a12 · · · a1n 1 0 0 · · · 0
 a21 a22 · · · a2n 0 1 0 · · · 0
 
 a31 a32 · · · a3n 0 0 1 · · · 0
 
 .. .. ... .. .. .. .. . . .. 
 . . . . . . . .
an1 an2 · · · ann 0 0 0 · · · 1

2. Se aplica la eliminación de Gauss Jordán a la matriz A, las operaciones que se le hagan a la


matriz A, también se le aplican a I. La matriz A se convierte en I. Se puede demostrar que
matriz I se convierte en A−1 .

3. Una vez hallada A−1 se procede a multiplicarla por bx = A−1 b.

3.3.2. Ejemplo de la inversa de una matriz


Encontrar la inversa de
 
1 3 4
2 8 6 
5 1 35
Primero se verifica que sea una matriz no singular, calculando el determinante de A, es igual a
0.5. Por lo tanto se procede a aumentar la matriz con la matriz identidad del mismo tamaño.
 
1 3 4 1 0 0
2 8 6 0 1 0
5 1 35 0 0 1

Se utiliza el método de Gauss-Jordan para encontrar la inversa de la matriz.


 
1 3 4 1 0 0
0 2 −2 −2 1 0
5 1 35 0 0 1
 
1 3 4 1 0 0
0 2 −2 −2 1 0
0 −14 15 −5 0 1
 
1 3 4 1 0 0
0 1 −1 −1 1 0
2
0 −14 15 −5 0 1

33
Dra. Dora-Luz Flores Métodos Numéricos

4 − 32 0
 
1 0 7
0 1 −1 −1 1 0
2
0 −14 15 −5 0 1
4 − 23 0
 
1 0 7
0 1 −1 −1 1
2
0
0 0 1 −19 7 1
1 0 7 4 − 23 0
 
0 1 0 −20 15 1
2
0 0 1 −19 7 1
1 0 0 137 − 101
 
2
−7
0 1 0 −20 15 1
2
0 0 1 −19 7 1
Por lo tanto, la matriz inversa es

137 − 101
 
2
−7
−20 15 1
2
−19 7 1

3.3.3. Ejercicios propuestos para inversión de matrices


 
1 4 6
1. 3
 1 2
4 5 8
 
1 3 −2 1
1 3 −1 2
2.  
0 1 −1 4
2 6 1 2

3.4. Método iterativo de Jacobi


Este método consiste en despejar una de las incógnitas de una ecuación dejándola en función de
las otras. La manera más sencilla es despejar a x1 de la primer ecuación, x2 de la segunda ecuación,
xi de la i-ésima ecuación, hasta xn de la n-ésima ecuación. Es necesario por razones obvias que
todos los elementos de la diagonal principal de la matriz de coeficientes del sistema lineal, sean
diferentes de cero.

34
Dra. Dora-Luz Flores Métodos Numéricos

Sea el sistema lineal:


a11 x1 + a12 x2 + a13 x3 + · · · + a1n xn = b1
a21 x1 + a22 x2 + a23 x3 + · · · + a2n xn = b2
a31 x1 + a32 x2 + a33 x3 + · · · + a3n xn = b3
.. ..
. .
an1 x1 + an2 x2 + an3 x3 + · · · + ann xn = bn

Al realizar los despejes propuestos se tiene x1 de la primera ecuación, x2 de la segunda ecuación,


etc., entonces:

b1 − (a12 x2 + a13 x3 + · · · + a1n xn )


x1 =
a11
b2 − (a21 x1 + a23 x3 + · · · + a2n xn )
x2 =
a22
b3 − (a31 x1 + a32 x2 + · · · + a3n xn ) (3.4)
x3 =
a33
..
.
bn − (an1 x1 + an2 x2 + · · · + a1(n−1) xn−1 )
xn =
ann
Para estimar la primera aproximación a la solución se debe partir de un vector inicial, el cual
puede ser el vector x0 = 0, o algún otro que se encuentre próximo al vector de solución x. En
general, el vector de aproximación a la solución después de las iteraciones se puede calcular de la
siguiente manera:
n
aij xkj
P
bi −
j=1
j6=i
xk+1
i = (3.5)
aii

3.4.1. Algoritmo del método de Jacobi


Para aplicar el método, se considera una primera aproximación al valor de las incógnitas x, se
sustituye esta primera aproximación en los segundos miembros de las ecuaciones 3.4, de manera
sucesiva, hasta llegar a la última ecuación y encontrar xn . Se obtiene de esta manera una nueva
aproximación a los valores de las incógnitas. 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 converja.
La condición suficiente para que el método de Jacobi converja es que la matriz de coeficientes
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.

35
Dra. Dora-Luz Flores Métodos Numéricos

X
|aii | > |aij | (3.6)
En la figura 3.4 se presenta un algoritmo para este método iterativo.

Figura 3.3: Algoritmo del método de Jacobi

36
Dra. Dora-Luz Flores Métodos Numéricos

3.4.2. Ejemplo del método de Jacobi


Resolver el siguiente sistema de ecuaciones utilizando el método de Jacobi. Emplear el vector
inicial de x0 = 0.

6x1 − x2 − x3 + 4x4 = 17
x1 − 10x2 + 2x3 − x4 = −17
3x1 − 2x2 + 8x3 − x4 = 19
x1 + x2 + x3 − 5x4 = −14
Al despejar las incógnitas correspondientes se tiene:

17 − (−x2 − x3 + 4x4 )
x1 =
6
−17 − (x1 + 2x3 − x4 )
x2 =
−10 (3.7)
19 − (3x1 − 2x2 − x4 )
x3 =
8
−14 − (x1 + x2 + x3 )
x4 =
−5
Y finalmente se llena la tabla 3.4.2 con los valores calculados.

3.5. Método de eliminación de Gauss-Seidel


Los métodos de Gauss y Gauss-Jordan forman parte de los métodos directos o finitos. Al cabo
de un número finito de operaciones, en ausencia de errores de redondeo, se obtiene x∗ solución del
sistema
Ax = b
El método de Gauss-Seidel forma parte de los métodos llamados indirectos o iterativos. En
ellos se comienza con x0 = (x01 ; x02 ; · · · ; x0n ), una aproximación inicial de la solución. A partir de
x0 se construye una nueva aproximación de la solución, x1 = (x11 ; x12 ; · · · ; x1n ). A partir de x1 se
construye x2 (aquı́ el superı́ndice indica la iteración y no indica una potencia). Ası́ sucesivamente
se construye una sucesión de vectores xk , con el objetivo, no siempre garantizado, de que

lı́m xk = x∗
k→∞

3.5.1. Algoritmo de eliminación de Gauss-Seidel


Generalmente los métodos indirectos son una buena opción cuando la matriz es muy grande y
dispersa, es decir, cuando el número de elementos no nulos es pequeño comparado con el número
total de elementos de A.

37
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 3.1: Valores generados a partir de las ecuaciones 3.7


i x1 x2 x3 x4 |Ex1 | |Ex2 | |Ex3 | |Ex4 |
0 0 0 0 0 - - - -
1 2.833333 1.700000 2.375000 2.800000 100 100 100 100
2 1.645833 2.178333 2.087500 4.181667 72.151899 21.958684 13.772455 33.041052
3 0.756528 1.863917 2.825104 3.982333 17.550945 16.868601 26.108919 5.005441
4 0.959948 1.942440 3.055073 3.889110 21.190747 4.042524 7.527439 2.397042
5 1.073512 2.018098 2.986768 3.991492 10.578776 3.748981 2.286907 2.565018
6 1.006483 2.005556 2.975894 4.015676 6.659766 0.625399 0.365414 0.602230
7 0.986458 1.994260 3.000917 3.997587 2.030015 0.566434 0.833855 0.452505
8 1.000805 1.999071 3.003342 3.996327 1.433584 0.240665 0.080719 0.031519
9 1.002851 2.001116 2.999007 4.000643 0.203982 0.102221 0.144546 0.107896
10 0.999591 2.000022 2.999290 4.000595 0.326059 0.054703 0.009464 0.001219
11 0.999489 1.999758 3.000233 3.999781 0.010259 0.013216 0.031418 0.020349
12 1.000145 2.000017 3.000104 3.999896 0.065556 0.012983 0.004312 0.002879
13 1.000090 2.000046 2.999937 4.000053 0.005505 0.001409 0.005552 0.003930
14 0.999962 1.999991 2.999984 4.000014 0.012786 0.002727 0.001578 0.000967
15 0.999986 1.999992 3.000014 3.999987 0.002459 0.000028 0.000983 0.000675
16 1.000009 2.000003 3.000001 3.999998 0.002301 0.000553 0.000415 0.000273
17 1.000002 2.000001 2.999997 4.000003 0.000752 0.000064 0.000150 0.000108
18 0.999998 1.999999 3.000000 4.000000 0.000384 0.000104 0.000101 0.000067
19 1.000000 2.000000 3.000001 3.999999 0.000193 0.000024 0.000020 0.000014
20 1.000000 2.000000 3.000000 4.000000 0.000056 0.000018 0.000022 0.000015

38
Dra. Dora-Luz Flores Métodos Numéricos

 
a11 a12 · · · a1m = b1
 a21 a22 · · · a2m = b2 
A =  .. (3.8)
 
.. ... .. .. 
 . . . = .
an1 an2 · · · anm = bn

En cada iteración del método de Gauss-Seidel, hay n subiteraciones. En la primera subiteración


se modifica únicamente x1 . Las demás coordenadas x2 , x3 , · · · , xn no se modifican. El cálculo de
x1 se hace de tal manera que se satisfaga la primera ecuación.

b1 − (a12 x02 + a13 x03 + · · · + a1n x0n


x11 =
a11
x1i = x0i , i = 2, · · · , n.
En la segunda subiteración se modifica únicamente x2 . Las demás coordenadas x1 , x3 , · · · , xn
no se modifican. El cálculo de x2 se hace de tal manera que se satisfaga la segunda ecuación.

b2 − (a21 x11 + a23 x13 + · · · + a2n x1n


x22 =
a22
x2i = x1i , i = 1, 3, 4, 5, · · · , n.
Ası́ sucesivamente, en la n-ésima subiteración se modifica únicamente xn . Las demás coorde-
nadas x1 , x2 , · · · , xn−1 no se modifican. El cálculo de xn se hace de tal manera que se satisfaga la
n-ésima ecuación.

bn − (an1 xn1 + an3 x3n−1 + · · · + ann xnn−1


xnn =
ann
xni = xn−1
i , i = 1, 2, · · · , n − 1.

3.5.2. Ejemplo de eliminación de Gauss-Seidel


Utilizar el método de Gauss-Seidel para obtener la solución del sistema de ecuaciones siguiente:

3x1 − 0.1x2 − 0.2x3 = 7.85


0.1x1 + 7x2 − 0.3x3 = −19.3
0.3x1 − 0.2x2 + 10x3 = 71.4
Primero se despeja la incógnita de la diagonal para cada una de las ecuaciones:

39
Dra. Dora-Luz Flores Métodos Numéricos

7.85 + 0.1x2 + 0.2x3


x1 =
3
−19.3 − 0.1x1 + 0.3x3
x2 =
7
71.4 − 0.3x1 + 0.2x2
x3 =
10

Primera iteración. Se supone que x2 y x3 son cero.

1. x2 = 0

x3 = 0

7.85 + 0.1(0) + 0.2(0)


x1 = = 2.616667
3

2. x1 = 2.616667

x3 = 0

−19.3 − 0.1(2.616667) + 0.3(0)


x2 = = −2.794524
7

3. x1 = 2.616667

x2 = −2.794524

71.4 − 0.3(2.616667) + 0.2(−2.794524)


x3 = = 7.005610
10

Segunda iteración.

1. x2 = −2.794524

x3 = 7.005610

40
Dra. Dora-Luz Flores Métodos Numéricos

7.85 + 0.1x1 + 0.2x2


x1 = = 2.990557
3

2.990557 − 2.616667
|erx1 | = ∗ 100 = 12.502353 %
2.990557

2. x1 = 2.990557

x3 = 7.005610
−19.3 − 0.1x1 + 0.3x3
x2 = = −2.499625
7

−2.499625 − (−2.794524)
|erx2 | =
∗ 100 = 11.797729 %
−2.499625

3. x1 = 2.990557

x2 = −2.499625
71.4 − 0.3x1 + 0.2x2
x3 = = 7.000291
10

7.000291 − (7.005610)
|erx3 | = ∗ 100 = 0.0759826 %
7.000291

3.5.3. Ejercicios propuestos


Resolver los sistemas de ecuaciones utilizando el método de Gauss-Seidel.

17c1 − 2c2 − 3c3 = 500


1. −5c1 + 21c2 − 2c3 = 200 con un error porcentual del 5 %.
−5c1 − 5c2 + 22c3 = 30
 
10 2 −1 0 26
1 20 −2 3 −15
2. 
−2
 con al menos tres iteraciones.
1 30 0 53 
1 2 3 20 47
 
−1 2 10 11
3.  11 −1 2 12 con al menos tres iteraciones.
1 5 2 8

41
Unidad 4

Aproximación 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ı́tica 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 aproxima-
ción y para ello existen métodos numéricos que por lo general utilizan funciones racionales enteras
(polinomios), de manera que la curva descrita por los mismos toque todos los puntos definidos.
Si no se requiere gran aproximación se deriva una curva simple que represente el comportamiento
general de los datos.

4.1. Método de interpolación


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.

4.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 4.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.
Se denominan primeras diferencias hacia adelante y se representan con ∆yi a las diferencias
entre dos valores consecutivos de y, es decir:

42
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 4.1: Valores para la función y = f (x)


x y
x0 y0
x1 = x0 + h y1
x2 = x0 + 2h y2
x3 = x0 + 3h y3
··· ···
xn = x0 + nh yn

a0 = y1 − y0
a1 = y2 − y1
a2 = y3 − y2 (4.1)
···
an−1 = yn − yn−1
Las diferencias de las primeras diferencias se llaman segundas diferencias hacia adelante, ∆2 yi :

b0 = a1 − a0
b1 = a2 − a1
b2 = a3 − a2 (4.2)
···
bn−2 = an−1 − an−2
Las diferencias de las segundas diferencias se llaman terceras diferencias hacia adelante, ∆3 yi .

c0 = b 1 − b0
c1 = b 2 − b1
c2 = b 3 − b2 (4.3)
···
cn−3 = an−2 − an−3
Siguiendo este proceso se definen las cuartas, quintas, etc., diferencias hacia adelante. Todas las
diferencias pueden arreglarse en una tabla de diferencias (ver tabla 4.2), en donde cada diferencia
se indica entre los dos elementos que la producen.

Tabla 4.2: Tabla de diferencias divididas


xi yi ∆yi ∆2 yi ∆3 yi
x0 y0 a0 b0 c0
x1 y1 a1 b1 c1
x2 y2 a2 b2 c2
··· ··· ··· ···

43
Dra. Dora-Luz Flores Métodos Numéricos

Si una de estas diferencias se vuelve constante (o aproximadamente constante), puede decir-


se que los valores tabulados pueden describirse por un polinomio de grado igual al orden de la
diferencia constante (o aproximadamente constante).

4.1.2. Diferencias divididas


Si se considera la función y = f (x) definida en forma tabular, parecida a la presentada en
la tabla 4.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
función tal como:

Pn (x) = a0 + (x − x0 )a1 + (x − x0 )(x − x1 )a2 + · · · + (x − x0 )(x − x1 ) · · · (x − xn−1 )an (4.4)

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 ), esta diferencia
se puede denotar también como yr o bien fr .
La diferencia de primer orden o diferencia de orden uno es igual a:

f (xs ) − f (xr )
f [xr , xs ] = (4.5)
xs − xr

4.2. Método de interpolación de Newton


La fórmula general para la interpolación de Newton para un polinomio de orden n es

fn (x) = b0 + b1 (x − x0 ) + b2 (x − x0 )(x − x1 ) + · · · + bn (x − x0 )(x − x1 ) · · · (x − xn−1 ) (4.6)

Para un polinomio de n-ésimo orden se requieren n + 1 puntos. Donde:

b0 = f (x0 )
b1 = f (x1 , x0 )
b2 = f (x2 , x1 , x0 )
b3 = f (x3 , x2 , x1 , x0 ) (4.7)
..
.
bn = f (xn , xn−1 , · · · , x1 , x0 )
Donde las evaluaciones puestas entre paréntesis son diferencias divididas finitas, es decir:

f (xi ) − f (xj )
f (xi , xj ) = (4.8)
xi − xj
La segunda diferencia dividida es:

44
Dra. Dora-Luz Flores Métodos Numéricos

f (xi , xj ) − f (xj , xk )
f (xi , xj , xk ) = (4.9)
xi − xk
La n-ésima diferencia dividida es:

f (xn , xn−1 , · · · , x2 , x1 ) − f (xn−1 , xn−2 , · · · , x1 , x0 )


f (xn , xn−1 , · · · , x1 , x0 ) = (4.10)
xn − x0
Estas diferencias se utilizan para evaluar los coeficientes b0 , b1 , · · · , bn , los cuales se utilizan
para obtener el polinomio de interpolación, el cual es conocido como el polinomio de interpolación
por diferencias divididas de Newton. De la ecuación 4.2 se obtiene:

fn (x) = f (x0 ) + f (x1 , x0 )(x − x0 ) + f (x2 , x1 , x0 )(x − x0 )(x − x1 ) + · · ·


(4.11)
+f (xn , xn−1 , · · · , x1 , x0 )(x − x0 )(x − x1 ) · · · (x − xn−1 )

De manera tabular se podrı́a expresar como se muestra en la tabla 4.2.

Tabla 4.3: Tabla que muestra las diferencias divididas del polinomio de Newton de orden 4
i xi f (xi ) Primero Segundo Tercero Cuarto
0 x0 f (x0 ) f (x1 , x0 ) f (x2 , x1 , x0 ) f (x3 , x2 , x1 , x0 ) f (x4 , x3 , x2 , x1 , x0 )
1 x1 f (x1 ) f (x2 , x1 ) f (x3 , x2 , x1 ) f (x4 , x3 , x2 , x1 )
2 x2 f (x2 ) f (x3 , x2 ) f (x4 , x3 , x2 )
3 x3 f (x3 ) f (x4 , x3 )
4 x4 f (x4 )

4.2.1. Ejemplo del polinomio de Newton


Utilizar la interpolación de polinomios de Newton para interpolar los puntos x0 = 1, x1 =
4, x2 = 6, x3 = 5 a un polinomio de tercer orden para estimar ln(2).

Tabla 4.4: Datos de las diferencias divididas obtenidas para estimar el valor de ln(2)
i xi f (xi ) f (xi , xi−1 ) f (xi , xi−1 , xi−2 ) f (xi , xi−1 , xi−2 , xi−3 )
0 1 0 0.46209812 −0.05187311 0.00786553
1 4 1.38629436 0.20273255 −0.02041100
2 6 1.79175947 0.18232156
3 5 1.60943791

Donde las diferencias divididas son:

b0 = f (x0 ) = 0
b1 = f (x1 , x0 ) = 0.46209812

45
Dra. Dora-Luz Flores Métodos Numéricos

b2 = f (x2 , x1 , x0 ) = −0.05187311
b3 = f (x3 , x2 , x1 , x0 ) = 0.00786553

Y el polinomio resultante queda:

f3 (x) = 0 + 0.46209812(x − 1) − 0.05187311(x − 1)(x − 4) + 0.00786553(x − 1)(x − 4)(x − 6)

Obteniendo el valor aproximado para x = 2, se tiene:

f3 (2) = 0.62876858

Donde el error es:



0.693147 − 0.62876858
|er | = ∗ 100 % = 9.287888 %
0.693147

4.2.2. Ejercicios
Usando los valores de la tabla siguiente, calcular el polinomio usando solo los 4 primeros puntos
y después calcular otro polinomio usando todos los puntos.

Tabla 4.5: Datos para calcular polinomios de Newton


i x f (x)
0 1 0
1 3 1.0986122
2 4 1.3862944
3 5 1.6094379
4 6 1.7917595

4.3. Método de interpolación de Lagrange de primer orden


El polinomio de interpolación de Lagrange se puede obtener de manera directa a partir de la
formulación del polinomio de Newton. Se hará esto únicamente para el caso del polinomio en primer
orden. Para obtener la forma de Lagrange, se reformulan las diferencias divididas. Por ejemplo, la
primera diferencia dividida se puede reformular como:

f (x1 ) − f (x0 )
f (x1 , x0 ) = (4.12)
x1 − x0
La cual es referida como:

f (x0 ) f (x1 )
f (x1 , x0 ) = + (4.13)
x0 − x1 x1 − x0

46
Dra. Dora-Luz Flores Métodos Numéricos

Por último, al agrupar términos similares y simplificar se tiene la forma del polinomio de
Lagrange:
x − x1 x − x0
f (x) = f (x0 ) + f (x1 ) (4.14)
x0 − x1 x1 − x 0
La interpolación de polinomios de Lagrange es simplemente una reformulación del polinomio de
Newton que evita el cálculo por diferencias divididas. Se puede expresar de manera concisa como:
n
X
f (x) = Li (x)f (xi ) (4.15)
i=0

donde:
n
Y x − xj
Li (x) = (4.16)
x − xj
j=0 i
j6=i
Q
Donde designa el producto de, por ejemplo, cuando n = 1 es
x − x1 x − x0
f (x) = f (x0 ) + f (x1 ) (4.17)
x0 − x1 x1 − x 0
Cuando n = 2 es

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )


f (x) = f (x0 ) + f (x1 ) + f (x2 ) (4.18)
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

y ası́ sucesivamente.
Para los casos en donde el orden del polinomio se desconozca, el método de Newton tiene
ventajas debido a que profundiza en el comportamiento de las diferentes fórmulas de orden superior.
En general puede integrarse fácilmente en los cálculos de Newton ya que la aproximación usa una
diferencia dividida. De esta forma, desde el punto de vista de cálculo, a menudo, se prefiere el
método de Newton.

4.3.1. Ejemplo del método de interpolación de Lagrange


Utilizar la interpolación de polinomios de Lagrange para interpolar los puntos x0 = 1, x1 =
4, x2 = 6 a un polinomio de primer orden para estimar ln(2).

x0 = 1 f (x0 ) = 0
x1 = 4 f (x1 ) = 1.386294
x2 = 6 f (x2 ) = 1.791759
Para n = 1 sustituyendo en la ecuación 4.17 los punto se tiene:

47
Dra. Dora-Luz Flores Métodos Numéricos

x − x1 x − x0
f (x) = f (x0 ) + f (x1 )
x0 − x1 x1 − x 0
x−4 x−1
f (x) = (0) + (1.386294)
1−4 4−1
para x = 2
2−4 2−1
f (2) = (0) + (1.386294) = 0.462098
1−4 4−1
Donde el error verdadero es:

Vv − Va 0.693147 − 0.462098
|Ev | = = × 100 = 33.333333 %
Vv 0.693147

4.4. Métodos de interpolación mediante polinomios de gra-


do n
Cuando n se va incrementando, se obtienen polinomios de mayor grado. Esto hace que en
muchas ocasiones el error se minimice y se obtenga una mejor aproximación a los datos que se
quieren mostrar mediante un polinomio.
Tomando el ejemplo anterior y usando n = 2, se tiene:

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )


f (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

(x − 4)(x − 6) (x − 1)(x − 6) (x − 1)(x − 4)


f (x) = (0) + (1.386294) + (1.791759)
(1 − 4)(1 − 6) (4 − 1)(4 − 6) (6 − 1)(6 − 4)
Para x = 2

f (2) = 0.565844
Donde el error verdadero es |Ev | = 18.3659 %.

4.5. Método de mı́nimos cuadrados


En este tipo de aproximación (también llamada aproximación funcional) 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 aproximadamente a esta curva.

48
Dra. Dora-Luz Flores Métodos Numéricos

4.5.1. Regresión lineal


El ejemplo más simple de aproximación por mı́nimos cuadrados es el ajuste de un conjunto de
datos a una lı́nea recta.
La expresión matemática de una recta es:

y = a0 + a1 x + e (4.19)
en donde a0 y a1 son coeficientes que representan la intersección con el eje y y la pendiente,
respectivamente, e es el error o diferencia entre el modelo y las observaciones. Reordenando, se
puede calcular el error como:

e = y − a0 − a1 x (4.20)
es decir, es la diferencia entre el valor real de y y el valor aproximado, a0 + a1 x 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
X n
X
Sr = e2i = (yi − a0 − a1 xi )2 (4.21)
i=1 i=1
Para encontrar los valores de a0 y a1 que minimicen la ecuación 4.21 se debe derivar esta
ecuación con respecto a los coeficientes indicados, es decir:
n
∂Sr X
= −2 (yi − a0 − a1 xi ) (4.22)
∂a0 i=1
n
∂Sr X
= −2 [(yi − a0 − a1 xi )xi ] (4.23)
∂a1 i=1
Para generar un mı́nimo, se igualan estas derivadas a cero y se expresan como un conjunto de
dos ecuaciones lineales con dos incógnitas a0 y a1 .
P P
P na0 + P x2i ai = P yi (4.24)
x i a0 + x i ai = xi yi
Si se resuelve este sistema se obtiene:
P P P
n xi y i − xi y i
a1 = (4.25)
n x2i − ( xi )2
P P

a0 = y − a1 x (4.26)
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:

49
Dra. Dora-Luz Flores Métodos Numéricos

r
Sr
Sy/x = (4.27)
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:
St − Sr
r2 = (4.28)
St
donde: St = (yi − y)2 es la cantidad de dispersión en la variable dependiente que existe antes de
P
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 el que se muestra en la figura 4.5.1.

4.5.2. Ejemplo de regresión lineal


Encontrar la ecuación de regresión lineal para los datos de la tabla siguiente que representan
las temperaturas (x) y las ventas (y) de bebidas refrescantes de cierta compañı́a.

x 5 7 10 12 16 20 23 27 19 14 9 6
y 9 11 15 16 20 24 27 29 22 20 14 9

Se formula la tabla 4.5.2 para obtener los valores de a0 y a1 .

x y xy x2
5 9 45 25
7 11 77 49
10 15 150 100
12 16 192 144
16 20 320 256
20 24 480 400
23 27 621 529
27 29 783 729
17 22 374 289
14 20 280 196
9 14 126 81
6 9 54 36
Suma 166 216 3502 2834
Media 13.833 18

50
Dra. Dora-Luz Flores Métodos Numéricos

Figura 4.1: Algoritmo para regresión lineal.

P P P
n xi yi − xi yi
a1 =
n x2i − ( xi )2
P P

(12)(3502) − (166)(216)
a1 =
(12)(2834) − (166)2
a1 = 0.9559826
Para obtener el valor de a0 , se realizan las siguientes operaciones:

a0 = y − a1 x

a0 = 18 − (0.9559826)(13.833)

51
Dra. Dora-Luz Flores Métodos Numéricos

a0 = 4.7755735
Para calcular el valor de R2 , primero se obtiene la tabla 4.5.2 con los valores de Sr y St .

e2i Sti
0.308565441 81
0.218511328 49
0.441693325 9
0.061189521 4
0.00508308 4
0.010977531 36
0.056086455 81
2.518901563 121
0.946187383 16
3.388064428 4
0.385122968 16
2.284539481 81
Sr = 10.6249225 St = 502

St − Sr
r2 =
St
502 − 10.6249225
r2 =
502
2
r = 0.9788348

4.5.3. Linealización de regresiones


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 = d1 eb1 x (4.29)
en donde d1 y b1 son constantes. Para linealizar este modelo se aplican logaritmos naturales, es
decir:

ln y = ln d1 + b1 ln e (4.30)
En una gráfica semilogarı́tmica de ln y con x se genera una lı́nea recta con pendiente b1 y
ordenada al origen: ln d1 .

52
Dra. Dora-Luz Flores Métodos Numéricos

Modelo de ecuación elevada a una potencia

y = d 2 xb 2 (4.31)
En donde d2 y b2 son coeficientes, puede linealizarse mediante logaritmos en base 10, o sea:

log y = b2 log x + log d2 (4.32)


de forma que en una gráfica logarı́tmica de log y y 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 (4.33)
b3 + x
Los coeficientes d3 y b3 son constantes y puede linealizarse si se invierte la ecuación 4.33, es
decir:
1 b3 1
= + (4.34)
y d3 x d3
1 1
Una gráfica de con será lineal, con pendiente b3 /d3 y ordenada al origen 1/d3 .
y x

4.5.4. 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 = a0 + a1 x + a2 x 2 + · · · + am x m (4.35)
En este caso la suma de cuadrados de los residuos es:
n
X
Sr = (yi − a0 − a1 xi − a2 x2i − · · · − am xm
i )
2
(4.36)
i=0

Siguiendo el procedimiento anterior, se deriva la ecuación con respecto a cada uno de los
coeficientes del polinomio, para obtener:
n
∂Sr X
= −2 (yi − a0 − a1 xi − a2 x2i − · · · − am xm
i ) (4.37)
∂a0 i=1
n
∂Sr X
= −2 (yi − a0 − a1 xi − a2 x2i − · · · − am xm
i )xi (4.38)
∂a1 i=1
···

53
Dra. Dora-Luz Flores Métodos Numéricos

n
∂Sr X
= −2 (yi − a0 − a1 xi − a2 x2i − · · · − am xm m
i )xi (4.39)
∂a1 i=1

Si estas ecuaciones se igualan a cero y se reordenan se obtiene un conjunto de ecuaciones


normales:

+ a2 P x2i + am P x m
P P P P
Pna0 + a1 P x i + ··· i = P yi
a0 P x i + a1 P x2i + a2 P x3i + ··· + am P xm+1
i = P xi yi
a0 x2i + a1 x3i + a2 x4i + ··· + am xm+2i = x2i yi (4.40)
···
xm xm+1 xm+2 x2m xm
P P P P P
a0 i + a1 i + a2 i + ··· + am i = i yi

El error de regresión polinomial se calcula con:


s
Sr
Sy/x = (4.41)
n − (m + 1)
Y el coeficiente de determinación con:
Sv − Sr
r2 = (4.42)
Sv

4.5.5. Ejemplo de regresión polinomial


Utilizar el mismo ejemplo anterior para encontrar el polinomio de orden 2.

54
Unidad 5

Integración numérica

La integración de una función dentro del ámbito de la ingenierı́a tiene tantas aplicaciones
que es una herramienta indispensable. Una integral representa un área bajo la curva sobre el eje
horizontal, acotada por un intervalo. La función a integrarse, en general deberá tener una de las
tres formas siguientes.

1. Una función simple y continua tal como un polinomio, una función exponencial o una función
trigonométrica.

2. Una función complicada y continua que es difı́cil o imposible de integrar directamente.

3. Una función tabulada en donde los valores de y se dan en un conjunto de puntos discretos,
como es el caso, a menudo, de datos experimentales.

5.1. Método analı́tico


Las fórmulas o ecuaciones de Newton-Cotes son esquemas de integración numérica donde se
reemplaza una función complicada con una función aproximada o fácil de integrar:
Z b Z b
I= ∼
f (x)dx = fn (x)dx (5.1)
a a

donde fn (x) es un polinomio de la forma:

fn (x) = a0 + a1 x + a2 x2 + · · · + an−1 xn−1 + an xn (5.2)


donde n es el orden del polinomio.
La integral se puede aproximar mediante una serie de polinomios aplicados por pedazos a la
función o datos sobre segmentos de longitud constante como se muestra en las figuras 5.1 y 5.2.
Se disponen de formas cerradas y abiertas de las ecuaciones de Newton-Cotes. Las formas
cerradas son aquellas en las que se conocen los datos al inicio y al final del intervalo de la integración
(figura 5.3). Las formas abiertas son aquellas en las cuales los lı́mites de integración se extienden
más allá del intervalo de los datos conocidos.

55
Dra. Dora-Luz Flores Métodos Numéricos

Figura 5.1: Estimación de una integral mediante una lı́nea recta.

Figura 5.2: Estimación de una integral mediante una parábola.

Figura 5.3: Forma cerrada de la estimación de una integral.

56
Dra. Dora-Luz Flores Métodos Numéricos

5.2. Método de la Regla del Trapecio


Es la primera forma o método de integración de Newton-Cotes. La integral aproximada es:
Z b
f (a) + f (b)
I= f (x)dx ∼
= (b − a) (5.3)
a 2
Geométricamente el método trapezoidal es un equivalente a aproximar gráficamente el área de
un trapezoide bajo la recta que una a f (a) y f (b), como se puede observar en la figura 5.4.

Figura 5.4: Método de la regla del trapecio.

El error aproximado está dado por:


b
(b − a)2
Z
Ea = − f 00 (x)dx (5.4)
12 a

5.2.1. Ejemplo del método de la regla del trapecio


Utilice el método de integración trapezoidal para integrar numéricamente la función

f (x) = 400x5 − 900x4 + 675x3 − 200x2 + 25x + 0.2

desde a = 0 y b = 0.8. El valor verdadero es 1.640533.

Evaluando a y b en f (x) se tiene:

a = 0, f (a) = 0.2

b = 0.8, f (b) = 0.232

Sustituyendo en los valores en la regla trapezoidal se tiene:

57
Dra. Dora-Luz Flores Métodos Numéricos

Figura 5.5: Regla del trapecio para dos segmentos.


Z b
f (a) + f (b) 0.2 + 0.232
I= f (x)dx ∼= (b − a) ∼
= (0.8 − 0) = 0.1728
a 2 2
(0.8 − 0)2 R 0.8 00
Ea = − 0
f (x)dx = 2.56
12
Z 0.8 Z 0.8
00
donde f (x)dx = (8000x3 − 10800x2 + 4050x − 400)dx = −48
0 0

El error verdadero es:



1.640533 − 0.1728
Ev = × 100 % = 89.47 %.
1.640533

5.2.2. Aplicación múltiple de la regla trapezoidal


Para mejorar la exactitud de la regla trapezoidal se divide el intervalo de integración de a a b
en un número n de segmentos y se aplica el método en cada uno de los nuevos segmentos, como
se observa en las figuras 5.5, 5.6 y 5.7.
Con lo cual se tiene que la regla trapezoidal múltiple es:
n−1
P
Z b f (x0 ) + 2 f (xi ) + f (xn )
I= f (x)dx ∼
=
i=2
(b − a) (5.5)
a 2n
De donde (b − a) es la anchura del intervalo de integración, y la división es la altura promedio
del trapecio.
b−a
Para calcular la anchura de los nuevos intervalos se tiene que h = .
n
Y el error aproximado está dado por:

(b − a)2 b 00
Z
Ea = − f (x)dx (5.6)
12n2 a

58
Dra. Dora-Luz Flores Métodos Numéricos

Figura 5.6: Regla del trapecio para tres segmentos.

Figura 5.7: Regla del trapecio para n segmentos.

59
Dra. Dora-Luz Flores Métodos Numéricos

5.2.3. Ejemplo de la aplicación múltiple de la regla del trapecio


Utilizar el método de integración trapezoidal múltiple para integrar numéricamente la función

f (x) = 400x5 − 900x4 + 675x3 − 200x2 + 25x + 0.2

desde a = 0 y b = 0.8. Utilizando dos y tres segmentos. El valor verdadero es 1.640533.

Para dos segmentos:


b−a 0.8 − 0
h= = = 0.4
n 2
x0 = 0, f (x0 ) = 0.2

x1 = 0.4, f (x1 ) = 2.456

x2 = 0.8, f (x2 ) = 0.232


n−1
P
Z b f (x0 ) + 2 f (xi ) + f (xn )
0.2 + 2(2.456) + 0.232
I= f (x)dx ∼
=
i=2
(b−a) = (0.8−0) = 1.0688
a 2n 2(2)
2
(0.8 − 0) R 0.8 00
Ea = − 0
f (x)dx = 0.64
12(2)2
Z 0.8 Z 0.8
00
donde f (x)dx = (8000x3 − 10800x2 + 4050x − 400)dx = −48
0 0

El error verdadero es:



1.640533 − 1.0688
Ev = × 100 % = 34.85 %.
1.640533

Para tres segmentos:


b−a 0.8 − 0
h= = = 0.266667
n 3
x0 = 0, f (x0 ) = 0.2

x1 = 0.266667, f (x1 ) = 1.432724

x2 = 0.533333, f (x2 ) = 3.487177

x3 = 0.8, f (x2 ) = 0.232

60
Dra. Dora-Luz Flores Métodos Numéricos

Z b
0.2 + 2(1.432724 + 3.487177) + 0.232
I= f (x)dx ∼= (0.8 − 0) = 1.369574
a 2(3)
(0.8 − 0)2 R 0.8 00
Ea = − 0
f (x)dx = 0.284444
12(3)2
Z 0.8 Z 0.8
00
donde f (x)dx = (8000x3 − 10800x2 + 4050x − 400)dx = −48
0 0

El error verdadero es:



1.640533 − 1.369574
Ev = × 100 % = 16.51 %.
1.640533

5.3. Método Simpson 1/3 y 3/8


Una manera de mejorar la exactitud del método trapezoidal es usar polinomios de mayor orden
para conectar los puntos. Por ejemplo, si existe un punto entre f (a) y f (b), a la mitad, estos puntos
se pueden conectar mediante una parábola.
Si hay dos puntos igualmente espaciados entre f (a) y f (b), los cuatro puntos se pueden conectar
mediante un polinomio de tercer orden.
A las ecuaciones que se utilizan para calcular las integrales bajo estos polinomios se conocen
como reglas de Simpson.

5.3.1. Regla de Simpson 1/3


Utilizando un polinomio de segundo orden se tiene que la aproximación del área bajo la curva
mediante tres puntos o una parábola esta dada por (ver figura 5.8):

Figura 5.8: Regla Simpson 1/3.

61
Dra. Dora-Luz Flores Métodos Numéricos

Z b
f (x0 ) + 4f (x1 ) + f (x2 )
I= f (x)dx ∼
= (b − a) (5.7)
a 6
Para calcular la anchura de los nuevos intervalos se tiene:
b−a
h= (5.8)
2
Y el error aproximado está dado por:
b
(b − a)4
Z
Ea = − f IV (x)dx (5.9)
2280 a

5.3.2. Ejemplo del método Simpson 1/3


Utilizar el método de Simpson 1/3 para integrar numéricamente la función

f (x) = 400x5 − 900x4 + 675x3 − 200x2 + 25x + 0.2

desde a = 0 a b = 0.8.

Obteniendo la anchura:
b−a 0.8 − 0
h= = = 0.4
n 2
Por lo tanto:

x0 = 0, f (x0 ) = 0.2

x1 = 0.4, f (x1 ) = 2.456

x2 = 0.8, f (x2 ) = 0.232


Z b
f (x0 ) + 4f (x1 ) + f (x2 ) 0.2 + 4(2.456) + 0.232
I= f (x)dx ∼ = (b − a) = (0.8 − 0) = 1.367467
a 6 6
El error verdadero es:

1.640533 − 1.367467
Ev = × 100 % = 16.64 %.
1.640533

62
Dra. Dora-Luz Flores Métodos Numéricos

5.3.3. Regla de Simpson 3/8


Utilizando un polinomio de tercer orden se tiene que la aproximación del área bajo la curva
mediante cuatro puntos o una ecuación cúbica esta dada por (figura 5.9):
Z b
f (x0 ) + 3f (x1 ) + 3f (x2 ) + f (x3 )
I= f (x)dx ∼
= (b − a) (5.10)
a 8

Figura 5.9: Regla Simpson 3/8.

Para calcular la anchura de los nuevos intervalos se tiene que:


b−a
h= (5.11)
3
Y el error aproximado está dado por:
b
(b − a)4
Z
Ea = − f IV (x)dx (5.12)
6480 a

5.3.4. Ejemplo del método Simpson 3/8


Utilizar el método de Simpson 3/8 para integrar numéricamente la función

f (x) = 400x5 − 900x4 + 675x3 − 200x2 + 25x + 0.2

desde a = 0 a b = 0.8.

Obteniendo la anchura:
b−a 0.8 − 0
h= = = 0.266667
n 3
Por lo tanto:

63
Dra. Dora-Luz Flores Métodos Numéricos

x0 = 0, f (x0 ) = 0.2

x1 = 0.266667, f (x1 ) = 1.432724

x2 = 0.533333, f (x2 ) = 3.487177

x3 = 0.8, f (x2 ) = 0.232


Z b
f (x0 ) + 4f (x1 ) + f (x2 )
I= f (x)dx ∼ = (b − a) =
a 6
0.2 + 3(1.432724) + 3(3.487177) + 0.232
(0.8 − 0) = 1.519170
8
El error verdadero es:

1.640533 − 1.519170
Ev = × 100 % = 21.99 %.
1.640533

5.4. Método de diferenciación


La diferenciación numérica, o aproximación numérica, es un método utilizado para evaluar
las derivadas de funciones por medio de valores funcionales de puntos de datos discretos. Si se
conocen los valores funcionales de dichos datos discretos, la función se puede expresar de una
forma aproximada por medio de una interpolación polinomial. Por lo que, al diferenciar dicho
polinomio, se pueden evaluar sus derivadas.
Por definición la derivada de una función esta dada por:

f (x + h) − f (x)
f 0 (x) = lı́m (5.13)
h→0 h

5.4.1. Diferenciación hacia adelante, hacia atrás y centrada


Se puede representar gráficamente como se muestra en las figuras 5.10, 5.11 y 5.12.
Lo que sigue es un resumen de las fórmulas de diferenciación que se pueden obtener a partir
de desarrollos en serie de Taylor.

Expresiones de primeras diferencias hacia adelante


f (x0 + h) − f (x0 )
f 0 (x0 ) =
h
f (x0 + 2h) − 2f (x0 + h) + f (x0 )
f 00 (x0 ) =
h2

64
Dra. Dora-Luz Flores Métodos Numéricos

Figura 5.10: Diferenciación hacia adelante.

Figura 5.11: Diferenciación hacia atrás.

f (x0 + 3h) − 3f (x0 + 2h) + 3f (x0 + h) − f (x0 )


f 000 (x0 ) =
2h3
f (x0 + 4h) − 4f (x0 + 3h) + 6f (x0 + 2h) − 4f (x0 + h) + f (x0 )
f IV (x0 ) =
h4

Expresiones de primeras diferencias hacia atrás


f (x0 ) − f (x0 − h)
f 0 (x0 ) =
h
f (x0 ) − 2f (x0 − h) + f (x0 − 2h)
f 00 (x0 ) =
h2

65
Dra. Dora-Luz Flores Métodos Numéricos

Figura 5.12: Diferenciación centrada.

f (x0 ) − 3f (x0 − h) + 3f (x0 − 2h) − f (x0 − 3h)


f 000 (x0 ) =
2h3
f (x0 ) − 4f (x0 − h) + 6f (x0 − 2h) − 4f (x0 − 3h) + f (x0 − 4h)
f IV (x0 ) =
h4

Expresiones de primeras diferencias centrales


f (x0 + h) − f (x0 − h)
f 0 (x0 ) =
2h
f (x0 + h) − 2f (x0 ) + f (x0 − h)
f 00 (x0 ) =
h2
f (x0 + 2h) − 2f (x0 + h) + 2f (x0 − h) − f (x0 − 2h)
f 000 (x0 ) =
2h3
f (x0 + 2h) − 4f (x0 + h) + 6f (x0 ) − 4f (x0 − h) + f (x0 − 2h)
f IV (x0 ) =
h4

5.4.2. Ejemplo de diferenciación numérica


Usar aproximaciones de diferencias finitas hacia adelante, hacia atrás y centradas para estimar
la primera derivada en x = 0.5 de:

−0.1x4 − 0.15x3 − 0.5x2 − 0.25x + 1.2

66
Unidad 6

Solución numérica de ecuaciones


diferenciales

Sea
dv c
=g− v (6.1)
dt m
una ecuación diferencial de primer orden, donde:

g, c y m son constantes,

v es una variable dependiente,

t es una variable independiente.

Las ecuaciones diferenciales se dividen en:

Ecuaciones diferenciales ordinarias (EDO), si poseen una sola variable independiente.

Ecuaciones diferenciales parciales (EDP), si poseen dos o más variables independientes.

La ecuación

d2 x dx
m 2
+ c + kx = 0 (6.2)
dt dt
es una ecuación diferencial de segundo orden, c y k son constantes. A estas ecuaciones se les conoce
como ecuación de segundo grado; en general a las ecuaciones de orden mayor a uno se les conoce
como ecuaciones de orden superior.
Las ecuaciones de orden superior se pueden reducir a un sistema de ecuaciones de primer orden
definiendo nuevas variables, si:
dx
y= (6.3)
dt

67
Dra. Dora-Luz Flores Métodos Numéricos

Entonces

dy d2 x
= 2 (6.4)
dt dt
Sustituyendo las ecuaciones anteriores se tiene que
dx
m + cy + kx (6.5)
dt
o
dy cy + kx
= (6.6)
dt m
La solución matemática de la EDO se puede obtener por separación de variables e integrando,
con lo que se obtiene:
c
dv = gdt −vdt (6.7)
m
Z h
c i
v= g − v dt (6.8)
m
Una EDO lineal es aquella que se ajusta a la forma general:

an (x)y n + an−1 (x)y n−1 + · · · + a1 (x)y + a0 (x) = f (x) (6.9)


donde: y n es la derivada n-ésima de y respecto a x y ai (x) y f (x) son funciones especı́ficas de
x.

6.1. Método de Euler


La primera derivada proporciona una estimación de la pendiente en tal como se observa en la
figura 6.1.
De la figura 6.1 se observa que

yi+1 = yi + φh (6.10)
Donde φ es la pendiente

φ = f (xi , yi ) (6.11)
Es decir, la función evaluada en los puntos (xi , yi ) rescribiendo la ecuación 6.11 se tiene que:

yi+1 = yi + f (xi , yi )h (6.12)


A esta ecuación se le conoce como método de Euler o punto medio de Euler.

68
Dra. Dora-Luz Flores Métodos Numéricos

Figura 6.1: Primera derivada para la estimación de una pendiente.

6.1.1. Ejemplo del método de Euler


dy
Utilizar el método de Euler para resolver la ecuación diferencial = −2x3 + 12x2 − 20x + 8.5,
dx
con condiciones iniciales x = 0, y = 1, hasta x = 4 con paso de iteración de 0.5.
Resolviendo la ecuación diferencial se tiene que la solución real es:

y = −0.5x4 + 4x3 − 10x2 + 8.5x + 1


Resolviendo usando Euler:
Calculando yi+1 se tiene que:

yi+1 = yi + f (xi , yi )h

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

y1 (0.5) = y0 + f (x0 , y0 )h = 1 + (8.5)(0.5) = 5.25


Se tiene que la aproximación para y1 es 5.25. Obteniendo el valor real:

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


El error verdadero relativo porcentual es:

3.21875 − 5.25
|Ev | = = 63.11 %
3.21875

69
Dra. Dora-Luz Flores Métodos Numéricos

Tabla 6.1: Valores obtenidos por el método de Euler para resolver una EDO.
i xi yi f (xi , yi ) vv |Ev |
0 0 1 8.50000 1.00000 0.00
1 0.5 5.25000 1.25000 3.21875 63.11
2 1.0 5.87500 -1.50000 3.00000 95.83
3 1.5 5.12500 -1.25000 2.21875 130.99
4 2.0 4.50000 0.50000 2.00000 125.00
5 2.5 4.75000 2.25000 2.71875 74.71
6 3.0 5.87500 2.50000 4.00000 46.88
7 3.5 7.12500 -0.25000 4.71875 50.99
8 4.0 7.00000 3.00000 133.33

Calculando y2 se tiene:

f (x1 , y1 ) = f (0.5, 5.25) = (−2)(0.5)3 + (12)(0.5)2 + (−20)(0.5) + 8.5 = 1.25

y2 (1) = y1 + f (x1 , y1 )h = 5.25 + (1.25)(0.5) = 5.875

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


El error verdadero relativo porcentual es:

3 − 5.875
|Ev | = = 95.83 %
3
Se continúa sucesivamente x = {1.5, 2, 2.5, 3, 3.5, 4}, de acuerdo al paso utilizado. Los valores
obtenidos son valores de la curva que resulta de la ecuación y se muestran en la tabla 6.1. Graficando
la solución aproximada y la verdadera se tiene lo que se observa en la grafica 6.2.

6.1.2. Análisis del error para el método de Euler


Los tipos de errores en las ecuaciones diferenciales ordinarias pueden ser de truncamiento y de
redondeo, los primeros son causados por las técnicas empleadas para aproximar, estas pueden ser
de dos tipos:

Error local: Es el resultado de aplicar el método en cuestión a un solo paso.

Error propagado: Resulta de las aproximaciones producidas en los pasos previos.

El error de redondeo es el resultado del número finito de cifras significativas que puede manejar
un procesador.
Se puede obtener información sobre la magnitud y propiedades del error de truncamiento al
derivar el método de Euler directamente de las series de Taylor:

70
Dra. Dora-Luz Flores Métodos Numéricos

Figura 6.2: Solución aproximada y verdadera de la EDO.

y 0 = f (x, y)
donde:
dy
y0 =
dx
Si la solución tiene derivadas continuas, entonces se puede representar mediante una expansión
en series de Taylor respecto a xi , yi ) de la forma:
(n)
yi00 2 y
yi+1 = yi + yi0 h
+ h + · · · + i hn + Rn (6.13)
2! n!
Donde h = xi + 1 y Rn es un término remanente definido por:

y (n+1) (ξ) n+1


Rn = h
(n + 1)!
Donde ξ es cualquier punto entre el intervalo definido entre xi y xi+1 .
Una forma alternativa se obtiene al sustituir las ecuaciones anteriores para obtener

f 0 (xi , yi ) 2 f (n−1) (xi , yi ) n


yi+1 = yi + f (xi , yi )h + h + ··· + h + Ohn+1 (6.14)
2! n!
Donde Ohn+1 especifica el error de truncamiento local, el cual es proporcional al tamaño del
paso elevado a la (n + 1)-ésima potencia.

71
Dra. Dora-Luz Flores Métodos Numéricos

Al comparar la ecuaciones se observa que el método de Euler corresponde a la serie de Taylor e


incluyendo el término f (xi , yi )h. De la comparación se observa que ocurre un error de truncamiento
debido a que se aproxima la solución verdadera mediante un número finito de términos de la serie
de Taylor.

f 0 (xi , yi ) 2
Ev = h + · · · + Ohn+1 (6.15)
2!
A lo cual se obtiene el error de truncamiento local verdadero. Si h es lo suficientemente pequeña,
los errores en la ecuación 6.15 disminuyen al elevar el orden de h y el resultado se presenta como:

f 0 (xi , yi ) 2
Ea = h
2!
o

Ea = Oh2
donde Ea es el error de truncamiento local aproximado.

6.1.3. Ejemplo del cálculo del error del método de Euler


f 0 (xi , yi ) 2
Usar la ecuación Ev = h + · · · + Ohn+1 para estimar el error del paso inicial del
2!
ejemplo se la sección 6.1.1. Utilizarla también para determinar el error debido a cada uno de los
términos de orden superior de la serie de Taylor.
Del ejemplo se tiene que:
dy
= −2x3 + 12x2 − 20x + 8.5
f (x, y) =
dx
Con condiciones iniciales x = 0, y = 1, hasta x = 4 con paso de iteración de 0.5. Para este
ejemplo se obtiene un Ev de:

f 0 (xi , yi ) 2 f 00 (xi , yi ) 3 f 000 (xi , yi ) 4


Ev = h + h + h
2! 3! 4!
Donde derivando f (x, y) se tiene:

f (x, y) = −2x3 + 12x2 − 20x + 8.5


f 0 (x, y) = −6x2 + 24x − 20
f 00 (x, y) = −12x + 24
f 000 (x, y) = −12
Tomado el punto inicial x = 0, se tiene el error para la primera derivada

f 0 (xi , yi ) 2 −6x2 + 24x − 20 (−6)(0) + (24)(0) − 20


h = (0.5)2 = (0.5)2 = −2.5
2! 2 2

72
Dra. Dora-Luz Flores Métodos Numéricos

Para la segunda derivada se tiene

f 00 (xi , yi ) 3 12x + 24 12(0) + 24


h = (0.5)3 = (0.5)3 = 0.5
3! 6 6
Para la tercera derivada se tiene

f 000 (xi , yi ) 4 −12


h = (0.5)4 = −0.03125
4! 24
El error verdadero se obtiene sumando los errores para cada término

Ev = −2.5 + 0.5 − 0.03125 = −2.03125


Recordando que se sabe el valor verdadero es 3.21875 y que la estimación obtenida con el
método de Euler es 5.25.
Valor verdadero = estimación de aproximación + error.
Valor verdadero = 5.25 − 1.03125 = 3.21875.

6.2. Método de Euler mejorado


Una fuente fundamental de error en el método de Euler es que la derivada al principio del
intervalo se supone que se aplica a través del intervalo entero. Existen dos modificaciones simples
para ayudar a evitar este inconveniente. Las modificaciones en realidad pertenecen a una clase
mayor de métodos de solución llamados métodos de Runge-Kutta. Sin embargo ya que tiene una
interpretación gráfica sencilla, se presenta antes de la derivación formal de los métodos de Runge-
Kutta. Para corregir estas deficiencias se plantean primero el método de Heun y posteriormente
los métodos de Runge-Kutta.

6.2.1. Método de Heun


Un método para mejorar la aproximación a la pendiente implica el cálculo de dos derivadas
del intervalo, en un punto inicial y otra en un punto final. En seguida se promedian las derivadas
y se obtiene una aproximación mejorada de la pendiente en el intervalo completo. Este esquema,
llamado método de Heun, se muestran en las figuras 6.3 y 6.4.
La pendiente al principio de un intervalo es

yi0 = f (xi , yi ) (6.16)


Se usa para extrapolar linealmente a yi+1 :

0
yi+1 = yi + f (xi , yi )h (6.17)
En el método estándar de Euler se pararı́a en este punto. Sin embargo en el método de Heun, la
0
yi+1 no es la respuesta final si no una predicción intermedia. Esto se debe a que se ha distinguido

73
Dra. Dora-Luz Flores Métodos Numéricos

Figura 6.3: Esquema gráfico del método de Heun. a) Predictor.

Figura 6.4: Esquema gráfico del método de Heun. b) Corrector.

0
a esta con el superı́ndice 0. La ecuación de yi+1 se llama ecuación predictora. Proporciona una
aproximación de yi+1 que permite el cálculo de una pendiente aproximada al final del intervalo:

0 0
yi+1 = f (xi+1 , yi+1 ) (6.18)
Por lo tanto, se pueden combinar las dos pendientes y obtener una pendiente promedio sobre
el intervalo:

yi0 + yi+1
0 0
f (xi , yi ) + f (xi+1 , yi+1 )
y0 = = (6.19)
2 2
Esta pendiente promedio se usa para extrapolar linealmente de yi a yi+1 usando el método de

74
Dra. Dora-Luz Flores Métodos Numéricos

Euler:
0
f (xi , yi ) + f (xi+1 , yi+1 )
yi+1 = yi + h (6.20)
2
Que se llama una ecuación correctora.
El método de Heun es un esquema predictor-corrector. Se puede expresar concisamente como:
0
Predictor yi+1 = yi + f (xi , yi )h

0
f (xi , yi ) + f (xi+1 , yi+1 )
Corrector yi+1 = yi + h
2
Nótese que debido a que la ecuación del corrector tiene yi+1 en ambos lados del signo igual, esta
puede aplicarse para “corregir” en un esquema iterativo. Esto es, se puede usar una aproximación
anterior varias veces para proporcionar una aproximación mejorada de yi+1 . Se debe entender que
este proceso no necesariamente converge a la respuesta correcta, sino converge a una aproximación
con un error de truncamiento finito.
El error aproximado esta dado por

Aproximacion actual − Aproximacion anterior
|Ea | =
× 100 %
Aproximacion actual

6.2.2. Ejemplo del método de Heun


Utilizar el método de Heun para resolver numéricamente la ecuación diferencial y 0 = 4e0.8x −0.5y
desde x = 0 a x = 4 con un tamaño de paso de 1. La condición inicial en (x, y) = (0, 2).

Donde la solución verdadera es: y = 3.076923 (e0.8x − e−0.5x ) + 2e−0.5x .

Resolviendo la pendiente f (xi , yi )

f (0, 2) = y 0 = 4e0.8x − 0.5y = 4e(0.8)(0) − (0.5)(2) = 3

0
Estimando el predictor yi+1

0
yi+1 = yi + f (xi , yi )h = 2 + (3)(1) = 5

0
Estimando la pendiente f (xi+1 , yi+1 )

f (1, 5) = 4e0.8x − 0.5y = 4e(0.8)(1) − (0.5)(5) = 6.40216

Obteniendo el corrector

75
Dra. Dora-Luz Flores Métodos Numéricos

0
f (xi , yi ) + f (xi+1 , yi+1 ) 3 + 6.40216
yi+1 = yi + h=2+ (1) = 6.70108
2 2
Obteniendo el error relativo porcentual:

6.19463 − 6.70108
|Ev | = × 100 = 8.18 %
6.19463

Este resultado proporciona un error relativo porcentual de 8.18 %, el cual comparado con el
método de Euler, es más pequeño aun cuando solo se usó un corrector. De la misma forma se
obtienen los valores mostrados en la tabla 6.2. En la figura 6.5 se observa una gráfica con el error
verdadero y la aproximación.

Tabla 6.2: Valores obtenidos por el método de Heun para resolver una EDO.
0 0
i xi yv yi f (xi , yi ) yi+1 f (xi+1 , yi+1 ) yi+1 |Ev |
0 0 2.00000 2.00000 3.00000 5.00000 6.40216 6.70108 0.00
1 1.0 6.19463 6.70108 5.55162 12.25270 13.68578 16.31978 8.18
2 2.0 14.84392 16.31978 11.65224 27.97202 30.10670 37.19925 9.94
3 3.0 33.67717 37.19925 25.49308 62.69233 66.78396 83.33777 10.46
4 4.0 75.33896 83.33777 10.62

Figura 6.5: Comparación del error verdadero con la aproximación.

6.3. Métodos de Runge-Kutta


Los métodos de Runge-Kutta logran la exactitud del procedimiento de las series de Taylor sin
requerir el uso de derivadas superiores. Existen diversas variantes, pero todas tienen la siguiente

76
Dra. Dora-Luz Flores Métodos Numéricos

forma:

yi+1 = yi + φ(xi , yi , h)h (6.21)


Donde φ(xi , yi , h) es conocida como la función incremento, la cual puede interpretarse como
una pendiente representativa en un intervalo. Esta función se escribe de forma general como:

φ = a1 k1 + a2 k2 + · · · + an kn (6.22)
Donde las a son constantes y las k son:

k1 = f (xi , yi )

k2 = f (xi + p1 h, yi + q11 k1 h)

k3 = f (xi + p2 h, yi + q21 k1 h + q22 k2 h)

···

kn = f (xi + pn−1 h, yi + qn−1,1 k1 h + qn−1,2 k2 h + · · · + qn−1,n−1 kn−1 h)


Obsérvese que las k son las relaciones de recurrencia; esto indica que k1 aparece en la ecuación
para k2 , y k2 aparece en la ecuación para k3 , etc.

6.3.1. Métodos de Runge-Kutta de segundo orden


La versión de segundo orden para yi+1 = yi + φ(xi , yi , h)h es:

yi+1 = yi + (a1 k1 + a2 k2 )h (6.23)


Donde

k1 = f (xi , yi )

k2 = f (xi + p1 h, yi + q11 k1 h)
Los valores para a1 , a2 , p1 y q11 son evaluados al igualar el término de segundo orden de
yi+1 = yi + (a1 k1 + a2 k2 )h con la expansión de la serie de Taylor. Para desarrollar esto se obtienen
tres ecuaciones con cuatro constantes desconocidas. Dichas ecuaciones son:
1 1
a1 + a2 = 1 a2 p 1 = a2 q11 = (6.24)
2 2
Debido a que se tienen cuatro incógnitas y tres ecuaciones, se propone el valor de una de estas
incógnitas para determinar las demás. Por ejemplo, si se propone un valor para a2 , se obtiene:

77
Dra. Dora-Luz Flores Métodos Numéricos

1
a1 = 1 − a2 p1 = q11 =
2a2
Debido a que se puede elegir un número infinito de valores para a2 , existen también un número
infinito de métodos o ecuaciones de Runge-Kutta de segundo orden. Las variantes más comunes
son: método de Heun de un solo corrector (donde a2 = 1/2), método del punto medio (a2 = 1) y
método de Ralston(a2 = 2/3).

6.3.2. Métodos de Runge-Kutta de tercer orden


Se puede llevar una derivación análoga a la del método de segundo orden, para n = 3. El resul-
tado de esta derivación es de seis ecuaciones con ocho incógnitas para determinar los parámetros
restantes. Una versión común que resulta es
 
1
yi+1 = yi + (k1 + 4k2 + k3 ) h (6.25)
6
Donde:

k1 = f (xi , yi )

1 1
k2 = f (xi + h, yi + k1 h)
2 2

k3 = f (xi + h, yi − k1 h + 2k2 h)

6.3.3. Métodos de Runge-Kutta de cuarto orden


Los Métodos de Runge-Kutta más populares son los de cuarto orden. Como sucede con los
métodos de segundo orden. Existe un número infinito de versiones. Se presentan las dos versiones
mas comunes de este método, la primera versión se basa en la regla de Simpson 1/3 y comúnmente
es llamada método clásico de Runge-Kutta como se describe continuación.
 
1
yi+1 = yi + (k1 + 2k2 + 2k3 + k4 ) h (6.26)
6
Donde:

k1 = f (xi , yi )

1 1
k2 = f (xi + h, yi + k1 h)
2 2
1 1
k3 = f (xi + h, yi + k2 h)
2 2

78
Dra. Dora-Luz Flores Métodos Numéricos

k4 = f (xi + h, yi + k3 h)
La segunda se basa en la regla de Simpson 3/8 y se escribe ası́
 
1
yi+1 = yi + (k1 + 3k2 + 3k3 + k4 ) h (6.27)
8
Donde:

k1 = f (xi , yi )

1 1
k2 = f (xi + h, yi + k1 h)
3 3
1 1
k3 = f (xi + h, yi + k2 h)
3 3

k4 = f (xi + h, yi + k1 h − k2 h)
Se pueden disponer de fórmulas de Runge-Kutta de orden superior, tales como el método de
Butcher, pero en general, la ganancia obtenida en exactitud por métodos de orden superior al
cuarto orden se contrapone con la complejidad y esfuerzo de cálculo.

6.3.4. Ejemplo del método de Runge-Kutta


Utilizar el método de Runge-Kutta de cuarto orden clásico para resolver numéricamente la
ecuación diferencial y 0 = 4e0.8x − 0.5y desde x = 0 a x = 4 con un tamaño de paso de 1. La
condición inicial en (x, y) = (0, 2).

Donde la solución verdadera es: y = 3.076923 (e0.8x − e−0.5x ) + 2e−0.5x .

Utilizando las ecuaciones de Runge-Kutta clásico de cuarto orden se calculan las kn

k1 = f (x0 , y0 ) = f (0, 2) = 4e(0.8)(0) − 0.5(2) = 3

1 1
k2 = f (x0 + h, y0 + k1 h) = f (0.5, 3.5) = 4e(0.8)(0.5) − 0.5(3.5) = 4.21730
2 2

1 1
k3 = f (x0 + h, y0 + k2 h) = f (0.5, 4.10865) = 4e(0.8)(0.5) − 0.5(4.10865) = 3.91297
2 2

k4 = f (x0 + h, y0 + k3 h) = f (1, 5.91297) = 4e(0.8)(1) − 0.5(5.91297) = 5.94568


Sustituyendo las kn obtenemos la siguiente aproximación

79
Dra. Dora-Luz Flores Métodos Numéricos

 
1
yi+1 = yi + (k1 + 2k2 + 2k3 + k4 ) h =
6
 
1
2 + (3 + 2(4.21730) + 2(3.91297) + 5.94568) (1) = 6.20104
6
Obteniendo le error verdadero

6.19463 − 6.20104
|Ev | =
× 100 = 0.10 %
6.19463
Y se obtienen los valores mostrados en la tabla 6.3.

Tabla 6.3: Valores obtenidos por el método de Runge-Kutta de cuarto orden para resolver una
EDO.
i xi yi Vv k1 k2 k3 k4 yi+1 |Ev |
0 0 2.00000 2.00000 3.00000 4.21730 3.91297 5.94568 6.20104 0.00
1 1.0 6.20104 6.19463 5.80165 8.72954 7.99756 12.71283 14.86248 0.10
2 2.0 14.86248 14.84392 12.38089 19.02976 17.36754 27.97769 33.72135 0.13
3 3.0 33.72135 33.67717 27.23203 42.10991 38.39044 62.07423 75.43917 0.13
4 4.0 75.43917 75.33896 0.13

80

También podría gustarte