• Author / Uploaded
• Paola

• Categories
• Análisis numérico
• Aproximación
• Física y matemáticas
• Matemáticas
• Análisis matemático

Citation preview

Universidad Católica Boliviana San Pablo Unidad Académica Tarija

DEPARTAMENTO DE CIENCIAS EXACTAS E INGENIERÍA CIVIL

Ing. Abel Barroso López Profesor de la Universidad Católica Boliviana San Pablo

Edición Preliminar Tarija, Invierno de 2010

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------CONTENIDO Numeral 1

Unidad Didáctica

Página

1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

Aproximaciones y errores Objetivo de la unidad didáctica Introducción Algoritmo ¿Qué es el análisis numérico? Objetivo del Análisis Numérico Introducción a los métodos numéricos Análisis de errores Las fórmula de Taylor y de Mc Laurin Problemas

3 3 4 4 4 4 8 11 16

2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8

Raíces de ecuaciones no lineales Introducción Método gráfico Método de la bisección Método de la regla falsa Método de Newton – Raphson Método de la secante Método de iteración del punto fijo Problemas

17 18 19 26 31 35 38 42

2

3

3.8 3.9

Sistemas de ecuaciones lineales Introducción Tipos de matrices Algunas operaciones elementales con matrices Eliminación Gaussiana Simple Eliminación Gaussiana con Pivoteo Parcial Método de Gauss – Jordán Aplicación del Método de Gauss-Jordan al cálculo de la matriz inversa Método de Gauss – Seidel Problemas

4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10

Interpolación Introducción Diferencias finitas divididas de Newton Polinomio de interpolación de Newton con diferencias divididas Polinomio de interpolación de Lagrange Interpolación de Splines Funciones de spline de grado 0 Funciones de Splines de grado 1 Funciones de Splines de grado 2 Funciones de Splines de grado 3 Problemas

3.1 3.2 3.3 3.4 3.5 3.6 3.7

4

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

44 46 47 50 54 56 59 62 68

70 75 76 78 81 83 83 87 93 95

1

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

Numeral 5

Unidad Didáctica

Página

Integración Numérica 5.1 Introducción 5.2 Teorema fundamental del cálculo integral 5.3 Fórmulas de integración de Newton – Cotes 5.4 Regla de integración del trapecio 5.5 Regla de Simpson de un tercio 5.6 Regla de Simpson de tres octavos 5.7 Integración en intervalos desiguales 5.8 Método de integración de Romberg 5.9 Algoritmo de integración de Romberg 5.10 Problemas

96 96 97 97 101 106 108 110 116 118

Ecuaciones diferenciales ordinarias Introducción Método de Euler Método de Taylor Método de Euler mejorado Método de Runge – Kutta Problemas

120 122 126 128 133 135

Bibliografía

136

6 6.1 6.2 6.3 6.4 6.5 6.6

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

2

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------UNIDAD DIDÁCTICA 1 APROXIMACIONES Y ERRORES

1.1 OBJETIVO DEL UNIDAD DIDÁCTICA Lograr que el estudiante de ingeniería, alcance una percepción clara de las estructuras generales que instituyen ésta área del conocimiento matemático y reconozca los diferentes tipos de errores y el efecto que los mismos causan en la determinación de soluciones de problemas de cálculo numérico. 1.2 INTRODUCCIÓN En razón de no existir una definición simple que detalle lo que es el análisis numérico, la mejor manera de entenderlo, parecería ser describirlo progresivamente a través de las necesidades que motivan su creación y las características que fue adquiriendo según esas necesidades. El análisis numérico aborda los mismos problemas que el análisis matemático cualitativo, pero manejando solo operaciones aritméticas. Es así como los métodos numéricos surgen de la necesidad de métodos alternativos de solución de problemas matemáticos que no son susceptibles de resolverse mediante métodos analíticos y, porque, no siempre se dispone del tiempo suficiente para investigar hasta hallar una solución analítica elegante. De este modo, en el análisis numérico se pierde la característica fundamental del análisis matemático, que es la solución, mediante lenguaje simbólica, de problemas matemáticos. Esta característica simbólica del análisis matemático – denominada también exacta – genera el poder del pensamiento que permite su generalización. Los métodos numéricos al eliminar la abstracción simbólica, pierden esos poderes, pero se ganan otros: la posibilidad de resolver, aunque sea aproximadamente, cualquier problema. En este contexto, para calcular valores numéricos que resuelven problemas científicos y tecnológicos, expresados analíticamente, mediante un modelo matemático, es necesario desarrollar una regla que permita calcular una aproximación al valor buscado. La tarea de construir y mejorar esa regla, es propia del análisis numérico, mientras que la regla en sí, es el método numérico. De este modo podemos inferir que el objeto de estudio del análisis numérico es la construcción y valoración de métodos numéricos para la solución de problemas del análisis matemático que tienen como resultado un valor numérico El desarrollo de métodos numéricos son procedimientos aproximados para resolver problemas científicos y tecnológicos que han dado lugar al nacimiento de una rama particular de las matemáticas, denominada matemáticas numéricas, o bien, Cálculo Numérico o Análisis Numérico.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

3

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Los métodos numéricos de esta rama del análisis son aproximados, ya que, en la práctica, cada cantidad se calcula con un cierto número de cifras significativas que, para las aplicaciones resulta suficiente, y no es importante conocer el valor exacto de una cantidad. 1.3 ALGORITMO Es un procedimiento que describe – de manera inequívoca – una serie finita de pasos a seguirse en un orden determinado. Su finalidad es implantar un procedimiento para resolver un problema o aproximar una posible solución. Los algoritmos son el objeto de estudio del análisis numérico. 1.4 ¿QUÉ ES EL ANÁLISIS NUMÉRICO? Por lo expuesto en 1.2, el análisis numérico y sus métodos constituyen un puente entre el análisis matemático cualitativo y el análisis matemático cuantitativo; puesto que, mientras el primero muestra que, bajo ciertas condiciones, se tiene “algo” (un modelo matemático), el segundo complementa al primero permitiendo calcular aproximadamente el valor numérico de aquello que existe. Así pues, se dice que el análisis numérico es “la teoría de los métodos constructivos en el Análisis Matemático” o también, que es “La rama de las Matemáticas que intenta obtener soluciones aproximadas, a los problemas del análisis matemático, por medio de operaciones elementales”. De manera general, el análisis numérico o cálculo numérico, puede entenderse como la rama de las matemáticas orientada a la búsqueda de algoritmos para su tratamiento aritmético. De ahí que, desde un punto de vista matemático, se lo entiende como un instrumento auxiliar que permite poner los algoritmos en práctica para llegar a resultados numéricos. 1.5 OBJETIVO DEL ANÁLISIS NUMÉRICO El propósito del Análisis Numérico es desarrollar algoritmos para hallar soluciones aproximadas a problemas matemáticos complejos, utilizando solamente operaciones simples de la aritmética. Es decir, se trata de resolver problemas difíciles mediante muchos pasos fáciles. 1.6 INTRODUCCIÓN A LOS MÉTODOS NUMÉRICOS Para comprender la forma de trabajar con métodos numéricos en la solución de problemas, resolveremos un problema común de la física elemental. Problema: Deseamos calcular la velocidad instantánea, cerca de la superficie terrestre, de un cuerpo, de masa m, que cae libremente, suponiendo que la velocidad inicial del cuerpo es igual a 0 y, que las únicas fuerzas que actúan sobre él son la fuerza de gravedad y la fuerza de resistencia del aire, la que suponemos es linealmente proporcional a la velocidad del cuerpo. Determinación del Modelo Matemático: Utilizamos la segunda ley de Newton: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

4

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------F = m a.

F=m

dv dt

Las fuerzas que actúan sobre el cuerpo indican que: F = FG + FR Donde: FG = m g

es la fuerza de gravedad (con g: aceleración de gravedad) y;

FR = – c v es la fuerza de resistencia del aire (con c: coeficiente de arrastre). dv =mg–cv dt

Por sustitución obtenemos:

m

De donde:

c dv =g– v m dt

… (1)

Que es el “modelo matemático” del problema que se trata de resolver. a) Solución Analítica: En este caso se identifica el modelo como una ecuación diferencial de primer orden de variables separables., por lo tanto, el camino se orienta a resolver la ecuación diferencial (1) Separando las variables:

dv = dt c g v m

Integrando ambos miembros de la ecuación:

dv = c g v m De lo cual obtenemos: Donde:

m ln g c

dt

c v =t+ k m

(2)

k = constante de integración

Para calcular la constante de integración, utilizamos la hipótesis dada: v = 0 si t = 0. Sustituyendo estos valores en (2) obtenemos: m ln g = k c

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

5

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------m ln g c

Finalmente tenemos:

m ln g c

c v =t m

De donde despejamos v en función de t, obteniendo:

v=

ct m

m g 1 e c

… (3)

Ecuación que resuelve el problema de manera exacta. b) Solución Numérica. Orientamos nuestra solución a buscar una expresión que aproxime el valor de v. Recordaremos la definición de la derivada de una función, teniéndose: dv ti dt

= lim ti

v ti

ti

ti

1

v ti

1

ti

1

Cuando ti+1 es un valor muy cercano a ti, se puede quitar el límite y obtener la siguiente aproximación: dv ti

v ti

dt

v ti

1

ti

1

ti

Expresión que, aplicada al modelo matemático, permite escribir: v ti 1 ti 1

v ti ti

g–

c v(ti) m

De donde se despeja v(ti+1) y se obtiene: v(ti+1) = v(ti) + [g –

c v(ti)] (ti+1 – ti) m

… (4)

Fórmula, permite calcular la velocidad aproximado de v(ti+1) si se conoce la velocidad v(ti) en el instante anterior. Ejemplo: Supongamos que se tenemos los siguientes datos: m = 70,000 g, c = 19,600 g/s, g = 980 cm/s2 Se calcularán los valores de v(0), v(1), v(2), … v(5). a) Solución analítica. Sustituimos los valores de m, c y g en la expresión (2) y, obtenemos:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

6

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------19,600t

70,000x980 v= 1 e 70,000 = 3500 (1 – e-0.28t) 19,600 Para los valores del tiempo, desde t = 0 hasta t = 5, determinamos los resultados que se ven en la siguiente tabla: 0 0

t (s) v (cm/s)

1 854.7569

2 3 4 5 1500.76828 1989.01317 2358.02072 2636.91063

Esta tabla presenta los valores exactos de las velocidades indicadas y que se han obtenido por un método analítico. b) Solución numérica: Veremos ahora, como se puede aproximar los datos dados utilizando un método numérico. Como la velocidad inicial es nula, es decir: v(0) = 0. Utilizamos la fórmula (4) y calculamos la velocidad que corresponde a los tiempos subsecuentes. Como estos cálculos son aproximaciones, entre más pequeños sean los intervalos de tiempo, más exactas serán dichas aproximaciones. Para nuestro ejemplo, tenemos: v(ti+1)

v(ti) + [980 -

19,600 v(ti)](ti+1 – ti) = v(ti) + [980 – 0.28v(ti)](ti+1 – ti) 70,000

Como expresamos anteriormente, comenzamos con v(0) = 0. Para aproximar v(1), tenemos dos opciones: aproximarla directamente, saltando del tiempo t = 0 al tiempo t = 1, o bien usar intervalos más pequeños de tiempo, como 0.2 segundos, para obtener una mejor aproximación. Para la primera opción tendremos v(1) siguientes resultados: 0 0

t (s) v (cm/s)

0.2 196

980, mientras que para la segunda tenemos los

0.4 381.024

0.6 555.686656

0.8 720.5682033

1 876.2163839

Evidentemente, con la segunda opción obtenemos una mejor aproximación para v(1), ya que los intervalos de tiempo son más pequeños, y lógicamente, si se reduce, aún más, la duración de estos intervalos de tiempo, se obtendrán mejores aproximaciones. Para obtener mejores aproximaciones, sin duda, se elegirá la segunda opción, con intervalos de tiempo de 0.2 s. Esto arroja los resultados que se muestran en la siguiente tabla: t (s) V (cm/s)

0 0

1 876.2163839

2 1533.074153

3 2025.489197

4 2394.629346

5 2671,356168

Tabla en la que se han omitido los datos intermedios para no hacer más larga la tabla. Si hacemos una comparación entre la tabla de valores exactos y esta última de valores aproximados, se verá que hay diferencias entre los valores obtenidos, es decir, en la ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

7

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------segunda tabla se han cometido ciertos errores que deben ser medidos y lo que es más importante, se hace necesario tener alguna forma de poder afirmar: "el resultado obtenido es lo suficientemente bueno". Observación: El precio que hay que pagar, en la solución aproximada, es el de efectuar cálculos cada vez más largos y a veces tediosos, y aquí es donde hace acto de presencia la poderosa herramienta computacional, que permite hacer cálculos, como los indicados, en poco tiempo y con mayor exactitud que si se los hiciera manualmente. Esta gran herramienta hace factible el camino de los métodos numéricos pues de otra forma, serían muy lentos los procesos y con muchos riesgos de cometer errores en cada paso. Es necesario destacar que, para poder elaborar un buen programa de computación, además de manejar algún lenguaje de programación, es preciso saber realizar el proceso "manualmente", ya que ello permitirá asimilar la lógica del proceso para implementar un mejor programa que contemple los posibles obstáculos que hay en el camino. 1.7 ANÁLISIS DE ERRORES El análisis del error implica el análisis de convergencia, es decir, el análisis de las características propias del método aplicado: “qué tan cerca está el valor calculado (valor aproximado) del valor buscado (valor verdadero) y, qué tan rápido se acerca el método elegido a ese valor. La pregunta que surge entonces es: ¿cómo sabemos cuándo se ha encontrado el valor verdadero? La respuesta está en el proceso del método, puesto que si el método numérico aplicado al problema proporciona, cada vez, un valor que se acerca más al valor verdadero, entonces la sucesión de valores aproximados obtenida tenderá a un número xr, y en este caso solo hay que comparar el nuevo valor aproximado xi con el valor aproximado anterior xi-1. Si la sucesión de valores converge en xr, la diferencia entre xi y xi-1, se irá haciendo más pequeña; cuando esta diferencia sea prácticamente cero (con una tolerancia prefija E), diremos que hemos encontrado una buena aproximación al valor verdadero, es decir si: xi-1 – xi-1

E

se ha alcanzado una buena aproximación

Uno de los objetivos principales del análisis numérico es la determinación de la precisión en los resultados de los cálculos. Fuentes básicas de los errores: a) Errores inherentes a la formulación del problema. b) Errores originados por el método empleado en la solución del problema Los errores del primer grupo se deben a que la formulación del problema es sólo una aproximación de la situación real, generalmente son errores despreciables. Los errores del segundo grupo tienen tres causas: a) Errores de datos o equivocaciones comunes en las operaciones, denominados errores instrumentales. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

8

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------b) Errores de truncamiento, Son los errores propios del método de cálculo y, resultan de representar aproximadamente un procedimiento matemático exacto. c) Errores de redondeo, cuando no se utilizan todas las cifras verdaderas de un número irracional. 1.7.1. Los errores de truncamiento, son los errores propios del método de cálculo y, resultan de representar aproximadamente un procedimiento matemático exacto. Por ejemplo, en la solución numérica al problema del objeto en caída libre, se usó una aproximación al proceso de derivación, la cual es una aproximación al procedimiento matemático exacto. dv v( t i 1 ) v( t i ) = ti 1 ti dt

Esto genera errores de truncamiento durante el procedimiento. 1.7.2. Los errores de redondeo, resultan de representar aproximadamente números que son exactos. Por ejemplo, aún en la "solución exacta" al problema del objeto en caída libre, los resultados impresos en la tabla de velocidades no son totalmente exactos puesto que el numero e es un número irracional y por lo tanto su extensión decimal es infinita y no periódica lo que impide escribirlo de forma completamente exacta. Usando 5 decimales, se tiene: e = 2.71828… Esto genera errores de redondeo durante los cálculos. En ambos casos se tiene que: valor verdadero = valor aproximado + error 1.7.3. Error absoluto. Se define el error absoluto a la diferencia entre el valor verdadero y el valor aproximado, es decir: Ev = valor verdadero – valor aproximado Esta definición de error tiene un pequeño defecto, como se verá en el siguiente: Ejemplo. Al medir la longitud de una varilla para construcción se obtiene el resultado aproximado de 19,999 cm mientras que al medir la longitud de un clavo, se obtiene el resultado de 9 cm Suponiendo que los valores verdaderos de la varilla y el clavo son de 20,000 cm y 10 cm respectivamente, calcular el error absoluto en ambos casos. Solución. Se tienen los siguientes resultados: Para el caso de la varilla, el error absoluto se calcula como: Ev = 20,000 – 19,999 = 1 cm Para el caso del clavo, el error absoluto se calcula como:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

9

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Ev = 10 – 9 = 1 cm En ambos casos, el error absoluto es igual, pero obviamente tiene mayor trascendencia el error en el caso del clavo que en el caso de la varilla, es decir, es necesario comparar el error absoluto contra el valor verdadero y esto da lugar a la siguiente definición. 1.7.4. Error relativo. Se define el error relativo al cociente del error absoluto sobre el valor verdadero, es decir: Er

error absoluto valor verdadero

Esto es: Er

valor verdadero valor aproximado valor verdadero

También se define el error relativo porcentual, como sigue: v=

Es decir:

v=

100

100 Er %

valor verdadero valor aproximado % valor verdadero

De hecho el error que más se usa es este último, ya que da una idea del porcentaje del error que se está cometiendo. Por ejemplo, en el caso de la varilla el error relativo porcentual es: v=

100

1 % = 0.005 % 20,000

Mientras que en el caso del clavo, el error relativo porcentual es: v=

100

1 % = 10 % 10

Se puede observar, que el error relativo porcentual refleja mejor la gravedad o no gravedad del error que se está cometiendo. Es claro, que en el caso de la varilla no es trascendente ya que representa solamente un 0.005% con respecto al valor verdadero, mientras que en el caso del clavo, el error si es importante ya que es del 10% del valor verdadero.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

10

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Finalmente, se mencionará que un proceso de aproximación puede detenerse cuando el valor absoluto del error relativo porcentual es menor que una cierta cota, fijada de antemano. Sin embargo, todavía se tiene un pequeño defecto en el análisis del error. Los métodos numéricos se aplican en realidad, a problemas que no se pueden resolver analíticamente; en el ejemplo del cuerpo en caída libre, en realidad no es necesario aplicar ningún método numérico, puesto que se conoce la solución exacta del problema. Por lo tanto, en una situación real, se desconocerá el valor verdadero de la solución al problema; luego entonces se estará imposibilitado de calcular el error relativo porcentual. La forma de resolver este problema es pensar que para obtener una cierta aproximación a un valor, se tuvo que haber obtenido una aproximación anterior al mismo valor. Una vez calculada la nueva aproximación se procede a calcular otra aproximación al mismo valor y así sucesivamente. Si el método realmente converge a un resultado (que se espera sea la solución del problema), todas estas aproximaciones se estarán aproximando entre sí y al valor al cual convergen. 1.7.5. Error aproximado porcentual. Se define el error aproximado porcentual, como sigue: a=

100

aproximación actual aproximación aproximación actual

previa

%

Como se mencionó anteriormente, el proceso se detiene cuando se ha logrado disminuir el valor absoluto del error aproximado porcentual hasta un cierto rango fijado de antemano; esto es, cuando: a

Se puede probar que si se toma:

< s

s

= 0.5 x 10 2-n (%),

Entonces se puede tener la seguridad de que la aproximación tiene, al menos, n cifras significativas, es decir, posee al menos n dígitos confiables. 1.8 LAS FÓRMULAS DE TAYLOR Y DE MC LAURIN En la solución de problemas de la ingeniería, utilizando los métodos numéricos, como ya se vio anteriormente, en primer lugar, es menester desarrollar un modelo matemático. Algunos de esos modelos matemáticos están determinados por funciones, cuya expresión algebraica, presenta algún grado de complejidad que obliga a buscar otro tipo de expresión más simple que, en un cierto intervalo, permita una aproximación aceptable con el modelo

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

11

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------matemático desarrollado. Las fórmulas de Taylor y de Mc Laurin son una ayuda muy práctica para este fin. Si f (x) es una función continua definida en una vecindad (entorno) de un punto a, con derivadas de todos los órdenes; es posible elegir un polinomio de grado n, pn(x), tal que en el punto x = a tenga el mismo valor que f (x) y las mismas n derivadas que f (x) en el mismo punto, este polinomio se aproximará más a f (x) en los puntos de abscisa x interiores a la vecindad del punto a, del dominio de f. Para estas condiciones, recordamos que la fórmula de Taylor muestra la siguiente forma: f (x) ≈ f (a) + (x – a)f ' (a) +

( x a ) n (n) ( x a) 2 f ”(a) + ... + f (a) n! 2!

Si f (x) es un polinomio algebraico de grado n, entonces la aproximación anterior se convierte en una verdadera igualdad. En general, si se tiene cualquier función, para aplicar la aproximación anterior y lograr que se torne en una verdadera igualdad, debemos añadir al segundo miembro un término más, llamado resto: f (x) = f (a)+ (x – a)f ' (a) +

( x a) 2 ( x a) n f“(a) + ... + f 2! n!

(n)

(a)+

( x a ) n 1 (n+1) f (c) ( n 1)!

Donde c es un punto, convenientemente elegido, al interior de la vecindad del punto a. La fórmula de Mc Laurin, es un caso particular de la fórmula de Taylor, en la que, en lugar de trabajar en la vecindad de un punto cualquiera a, se trabaja en la vecindad del origen de coordenadas, es decir, se hace a = 0: f (x) = f (0) + x f '(0) +

xn 1 x2 x n (n) f“(0) + ... + f (0) + f (n+1)(c) ( n 1)! 2! n!

Ejemplo 1. Utilicemos la serie de Mc Laurin para aproximar el número irracional e hasta un valor que tenga 4 cifras significativas confiables. Solución. Primero determinamos el valor de la cota del error s utilizando la expresión: s

= 5 x 102-n (%) = 5 x 102-4 = 0.05 %.

Ahora, implementamos la lógica de nuestra resolución: nos servimos de la serie de Mc Laurin, para desarrollar la aproximación de e y, luego iremos por pasos, empezando por el primer término de esta aproximación y, luego, agregamos un término en cada aproximación, para obtener nuevas valores hasta lograr que se cumpla: a < 0.005%. Partimos de: f(x) = ex y con ello obtenemos: f ’(x) = f ”(x) = f’’’(x) = … = f (n)(x) = ex,

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

12

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

ex

Lluego:

e0

x2 0 xn 0 e  e 2! n!

xe0

xn 1 c e n 1!

12 1n   Como e = e , hacemos x = 1 y, como e = 1, entonces resulta: e 1 1 2! n! 1 Y, en forma concentradas: e n 0 n! 1

0

Aplicando la lógica anterior, en un primer paso simplemente tendremos un término: 1 e =1 0! En el segundo paso, tenemos la suma de dos términos: e

1 0!

1 =2 1!

Aquí, calculamos el primer error de la aproximación: a

= 100

actual

previa 2 1 = 100 = 50 % actual 2

Puesto que no se ha cumplido el objetivo, seguimos agregando un nuevo término a la serie: e

1 0!

1 1!

1 = 2.5 2!

Y nuevamente calculamos el error aproximado correspondiente, es decir: a=

100

2.5 2 % = 20 % 2.5

El proceso continua hasta lograr la meta. Los resultados se resumen en la siguiente tabla: Nº términos 1 2 3 4 5 6 7 8 9

Aprox. al valor e 1 2 2.5 2.666666667 2.708333333 2.716666667 2.718055556 2.718253968 2.718278770

De este modo, el resultado que obtenemos es:

Error aproximado 50% 20% 6.25% 1.54% 0.307% 0.051% 0.007% 0.0009% e

2.718281526

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

13

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Que en realidad tiene 8 cifras significativas. Sin embargo, la cota definida por s , asegura que, al menos, se tendrá n = 4 cifras significativas sin error. Ejemplo 2: En la vecindad de a = 0, analicemos gráficamente las diferencias entre la función sen(x) y su desarrollo con la fórmula de Mc Laurin con dos y con cinco términos. Solución. El análisis lo realizaremos en la vecindad de 0. Sen (0) = 0 d sen 0 = cos (0) = 1 dx d2 sen 0 = – sen (0) = 0 dx 2 d3 sen 0 = – cos (0) = – 1 dx3 d4 sen 0 = sen (0) = 0 dx 4 d5 sen 0 = cos (0) = 1 dx5 d6 sen 0 = – sen (0) = 0 dx 6 d7 sen 0 = – cos (0) = – 1 dx 7 d8 sen 0 = sen (0) = 0 dx8 d9 sen 0 = cos (0) = 1 dx9 Para el primer caso, desarrollamos la formula de Mc Laurin: x2 d 2 x3 d 3 d sen 0 sen 0 sen (x) sen(0) + x + sen 0 + 3 ! dx3 2! dx 2 dx O bién: x2 x3 x3 sen (x) 0 + x (1) + (0) + (– 1) = x – 3! 3! 2! Graficando esta expresión y la de la función sen (x), observamos que la fórmula de Mc Laurin con dos términos, es congruente con la función seno en el intervalo [–1; 1]. Lo que nos permite concluir que: x3 sen (x) = x – x [–1; 1] 3! La figura 1.1 muestra claramente el análisis de esta aproximación.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

14

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

Figura 1.1 Pasemos ahora al análisis del segundo caso. Desarrollamos la formula de Mc Laurin: x2 d 2 x3 d 3 x4 d 4 d sen 0 + sen 0 + sen 0 + sen(x) sen(0)+x sen 0 + 3 ! dx3 4 ! dx 4 2! dx 2 dx x5 d 3 x6 d 6 x7 d 7 x8 d 8 x9 d 9 sen 0 sen 0 sen 0 sen 0 sen 0 + + + + + 5 ! dx3 6 ! dx 6 7 ! dx 7 8 ! dx8 9 ! dx9 x3 x4 x5 x6 x7 x8 x9 x2 sen(x) 0+x(1)+ (0)+ (–1)+ (0)+ (1)+ (0)+ (–1)+ (0)+ (1) 3! 4! 5! 6! 7! 8! 9! 2! Quedando finalmente: x3 x5 x7 x9 sen (x) x – + – + 3! 5! 7! 9!

Figura 1.2

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

15

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------La figura 1.2 muestra esta expresión y la de la función sen (x), donde observamos que la fórmula de Mc Laurin con cinco términos, es congruente con la función seno en el intervalo [–3; 3]. Lo que nos permite concluir que: x3 x5 x7 x9 sen (x) = x – + – + x [–3; 3] 3! 5 ! 7 ! 9 ! La figura siguiente muestra claramente el análisis de esta aproximación. 1.9. PROBLEMAS Calcula el error absoluto y el relativo en las aproximaciones de p por p*. a) p = ,

p* = 22/7

b)

p= ,

p* = 3.1416

c) p = e,

p*=2.718

d)

p=

e) p = e10,

p*=22000

f)

p=10 ,

p*=1400

g) p = 8!,

p*=39900

h)

p = 9!,

p*= 18 ( 9 / e )9

2,

p*=1. 414

i) Determina el valor aproximado de la función sen( /4) hasta la cuarta cifra significativa. j) Determina el valor aproximado de la función cos( /6) hasta la cuarta cifra significativa. k) Desarrolla en serie de Mc Laurin, para x = , las siguientes funciones: f (x) = sen (x) f (x) = cos (x) l) Desarrolla en serie de Mc Laurin las siguientes funciones: f (x) = log(1 + x) f (x) = ln (1 + x)

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

16

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------UNIDAD DIDÁCTICA 2 RAÍCES DE ECUACIONES NO LINEALES 2.1. INTRODUCCIÓN En la presente unidad didáctica estudiaremos uno de los problemas básicos de la aproximación numérica: el problema de la búsqueda de raíces. El cero de una función f (x) es un número xo para el que se cumple con f (x0) = 0. Es decir, la búsqueda de raíces consiste en determinar ese número x0 que satisface a la ecuación f (x) = 0. Muchas veces la solución de ecuaciones no puede obtenerse en forma analítica, por ello se hace necesario determinarlas en forma numérica, las que constituyen una aproximación al valor verdadero. En esta unidad didáctica, solo se considerarán raíces reales. Geométricamente, sabemos que una función y = f(x) es una curva en el plano E2, de modo que la raíz de dicha función representa el punto donde la gráfica de f (x) corta al eje x, tal como lo podemos apreciar en la figura 2.1.

Figura 2.1 Analicemos algunos ejemplos de ecuaciones y sus raíces: 1. Las raíces de f (x) = x2 – 9 son x =

3

2. La función f (x) = x4 + x2 + 1 no tiene raíces reales. 3. La función f (x) = 5 – sen x no tiene raíces. 4. Las raíces de f (x) = (x + 1)(x – 3)(x + 7) son x = – 1, x = 3 y x = – 7 El la presente unidad didáctica estudiaremos algunos métodos numéricos para aproximar las raíces de una ecuación de variable real.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

17

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------2.2 MÉTODO GRÁFICO Este método básicamente se usa para determinar un intervalo en el cual una función dada tiene alguna raíz. Ejemplo 1. Localicemos un intervalo donde la función: f (x) = e-x – ln x tenga una raíz. Solución Para calcular la raíz de f (x) hacemos f (x) = 0, Es decir:

e-x – ln x = 0

De donde:

e-x = ln x

Esta expresión implica que, el problema equivale a encontrar el punto de intersección de dos funciones: g (x) = e-x y h(x) = ln x. La figura 2.2 muestra el comportamiento de estas funciones.

Figura 2.2 En la figura 2.2 observamos que la única raíz de la función dada, corresponde al punto de intersección de los gráficos de las funciones ln(x) y e-x, se halla en el intervalo [1; 1.5]. En realidad, no interesa ser más precisos en la búsqueda del intervalo, ya que posteriormente aplicaremos métodos más sistemáticos para aproximar la raíz. Podemos afirmar que el interés del método gráfico radica en proporcionar un intervalo en el cual se comienza a trabajar con un método numérico. Ejemplo 2. Localicemos un intervalo donde la función f (x) = arc tg (x) + x – 1 tenga una raíz. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

18

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Solución. Nuevamente, para calcular la raíz de f (x) = 0. Hacemos:

arc tg x + x – 1 = 0

De donde obtenemos:

arc tg x = 1 – x.

Ahora, el problema equivale a encontrar el punto de intersección de las gráficas de las funciones: g(x) = arc tg x

h(x) = 1 – x.

y

Graficando ambas funciones podemos observar que la raíz de f se encuentra en intervalo [0; 1] donde tenemos el punto de intersección de las gráficas de ambas funciones, tal como se muestra en la figura 2.3.

Figura 2.3 2.3 MÉTODO DE LA BISECCIÓN El método de la bisección se basa en el Teorema del Valor Intermedio (TVI) Teorema del Valor Intermedio Si f (x) es una función continua en un intervalo [a, b] y las ordenadas f (a) y f (b) son valores diferentes, entonces para cada valor real k, comprendido entre los valores f (a) y f (b), existe, al menos, una abscisa x0 (a, b) tal que f (x0) = k. Esta conclusión se vale para cualquiera de los casos: f (a) > f (b)

o

f (a) < f (b)

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

19

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Básicamente el TVI expresa que toda función continua en un intervalo cerrado, una vez que alcanzó ciertos valores en los extremos del intervalo, entonces puede alcanzar todos los valores intermedios. El Teorema de Bolzano es un caso particular del TVI y expresa: Si f (x) es una función continua en un intervalo [a, b] y las ordenadas f (a) y f (b) tienen signos opuestos, entonces debe existir, al memos, un valor x0 (a, b) tal que f (x0) = 0, es decir, debe haber por lo menos una raíz de f (x) en el intervalo abierto (a, b). Pasemos ahora a explicar la lógica del método de la bisección para una cota del error de aproximación dada s: i) En un primer paso, siendo f (x) una función continua, elegimos – por el método gráfico explicado en el numeral anterior – el intervalo [a, b], de modo que f (a) y f (b) tomen signos opuestos es decir: f (xa) f (xb) < 0 ii) Seguidamente, en una primera aproximación, para hallar la raíz, ésta se asume igual al a b punto medio entre a y b, es decir: xr = 2 iii) Se calcula el valor de la ordenada f (xr), con lo que, forzosamente se caerá en alguno de los siguientes casos: a)

f (a)

f (xr) < 0

En este caso, se tiene que f (xa) y f (xr) tienen signos opuestos, y por lo tanto la raíz se encuentra en el intervalo (a, xr). b)

f (a)

f (xr) > 0

En este caso, se tiene que f (a) y f (xr) tienen el mismo signo, resultando por lo tanto que f (xr) y f (b) tienen signos opuestos; por lo que, la raíz se encuentra en el intervalo (xr, b). c)

f (xa)

f (xr) = 0

En este caso se tiene que f (xr) = 0 y por lo tanto ya se determinó la raíz. Si se tiene, ya sea el caso a) o b), el proceso se vuelve a repetir con el nuevo intervalo, calculando la nueva aproximación y verificando que:

a

= 100

xactual

x previa

xactual

%
0 Calculamos la primera aproximación a la raíz: xr1 = Calculamos:

1 1.5 = 1.25 2

f (1.25) = e-1.25 + 4 1.252– 5 = 1.536505 > 0

Para identificar mejor en que nuevo intervalo se encuentra la raíz, hacemos la siguiente tabla: f (1) –

f (1.25) +

f (1.5) +

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

24

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Por lo tanto, la raíz se encuentra en el intervalo [1; 1.25]. Calculamos el punto medio (que es una segunda aproximación a la raíz): xr2 =

1 1.25 = 1.125 2

Ahora si podemos calcular el primer error de la aproximación, puesto que se cuenta ya con la aproximación actual y la aproximación previa:

a

xr2

= 100

xr1 xr2

% = 1.11 %

Puesto que no hemos logrado el objetivo, continuamos con el proceso. f (1.125) = e-1.125 +4G1.1252 – 5 = 0.387152 > 0

Evaluamos: Construimos la tabla:

f (1) –

f (1.125) f (1.25) + +

Observamos que la raíz se encuentra en el intervalo [1; 1.125] Calculamos el punto medio:

xr3 =

1 1.125 = 1.0625 2

Y calculamos el nuevo error aproximado:

a

= 100

xr3

xr2 xr3

% = 5,88 %

Seguimos hasta cumplir la octava aproximación pedida, los valores así obtenidos son: N de Iteración 1 2 3 4 5 6 7

xa 1.00 1.00 1.00 1.0625 1.0625 1.078125 1.078125

xb 1.5 1.25 1.125 1.125 1.09375 1.09375 1,0859375

xri 1.25 1.125 1.0625 1,09375 1,078125 1,0859375 1,08203125

De lo cual, se ve que la aproximación buscada es x r 7

f(xri) 1.536505 0.387152 –0.138784 0.120114 –0.010353 0.05463 0.022071

Error % 11.11 5,88 2.86 1.45 0.72 0.36

0.022071

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

25

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Observación: El método de bisección por lo general es lento, y en los casos en los que la función tiene su raíz muy próximo a uno de los extremos del intervalo de análisis, como el de la siguiente gráfica, puede ser exageradamente lento.

Figura 2.5 En tal caso, el proceso de bisección comienza a acercarse a la raíz de forma muy lenta, ya que el método solamente toma en cuenta que la raíz se encuentra en el centro del intervalo, sin importar si se halla más cerca de alguno de los extremos del intervalo. 2.4 MÉTODO DE LA REGLA FALSA En este método, se considera si la raíz de la ecuación se localiza más cerca de alguno de los extremos del intervalo. Por tanto, la lógica del Método de la Regla Falsa aprovecha la forma de la función y, para explicarlo podemos utilizar la figura 2.6.

Figura 2.6 Tracemos, en esta gráfica, la línea recta que una los puntos definidos por las coordenadas (a, f(a) y (b, f(b)). Es claro que si en lugar de considerar el punto medio del intervalo, elegimos el punto x1 donde esta recta intercepta al eje x, se determinará mucho más rápidez la aproximación de la raíz. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

26

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Ésta es en sí, la idea central del método de la regla falsa y ésta es realmente la única diferencia con el método de bisección, puesto que en todo lo demás, ambos dos métodos son similares. Si f (x) una función continua en el intervalo [a, b] tal que f (a) y f (b) tienen signos opuestos; determinamos, primero, la ecuación de la línea recta que une los puntos (a, f (a)) y (b, f (b)), cuya pendiente está definida por la expresión: m=

f b

b a

y – f (a) =

Y, la ecuación de esta recta es:

f a

f b

f a

(x – a)

b a

Para obtener su intersección con el eje x, se hace y = 0, resultando: – f (a) =

f b

f a

b a

(x – a)

De donde despejamos x, obteniendo: x=a–

f a b a f b

f a

Que es el valor que es una aproximación a xr Sobre la base de este análisis, podemos aproximar la raíz de f(x) en el intervalo [a, b], siguiendo los siguientes pasos: i. ii.

Se determinan valores iniciales de a y b de modo que f (a) y f (b) tienen signos opuestos, es decir: f (a).f (b) < 0 La primera aproximación a la raíz se calcula con la expresión: x1 = a –

iii.

f a b a f b

f a

Se calcula la ordenada f (x1) y, forzosamente este valor deberá caer en uno de los siguientes casos: a) Si f (a).f (x1) < 0, entonces, tenemos que f (a) y f (xr) tienen signos opuestos, y por lo tanto la raíz de f se encuentra en el intervalo [a, x1]. b) Si f (a).f (x1) > 0, en este caso, f (a) y f (x1) tienen el mismo signo, y por tanto f (x1) y f(b) tienen signos opuestos; por lo que la raíz de f se encuentra en el intervalo [x1, b]. c) Si f (a).f (x1) = 0, entonces se tiene que f (x1) = 0 y, por lo tanto la raíz es x1.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

27

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Par cualquiera de las situaciones a) o b), l proceso se vuelve a repetir en el nuevo intervalo, verificando luego que se cumpla: a < s . De no cumplirse con este objetivo, deberá repetirse este proceso hasta que dicho objetivo se cumpla. Ejemplo 1. Utilicemos el método de la regla falsa para aproximar la raíz de la función definida por f (x) = e-x – ln x, hasta que a < 1 %. Comencemos en el intervalo [1, 2]. Solución. Este es el mismo ejemplo 1 del método de la bisección, de modo que ya sabemos que f (x) es continua en el intervalo dado y que toma signos opuestos en los extremos de dicho intervalo. Por lo tanto podemos aplicar el método de la regla falsa y procedemos a calcular la primera aproximación: Calculemos: f (1) = e-1 – ln 1 = e-1 = 0.367879441 y f (2) = e-2 – ln 2 = – 0.557811897 x1 = b –

f b a b f a

f b

=2–

f 2 1 2 = 1.397410482 f 1 f 2

Como solamente tenemos la primera aproximación, seguimos con el proceso para calcular la segunda aproximación. Evaluamos: f (1.397410482) = e-1.397410482 – ln 1.397410482 = – 0.087384509 < 0 Hacemos la tabla de signos: f (1)

f (1.397410482)

f (2)

+

–

–

Observamos que la raíz se encuentra en el intervalo [1, 1.397410482]. Para este nuevo intervalo, calculamos la nueva aproximación: x2

f 1.397410482 1 1.397410482 f 1 f 1.397410482

1.397410482

1.321130513 3

Ahora, podemos calcular el primer error de la aproximación: a = 100

1.321130513 1.397410482 1.321130513

% = 5.77%

Por no se cumplirse el objetivo proseguimos con el proceso. Evaluamos: f (x2) = f (1.321130513) = –0.011654346 Luego: f (1)

f (1.397410482)

f (1.321130513)

+

–

–

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

28

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------De donde se ve que la raíz se encuentra en el intervalo [1, 1.321130513], en el que debemos calcular la nueva aproximación: x3= 1.321130513 –

Y el error aproximado:

f 1.321130513 1 1.321130513 = 1.311269556 f 1 f 1.321130513 a = 100

1.311269556 1.321130513 1.311269553

% = 0.075%

Como se ha alcanzado el objetivo, concluimos con el proceso de aproximación obteniendo: xr

x3 = 1.311269556

Es interesante observar la rapidez con la cual se converge a la raíz con el método de la regla falsa, a diferencia de la lentitud que presenta el método de la bisección. Ejemplo 2. Con el método de la regla falsa aproximemos la raíz de f (x) = arc tg x + x – 1 hasta que a < 1 %, comenzando en el intervalo [0, 1]. Solución. Este es el mismo ejemplo 2 resuelto por el método de la bisección, de modo que ya sabemos que se cumplen las condiciones necesarias para poder aplicar este método, es decir, que f (x) sea continua en el intervalo dado y que asume signos opuestos en los extremos de dicho intervalo. Calculemos la primera aproximación: f (0) = arc tg 0 + 0 – 1 = – 1 f (1) = arc tg 1 + 1 – 1 = 0.785398163 x1 = xb –

f xb x a xb f 1 0 1 =1– = 0.5600991535 f xa f xb f 0 f 1

Evaluamos: f ( xr1 ) = arc tg 0.5600991535 + 0.5600991535 – 1 = 0.070662953 > 0 Y hacemos la tabla de signos: f (0)

f (05600991535)

f (1)

–

+

+

De lo cual deducimos que la raíz se localiza en el intervalo [0, 0.5600991535]. Reiniciamos el proceso calculando la nueva aproximación: x2 = 0.5600991535 –

f 0.56009915 35 0 0.56009915 35 = 0.5231330281 f 0 f 0.56009915 35

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

29

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Y determinamos el error aproximado: a = 100

0.52313302 81 0.5600991535 0.52313302 81

% = 7.06%

Como no se cumple el objetivo, se prosigue con el proceso: Evaluando: f (x2) = arc tg 0.5231330281 + 0.5231330281 – 1 = 0.00511533 Construimos la tabla de signos: f (0)

F (0.5231330281)

f (05600991535)

–

+

+

De los cual tenemos que la raíz se localiza en el intervalo [0, 0.5231330281], luego podemos calcular la siguiente aproximación: x3 = 0.5231330281 –

f 0.52313302 81 0 0.52313302 81 =0.5204706484 f 0 f 0.52313302 81

Y el siguiente error aproximado: a = 100

0.52047064 84 0.523133081 0.52047064 84

% = 0.51 %

Como se ha cumplido el objetivo, concluimos que la aproximación buscada es: xr

x3 = 0.5204706484

Nuevamente observamos el contraste entre la rapidez del método de la regla falsa contra la lentitud del método de la bisección. No debemos descartar la posibilidad de que el método de la regla falsa encuentre la aproximación a la raíz de forma más lenta que el método de la bisección. Ejercicios para que resuelva el estudiante en aula. 1) La función f(x) = e–x + 4x2 – 5, tiene una raíz en el intervalo [1; 1.5], usa ocho iteraciones del método de la regla falsa para aproximar la raíz de f(x). Tabula el error después de cada iteración. 2) Dada la función f (x) = x6 – 1, determina su raíz en el intervalo [0, 1.5], por el método de bisección hasta que a < 1 %, y compara el número de aproximaciones que requerirás con el método de la regla falsa.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

30

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------2.5 MÉTODO DE NEWTON-RAPHSON El Método de Newton - Raphson es un método iterativo que, por su rapidez, es uno de los más usados y efectivos. A diferencia de los métodos anteriores, el método de NewtonRaphson no trabaja sobre un intervalo sino que basa su fórmula en un proceso iterativo. Supongamos que tenemos la aproximación x1 a la raíz xr de f (x), como se muestra en la figura 2.7. Para determinar la primera aproximación x1 a la raíz de f (x), por un punto (x0, f (x0)) asumido sobre la gráfica de f, se traza la recta tangente a dicha gráfica, esta tangente corta al eje x en un punto de abscisa x1 que, en este método, será la aproximación a la raíz xr. Para calcular esta abscisa x1, se determina primero la ecuación de la recta tangente, cuya pendiente es: m = f ´(x0). Luego determinamos la ecuación de la recta tangente que resulta: y – f (x0) = f ´(x0)(x – x0)

Figura 2.7 Haciendo y = 0 Despejando x1:

x = x1, luego, se tiene: x1 = x0 –

– f (x0) = f ´(x0)(x1 – x0). f x0 f ' x0

El siguiente paso usa el nuevo valor x1 para determinar de forma similar la aproximación x2 y así sucesivamente, hasta calcular el valor para el cual se cumple el objetivo a < s

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

31

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------xi+1 = xi –

f xi f ' xi

Expresión que es la fórmula iterativa del método de Newton-Raphson para calcular la siguiente aproximación xi+1 a la raíz de f. Observaciones: Notamos que el método de Newton-Raphson no trabaja en un intervalo en el que se encuentre la raíz de la ecuación; sin embargo, hay que conocer ese intervalo, para poder elegir un valor inicial que permita agilizar el proceso de cálculo propuesto. Existen ejemplos en los que este método no converge a la raíz, en cuyo caso se dice que el método diverge. Sin embargo, en los casos donde si converge a la raíz, lo hace con una rapidez impresionante, por lo cual es uno de los métodos preferidos por excelencia. También podemos observar que en el caso de que f ´(xi) = 0, el método no se puede aplicar. En tal caso, geométricamente esto significa que la recta tangente es horizontal y por lo tanto no corta al eje x en ningún punto, a menos que la tangente coincida con éste eje, en cuyo caso xi es la misma raíz de f (x). Ejemplo 1 Utilicemos el método de Newton – Raphson, para aproximar la raíz de la función f (x) = e-x – ln x, hasta que a < 1 %, comenzando con xo = 1. Solución: La derivada de f (x) es: f´(x) = – e-x –

1 x

La fórmula del método resulta:

xi 1

xi

e xi

ln xi = xi 1 xi e xi

e xi

ln xi 1 e xi xi

Comenzando el proceso iterativo con xo = 1, y obtenemos: x1 = 1

El error aproximado es:

e 1 ln 1 = 1.268941421 1 e 1 1 a = 100

1.268941421 1 % = 21.19% 1.268941421

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

32

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Continuamos con el proceso hasta lograr el objetivo de la aproximación de Resumimos los resultados en la siguiente tabla: i 1 2 3

xi 1.000 1.268941421 1.309108403

f (xi) 0.367879 0.042947 0.000715

f ' (xi) –0632121 –0.50693 –0.493818

a < 1 %.

xi+1 1.268941421 1.309108403 1.309799389

De lo que concluimos que la aproximación obtenida es: xr

a

21.19 3.06 0.052

x3 = 1.309799389.

Ejemplo 2 Utilizando el método de Newton-Raphson para aproximar la raíz de la función: f (x) = arc tg x + x – 1, hasta que a < 1 %, comenzando con x0 = 0. Solución: Calculemos la derivada: f ´(x) =

1 1 x2

+1

Por sustitución en la fórmula de Newton-Raphson obtenemos: xi+1 = xi –

arctg xi 1

xi 1

1 xi2 Comenzamos con x0 = 0 y obtenemos:

x1 = 0 –

1

arctg 0 1 1 02

En este caso tenemos un error aproximado de

a = 100

0 1

= 0.5

1

0 .5 0 % = 100 % 0.5

Continuamos con el proceso hasta lograr el objetivo. Resumimos los resultados en la siguiente tabla: i 1 2 3

xi 0 0.5 0.5201957728

f (xi) –1 –0.036352 –0.00013

Concluimos con la aproximación obtenida de: xr

f ' (xi) 2 1.8 1.787027

xi+1 0.5 0.5201957728 0.5202689918

a

100 3.88 0.01

x3 = 0.5202689918

Ejemplo 3. Usemos el método de Newton-Raphson para aproximar las raíces cuadradas de números reales positivos. Solución Si R un número real positivo. Queremos calcular el número real x de tal manera que x = R . ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

33

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Elevando al cuadrado la anterior expresión, obtenemos: x2 = R, o bien: x2 – R = 0 Expresión que define la función: f (x) = x2 – R. Derivándola: f´(x) = 2x. Al sustituir estas expresiones en la fórmula de Newton-Raphson tenemos: x2 R xi+1 = xi – i 2 xi

De donde resulta:

xi+1 =

1 xi 2

R xi

Fórmula que ya era conocida por los antiguos griegos (Herón). Para tomar un ejemplo de su uso, asumiremos: R =26. Apliquemos la fórmula anterior, comenzando con xo = 5 Resumiendo los resultados del proceso en la siguiente tabla: i 1 2 3

xi 5 5.1 5.099019608

Concluimos que

26

f (xi) –1 0.01 0.000005

f ' (xi) 10 10.2 10.09804

xi+1 5.1 5.099019608 5.099019514

a

1.96 0.019 0.0000018

5.099019514, valor que es correcto en todos sus dígitos.

La misma idea puede aplicarse para crear algoritmos que aproximen raíces enésimas de números reales positivos. Ejemplo 4. La función f(x) = e–x + 4x2 – 5, tiene una raíz próxima al punto x = 1, usa ocho iteraciones del método de Newton - Raphson para aproximar la raíz de f(x). Tabula el error después de cada iteración. Resuélvelo por tu cuenta. Observación: Cuando el método de Newton-Raphson converge a la raíz, lo hace de una forma muy rápida y de hecho, se observa que el error aproximado disminuye rápidamente en cada paso del proceso. Aunque no es el objetivo establecer formalmente las cotas para los errores en cada uno de los métodos que se han estudiado, debemos mencionar que evidentemente, existen estas cotas que miden con mayor precisión la rapidez ó lentitud del método en estudio. 2.6 MÉTODO DE LA SECANTE ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

34

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Este método se basa en la fórmula de Newton-Raphson, pero evita el cálculo de la derivada usando en su lugar la siguiente aproximación:

f xi

f ´(xi)

f xi

xi

xi

1

1

(Recordemos la fórmula de la derivada usada en la solución numérica al problema del cuerpo en caída libre). Sustituimos esta expresión en la fórmula de Newton-Raphson y obtenemos: xi+1 = xi –

f xi f ' xi

xi –

f xi f xi

f xi

xi

Luego, tenemos:

xi+1

xi –

f xi f xi

xi

xi

1

f xi

1

xi

1

1

Que es la fórmula del método de la secante. Notamos que para poder calcular el valor de xi+1, es necesario conocer los dos valores anteriores xi y xi-1. La lógica de este método, gráficamente se explica en la figura 2.8.

Figura 2.8 Observación: Existe un gran parecido con la fórmula del método de la regla falsa. La diferencia entre una y otra es que mientras el método de la regla falsa trabaja sobre intervalos cerrados, el método de la secante es un proceso iterativo y por lo mismo, encuentra la aproximación casi con la misma rapidez que el método de Newton-Raphson. Claro que se corre el mismo riesgo de éste último de no converger a la raíz, mientras que el método de la regla falsa si es seguro. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

35

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------2 Ejemplo 1. Utilicemos el método de la secante para aproximar la raíz de f (x) = e x – x, comenzando con x0 = 0 y x1 = 1, hasta que a < 1%.

Solución Tenemos que f (x0) = 1 y f (x1) = – 0.63210558; los que sustituimos en la fórmula de la secante para calcular la aproximación x2: x2 = x1 –

f x1 x0 x1 0.63210558 0 1 =1– = 0.612699837 f x0 f x1 1 0.63210558

El error aproximado es: a = 100

0612699837 1 % = 63.2 % 0.612699837

Como no alzamos el objetivo, continuamos con el proceso. Resumimos los resultados en la siguiente tabla: i 1 2 3

xi-1 xi 0 1 1 0.612699837 0.612699837 0.653442133

f (xi-1) 1 – 0.632121 0.074314

Concluimos que la aproximación a la raíz es: xr

f (xi) – 0.632121 0.074314 0.00097

xi+1 0.612699837 0.653442133 0.652917265

x4 = 0.652917265.

Ejemplo 2. Con el método de la secante aproximemos la raíz de f (x) = arc tg x – 2x + 1, comenzando con x0 = 0 y x1 = 1, hasta que a < 1 %. Solución: Calculamos los valores: f (xo) = arc tg 0 – 2 x 0 + 1 = 1 f (x1) = arc tg 1 – 2 x 1 + 1 = 0.214601836 Los que se sustituyen en la fórmula de la secante para obtener la aproximación x2: x2 = x1 –

f x1 x0 x1 0.214601836 0 1 =1– = 0.823315073 f x0 f x1 1 0.214601836

Con un error aproximado de:

a = 100

0.823315073 1 % = 21.46 % 0.823315073

Como aún no se logra el objetivo, continuamos con el proceso. Resumiendo los resultados en la siguiente tabla: I

xi-1

xi

f (xi-1)

f (xi)

xi+1

a

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

36

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------1 2 3

0 1 1 0.823315073 0.823315073 0.852330280

1 0.214602 0.042167

0.214602 0.042167 .001185

De lo cual se concluye que la aproximación a la raíz es: xr

0.823315073 0.852330280 0.853169121

21.4 3.40 0.09

x4 = 0.853169121

2.7 MÉTODO DE ITERACIÓN DE PUNTO FIJO Un punto fijo a una función g(x) es un número real x tal que: x = g(x) Este concepto aplicaremos para desarrollar un método que permita resolver numéricamente ecuaciones de la forma: f (x) = 0. En primer lugar, transformamos la ecuación f (x) = 0, de alguna manera, en otra de la forma x = g(x); esto se consigue, ya sea despejando x de alguno de los términos de la ecuación, o bien sumando x a ambos miembros de la ecuación. Veamos algunas formas de despejar x:. 1) Para la ecuación: f (x) = 2 x2 – x – 5 = 0 , algunas de las posibles formas de x = g(x) son: x = 2x2 – 5 x 5 x= 2 5 x= 2x 1 2x 2 x 5 x=x– 4x 1

Despejando x del segundo término: Despejando x del primer término: Factorizando x y despejándola: Aplicando Newton – Raphson: 2)

Para la ecuación cos x – x = 0 se puede obtener la forma: x = cos x

3)

Para la ecuación tg x – e-x = 0 se puede sumar de x en ambos miembros, obteniéndose: x = tg x – e-x + x.

Para continuar con el desarrollo del método, luego optamos por un valor x0, el que puede ser elegido próximo a la raíz xr utilizando el método gráfico. Asumiendo éste como una primera aproximación de xr. El método iterativo surge de seguir los siguientes pasos: Calculamos: x1 = g(x0)

y, verificamos si

a

= 100 ( x1 – x0)/ x1 < r. Si no es así:

Calculamos: x2 = g(x1)

y, verificamos si

a

= 100 ( x2 – x1)/ x2 < r. Si no es así:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

37

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------… y así sucesivamente, hasta calcular el valor: xn+1 = g(xn) que llegue a satisfacer la condición dada por el error de aproximación, es decir: a < r Convergencia: En general, dada la aproximación xi, la siguiente iteración se calcula con la fórmula: xi+1 = g (xi) Supongamos que la raíz verdadera es xr, es decir:

xr = g (xr)

Restando las últimas ecuaciones se obtiene:

xr – xi+1 = g (xr) – g (xi)

Por el Teorema del Valor Medio para derivadas, se sabe que si g (x) es continua en [a, b] y diferenciable en (a, b) entonces existe (a, b) tal que: g´ ( ) = En el presente caso, existe, al menos, un g´ ( ) =

gb ga b a

en el intervalo determinado por xi y xr tal que: g xr xr

g xi xi

De aquí se tiene que:

g (xr) – g (xi) = g´ ( ).(xr – xi)

O bien:

xr – xi+1 = g´( ).(xr – xi)

Tomando valor absoluto en ambos lados, resulta: │xr – xi+1│ = │g´( )│.│(xr – xi)│ Observemos que el término │xr – xi+1│es precisamente el error absoluto en la (i + 1)ésima iteración, mientras que el término │(xr – xi)│ corresponde al error absoluto en la iésima iteración. Por lo tanto, solamente si │g´( )│< 1, se disminuirá el error en la siguiente iteración. En caso contrario, el error irá en aumento. En resumen, el método de iteración del punto fijo converge a la raíz si │g´( )│< 1 para todo [a, b] que contiene a la raíz de f(x) y donde g (x) es continua y diferenciable, pero diverge si │g´( )│> 1 en dicho intervalo. Analicemos ejemplos anteriores: En el ejemplo: g (x) = cos x claramente se cumple la condición de que │g´(x)│< 1; por lo tanto el método sí converge a la raíz. En el ejemplo: g (x) = x + tg x – e-x, se tiene: │g´(x)│= │1 + sec2x + e-x│> 1. Por lo tanto, el método no converge a la raíz. Gráficamente, el método de punto fijo, es un método de dos funciones: una formada por la recta y = x y la otra por y =g(x). La raíz es el punto de intersección de ambas curvas, tal como lo podemos observar en la figura 2.9.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

38

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

Figura 2.9 Para entender el método veamos los ejemplos siguientes: Ejemplo 1. Usemos el método de iteración del punto fijo para aproximar la raíz de f (x) = x2 – 5x – ex hasta que a < 1. Solución. Utilizamos el método gráfico (ver figura 2.10) para determinar el intervalo donde se halla la raíz de esta función, observándose que está en [– 0.5; 0]

Figura 2.10 Como primer paso del método, despejamos x del término lineal, observamos que la ecuación equivalente tiene la forma: x=

x2

ex 5

g (x) =

x2

ex 5

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

39

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Luego analizamos la convergencia, calculando la derivada de g(x): g´(x) =

2x e x . 5

Figura 2.11 Graficamos la gráfica de esta derivada (ver figura 2.11) y, vemos que │g´(x)│< 1 para todo x [–0.5, 0] , por lo que deducimos que el método sí converge a la raíz buscada. Aplicando la fórmula iterativa, empezando con x0 = 0, y obtenemos: x1 = g (xo) = – 0.2 Con lo que el error de la aproximación es de 100%. Aplicando nuevamente la fórmula iterativa, obtenemos: x2 = g (x1) = – 0.1557461506 Y el error de la aproximación es igual a 28.41%. Continuamos con el método hasta la quinta iteración, para reducir el error menor al 1%. Los resultados los resumimos en la siguiente tabla: Aproximación a la raíz 0 –0.2 –0.1557461506 –0.1663039075 –0.163826372 –0.164410064

Error aproximada 100% 28.41% 6.34% 1.51% 0.35%

Resultando que la aproximación buscada es: x5 = – 0.164410064 Ejemplo 2. Usemos el método de iteración del punto fijo para aproximar la raíz de la función: f (x) = cos x – x, comenzando con xo = 0 hasta que a < 1 ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

40

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Solución Como ya vimos anteriormente, el método sí converge a la raíz, por tanto, aplicando la fórmula iterativa se tiene: x1 = cos xo = cos 0 = 1 Con el error de la aproximación es de 100%. Aplicando nuevamente la fórmula iterativa y tenemos: x2 = cos x1 = cos 1 = 0.510302305 En este caso, el error aproximado de 85.08 %. Se puede intuir que el error aproximado se irá reduciendo muy lentamente. En efecto, necesitaremos de 13 iteraciones para lograr reducir el error aproximado a menos del 1%. El resultado final que se obtiene es: x13 = 0.7414250866 Siendo el error de la aproximación igual al 0.78 %

2.8 PROBLEMAS 3 1. Usa el método de bisección para aproximar la raíz de f (x) = e x – 2x comenzando en el intervalo [0.75, 1] y hasta que a < 1

2. Usa el método de bisección para aproximar la raíz de f (x) = comenzando en el intervalo [0.5, 1] y hasta que a < 1

+ 1

x 2 1 – tg x

3. Usa el método de la regla falsa para aproximar la raíz de f (x) = 4 – x2 – x3 comenzando en el intervalo [1, 2] y hasta que a < 1. 4. Usa el método de la regla falsa para aproximar la raíz de f (x) = ln x + x2 – 4 comenzando en el intervalo [1, 2] y hasta que a < 1. 5. Usa el método de Newton-Raphson para aproximar la raíz de f (x) = 1 – x2 – arc tg x comenzando con xo = 0.5 y hasta que a < 1. 6. Usa el método de Newton-Raphson para aproximar la raíz de f (x) = cos x – x, comenzando con xo = 1 y hasta que a < 1. 7. Usa el método de la secante para aproximar la raíz de f (x) = arc sen x – e-2x comenzando con xo = 0, x1 = 0.5 y hasta que a < 1. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

41

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------8. Usa el método de la secante para aproximar la raíz de f (x) = e-x – x comenzando con xo = 0, x1 = 1 y hasta que a < 1. 9. Usa el método de iteración del punto fijo para aproximar la raíz de f (x) = sen x + x – 1 comenzando con xo = 0.52 y hasta que a < 1. 10. Usa el método de iteración del punto fijo para aproximar la raíz de f (x) = ln x + 2x – 4 comenzando con xo = 1.5 y hasta que a < 1. 11. La función f(x) = e–x + 4x2 – 5, tiene una raíz en el intervalo [1; 1.5], usa ocho iteraciones del método de la bisección para aproximar la raíz de f(x). Tabula el error después de cada iteración. 12. Encontrar la raíz cerca de x = 0.75 de f (x) = ex–1 – 5x3 empezando con x0 = 0.75. ¿Cuán exacta es la estimación después de cuatro iteraciones del método de Newton? Tabula el número de dígitos correctos en cada iteración del método de Newton.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

42

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------UNIDAD DIDÁCTICA 3 SISTEMAS DE ECUACIONES LINEALES

3.1 INTRODUCCIÓN. En diversos estudios de ingeniería, con frecuencia encontramos problemas que involucran sistemas de ecuaciones lineales que se estudian en la asignatura de álgebra lineal. Se han desarrollado muchos métodos para resolver sistemas de ecuaciones, lo que demuestra que es aparente la sencillez de su solución y, que existen una serie de fallas en los diferentes algoritmos. Las soluciones pueden ser exactas o por aproximaciones sucesivas. Recordemos, primeramente que un sistema de ecuaciones lineales, es un conjunto de ecuaciones de primer grado que deben resolverse simultáneamente, cuya forma general es: a11x1

a12 x 2



a1n x n

b1

a 21x1

a 22 x 2



a 2n x n

b2









a n1 x1

an2 x2

a nn x n

bn



Donde los coeficientes aij (i = 1, 2, 3, ..., n y j = 1, 2, 3, ..., n) y los términos independientes bi son constantes y xj son las variables o incógnitas del sistema. El problema consiste en encontrar los valores desconocidos de las variables x1, x2 … xn que satisfacen las n ecuaciones. 3.1.1. Tipos de sistemas de ecuaciones lineales Recordemos que los sistemas de ecuaciones se pueden clasificar según el número de soluciones que pueden presentar. De acuerdo con ese caso se pueden presentar los siguientes casos: Sistema incompatible si no tiene ninguna solución. Sistema compatible si tiene alguna solución, en este caso además puede distinguirse entre: o Sistema compatible determinado cuando tiene un número finito de soluciones. o Sistema compatible indeterminado cuando admite un conjunto infinito de soluciones. 3.1.2. Expresión matricial de un sistema de ecuaciones lineales Recordemos asimismo que un sistema de n ecuaciones lineales con n incógnitas, también puede escribirse en la siguiente forma matricial: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

43

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

a11x1

a12 x2

a21x1

a22 x2 

an1 x1

b1

... a1n xn ... a2n xn

an2 2 x2

=

... ann xn

b2  bn

Si aplicamos la definición de producto entre matrices, esta ecuación matricial, podemos escribirla en forma de producto matricial:

a11

a12

a21 a22   an1

an2

x1 ... a1n x2 ... a2n x   xn ... ann

b1 b2 =  bn

En cuyos los elementos aij , el subíndice i se refiere a la fila, y el subíndice j se refiere a la columna donde está ubicado el elemento correspondiente. De esta manera definimos tres matrices:

a11 a) La matriz de coeficientes:

b) La matriz de incógnitas:

A=

a12

... a1n

a21 a22  

... a2n 

an1

... ann

an2

x1 x2 X=  xn

c) La matriz de términos independientes:

b1 b2 B=  bn

De este modo el sistema es equivalente a la ecuación matricial: A

X=B

De los variados y diferentes métodos para resolver sistemas de ecuaciones lineales, solo estudiaremos algunos de ellos en la presente unidad didáctica. Entre los objetivos de la presente unidad didáctica no está previsto repetir el curso de álgebra lineal, por lo que únicamente recordaremos algunas condiciones básicas para la solución de sistemas de ecuaciones lineales.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

44

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------3.1.3. Teorema fundamental de equivalencia. Si en un sistema de ecuaciones lineales se sustituye una de ellas por una combinación lineal de las ecuaciones, se obtiene un nuevo sistema equivalente que sustituye al anterior. 3.2 TIPOS DE MATRICES Las matrices se utilizan en el cálculo numérico, en la resolución de sistemas de ecuaciones lineales, de las ecuaciones diferenciales y de las derivadas parciales. Además de su utilidad para el estudio de sistemas de ecuaciones lineales, las matrices aparecen de forma natural en geometría, estadística, análisis estructural, economía, informática, física, etc. Hay algunas matrices que aparecen frecuentemente y que según su forma, sus elementos, y otras características, reciben nombres diferentes: a) Matriz cuadrada de orden n x n: Es una matriz de n fila y n columnas, por ejemplo, la matriz A de los coeficientes del sistema de ecuaciones lineales simultáneas. b) Matriz Rectangular de orden n x m: Es una matriz de n fila y m columnas, donde n m c) Vector Fila de orden de orden 1 x n: Es una matriz de una sola fila y n columnas. d) Vector Columna de orden n x 1: Es una matriz de n filas y una sola columna. e) Matriz escalonada: Es una matriz de orden n x m, m > n, en la cual los elementos aij = 0, para i > j Ejemplos: La matriz A es escalonada, y la matriz B no es escalonada:

A=

6 12

3

0

0

5

7 10

0

0

7

14

0

0

0

5

2 3

B= 0 5 0 1

f) Matriz Diagonal: Es una matriz cuadrada en la que todos los elementos aij = 0 para i j, es decir, es una matriz en la que solo los elementos de la diagonal principal no son ceros. Un caso particular de la matriz diagonal es la matriz unitaria o matriz identidad In, en la que aij = 1 para i = j, es decir:

In =

1

0 ... 0

0

1 ... 0

... ... ... ... 0

0

0

1

g) Matriz aumentada de la matriz A de orden n x m y la matriz B de orden n x p, es otra matriz, C, de orden n x (m + p) en la que, a las m columnas de la matriz A se unen a las p columnas de la matriz B y, se denota: C = (A  B) ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

45

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Ejemplo: Dadas las matriz A y B, hallar la matriza aumentada (A  B). 1

A=

2

2 5

0.5 1.5

0.2 1.75

B=

1 1

Solución:

5

(A  B) =

2

0 10

2

0.5

5

0.2 1.75

5

1.5 0 1 10

h) Matriz inversa: Si A es una matriz de orden n x n; se denomina matriz inversa de A, a una matriz B del mismo orden n x n tal que: A x B = In B x A = In

Y también:

Y se denota B = A-1 para representar la matriz inversa; luego resulta: A x A-1 = A-1 x A = In Cuando la matriz inversa, de una matriz dada, existe, ésta es única, sin embargo debe destacarse que no siempre existe la matriz inversa. El álgebra lineal prueba que la matriz inversa A-1 existe si y solo si el determinante de dicha matriz A es distinto de cero. 3.3 ALGUNAS OPERACIONES ELEMENTALES CON MATRICES Ciertas operaciones sobre filas y/o columnas de una matriz tienen las virtudes de conservar, o modificar de manera controlada, ciertos objetos asociados: determinante, rango, conjunto de las soluciones del sistema de ecuaciones asociado. Presentamos estas operaciones (“operaciones elementales”), empezando por su utilización en un calculo de determinante. Seguimos por su utilización en la inversión de una matriz cuadrada regular. Se plantea el problema de usar estas operaciones de manera ordenada: nos conduce a presentar las matrices escalonadas y el algoritmo de reducción de Gauss. Seguimos presentando como usar las operaciones elementales para calcular el rango. Acabamos remarcando que las operaciones elementales pueden ser realizadas por multiplicación por ciertas matrices. Para una matriz A recordaremos tres operaciones elementales por filas. Cuando sobre una matriz A se aplican estas operaciones elementales se obtiene una matriz equivalente, y utilizaremos el símbolo de equivalencia “ ”.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

46

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------I.- Intercambio de dos filas: 4

Ejemplo: Si se intercambian las filas 1 y 3:

5

8

1 2

3

6

7

6

~

7

10

1 2

10

4

3

5

8

II.- Multiplicación de una fila por una constante distinta de cero (escalar)

Ejemplo: Si se multiplica la fila 3 por la constante 2:

2

3

0

3

1

2

3

1

4 ~ 0 3 4 7 2 6 14

1 3

III.- Suma de una fila a otra fila: Se aplica el teorema que dice: “Si en un sistema de ecuaciones lineales se sustituye una de ellas por una combinación lineal de las ecuaciones, se obtiene un nuevo sistema equivalente que sustituye al anterior” 4

3

Ejemplo: Si se suma la fila 3 a la fila 2: 6 1

1

7

4

3

2 ~ 7 1 1

1

1

6

1

1

1

Las operaciones II y III se combinan para sumar un múltiplo de una fila a otra fila. 1

Ejemplo: (i) Se comienza con la matriz:

2 3

1 5 6 3

0 7 1

2 3 2 2 4 6 1 5 6 ~ 1 5 6

(ii) Se multiplica la fila 1 por la constante 2:

3 2

(iii) Se suma la fila 1 a la fila 2:

4 6

0 7 2 4

3

0 7

6

1 5 6 ~ 1 9 12 3 0 7 3 0 7 2 4

6

1 (iv) Finalmente se multiplica por la fila 1: 1 9 12 2 3 0 7

1

Ahorrando pasos se puede escribir simplemente:

2 3 2,

1 5 6 3

0 7

1 2

1 2

1 2

3

~ 1 9 12 3 0 7

1 2

3

~ 1 9 12 3 0 7

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

47

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Estas operaciones elementales se utilizan para “hacer ceros” debajo de algún elemento aij 0 con el objetivo final de transformar una matriz A en una matriz escalonada. Ejemplo: Hagamos ceros debajo del elemento a11 en la siguiente matriz: 1

3

5

2 5

A=

3

6

0

1

Solución. Observamos que para lograr el objetivo, debemos multiplicar la fila 1 por 2, y sumarla a la fila 2. Del mismo modo podemos multiplicar la misma fila 1 por –3 para luego sumarla a la fila 3, obteniendo: 1

3

5 2, 6

2 5 3

0

1

3

5

~ 0 0

11

4

9

14

3

1

Ejemplos: 1) Usando operaciones elementales, escalonemos la siguiente matriz:

A=

6

2

2

12

8

6 10

3

4

13 9

6

4

3

1

8

Solución. La notación se explica por sí sola: 6

2

2

12

8

6 10

3

4

13 9

6

4

1 , 1 2

2,

~

3

1

8

6

2

2

4

0

4

2

2

0

12 8

1

0

2

3

4

6 1 0 3, 2~ 0 0

6

2

2

4

0

4

2

2

0

12 8

1

0

2

3

2 2

4

4 2

2

0

2

5

0

4

3

4

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

48

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------6

2 2

4

0

4 2

2

0

0

2

5

0

0

4

3

~

2

6

2 2

4

0

4 2

2

0

0

2

0

0

0

1 2

2)

5 7

3

Escalonemos la siguiente matriz: A = 4 5 6 7 8 10

Solución. Procedemos como en el ejercicio anterior: 1 2

3

4 5

6

7 8 10

4,

1

7

~ 0 0

2

3

3

6

6

11

1

2 ~ 0 0

2

3

3 0

6 1

Ahora, tenemos todas las herramientas para estudiar los dos primeros métodos numéricos de solución a sistemas de ecuaciones lineales. 3.4 ELIMINACIÓN GAUSSIANA SIMPLE Este método es aplicado para resolver sistemas lineales de la forma: A .X = B. El método de eliminación Gaussiana (simple), consiste en escalonar la matriz aumentada del sistema: (A  B), para obtener un sistema equivalente de la forma: a11x1

a12 x 2



a1n x n

b1

a´ 22 x 2



a´ 2n x n

b´ 2





a´ nn x n

b´ n

Donde la notación a´ij se usa para denotar que los elementos aij cambiaron. De este modo, se despejan las incógnitas comenzando con la última ecuación y sustituyendo en la penúltima y despejando hacia arriba. Por esta razón, muchas veces se dice que el método de eliminación Gaussiana consiste en la eliminación hacia adelante y sustitución hacia atrás. Ejemplo 1: Usando el método de eliminación Gaussiana (simple), resolver el siguiente sistema de ecuaciones:

x1

2 x2

3 x3

1

4 x1 7 x1

5 x2 8 x2

6 x3 10 x3

2 5

Solución: Escalonando la matriz aumentada del sistema, se obtiene: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

49

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------1

2

3

4

5

6

7

8

10

1

4,

1

7

2

~

5

2

3

1

1

0

3

6

6

0

6

11

2

2

2~ 0 0

3 0

3

1

6 1

6 10

Y multiplicando la segunda fila entre –1/3, se tiene la matriz equivalente: 1

2

0

3

3

0

1

6

0

1

1 1 ~ 0 3 0

6 10

2

3

1 0

1

2

2

1

10

Por lo tanto, el sistema dado equivale a:

x1

2 x2

3 x3

1

x2

2 x3 x3

2 10

De esta ecuación obtenemos x3 = 10; valor que sustituido en la ecuación de arriba nos permite obtener x2 = – 18; sustituyendo estos valores en la ecuación de arriba obtenemos x1 = 7. Por lo tanto, la solución del sistema es: x1

7, x2

18 , x3

10

Ejemplo 2: Usando eliminación Gaussiana (simple), resolvamos el sistema siguiente:

x1

3x2

2 x3

2 x1 3x1

5 x2 x2

x3 6 x3

12 10 4

Solución. Escalonando la matriz aumentada del sistema, obtenemos: 1 2 3

3 5 1

2 1 6

12

2, 3

10

1

~ 0 0

4

3

2 12

1

5 14

10 0 40

1 3

10 ~ 0 0

2

12

1

5

14

0

50

100

Por lo tanto, el sistema equivalente resulta:

x1

3x2

2 x3

12

x2

5 x3 50 x3

14 100

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

50

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------De la tercera ecuación se obtiene x3 = 2; valor que sustituido arriba da x2 = 4; y sustituyendo arriba para obtener x1 = 4. Por lo tanto la solución del sistema es: x1

4, x2

4 , x3

2

El método de eliminación Gaussiana (simple) puede presentar un problema cuando uno de los elementos que se usan para hacer ceros, es cero. Así por ejemplo, en la matriz: 1

2

3

7 6

0

0

10

0

5

8

3

El elemento a22 = 0 no puede usarse para hacer ceros. En este caso el problema se puede resolver intercambiando las filas 2 y 3. En este caso el resultado que se obtiene es la matriz escalonada: 1

2

3

0

5

8

0

0

10

7 3 6

El problema puede presentarse también cuando el elemento aquel es muy cercano a cero. Ejemplo: Resolver el siguiente sistema, usando eliminación Gaussiana (simple)

Solución.

0.00005x1

5 x2

3.00005

x1

x2

1

Usando eliminación Gaussiana (simple) se obtiene:

0.00005 5 3.33335 1

1

1

1 0.00005 0.00005 ~ 0

5 99999

3.33335 66666

Que da el sistema equivalente:

0.00005x1

5 x2 99999x2

De donde se obtiene: x2 =

3.00005 66666

2 ; valor que sustituido arriba, permite obtener: 3

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

51

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------3.33335

x1 =

5

2 3

0.00005

El resultado cambia de acuerdo al número de cifras significativas que se usen: Nº de cifras significativas x2 3 0.667 4 0.6667 5 0.66667 6 0.666667 7 0.6666667

Error relativo porcentual (*) 10,000 % 1,000 % 100 % 10 % 1%

x1 -33 -3 0 .3 0.33

*) Para calcular este error se tomó el valor verdadero de x1 = 1/3. Resolviendo el mismo sistema, pero intercambiando las filas 1 y 2, se obtiene:

1

1

0.00005

1

0.00005 5 3.33335

~

1

1

1

0 4.99995 3.3333

Lo que da el sistema equivalente:

x1

De donde se obtiene x2 =

x2

1

4.99995x2

3.3333

2 ; sustituyendo arriba da: 3

x1 = 1 –

2 1 = 3 3

Tomando distintas cifras significativas, los resultados se muestran en la siguiente tabla: Número de cifras significativas

x2

x1

Error relativo porcentual

3 4 5 6 7

0.667 0.6667 0.66667 0.666667 0.6666667

0.333 0.3333 0.33333 0.333333 0.3333333

0.1 % 0.01 % 0.001 % 0.0001 % 0.00001 %

En este último caso, se ve que el error relativo porcentual no varía como en la solución anterior. De este modo, se deduce que los elementos que son cercanos a cero, son elementos malos para hacer ceros. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

52

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------3.5 ELIMINACIÓN GAUSSIANA CON PIVOTEO PARCIAL En general, para evitar este problema se elige como elemento para hacer ceros (elemento que recibe el nombre de elemento pivotal o simplemente pivote), aquel elemento de mayor valor absoluto de entre todos los candidatos. A este procedimiento se le llama método de eliminación Gaussiana con pivoteo parcial. El pivoteo parcial se puede resumir como sigue: Para elegir el elemento pivote en la primera columna, se selecciona aquel elemento de mayor valor absoluto de dicha primera columna y, toda la fila que contiene a este elemento, se la traslada como primera fila. Luego se procede a hacer ceros debajo de este pivote. Para elegir el elemento pivote en la segunda columna, exceptuando el elemento de la primera fila, nos olvidamos de toda la fila 1 y de toda la columna 1 y, procedemos a escoger el elemento de mayor valor absoluto de toda la segunda columna, luego, toda la fila que contiene a este elemento, se la relocaliza como segunda fila. De este modo procedemos a hacer ceros debajo de este pivote. Para la tercera columna, nos olvidándonos de las filas 1 y 2 y de las columnas 1 y 2 y, elegimos como su elemento pivote el de mayor valor absoluto de esta columna, luego, la fila que contiene a este elemento, se la localiza como tercera fila. Luego se procede a hacer ceros debajo de este pivote. Y, así sucesivamente. En el siguiente diagrama matricial, se muestran los elementos pivotes elegibles de cada columna, siguiendo el procedimiento anterior:

a11 a 21 a 22 a31 a32   a n1 a n 2

a33  a n3  a nn

Ejemplo 1: Usar el método de eliminación Gaussiana con pivoteo para resolver el siguiente sistema:

x1 2 x1 0.2 x1

2 x2

0.5 x3

5

5 x2 1.75x2

1.5 x3 x3

0 10

Solución. Escribimos la matriz aumentada del sistema: 1 2

2 5

0.2 1.75

0.5

5

1.5

0

1

10

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

53

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Para escoger el primer elemento pivote en la columna 1, se elige el elemento de mayor valor absoluto entre –1 , – 2 y – 0.2, siendo, por supuesto, el – 2 ; por lo tanto se intercambian las filas 1 y 2 (éste es el primer pivoteo realizado): 2

5

1

1.5

0

2

0.5

5

0.2 1.75

1

10

Luego, se procede a hacer ceros debajo del pivote. Para ello, se multiplica la fila 1 por 1 2

0 .2 2

y se la suma a la fila 2. También, se multiplicamos la fila 1 por

y se la suma a la

fila 3. Es decir: 2 1

5

1.5

0

2

0.5

5

0.2 1.75

1

10

1 , 2

0 .2 2

2

~ 0 0

5

1.5

0

0.5

0.25

1.25

0.85 10

5

Olvidándonos, ahora, de la fila 1 y de la columna 1, procedemos a escoger el pivote de la columna 2, es decir únicamente entre 0.5 y 1.25 que, como es obvio, resulta ser 1.25. Por lo tanto intercambiando las filas 2 y 3 (éste es el segundo pivoteo realizado): 2

5

1.5

0

0

1.25

0.85 10

0

0.5

0.25

5

Luego, procedemos a hacer ceros debajo del elemento pivote elegido. multiplicamos la fila 2 por

Para ello

0.5 y se suma a la fila 3 para obtener: 1.25

2

5

0

1.25

0

0

1.5

0

0.85 10 0.09

9

La que resulta ser una matriz escalonada. El sistema equivalente es:

2 x1

5 x2

1.5 x3

0

1.25x2

0.85x3 0.09 x3

10 9

Realizando la sustitución hacia arriba, se obtiene la siguiente solución del sistema: xx x1 = – 100,

x2 = – 60

y

x3 = – 75

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

54

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Ejemplo 2. Usar eliminación Gaussiana con pivoteo para resolver el siguiente sistema de ecuaciones:

0.4 x1

1.5 x2

0.75x3

20

15x2 9 x2

10 x3 2.5 x3

10 30

0.5 x1 10 x1

Solución. La matriz aumentada del sistema es: 0.4

1.5 0.75

20

0 .5

15

10

10

10

9

2.5

30

Seleccionamos como elemento pivote en la columna 1 al –10, intercambiar las filas 1 y 3 y, hacemos ceros debajo de dicho pivote: 10

9

2.5

0 .5

15

10

0.4

0 .5 0 .4 , 10 10

30 10

1.5 0.75

10

~

20

9

lo que obliga a

2.5

30

0

14.55 9.875

11.5

0

1.86

18.8

0.85

Ahora el elemento pivote en la columna 2 resulta ser el -14.55, el mismo que está bien colocado, y no hay necesidad de intercambiar filas. Procedemos a hacer ceros debajo del pivote: 10

9

2.5

30

0

14.55 9.875

11.5

0

1.86

18.8

0.85

1.86 ~ 14.55

10 0

9

2.5

14.55

9.875

0

0

0.412371

30 11.5 17.3299

Escribiendo el sistema equivalente, y resolviéndolo con la sustitución hacia arriba, se obtiene la solución del sistema: x1 = – 18.875;

x2 = 29.3125;

x3 = 42.0250

3.6 MÉTODO DE GAUSS - JORDÁN Este método utiliza las mismas técnicas que el de eliminación Gaussiana, incluyendo el pivoteo, pero con el objetivo de finalizar con una matriz de la siguiente forma: (A

B)

(In

Bˆ )

Donde In es la matriz identidad del mismo orden de la matriz A. Para lograr esto, se sigue la técnica del pivoteo buscando hacer ceros hacia abajo y hacia arriba de los términos de la diagonal principal. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

55

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Ejemplo 1:

Usar el método de Gauss-Jordán para resolver el siguiente sistema:

x1

2 x2

2 x1 0.2 x1

0.5 x3

5 x2 1.75x2

5

1.5 x3 x3

0 10

1

Solución.

2

2

Comenzamos con la matriz aumentada:

0.5

5

0.2 1.75

5

1.5

0

1

10

Intercambiamos las filas 1 y 2 y procede a hacer el primer pivoteo: 2

5

1

1.5

1 , 2

0

2

0.5

5

0.2 1.75

1

10

0 .2 2

2

~ 0 0

5

1.5

0

0.5

0.25

1.25

0.85 10

5

Colocamos adecuadamente el segundo pivote intercambiando las filas 2 y 3. Luego, 5 hacemos ceros arriba del pivote 1.25, multiplicando la fila 2 por y sumándola a la 1.25 fila 1. A continuación hacemos ceros debajo del mismo pivote, para lo cual se multiplica 0 .5 la misma fila 2 por y se la suma al fila 3. Lo que arroja el siguiente resultado: 1.25 2

5

1.5

0

0

1.25

0.85 10

0

0.5

0.25

5 , 1.25

5

2 0 0.5 ~ 0 1.25 1.25 0 0

1.9 0.85 0.09

40 10 9

Luego, procedemos a hacer ceros arriba del pivote 0.09, para lo cual se multiplica la fila 3 0.85 1 .9 por y se la suma a la fila 2; igualmente multiplicamos la fila 3 por y la 0.09 0.09 sumamos a la fila 1. Todo ese conjunto de operaciones da los siguientes resultados: 2

0

0

1.25

0

0

1.9 0.85 0.09

40

2

10 9

0.85 , 0.09

~ 0 1.9 0 0.09

0

0

150

1.25

0

75

0

0.09

9

Finalmente para obtener los unos en la diagonal principal, multiplicamos las filas 1, 2, y 1 1 1 3 por , y , respectivamente. Obteniéndose: 2 1.25 0.09

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

56

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------1 2 2 0 0 150 1 0 0 1 0 1.25 0 75 ~ 0 1 0 1.25 0 0 0.09 9 0 0 1 1 0.09

75 60 100

Con lo que obtenemos la siguiente solución del sistema de ecuaciones: x1 = – 75;

x2 = – 60;

x3 = – 100

Ejemplo 2. Usar el método de Gauss-Jordán para resolver el siguiente sistema:

x1

2 x2

3x3

1

0.4 x1 0.5 x1

2 x2 3x2

x3 x3

10 15

Solución. La matriz aumentada del sistema es:

1

2

0.4

2

0.5

3

1

1 10

3

1

15

Se observa que el primer elemento pivote está bien colocado y por lo tanto no hay necesidad de intercambiar filas. Por lo tanto hacemos ceros debajo del pivote a11 = 1; para ello, multiplicamos la fila 1 por 0.4 y la sumamos a la fila 2, y luego multiplicamos la misma fila 1 por –0.5 y luego la sumamos a la fila 3. Es decir: 1

2

0. 4

2

0.5

3

1 0.4, 1 10

3

1

15

1

0.5

2

~ 0 2.8 0 4

3

1

0.2

10.4

0.5 14.5

Para elegir el segundo pivote, debemos elegir el elemento de mayor valor absoluto entre a22 = 2.8 y a32 = – 4, resultando este último. Por lo tanto, intercambiamos la fila 2 y la fila 3, y luego, procedemos a hacer ceros arriba y abajo del segundo elemento pivote; para lo cual, multiplicamos la fila 2 por 0.5 y la sumamos a la fila 1, y de igual modo, 2.8 multiplicamos la misma fila 2 por y la sumamos a la fila 3. Esto es: 4 1 0

2 4

0 2.8

3

1

1 2.8 0.5 14.5 0.5, ~ 0 4 0.2 10.4 0

0 4 0

2.75

8.25

0.5

14.5

0.15 20.55

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

57

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------El tercer elemento pivote es a33 = – 0.15. Para hacer ceros arriba de este elemento, 0 .5 multiplicamos la fila 3 por y la sumamos a la fila 2, e igualmente multiplicamos la 0.15 2.75 misma fila 3 por y la sumamos a la fila 1. De todo esto resulta: 0.15 1 0 0

0 4 0

2.75

8.25

1

0.5

14.5

~ 0 0.5 2.75 0 , 0.15 0.15

0.15 20.55

0 4 0

0

385

0

54

0.15 20.55

Finalmente, obtenemos los unos de la diagonal, multiplicando la fila 2 por por

1 y la fila 3 4

1 . Esto es: 0.15 1 0 0

0 4 0

0

385

0

54

1 0 0 1 ~ 0 1 0 4 0 0 1 1 0.15

0.15 20.55

385 13.5 137

Por lo tanto, la solución del sistema de ecuaciones es: x1 = 385;

x2 = 13.5;

x3 = –137

3.7 APLICACIÓN DEL MÉTODO DE GAUSS – JORDAN AL CÁLCULO DE LA MATRIZ INVERSA. El método de Gauss-Jordán se basa en la siguiente lógica: (A  In)

(In  A-1)

Es decir, se comienza por escribir la matriz A y, a su derecha se aumenta la matriz identidad In del mismo orden que la matriz A; luego se aplica el método de Gauss-Jordán para hacer los ceros y unos y obtener del lado izquierdo la matriz identidad In. En el lado derecho lo que se obtiene es la matriz inversa de A-1. Ejemplo 1.

Usar el método de Gauss-Jordán para calcular la matriz inversa de la 4 11 siguiente matriz: A= 1 3 Solución. A esta matriz A se aumenta la matriz identidad I2, obteniéndose:

4 11 1 0 1

3

0 1

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

58

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------El primer elemento pivote a11 = 4 está bien colocado y se procede a hacer ceros debajo de 1 este elemento. Para ello, multiplicamos la fila 1 por y la sumamos a la fila 2. Esto es: 4

1 4 11 4~ 0 0.25

4 11 1 0 1

3

0 1

1

0

0.25 1

El siguiente elemento pivote es a22 = 0.25. Para hacer ceros arriba de este elemento, 11 multiplicamos la fila 2 por y la sumamos a la fila 1. Esto da: 0.25

4

11

1

0 0.25

0

4 0 11 ~ 0 0.25 0.25

0.25 1

12

44

0.25

1

Finalmente, determinamos los unos de la diagonal principal, se multiplicando la fila 1 por 1 1 y el fila 2 por . Lo que da la matriz final: 4 0.25

4

0

0 0.25

12

44

0.25

1

1 4 ~ 1 0 1 0 1 0.25

Por lo tanto, se concluye que la matriz inversa de A es:

3

11

1

4

A-1 =

3 1

11 4

Ejemplo 2. Usar el método de Gauss-Jordán para calcular la matriz inversa de: 2

A=

0.5 0.3125

4 1.2

0 0.1

0.625 3.125

Solución. Aumentamos a la derecha de la matriz A la matriz identidad y, observamos que el primer elemento pivote a11 = 2 está bien colocado, procediendo a hacer ceros debajo de este elemento. Para lo cual multiplicamos la fila 1 por multiplica la misma fila 1 por

0.5 y la sumamos a la fila 2; igualmente, 2

0.3125 y la sumamos a la fila 3. Esto es: 2

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

59

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------2 0.5 0.3125

4

0

1.2

1 0 0

0.5 0.3125 , 2 2

0.1 0 1 0

0.625 3.125 0 0 1

2

4

~ 0 0

0.2

0 0.1

1

0 0

0.25

1 0

1.25 3.125 0.15625 0 1

Para el segundo elemento pivote, escogemos el elemento con mayor valor absoluto entre a22 = 0.2 y a32 = –1.25, siendo éste último. Por lo que intercambiamos la fila 2 y la fila 3. Luego procedemos a hacer ceros arriba y 4 abajo del segundo elemento pivote; para ello multiplicamos la fila 2 por y la 1.25 0.2 sumamos a la fila 1, y del mismo modo multiplicamos la misma fila 2 por y la 1.25 sumamos a la fila 3. Es decir: 2

4

0

1

0 0

0

1.25 3.125 0.15625 0 1

0

0.2

4 0.2 ~ , 1.25 1.25

3.125 0.15625 1 0

2

0

10

0

1.25 3.125 0.15625 0

0

0

0.4

0.5

0

0.275

1

3.2 1 0.16

El tercer elemento pivote resulta a33 = 0.4. Para hacer ceros arriba de este elemento, 3.125 multiplicamos la fila 3 por y la sumamos a la fila 2, e igualmente se multiplica la 0.4 10 misma fila 3 por y se la suma a la fila 1. Este proceso tiene el siguiente resultado: 0 .4 2

0

0

1.25 3.125 0.15625 0

0

10

0

2

~ 0 0

0.4

0.275

0

0

1.25

0

0

0.5

0.4

0 1

7.375 1.9921875 0.275

3.2 1

~ 3.125 10 , 0.4 0.4

0.16

25

0.8

7.8125 1

0.25 0.16

Finalmente, determinamos los unos de la diagonal principal. Para ello multiplicamos las 1 1 1 filas 1, 2 y 3 por , y , respectivamente. Esto da la matriz final: 2 1.25 0.4

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

60

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------1 0 0

3.6875

12.5 0.4

0 1 0 1.59375 6.25 0.2 0 0 1

0.6875

2.5

0.4

Por lo tanto, la matriz inversa de A es: 3.6875

12.5 0.4

-1

A = 1.59375 6.25 0.2 06872 2.5 0.4 3.8 MÉTODO DE GAUSS-SEIDEL El método de Gauss-Seidel, es un método iterativo y por lo mismo, resulta ser un método bastante eficiente. Se considera el siguiente sistema de ecuaciones: a11x1

a12 x 2



a1n x n

b1

a 21x1

a 22 x 2



a 2n x n

b2









a n1 x1

an2 x2

a nn x n

bn



Despejando de las n ecuaciones las n variables: x1, x2, …, xn, se obtiene: b x1 = 1

a12 x 2 ... a1n x n a11

(1)

b x2 = 2

a 21 x1 ... a 2n x n a 22

(2)

…. b xn = n

a n1 x1

... a nn 1 x n 1 a nn

(n)

Estas ecuaciones dan origen a las fórmulas iterativas. Para comenzar este proceso iterativo, en la ecuación (1), se da el valor cero a las variables x2, …, xn; lo que permite obtener un primer valor para x1, resultando: b x1 = 1 a11

Seguidamente se sustituye este valor de x1 en la ecuación (2), y las restantes variables x3, …, xn siguen manteniendo el valor de cero. Esto permite obtener el siguiente valor para x2:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

61

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------b2

x2 =

a 21

b1 a11

a 22

Estos últimos valores de x1 y x2, son sustituidos en la ecuación 3, mientras que x4, …, xn siguen manteniendo el valor de cero; y así sucesivamente hasta llegar a la última ecuación. Todos estos pasos, proporcionarán una lista de primeros valores para las incógnitas, los que conforman el primer paso en el proceso iterativo. Obteniéndose los primeros valore: x1 = 1 ,

x2 =

…,

2,

xn =

n

Se repite este proceso, pero ahora sustituyendo estos últimos valores en vez de los ceros inicialmente asignados. Con ello se obtiene una segunda lista de valores para cada una de las incógnitas. Teniéndose: x1 = 1 ,

x2 =

…,

2,

xn =

n

A partir de estos valores se pueden calcular los errores aproximados relativos para cada una de las incógnitas. Obteniéndose la siguiente lista de errores: 1 ,1 = 100

1 % 1

2 ,2 = 100

2 % 2

…

n

n

,n = 100

%

n

Se repite este proceso hasta que:

,i
xj, y si la función de interpolación está dada por la expresión y = f (x), entonces la diferencia finita dividida, para estos dos puntos, está definida por: f xi

f [xi, xj] =

f xj

xi

xj

b) Si se tienen tres puntos xi, xj y xk, siendo xi > xj > xk, y si la función de interpolación está dada por y = f (x), entonces la diferencia finita dividida estará definida por:

f [xi, xj, xk] =

f xi ,x j xi

f x j ,xk xk

c) Y, en general, si se tienen n + 1 puntos xn, xn-i … x1, xo, siendo xn > xn-i > … > x1 > xo, y si la función de interpolación está definida por y = f (x), entonces la diferencia finita dividida está dada por: f [xn, xn-1, …, x1, xo] =

f x n , , x1 f x n 1 , x0 x n x0

Ejemplo: Se tienen los puntos x3, x2, x1 y x0, de modo que x3 > x2 > x1 > x0. Para la función de interpolación f (x), se tienen las siguientes diferencias finitas divididas:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

74

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------f [x3, x2, x1, x0] =

f x3 , x 2 , x1 x3

f x 2 , x1 , x0 x0

f [x3, x2, x1] =

f x3 , x 2 x3

f x 2 , x1 x1

f [x2, x1, x0] =

f x 2 , x1 x2

f x1 , x0 x0

f [x3, x2] =

f x3 x3

f x2 x2

f x2 , x1 =

f ( x 2 ) f ( x1 ) x 2 x1

f x1 , x0 =

f ( x1 ) f ( x0 ) x1 x0

4.3 POLINOMIO DE INTERPOLACIÓN DE NEWTON CON DIFERENCIAS FINITAS DIVIDIDAS Dados los “n + 1” pares ordenados siguientes: x y

xo yo

x1 y1

… …

xn yn

El polinomio de interpolación de Newton se define como sigue: f (x) = bo + b1(x – xo) + b2 (x – xo)(x – x1) + ... + bn (x – xo)(x – x1) ... (x – xn-1) Donde:

bo = f (xo) f ( x1 ) f ( x0 ) x1 x0

b1 = f x1 , x0 =

b2 = f [x2, x1, x0] =

f x 2 , x1 x2

f x1 , x0 x0

… bn = f [xn, xn-1, …, x1, xo] =

f x n , , x1 f x n 1 , x0 x n x0

Reemplazando estos valores en el polinomio de Newton, finalmente obtenemos: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

75

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------f(x) = f(x0) + f[x1, x0](x – x0) + f[x2, x1, x0](x – x0)(x – x1)+ … + + f[xn, …x1, x0](x – x0) (x – x1)… (x – xn-1) Para calcular los coeficientes b0, b1,…, bn, es conveniente construir una tabla de diferencias finitas divididas como la siguiente: b0

b1

b2

b3

b4

x0

f (x 0)

f [x 1, x 0]

f [x 2, x 1, x 0]

f [x 3, x 2, x 1, x 0]

x1

f (x 1)

f [x 2, x 1]

f [x 3, x 2, x 1]

x2

f (x 2)

f [x 3, x 2]

x3

f (x 3)

x4 M

f (x 4) M

f [x 4, x 3] M

f [x 4, x 3, x 2] M

f [x 4, x 3, x 2, x 1] M

f [x 4, x 3, x 2, x 1, x 0] M

M

Observamos que los coeficientes del polinomio de interpolación de Newton, se encuentran en la parte superior de la tabla de diferencias finitas divididas. Ejemplo 1. Calculemos el polinomio de diferencias divididas finitas de Newton para: x y

–2 4

–1 6

2 9

3 3

Solución. Procedemos con la siguiente lógica:

Por lo tanto el polinomio de interpolación de Newton con diferencias finitas divididas es: f(x) = 4 + 2(x + 2) – 0.25(x + 2)(x + 1) – 0.3(x + 2)(x + 1) (x – 2) Ejemplo 2. Calculemos la tabla de diferencias divididas finitas y usemos la información en la misma, para construir el polinomio de interpolación de Newton con los siguientes datos: x y

–3 5

–2 8

0 4

4 2

Solución. Con los datos dados, procedemos a construir la tabla de diferencias divididas como sigue:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

76

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

Luego, el polinomio de interpolación de Newton tiene la siguiente forma: f(x) = 5 + 3(x + 3) – 1.66667(x + 3)(x + 2) – 0.20238(x + 3)(x + 2) (x) TEOREMA. Si x0, x1, …, xn son números reales distintos, entonces para valores arbitrarios y0, y1, … yn existe un única polinomio de Newton fn(x) de máximo grado n, tal que: fn (xi) = yi ,

: i = 0, 1, 2, …, n

DEMOSTRACIÓN. No se probará formalmente la existencia de un polinomio de interpolación, aunque se acepta que dada cualquier tabla de datos, el polinomio de Newton siempre existe. Sin embargo, se probará la unicidad del polinomio de interpolación. Supóngase que gn(x) es otro polinomio de interpolación de máximo grado n, y sea: hn(x) = fn (x) – gn(x). Luego: hn(xi) = fn (xi) – gn(xi) = yi – yi = 0, para todo i = 1, 2, …, n. Por lo tanto, hn(x) tiene n + 1 raíces distintas, y es un polinomio de grado a lo más n, esto solamente es posible si hn(x) = 0. Luego:

fn (x) = gn(x)

Que es lo que se quería probar. Sin embargo, aunque el polinomio de interpolación es único, pueden existir diversas formas de encontrarlo: Una, es mediante el polinomio de Newton, otra mediante el polinomio de Lagrange. 4.4 POLINOMIO DE INTERPOLACIÓN DE LAGRANGE Consideremos los siguientes pares ordenados: x y

x0 y0

x1 y1

… …

xn yn

El polinomio de interpolación de Lagrange se define como sigue: n

P (x) = y0l0(x) + y1l1(x) + … + ynln(x) =

y jl j x j 0

Donde los polinomios lj(x), para todo j = 0, 1, 2, …, n, se denominan funciones cardinales del polinomio de Lagrange, correspondientes a los pares dados. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

77

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Un análisis lógico de esta expresión sugiere que, si x = x0, entonces debe satisfacerse que P(x0) = y0 , lo cual se cumple si l0(x0) = 1 y también si li(x0) = 0 para toda i 0. Del mismo modo, si x = x1 , entonces debe satisfacerse que: P(x1) = y1 , lo cual se cumple si l1(x1) = 1 y si li(x1) = 0 , para toda i 1. Y así sucesivamente, si x = xn , entonces debe cumplirse que P(xn) = yn , lo cual es evidente si ln(xn) = 1 y si li(xn) = 0 para toda i n. Este análisis lógico sugiere como determinar las funciones cardinales del polinomio de Lagrange. Para ser más claros, analicemos detenidamente la lógica para determinar la función l0(x0). De acuerdo al análisis anterior es evidente que la función cardinal l0(x) debe cumplir con las siguientes condiciones: o o

Primera condición: l0(x0) = 1. Segunda condición: l0(xj) = 0, para todo j

0.

Por lo tanto, para que se cumplan ambas condiciones, el polinomio l0(x) debe tener la forma que sigue: l0(x) = c (x – x1)( x – x2) … (x – xn) El valor de la constante c se determina aplicando la primera condición a la anterior expresión, es decir: l0(x0) = 1 = c (x0 – x1)( x0 – x2) … (x0 – xn) Luego, resulta:

c=

1 x1 x0 x 2 ... x0

x0

xn

Por lo tanto el polinomio l0(x) queda definido como: x - x1 x - x 2 ... x - x n x0 - x1 x0 - x 2 ... x0 - x n

l0(x) = Análogamente se podemos deducir: n

x

xi

i 0 i j

xj

xi

lj(x) =

, con j = 1, 2, …, n

Para el conjunto de nodos x0, x1, x2, …, xn. Ejemplo 1: Construyamos las funciones cardinales y el polinomio de interpolación de Lagrange correspondiente. para los siguientes pares: x y

1 –2

3 1

5 2

7 –3

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

78

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Solución. Se tiene que: P (x) = y0l0(x) + y1l1(x) + y2l2(x) + y3l3(x) P (x) = – 2l0(x) + l1(x) + 2l2(x) – 3l3(x) x 3 x 5 x 7 x 3 x 5 x 7 = 2 4 6 48 x 1 x 5 x 7 x 1 x 5 x 7 l1(x) = = 2 2 4 16 x 1 x 3 x 7 x 1 x 3 x 7 l2(x) = = 4 2 2 16 x 1 x 3 x 5 x 1 x 3 x 5 l3(x) = = 6 4 2 48 Sustituyendo en el polinomio de Lagrange definido se tiene:

Donde:

l0(x) =

x 3 x 5 x 7 x 1 x 5 x 7 x 1 x 3 x 7 + +2 – 48 16 16 x 1 x 3 x 5 3 48 x 3 x 5 x 7 x 1 x 5 x 7 x 1 x 3 x 7 P (x) = + – – 24 16 8 x 1 x 3 x 5 16

P (x) = – 2

Ejemplo 2. Construyamos las funciones cardinales y el polinomio de interpolación de Lagrange correspondiente, usando los siguientes datos: x y Solución. Se tiene:

–2 1

0 –1

2 3

4 –2

P (x) = y0l0(x) + y1l1(x) + y2l2(x) + y3l3(x) P (x) = y0l0(x) – l1(x) + 3l2(x) – 2l3(x)

Donde:

x 0 x 2 x 4 xx = 2 4 6 x 2 x 2 x 4 x 2 l1(x) = = 2 2 4 x 2 x 0 x 4 xx l2(x) = = 4 2 2 x 2 x 0 x 2 x 2 l3(x) = = 6 4 2

l0(x) =

2 x 4 48 x 2 x 4 16 2 x 4 16 x 0 x 2 48

Sustituyendo en el polinomio definido de Lagrange queda: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

79

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

xx 2 x 4 x 2 x 2 x 4 xx 2 x 4 x 2 x 0 x 2 – +3 – 2 48 16 16 48 xx 2 x 4 x 2 x 2 x 4 xx 2 x 4 x 2 x 0 x 2 P (x) = – +3 – 48 16 16 24

P (x) =

Ejemplo 3: Construyamos las funciones cardinales y el polinomio de interpolación de Lagrange correspondiente, para el conjunto de pares dado en la siguiente tabla: x y

–7 –23

5 1

–6 –54

0 –954

Solución: Las funciones cardinales, resultan ser: l0(x) =

5 7 5 6 5 0

7 5

7 6

7 0

x 5 x 7 x 0

l2(x) =

6 5

6 7

6 0

x 5 x 7 x 6 0 5 0 7 0 6

x 7 x 6 x

=

x 5 x 6 x 0

l1(x) =

l3(x) =

x 7 x 6 x 0

=

660

;

x 5 x 6 x

=

84

=

x 5 x 7 x 72

;

;

x 5 x 7 x 6 210

El polinomio de interpolación de Lagrange resulta: P(x) = l0(x) – 23 l1(x) – 54 l2(x) – 954 l3(x) P(x)=

x 7 x 6 x 660

+23

x 5 x 6 x 84

–54

x 5 x 7 x 72

+954

x 5 x 7 x 6 210

4.5 INTERPOLACIÓN DE SPLINES En esta tipo de interpolación – denominada también interpolación segmentaria o interpolación por splines – la idea central radica en que, en vez de usar un solo polinomio para interpolar todos los pares ordenados, se utilizan segmentos de polinomios para cada subintervalo [xi-1, xi], para luego unirlos adecuadamente entre si bajo ciertas condiciones de continuidad y derivabilidad y formar la función de interpolación polinómica suave. De este modo, resulta que, un spline es una curva definida en fragmentos mediante polinomios, cada uno definido entre los dos puntos extremos de cada subintervalo.

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

80

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------En los problemas de interpolación, a menudo se utiliza la interpolación splines porque da lugar a resultados similares requiriendo solamente el uso de polinomios de bajo grado a la vez que se evitan las oscilaciones, que en la mayoría de las aplicaciones resultan indeseables, que aparecen al interpolar mediante polinomios de grado elevado. 4.5.1. Definición de un Splines de grado k. Para una partición en n subintervalos del intervalo [x0, xn] y los correspondientes valores yi (con i = 0, 1, 2, …, n) que se presentan en una tabla de datos de la siguiente forma: x y

xo yo

… …

x1 y1

xn yn

Y, dado un número entero positivo k, se define como una función de interpolación spline de grado k, para el conjunto de n + 1 pares ordenados dados, a aquella función S(x) compuesta de n polinomios sj(x) tal que: s(xi) = yi, i = 0, 1, 2, …, n. sj(x), (con j = 1, 2, 3, … n) es un polinomio de grado m k en cada subintervalo [xj-1, xj]. iii ) s(x), es derivable, al menos, hasta el orden k-1 sobre el intervalo [x0, xn]. i) ii)

4.6 FUNCIONES SPLINES DE GRADO CERO Los splines de grado 0 son funciones constantes por subintervalos. Una forma explícita de presentar un spline de grado 0 es la siguiente:

S(x) =

s1 x

c1

x

x0 ,x1

s2 x

c2

x

x1 ,x2  xn 1 ,xn



 sn x

cn

x

Y, su gráfica tiene la forma semejante a la que se muestra en la siguiente figura:

Figura 4.5. Observación: Los polinomios de grado cero no cumplen con la condiciones de continuidad y derivabilidad en los puntos extremos que definen cada subintervalo [xi, xi+1] (con i = 0, 1, 2, 3, …, n)

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

81

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------4.7. FUNCIONES SPLINES DE GRADO UNO Para los “n + 1” pares ordenados: x y

xo yo

… …

x1 y1

xn yn

Una función spline de grado 1 que interpole estos datos se define al unir cada uno de los puntos dados mediante segmentos de recta y, su expresión analítica tiene la forma:

S(x) =

s1 x

a1 x b1

x

x0 ,x1

s2 x

a2 x b2 

x

x1 ,x2  xn 1 ,xn

sn x



an x bn

x

La misma que cumple con las condiciones que definen la spline de grado 1, es decir: 1. sj(x), con j = 1, 2, 3, …, n, es un polinomio de grado m = 1 2. sj(x) tiene, al menos, una derivada de orden k – 1 = 0. 3. s(xi) = yi, i = 0, 1, 2, …, n. Y, su gráfica tiene asume una forma semejante a la que se muestra en la siguiente figura 4.6

Figura 4.6. Observación: Los polinomios de grado uno cumplen con la condición de continuidad, pero no cumplen con la condición de derivabilidad en los puntos extremos que definen cada subintervalo [xi, xi+1] (con i = 0, 1, 2, 3, …, n) 4.8 FUNCIONES SPLINES DE GRADO 2 En este caso, los polinomios sj(x), a través de los que construimos el spline tienen grado 2. Esto significa que tendrán la forma: sj(x) = ajx² + bjx + cj La interpolación cuadrática nos va a asegurar que la función que generemos a fragmentos con los distintos sj(x) va a ser continua y derivable, ya que para obtener las expresiones que ajusten el polinomio y definan el spline, tenemos como condiciones: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

82

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------La condición de continuidad: Que establece que los polinomios sj(x) y sj-1(x) pasen por un mismo punto xi, es decir, que las valores sj-1(xi) sj(xi) deben ser iguales en cada uno de los puntos xi. La condición de derivabilidad: Que define que las derivadas s’j-1(xi) y s’j(xi) sean iguales, es decir la función s(x) debe ser derivable en xi. Dado que estas condiciones son las que nos permitirán determinar los parámetros aj, bj, y cj (con j = 1, 2, 3, …, n) para los n polinomios sj, las mismas, no son suficientes, y necesitamos alguna otra condición más; en razón a que, por cada spline sj(x) que pasa por dos puntos, de cada subintervalo [xi-1, xi] del intervalo [x0, xn], tenemos 3 incógnitas por lo que resultan 3n incógnitas, teniendo por otro lado que, por la condición de continuidad, se planten dos ecuaciones por cada polinomio sj(x) obteniéndose 2n ecuaciones para los n polinomios sj(x) y, por la condición de derivabilidad, podemos obtener una ecuación por cada punto de contacto entre dos polinomios consecutivos, con lo que obtenemos n – 1 ecuaciones, lo que, en total, hacen 3n – 1 ecuaciones, necesitándose un sistema de 3n ecuaciones, por lo que, en este caso, el sistema tiene un grado de libertad. Para tener las 3n, necesitamos una ecuación más, la que puede obtenerse asignando un valor a la derivada s’j(x) en algún punto, forzándose a uno de los sj(x), o alternativamente, suministrando un valor a alguna de las incógnitas. Resolvamos un ejemplo concreto para obtener una mejor comprensión de esta lógica. Ejemplo 1. Consideremos los siguientes datos: xi yi

3 2.5

4.5 1

7 2.5

9 0.5

Con los datos dados, se forman tres subintervalos: [3, 4.5], [4.5, 7] y [7, 9], por cada uno de los cuales debemos definir una función polinomial de grado 2, que definirán el spline S(x) de la forma siguiente: s1 x

a1 x 2

b1 x c1

x

3, 4.5

S(x) = s2 x

a2 x 2

b2 x c2

x

4.5, 7

s3 x

a3 x 2

b3 x c3

x

7 ,9

Lo que nos indica que tenemos 9 incógnitas, cuya solución nos obliga a determinar 9 ecuaciones. Para cumplir con la condición de continuidad, debemos hacer que cada polinomio pase por los puntos dados de la tabla de datos. Es decir: s1(3) = 2.5,

s1(4.5) = 1,

s2(4.5) = 1,

s2(7) = 2.5,

s3(9) = 2.5,

s3(9) = 0.5

De este modo, conformamos las siguientes seis ecuaciones: s(3) = 2.5 = 32a1 + 3b1 + c1

(1)

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

83

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------s(4.5) = 1 =

4.5 2 a1 4.5b1 4.5 2 a 2 4.5b2 7 2 a2

2

c1 c2

3

c2

4

7 a3 7b3 c3 s(9) = 0.5 = 92a3 + 9b3 + c3

5

s(7) = 2.5 =

7b2

2

(6)

Ahora nos restan determinar 3 ecuaciones más, las que debemos deducir a partir de la condición de derivabilidad de la función spline. Los polinomios de grado 2, deben tener derivadas de orden k – 1 = 1, lo que implica una derivada por polinomio, es decir: 2a1 x b1

x

3,4.5

s´(x) = 2a 2 x b2 2a3 x b3

x

4.5,7

x

7 ,9

La condición de derivabilidad de una función en un punto estipula que las derivadas laterales en dicho punto deben ser iguales, es decir que, en los puntos de abscisas x = 4.5 y x = 7, se debe cumplir: 2a1(4.5) + b1 = 2a2(4.5) + b2 2a2(7) + b2 = 2a3(7) + b3

9a1 + b1 = 9a2 + b2 14a2 + b2 = 14a3 + b3

(7) (8)

Con lo que obtenemos 8 ecuaciones para las 9 incógnitas; lo que nos da un grado de libertad para elegir alguna condición que podríamos aplicarla sobre la primera derivada o asignando un valor para alguna de las incógnitas. Elegimos por simple conveniencia la segunda opción, haciendo: a1 = 0. De este modo eliminamos una incógnita y ahora tenemos un sistema de 8 ecuaciones con 8 incógnitas que resulta: 3b1 + c1 = 2.5 4.5b1 + c1 = 1 20.25a2 + 4.5b2 + c2 = 1 49a2 + 7b2 + c2 = 2.5 49a3 + 7b3 + c3 = 2.5 81a3 + 9b3 + c3 = 0.5 b1 – 9a2 – b2 = 0 14a2 + b2 – 14a3 – b3= 0 Cuya forma matricial se expresa de la siguiente manera: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

84

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

1

0

0

0

0

0

0

b1

2 .5

4.5 1

0

0

0

0

0

0

1

0

0 20.25 4.5 1

0

0

0

0

0

49

7

1

0

0

0

c1 a2 b2

0

0

0

0

0

49

7

1

0

0

0

0

0

81

9

1

1

0

9

1 0

0

0

0

c2 a3 b3

0

0

14

1 0

c3

3

1

0

14

1 =

2 .5 2 .5 0 .5 0 0

De cuya resolución obtenemos los siguientes valores: b1 = – 1, c1 = 5.5, b2 = – 6.75, c2 = 18.46, b3 = 24.6, c3 = – 91.3.

a1 = 0, a2 = 0.64, a3 = – 1.6,

Sustituimos estos valores, en la función s(x), obtenemos la función spline cuadrática que interpola la tabla de datos dada: x 5.5c1

s(x) = 0.64 x 2 6.75 x 18.46 1.6 x 2 24.6 x 91.3

x

3,4.5

x

4.5,7

x

7 ,9

La figura 4.7, contiene tanto los puntos iniciales de la tabla de datos, así como la spline cuadrática determinada.

Figura 4.7. Ejemplo 2. Este ejemplo es para que lo resuelvas tú: Utilizando un spline de segundo grado interpola los siguientes datos: x y

–1 –1

1 1

2 5

4 –2

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

85

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------4.9 FUNCIONES SPLINES CÚBICOS El spline cúbico (k = 3) es el spline más empleado, debido a que proporciona un excelente ajuste a los puntos tabulados y su cálculo no es excesivamente complejo. Si tenemos n + 1 puntos definidos por los datos de la tabla: xi yi

x0 y0

x1 y1

x2 y2

xn yn

Cada uno de los n subintervalos [x0, x1], [x1, x2], …, [xn-1, xn] contiene un polinomio Sj (con j = 1, 2, 3, …, n) de tercer grado de la forma: Sj(x)= ajx3 + bjx2 + cjx + dj Por lo que el spline S(x) queda definido por n polinomios diferentes de tercer grado Sj(x), cuya forma toma la representación siguiente:

S(x) =

s1 x

a1 x 3 + b1 x 2 + c1 x + d1

x

x0 ,x1

s2 x

a2 x3 + b2 x 2 + c2 x + d 2 

x

x1 ,x2

sn x





an x3 + bn x 2 + cn x + d n

x

xn 1 ,xn

Donde aj, bj, cj y dj (con j = 1, 2, 3, …, n) son los 4 parámetros que hay que determinar para cada uno de los n subintervalos, es decir que, en todo el spline tenemos 4n parámetros desconocidos, los que deberán calcularse sobre la base de las condiciones de continuidad y derivabilidad para la primer y segunda derivadas, que deben generar también 4n ecuaciones para que la spline quede perfectamente definida como una curva suave. Condición de continuidad: Cada uno de los polinomios Sj (con j = 1, 2, 3, …, n) contiene a los límites de sus respectivos intervalos [xi-1, xi] con i = j, cumpliéndose en consecuencia que: Para el intervalo [x0, x1]:

y0 = S1(x0), y1 = S1 (x1)

Para el intervalo [x1, x2]:

y1 = S2(x1) y2= S2(x2)

Y así sucesivamente, para el intervalo [xn-1, xn]:

yn-1 = Sn(xn-1) yn= Sn(xn)

De este modo obtenemos dos ecuaciones por cada subintervalo, creándose por esta condición, un total de 2n ecuaciones

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

86

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Condición de derivabilidad para la primera y segunda derivada: Para cada Sj(x) resultan las siguientes derivadas: Sj'(x) = 3ajx2 + 2bjx + cj Sj''(x) = 6ajx + 2bj Las que, en cada punto de contacto de cada Sj-1(x) y Sj(x), deben cumplir con la condición de existencia de la primera deriva, generándose n – 1 ecuaciones de la forma: S’ j(xi) = S’j(xi) con i = j = 1, 2, 3, …, (n – 1) Y con la condición de existencia de la segunda deriva, generándose n – 1 ecuaciones de la forma: S”j(xi) = S”j(xi) con i = j = 1, 2, 3, …, (n – 1) Con lo que resulta que el número de ecuaciones que se forman por estas condiciones es: (2n) + (n +1) + (n + 1) = 4n – 2 Con lo que construimos un sistema de 4n ecuaciones lineales con 4n – 2 incógnitas, es decir, se tiene un sistema con cos grados de libertad, los que definimos con las siguientes condiciones: S”1(x0) = 0

S”n(xn) = 0

y

Completándose de esta manera el sistema de 4n ecuaciones con 4n incógnitas, el mismo que puede ser resuelto por cualquiera de los métodos desarrollados en la unidad didáctica anterior. Ejemplo 1: Determinemos la spline de grado 3 para los datos de la siguiente tabla: –1 –1

x y

1 1

2 5

4 –2

Solución.- Como los datos dados presentan tres subintervalos, entonces, el número de polinomios cúbicos del spline es n = 3, por lo que el número de incógnitas será: 4

n=4

3 = 12

Luego, el spline cúbico estará definido de la siguiente manera:

(A)

a1 x3

b1 x 2

c1 x d1

x

s(x) = a2 x3

b2 x 2

c2 x d 2

x

1; 2

a3 x3

b3 x 2

c3 x d3

x

2; 4

1;1

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

87

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Y su primera y segunda derivada:

(B)

3a1 x 2

2b1 x c1

x

2

2b2 x c 2 2b3 x c3

x

1,2

x

2,4

S´(x) = 3a 2 x 3a3 x 2

1,1

6a1 x 2b1 x 1,1 x 1,2 (C) S”(x) = 6a 2 x 2b2 6a3 x 2b3 x 2 ,4 i) Por la condición de continuidad, partimos de la expresión (A) y obtendremos:

2n=2

3 = 6 ecuaciones, es decir:

s(– 1) = – a1 + b1 – c1 + d1 = – 1

… (1)

s(1) = a1 + b1 + c1 + d1 = 1

… (2)

s(1) = a2 + b2 + c2 + d2 = 1

… (3)

s(2) = 8a2 + 4b2 + 2c2 + d2 = 5 … (4) s(2) = 8a3 + 4b3 + 2c3 + d3 = 5 … (5) s(4) = 64a3 + 16b3 + 4c3 + d3 = – 2 … (6) ii) Por la condición de derivabilidad para la primera derivada, serán: n – 1 = 3 – 1 2 ecuaciones Para que exista S´(x) en x = 1 y en x = 2, en la expresión (B) debemos igualar las ecuaciones correspondientes en ambos puntos, obteniendo: 3a1 + 2b1 + c1 – 3a2 – 2b2 – c2 = 0

… (7)

12a2 + 4b2 + c2 – 12a3 – 4b3 – c3 = 0 … (8) De igual manera, para que exista S”(x) en x = 1 y en x = 2, en la expresión (C) debemos igualar las respectivas ecuaciones en ambos puntos, obteniendo: 6a1 + 2b1 = 6a2 + 2b2 2a2 + 2b2 = 12a3 + 2b3

3a1 + b1 – 3a2 – b2 = 0 1 6a2 + b2 – 6a3 – b3 = 0

… (9) … (10)

Finalmente, agregamos las condiciones para que la segunda derivada se anule en los puntos inicial y final del intervalo [–1, 4]; obteniendo: s”( –1) = – 6a1 + 2b1 = 0 s”( 4) = 24a3 + 2b3 = 0

– 3a1 + b1 = 0 … (11) 12a3 + b3 = 0

… (12)

Con lo que completamos el sistema de doce ecuaciones con doce incógnitas, cuya forma matricial es la siguiente:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

88

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

1 1 1 1

1 1 1 1

0 0

0 0

0 0

0 0

0 0

0 0

0 0

0 0

a1 b1

1 1

0

0

0

0

1

1

1

1

0

0

0

0

c1

1

0

0

0

0

8

4

2

1

0

0

0

0

d1

5

0 0

0 0

0 0

0 0

0 0

0 0

0 0

0 0

8 64

4 16

2 4

1 1

a2

5 2

3

2

1

0

1 0

0

0

0

0

c2

0

0

0

0 12

4

1

0

12

4

1 0

d2

0

3 0

1 0

0 0

0 0

3 6

1 1

0 0

0 0

0 6

0 1

0 0

0 0

a3

0 0

3 1 0 0

0 0

0 0

0 0

0 0

0 0

0 0

0 12

0 1

0 0

0 0

c3

3

2

b2

b3 d3

=

0

0 0

Resolviendo esta ecuación matricial, obtenemos la siguiente solución: 51 , 140 153 , b1 140 89 , c1 140 153 , d1 140 a1

a2 b2 c2 d2

21 24 , a3 , 10 35 297 288 , b3 , 35 35 473 1867 , c3 , 70 70 48 732 , . d3 35 35

Figura 4.8 Por lo tanto, la spline cúbica tiene la siguiente expresión analítica:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

89

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

51 3 x 140 21 3 s(x) = x 10 24 3 x 35

153 2 89 x x 140 140 297 2 473 x x 35 70 288 2 1867 x x 35 70

153 140 48 35 732 35

x

1,1

x

1, 2

x

2, 4

Cuya gráfica correspondiente se muestra en la figura 4.8 Ejemplo 2. Interpolemos los siguientes datos mediante una spline cúbica: x y

2 –1

3 2

5 –7

Solución. Se tendrá un polinomio cúbico en cada uno de los intervalos que se forman: a1 x 3 b1 x 2 s (x) = a 2 x 3 b2 x 2

c1 x d1 c2 x d 2

x x

2 ,3 3,5

Como tenemos dos subintervalos (n = 2), entonces tendremos que determinar (4n) = 8 ecuaciones para determinar las 8 incógnitas a1, b1, c1, d1, a2, b2, c2 y d2 i) Por la condición de continuidad se tiene: s(2) = 8a1 + 4b1 + 2c1 + d1= – 1

(1)

s(3) = 27a1 + 9b1 + 3c1 + d1 = 2

(2)

s(3) = 27a2 + 9b2 +3c2 +d2 = 2 s(5) = 125a2 + 25b2 +5c2 +d2 = – 7

(3) (4)

ii) Ahora trabajamos con la condición de derivabilidad, por lo que calculamos la primera derivada s´(x): s´(x) =

3a1 x 2 3a 2 x 2

2b1 x c1 2b2 x c 2

x x

2 ,3 3,5

En este caso, el punto donde se cambia de intervalo es x = 3. Entonces, para que exista s´(x), en x = 3 se reemplaza este valor en los dos polinomios anteriores y luego se los iguala, es decir: 3a1(3)2+2b1(3)+c1=3a2(3)2+2b2(3)+c2

27a1+6b1+c1=27a2+6b2+c2

(5)

iii) De igual modo procedemos con la segunda derivada: s”(x) =

6a1 x 2b1

x

2,3

6a2 x 2b2

x

3,5

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

90

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------En este caso, para lograr que s”(x) exista en x = 3, se debe cumplir: 6a1(3)+2b1=6a2(3)+2b2

18a1 + 2b1 = 18a2 + 2b2

(6)

Para los 2 grados de libertad que, tenemos las siguientes 2 condiciones: s”(x0) = 0

y

s”(xn) = 0

De este modo obtenemos: s”(2) = 6a1(2) + 2b1 = 0 s”(5) = 6a2(5) + 2b2 = 0

12a1 + 2b1 = 0 30a2 + 2b2 = 0

(7) (8)

Con lo cual, hemos completado el sistema de 8 ecuaciones con 8 incógnitas, el mismo puede escribirse en forma matricial de la manera siguiente:

4 2 1

0

0

0

0

a1

27 9 3 1

0

0

0

0

0

0 0 0

27

9

3

1

b1 c1

0

0 0 0

125

25

5

1

27 6 1 0

27

6

18 2 0 0

18

2

12 2 0 0

0

0

30

8

0 0 0

d1 a2

1 0

1 2 2 =

7 0 0

0

b2 c2

0

d2

0

0

0

0

0

2

0

0

De donde obtenemos la siguiente solución: a1 = – 1.25, b1 = 7.5, a2 = 0.625, b2 = – 9.75,

c1 = – 10.75, d1 =0.5, c2 = 39.875, d2 = –50.125,

Figura 4.9 Sustituyendo estos valores en la expresión inicial (A), obtenemos: s (x) =

1.25x 3 0.625x

3

7.5 x 2 10.75x 0.5

9.375x

2

39.875x 50.125

x

2,3

x

3,5

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

91

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Cuya gráfica correspondiente es la que se muestra en la figura 4.9. Ejercicio 3 Dado que ln(2) = 0.69315, ln(3) = 1.0986 y ln(6) = 1.7918, interpolemos con un polinomio spline cúbico el logaritmo natural de cada entero desde 2 hasta 7 y tabulemos los valores de s(x), los valores verdaderos de ln(x) y el error en cada punto. Solución.- Organizamos la tabla de valores: xi f(xi)

2 0.69315

3 1.09860

6 1.79180

Tenemos n = 2 subintervalos, por lo que tendremos que plantear un sistema de 8 ecuaciones para las 8 incógnitas que resultan. Resultando un spline de la forma:

s (x) =

a1 x3

b1 x 2

c1 x d1

x

2,3

a2 x3

b2 x 2

c2 x d2

x

3, 6

i) Por la condición de continuidad tenemos las siguientes ecuaciones: s(2) = 8a1 + 4b1 + 2c1 + d1= 0.69315

(1)

s(3) = 27a1 + 9b1 + 3c1 + d1 = 1.09860

(2)

s(3) = 27a2 + 9b2 +3c2 +d2 = 1.09860 s(6) = 216a2 + 36b2 +6c2 +d2 = 1.79180

(3) (4)

ii) Con la condición de derivabilidad para la primera derivada s´(x), tenemos:

s´(x) =

3a1 x 2

2b1 x c1

x

2,3

3a2 x 2

2b2 x c2

x

3,6

En este caso, el punto donde se cambia de intervalo es x = 3. Entonces, para que exista s´(x), en x = 3 se reemplaza este valor en los dos polinomios anteriores y luego se los iguala, es decir: 3a1(3)2+2b1(3)+c1=3a2(3)2+2b2(3)+c2

27a1+6b1+c1–27a2–6b2–c2 = 0

(5)

iii) Del mismo modo para la segunda derivada:

6a1 x 2b1 x 2,3 6a2 x 2b2 x 3,5 En este caso, para lograr que s”(x) exista en x = 3, se debe cumplir: s”(x) =

6a1(3)+2b1=6a2(3)+2b2

9a1 + b1 – 9a2 –b2 = 0

(6)

Para los 2 grados de libertad que, tenemos las siguientes 2 condiciones:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

92

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------s”(x0) = 0

y

s”(xn) = 0

De este modo obtenemos: s”(2) = 6a1(2) + 2b1 = 0 s”(6) = 6a2(6) + 2b2 = 0

6a1 + b1 = 0 18a2 + b2 = 0

(7) (8)

Con lo cual, hemos completado el sistema de 8 ecuaciones con 8 incógnitas, el mismo puede escribirse en forma matricial de la manera siguiente:

4 2 1

0

0

0

0

a1

0.69315

27 9 3 1

0

0

0

0

1.09860

0

0 0 0

27

9

3

1

b1 c1

0

0 0 0 216 36

6

1

8

27 6 1 0

27

6

9

1 0 0

9

1

6

1 0 0

0

0

0 0 0

18

1 0

d1 a2

1.09860 =

1.79180 0 0

0

b2 c2

0

d2

0

0

0

0

0

1

0

0

De donde obtenemos la siguiente solución: a1 = – 0.02179, a2 = 0.0072659,

b1 = 0.13078, c1 = 0.16567, d1 =0.01303, b2 = – 0.1307, c2 = 0.95039, d2 = –0.7716

Finalmente, el spline queda definido en la siguiente expresión analítica

s (x) =

0.02179 x3

0.13078 x 2

0.0072659 x3

0.1307 x 2

0.16567 x 0.01303

x

2,3

095039 x 0.7716

x

3,6

Luego, el cálculo de los valores de interpolación con un polinomio spline cúbico de cada entero desde 2 hasta 7 se presentan junto a los valores verdaderos de ln(x) y el error en cada punto en la siguiente tabla: x 2 3 4 5 6 7

s(x) 0,69317 1,09873 1,4037776 1,6210875 1,7949744 1,9690337

ln(x) 0,6931472 1,0986123 1,3862944 1,6094379 1,7917595 1,9459101

0,003292149 0,010714547 1,261149102 0,723829573 0,179428703 1,18831545

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

93

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------4.10 PROBLEMAS NOTA: Cuando sea necesario, redondear a cinco decimales. 1. Calcular el polinomio de interpolación de Newton para los siguientes datos: x y x y

–2 –3

2 0.5

0.3 –3

0.6 0

1 2.4

4 7.8

0.9 –6

1.2 9

1.5 –12

2. Calcular el polinomio de Lagrange para los siguientes datos: x y

1 1.56

–2 3.54

–1.5 –0.5 9 –2

x y

3 –2.57 –2 33

1 5

–5 –8.9 –4 0

3. Dado que ln(2) = 0.69315, ln(3) = 1.0986 y ln(6) = 1.7918, interpola con un polinomio de Lagrange el logaritmo natural de cada entero en el intervalo [1, 10]. Tabula lo anterior junto con el error en cada punto. 4. Calcular las splines cúbicas para los siguientes datos: x y x y

–5 20

–2 40 –2 4

1 –5

3 –20 3 –6

7 40

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

94

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------UNIDAD DIDÁCTICA 5 INTEGRACIÓN NUMÉRICA 5.1 INTRODUCCIÓN El objetivo de la presente unidad didáctica se orienta a la solución del problema del cálculo de un área subtendida por la gráfica de una función (continua o discreta) en un intervalo, definido por dos límites conocidos. Específicamente se buscamos: 1. Comprender las bases conceptuales de la integración aproximada. 2. Comprender los procedimientos generales de la integración aproximada. b

En general, la evaluación de una integral definida de la forma

f x dx , utilizando a

métodos exactos, es con frecuencia difícil o imposible. Por tanto, se hacen necesarios otros métodos de solución, tanto para esos casos como para el problema más general de la integración, en el que solamente se conocen valores discretos de f(x) para puntos de abscisas conocidas x0, x1, x2, , xn; una evidente alternativa, en este caso, consiste en encontrar una función polinómica p(x) que interpole dichos puntos y, para el primer caso, sea una adecuada aproximación de f(x) que, además, sea fácil de integrar analíticamente. Es decir, se asume: b

b

f x dx a

p x dx a

Debemos recordar que los polinomios de interpolación p(x) desarrollados en la unidad didáctica anterior, generan aproximaciones adecuadas y son fáciles de integrar. Esta combinación de características es una razón más por la cual estudiamos estos polinomios dentro del análisis numérico. 5.2 TEOREMA FUNDAMENTAL DEL CÁLCULO Si f(x) es una función continua en el intervalo [a, b] y si F(x) es la antiderivada de f(x). El teorema fundamental del cálculo integral se expresa en la siguiente expresión: b

f x dx = F(b) – F(a)

F ’(x) = f (x)

a b

Cuya propiedad básica es:

c

a

b

f ( x)dx +

f x dx = a

f ( x)dx

c [a, b]

c

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

95

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------En la práctica, el problema que el análisis numérico resuelve, es cuando se tiene que evaluar la antiderivada de una función que es imposible de determinar con el teorema fundamental del cálculo, aún para integrales aparentemente sencillas como la siguiente: 1

2

e x dx 0

En esta unidad didáctica estudiaremos algunos métodos numéricos que nos permitirán obtener aproximaciones bastante exactas a integrales imposibles de resolver con el teorema fundamental, especialmente se verán dos tipos de integración numérica: las fórmulas de Newton-Cotes y el algoritmo de Romberg. 5.3 FORMULAS DE INTEGRACIÓN DE NEWTON-COTES Las fórmulas de Newton-Cotes están conformadas por las conocidas reglas del trapecio y de Simpson (regla de un tercio y de tres octavos). En general, estas fórmulas se basan en la integración de una función polinómica pn(x) como una aproximación de la función original f(x) en el intervalo [a, b], es decir: b

b

f x dx a

pn x dx a

Donde pn (x) = a0 + a1x + a2x2 +…+ anxn es un polinomio de grado que interpola adecuadamente ciertos puntos de f(x) escogidos apropiadamente o que se presentan en una tabla de datos, los que se pueden ajustar a dicho polinomio de interpolación. Dentro de las fórmulas de Newton-Cotes, existen las formas cerradas y abiertas. En las formas cerradas se conocen los valores de f(a) y f(b), en el caso opuesto, se llaman formas abiertas. En la presente unidad didáctica estudiaremos únicamente las formas cerradas, por lo que, siempre se suponen conocidos los valores f (a) y f (b). 5.4 REGLA INTEGRACIÓN DEL TRAPECIO En esta regla, el polinomio de interpolación es un polinomio lineal, es decir: b

b

f x dx a

p1 x dx a

Donde p1 (x) es un polinomio línea que se ajusta a los datos: x y

a b f (a) f (b)

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

96

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Como ya vimos en la unidad didáctica anterior, este polinomio esta dado por la expresión: f (b) f (a) (x – a) b a

p1(x) = f (a) +

Integrando entre los límites a y b, obtenemos: b

a

f (b) f (a) ( x a ) 2 p1 ( x )dx = f (a) x + 2 b a

b

= a

f (b) f (a) (b a) 2 = f (a) (b – a) + = b a 2 = (b – a)

f (a)

f (b)

f (a) 2

=

b a [f(b) + f(b)] 2

Figura 5.1 Finalmente, obtenemos la siguiente expresión: b

f x dx a

b a [f(a) + f(b)] 2

Que define a la Regla del Trapecio. El nombre que recibe obedece a la interpretación geométrica que se da a la fórmula; toda vez que, la integral obtenida, corresponde al valor del área subtendida por la línea recta y el eje de las abscisas en el intervalo [a, b], que es precisamente el área del trapecio que se forma, tal como lo muestra la figura 5.1. 1

2

e x dx

Ejemplo 1: Utilicemos la regla del trapecio para aproximar la integral: 0

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

97

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Solución. Determinamos los datos requeridos para aplicar la regla: 2 f (x) = e x , a = 0 ,

b=1

2 2 f (0) = e 0 = 1, f (1) = e1 = e

1

1 0 f 0 2

2

Aplicando la fórmula tenemos: e x dx 0

=

f 1

1 e = 1.85914 2 4

ex Ejemplo 2. Usemos la regla del trapecio para aproximar la integral: dx x 2 Solución. Igual que en el ejemplo anterior, sustituimos los datos de manera directa en la fórmula del trapecio. En este caso, tenemos:

f (2) =

e2 , 2

4 x

Por lo tanto, se tiene: 2

f (x) =

f (4) =

4 2 f 2 2

e dx x

ex , a = 2, x

b=4

e4 4

f 4

e2 2

=

e4 = 17.3441 4

4 x 2

e dx x

17.3441

5.4.1. Mejoramiento de la aproximación de la Regla de Trapecio La regla del trapecio puede mejorar su aproximación, si se subdivide el intervalo [a, b] en b a n subintervalos, todos de la misma longitud h = . n Sea P = {x0, x1, …, nn}el conjunto de los límites de los subintervalos que se forman con la partición del intervalo [a, b], con x0 = a y xn = b. De acuerdo con las propiedades de la integral definida, podemos escribir: b

a

b

f x dx + … +

f x dx +

f x dx = a

x2

x1

x1

f x dx xn

1

Aplicando la regla del trapecio en cada una de las integrales, obtenemos:

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

98

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------b

x1

f x dx

x0 2

a

[f(x0) +f(x1)] + … +

xn

xn 2

1

[f(xn-1) + f(xn)]

Como todos los subintervalos tienen la misma longitud h, se tiene: b

h [f (x0) + 2 f (x1) + 2 f (x2) +… + 2 f (xn-1) + f (xn)] 2

f x dx a

Usando la notación de la sumatoria, finalmente se tiene: b

f x dx a

n 1 b a [f (a) + 2 f xi + f (b)] 2n i 1

Esta es la regla del trapecio para n subintervalos. Es claro que entre más subintervalos se utilicen, mejor será la aproximación a la integral. 1

2

e x dx ,

Ejemplo 1: Apliquemos la regla del trapecio para aproximar la integral 0

subdividiendo [0, 1] en 5 subintervalos. Solución. En este caso tenemos n = 5, y la partición generada es: P = {0, 0.2, 0.4, 0.6, 0.8, 1}, resultando: 2 2 2 f (0) = e 0 = 1, f (0.2) = e 0.2 , f (0.4) = e 0.4 , 2 2 f (0.6) = e 0.6 , f (0.8) = e 0.8 , f (1) = e

Al aplicar la fórmula tenemos: 1

2

e x dx 0

1

2

e x dx 0

1 0 [ f (0) + 2 f (0.2) + 2 f (0.4) +2 f (0.6) + 2 f (0.8) + f (1)] 25

2 2 2 2 1 [1 + 2 e 0.2 + 2 e 0.4 + 2 e 0.6 + 2 e 0.8 + e] = 1.48065 10

Es menester mencionar que el valor verdadero de esta integral es de 1.4626… Así podemos ver que con 5 subintervalos, la aproximación no es tan mala. Para hacer cálculos con más subintervalos, sugerimos elaborar un programa que aplique la fórmula con el número de subintervalos que se desee. ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

99

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------5.5 REGLA DE SIMPSON DE UN TERCIO Supongamos que para una función continua y = f (x), tenemos los datos: a f(a)

xm f(xm)

b f(b)

Siendo xm es el punto medio del intervalo definido por los puntos a y b. Para este caso tenemos: b

b

f x dx a

f 2 x dx a

Donde f2(x) es el polinomio de interpolación de Lagrange para los datos en la tabla anterior. Si aplicamos la formula de este polinomio, obtenemos: f2(x) = f (a)

x a

xm x b x a x xm x a x b + f (xm) + f (b) xm a b b a b xm xm a xm b

Si denotamos con: h =

b a = xm – a = b – xm, entonces, la anterior expresión se puede 2

escribir: f2(x) = f (a)

x

xm x b x a x xm x a x b + f (xm) + f (b) h 2h 2h h h h

Simplificando términos: f2(x) =

f a 2h 2

(x – xm)(x – b) –

f xm f b (x – a)(x – b) + (x – a)(x – xm) h2 2h 2

En esta expresión observamos que cada uno de los términos tienen básicamente la misma forma algebraica, es decir, todos tienen el factor general: (x – )(x – ); de modo que podemos calcular una integral general que tenga la siguiente forma:

x

x

dx

Para ello utilizamos el método de integración por partes, haciendo los siguientes cambios: u=x–

du = dx y, dv = (x – )dx

v=

x

dx =

2

x 2

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

100

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas -----------------------------------------------------------------------------------------

Resultando:

x

dx = (x –

x

x

x

)

dx = (x –

)

2

x

–

2 2

x

–

2

2

x

dx

2 3

x 6

Utilizamos esta fórmula para calcular la integral de cada uno de los tres términos de f2(x), obteniendo: b

f 2 x dx a

f a

b

2h 2 a

b

x xm x b dx

b

Hacemos:

I1 =

x a

I1 = – (a – xm)

f xm x a x b dx h2 a

I2 = a

2h 2 a

a b2 a b3 + =– ( – h) 2 6

2h 2 + 2

a

2h 3 = 6

4 3 2 3 h = h 3 3

a

I3 =

x a x xm dx

a

b x b2 a b3 x b3 – = x a x b dx = (x – a) 2 6 6

b

b

b x b2 x b3 – = xm x b dx = (x – xm) 2 6

I1 = 2h3 – b

f b

2h 3 = 6

4 3 h 3

b x xm 3 x xm 2 – x a x xm dx = (x – a) 6 2

a

b xm 2 b xm 3 a xm 3 I3 = (b – a) – + 2 6 6 I3 = (2h)

h2 h3 – + 2 6

b

Luego tenemos:

f 2 x dx = a

h3 h 3 2h 3 = h3 – = 6 3 3

f a 2h

2

I1 –

f xm f b I2 + I3 2 h 2h 2

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

101

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------b

f 2 x dx = a

f a 2h 2

[

b

h 4h h + f (xm) + f (b) 3 3 3

f 2 x dx = f (a) a b

h f ( a ) 4 f ( xm ) 3

f 2 ( x )dx

Y, finalmente tenemos:

f b 2h 3 4 3 h]+ [ ] 3 3 2h 2

f xm 2 3 h]– [ 3 h2

a

f (b)

Debido al factor h/3 este método recibe el nombre de la regla de Simpson de un tercio. En b a la práctica, se sustituye el valor de h = para obtener la siguiente fórmula final: 2 b

f x dx a

b a [f (a) + 4f (xm) + f (b)] 6

Ejemplo 1. Usemos la regla de Simpson de 1/3 para aproximar la siguiente integral: 1

2

e x dx 0

Solución. Se aplica directamente la fórmula con los siguientes datos: a = 0,

b = 1,

xm = 0.5,

y

f(x) = e x

2

Utilizamos la fórmula obtenida anteriormente, obtenemos: 1

2

e x dx 0

1 0 6

[ e0

2

+ 4 e 0.5

2

2 + e1 ] = 1.47157

Ejemplo 2. Usemos la regla de Simpson de 1/3, para aproximar la siguiente integral: 4 x 2

e dx x

Solución. Al igual que en el ejercicio anterior, sustituimos los datos adecuadamente: 4 x 2

e dx x

4 2 6

[e

2

2

+4

e3 e4 + 3 4

] = 1.47082

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

102

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------Del mismo modo que con la regla del trapecio, se puede mejorar la regla de Simpson de b a 1/3, si se parte el intervalo [a, b] en n subintervalos de la misma longitud h = n Sea P = {x0, x1, …, xn} el conjunto de los límites de los subintervalos de la partición que se aplica al intervalo [a, b] y, xmi el punto medio en cada subintervalo [xi-1, xi]. Ahora, aplicamos la propiedad básica de la integral definida obtenemos: x1

b

f x dx + … +

f x dx +

f x dx = x0

a

xn

x2 x1

f x dx xn 1

Luego aplicamos la regla de Simpson de 1/3, en cada una de las integrales anteriores, obteniendo: b

x1

f x dx a

+

xn

Sustituimos x0 = a, xn = b y h =

x0 6

[f (x0) + 4f ( xm1 )

+ f (x1)] + … +

xn 1 6

[f (xn-1) + 4f ( xmn )

+ f (xn)]

b a y usando la notación de sumatoria, podemos n

escribir: b

f x dx a

n

b a f a 6n

4

n 1

f xmi i 1

2

f xi

f b

i 1

Ejemplo 1. Apliquemos la regla de Simpson de 1/3 y subdividamos en 5 intervalos, para 1

2

aproximar la siguiente integral: e x dx 0

Solución. Como n = 5, la partición que se genera es: P = {0, 0.2, 0.4, 0.6, 0.8, 1} Y los puntos medios de cada subintervalo son: Pm = {0.1, 0.3, 0.5, 0.7, 0.9} Ahora, sustituimos los datos en la anterior fórmula y obtenemos: 1

2

e x dx 0

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

103

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------1 0 f 0 6 5

1

2

e x dx

4 f 0.1

1 4 e

f 0.3

0.1 2

e

...

0.3 2

f 0.9

... e

2 f 0.1

0.9 2

2 e

0.1 2

f 0.4

e

0.4 2

...

... e

f 0.8

0.8 2

f 1

e

30

0

1

2

e x dx

1.4626

0

Notamos que esta aproximación ya es más exacta hasta en el cuarto decimal. Ejemplo 2. Utilizando la regla de Simpson de 1/3 y subdividiendo en 4 subintervalos, 4 x

aproximemos la siguiente integral: 2

e dx x

Solución. Como n = 4, y la partición que se genera es: P = {2, 2.5, 3, 3.5, 4} Y, los puntos medios de cada subintervalo son: Pm = {2.25, 2.75, 3.25, 3.75} Sustituimos todos estos datos en la fórmula y obtenemos la siguiente aproximación: 4 x

e dx y x

2

y

1 e2 12 2

4

e 2.25 2.25

e 2.75 2.75

e 3.25 3.25

e 3.75 3.75

e 2.5 2.5

2

e3 3

e 3.5 3.5

e4 = 14.6767 4

5.6 REGLA DE SIMPSON DE TRES OCTAVOS b

Corresponde a este caso n = 3, es decir:

b

f 3 x dx ;donde f3(x) es un polinomio

f x dx a

a

de interpolación de tercer grado para los siguientes datos: a f(a)

p f(p)

q f(q)

b f(b)

Donde p y q son los puntos que dividen en tres partes iguales al intervalo [a, b]. De igual modo que en el caso anterior, usamos el polinomio de interpolación de Lagrange, y el método de integración por partes llegando a la siguiente fórmula: ---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

104

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------b

f x dx a

3 h[f(a) + 3 f(p) + 3 f(q) + f(b)] 8

Siendo h = (b – a)/3. Debido al factor (3/8)h que aparece en ésta formula, es que se le dio el nombre de Regla de Simpson de 3/8. En la práctica, se sustituye el valor de h para obtener: b

b a [f(a) + 3 f(p) + 3 f(q) + f(b)] 8

f x dx a

Ejemplo 1. Utilicemos la regla de Simpson de 3/8 para aproximar la siguiente integral: 4

e x ln xdx 1

Solución. En este caso, tenemos los siguientes datos: a = 1,

p = 2,

q =3,

b=4

y,

f(x) = ex ln x

Los que sustituimos en la fórmula de Simpson 3/8 y obtenemos: 4

e x ln xdx 1 4

e x ln xdx 1

4 1 [f(1) + 3 f(2) + 3 f(3) + f(4)] 8

3 [e ln(1) + 3 e2ln(2) + 3 e3ln(3) + e4ln(4)] = 58.9698 8

Mejoramiento de la aproximación de la regla de Simpson de 3/8. Del mismo modo que en los dos casos anteriores, podemos mejor la aproximación del valor encontrado por este método, si se subdividimos el intervalo [a, b] en n subintervalos de la misma longitud, b a haciendo: h= n Consideremos que x0, x1, …, xn son los límites de los subintervalos de la partición determinada y, dividimos cada uno de estos subintervalos [xi-1, xi], en tres partes iguales, siendo pi y qi los puntos que determinan esta división y aplicamos la regla de Simpson 3/8 en cada uno de los subintervalos, obtenemos: x1

b a

f x dx + … +

f x dx +

f x dx = x0

xn

x2 x1

f x dx xn 1

---------------------------------------------------------------------------------------------------------------Apuntes universitarios: MAT – 362: Análisis NuméricoI

Profesor: Ing. Abel Barroso López

105

Universidad Católica Boliviana San Pablo – Unidad Académica Regional Tarija Departamento de Ingenierías y Ciencias Exactas ----------------------------------------------------------------------------------------3h f ( x0 ) 3 f ( p1 ) 3 f ( q1 ) 8

…+ b

f ( x)dx a

f ( x1 ) +

3h f ( x1 ) 3 f ( p2 ) 3 f ( q2 ) 8

3h f ( xn 1 ) 3 f ( pn ) 3 f ( qn ) 8

3h f ( x0 ) 3 8

f ( x2 ) + …

f ( xn )

n

n 1

f ( pi )

f ( qi )

i 1

2

f ( xi )

f ( xn )

i 1

Sustituyendo, en la anterior expresión: x0 = a, xn = b y h = (b – a)/3n, resulta: b

f ( x)dx a

n

b a f a 8n

3

n 1

f pi

f qi

2

i 1

f xi

f b

i 1

Esta última, es la regla de Simpson de 3/8 para n subintervalos todos de la misma longitud. Ejemplo 2. Aplicando la regla de Simpson de 3/8 y dividiendo en 3 los subintervalos, 4

aproximemos la siguiente integral: e x ln xdx 1

Solución. Como tenemos n = 3, la partición correspondiente resulta: P = {1, 2, 3, 4}. Al considerar los puntos que dividen en tres partes iguales a cada subintervalo, obtenemos los siguientes datos: 1