FUNDAMENTOS DE CLASE
f(b)
a
c
r1
b
f(a)
c=
x1f(x 2 ) − x 2f(x1 ) f(x 2 ) − f(x1 )
Rodrigo A. Ocoró
Métodos Numéricos
Rodrigo A. Ocoró
1
UNIVERSIDAD SANTIAGO DE CALI FACULTAD: PROGRAMA: ASIGNATURA: PRERREQUISITO: PROFESOR: E-mail : SEMESTRE:
INGENIERIAS INGENIERÍA DE SISTEMAS METODOS NUMERICOS CALCULO 2 Rodrigo A. Ocoró
[email protected] QUINTO
OBJETIVO GENERAL: Al finalizar el curso, el estudiante tendrá las capacidades para analizar, calcular y argumentar conceptos matemáticos que le permitan: permitan: 1. Resolver problemas matemáticos matemáticos que no tienen solución por métodos analíticos 2. Construir un asistente matemático en un lenguaje de programación que le permitan realizar cálculos por métodos numéricos. INTRODUCCION El análisis numérico y sus métodos son una reflexión sobre los cursos tradicionales de cálculo, álgebra lineal, ecuaciones diferenciales, etc., desde el punto de vista numérico concretado en una serie de métodos o algoritmos algoritmos con aplicaciones en diferentes áreas de la ingeniería y de las ciencias. Por lo tanto, consiste en procedimientos que resuelven problemas y realizan cálculos puramente aritméticos ya que tiene que ver con el diseño, desarrollo y evaluación de métodos como parte del procesamiento de información donde los datos constituyen las entradas, los métodos o algoritmos los procedimientos y los resultados las salidas. Dado que cada algoritmo implica numerosas operaciones lógicas, aritméticas y en algunos casos graficaciones, el computador es fundamental para el estudio y uso de ellos aunque el curso de métodos numéricos no incluye detalles sobre programación de computadores. Los algoritmos se presentan en pseudocódigos para que el estudiante los programe en el lenguaje que el considere apropiado y presente los trabajos propuestos por el profesor. ESTRUCTURA TEMATICA El curso se desarrollará en 16 secciones seccio nes de clases (teorica-prácticas) así: Modulo 1: APROXIMACIONES Y ERRORES (2 secciones) Objetivo: - Manejar tipos de errores y algunas formas de evitarlos. Contenido detallado: 1.1.introducción
Métodos Numéricos
Rodrigo A. Ocoró
2
1.2.Tipos de errores 1.3.Formas de medir el error 1.4.Cifras significativas 1.5.Reglas de redondeo manual 1.6.Números en la máquina 1.7.Serie de Taylor 1.8.Actividad: desarrollo de ejercicios propuestos y lecturas complementarias. complem entarias. Modulo 2: SOLUCION DE ECUACIONES NO LINEALES (3 secciones) Objetivo: Presentar métodos analíticos para solución de ecuaciones ec uaciones no lineales Contenido detallado: 2.1.Método de bisección 2.2.Método de regula falsi 2.3.Método de la secante 2.4.Método de Newton Raphson 2.5.Método de punto fijo 2.7.Actividad: desarrollo de algoritmos para diversos métodos y el programa de los mismos en el computador aplicándolos a varios ejercicios. Modulo3: ÁLGEBRA LINEAL NUMÉRICA (4 secciones) Objetivo: Solucionar sistemas de ecuaciones lineales mediante la utilización de diferentes métodos. Contenido detallado: 3.1.Fundamentos preliminares 3.2.Solución de sistemas de ecuaciones lineales 3.3.Eliminación Gaussiana 3.4.Efecto de los errores de redondeo 3.5.Método de Gauss-Jordan 3.6.Método de Jacobi 3.7.Método de Gauss-Seidel 3.5.Actividad: Desarrollo de ejercicios propuestos y creación de programas para la solución de ejercicios. Modulo 4: AJUSTE DE CURVAS (3 secciones) Objetivo: Aproximar funciones disponibles en forma discreta con funciones analíticas sencillas. Contenido detallado: 4.1.Introducción 4.2.Ajuste lineal por mínimos cuadrados 4.3.Ajuste exponencial, potencial y logarítmico 4.4.Ajuste lineal múltiple 4.5.Aproximación polinomial con mínimos cuadrados 4.6. Polinomios de Lagrange 4.7.Diferencias divididas
Métodos Numéricos
Rodrigo A. Ocoró
3
4.5.Actividad: desarrollo de ejercicios propuestos y creación de diferentes algoritmos y sus respectivos programas para solución de ejercicios. Modulo 5: INTEGRACION NUMERICA (2 secciones) Objetivo: Aplicar algunos métodos iterativos para el cálculo de integrales. Contenido detallado: 5.1.Métodos de Newton-Cotes 5.2.Regla trapezoidal simple 5.3.Regla trapezoidal compuesta 5.4.Reglas de Simpson 5.5.Cuadratura de Gauss 5.4.Actividad: desarrollo de ejercicios propuestos. Trabajo complementario. Modulo 6. ECUACIONES DIFERENCIALES NUMÉRICAS (2 secciones) Objetivo: Utilizar algunos métodos numéricos en la solución de ecuaciones diferenciales ordinarias. Contenido detallado: 6.1.Método de Euler 6.2.Análisis de error mediante la serie de Taylor 6.3.Método de Euler modificado 6.4.Métodos de Runge-Kutta 6.5. Diseño y aplicación de Softwares de los métodos estudiados.
METODOLOGIA Con el constante apoyo y orientación del profesor, el alumno se irá apropiando del conocimiento necesario para cumplir con los propósitos. Cada sección contemplará desarrollo teórico y practica sobre ejercicios propuestos y verificación de software en el computador además del trabajo extraclase que debe realizar cada estudiante. En tal sentido, se realizarán las siguientes actividades: • Exposición del profesor sobre el tema a tratar • Ejercicios de ejemplificación, motivación, y apoyo para la orientación del estudiante. • Talleres de practicas para controlar el avance en el aprendizaje de los temas. • Desarrollo por parte de los estudiantes de las actividades propuestas • Creación de algoritmos y desarrollo de estos en programas para el computador. • Sustentación de los trabajos realizados por los estudiantes. • Análisis, discusión y evaluación de los trabajos realizados. RECURSOS -Aula de clases -Módulos proporcionados por el docente -Bibliografía especificada -Computador personal -Salas de cómputo de la Universidad
Métodos Numéricos
Rodrigo A. Ocoró
4
-Plataforma virtual Dokeos -Calculadora EVALUACION Primer parcial: Segundo parcial Evaluación escrita: Talleres y software: Evaluación final: Evaluación Trabajo final
30% 15% 15% 30% 10%
BIBLIOGRAFIA •
AKAI, Métodos Numéricos aplicados a la Ingeniería. Ed. Limusa Wiley
•
CHAPRA S.C, CANALES R.P., Métodos numéricos para ingenieros Ed. Mc Graw Hill.
•
CONTE S.D., Análisis numérico. Ed. Mc.Graw Hill.
•
CURTIS F. G.,Análisis numérico. Ed. Alfaomega.
• •
NIEVES A., DOMINGUEZ F., Métodos numéricos aplicados a ingenieros Ed. Cecsa. RINCON L.A., Análisis Numérico. Universidad Nal. Palmira.
•
SCHEID F., Análisis numérico. Colección Shawm.
•
SHOICHIRO N., Métodos numéricos aplicados con Software Ed. Prentice Hall.
•
EN INTERNET: 1. METODOS NUMERICOS I Métodos Numéricos I. Temario: INTRODUCCION AL ANALISIS NUMERICO. Problemas numéricos y algoritmos. Fuentes de error. Errores absolutos y relativos.... URL: www.csi.ull.es/docencia/asignaturas/311.html 2. Laboratorio de Métodos Numéricos Laboratorio de Métodos Numéricos. Programas para calcular: Las raíces de ecuaciones [ f(x)=0 ] por el método de Newton-Raphson. La integral de una función. URL: gaia.fim.utp.ac.pa/Laboratorios/LabMetNum/
Métodos Numéricos
Rodrigo A. Ocoró
5
4. METODOS NUMERICOS II Métodos Numéricos II. Temario: OBSERVACIONES GENERALES. INTERGRACION DE E.D.O. POR EL PROCEDIMIENTO DE LAS SERIES DE POTENCIAS. METODOS DE EULER PARA LA... URL: filemon.csi.ull.es/docencia/asignaturas/331.html 5. Métodos Numéricos de Resolución de E.D.P.- Ingeniería Técnica Informática URL: www.etsimo.uniovi.es/~cesarm/DOCENCIA/mnredp_eutio.html 6. Métodos Numéricos y álgebra lineal Panorama de las aplicaciones de los métodos numéricos. Fundamentos de los métodos numéricos. Secuencias y... URL: webdia.cem.itesm.mx/dia/ac/nbatres/cb851.html ___________________________________________________________ CRONOGRAMA DESARROLLO TEMATICO METODOS NUMERICOS INGENIERIA DE SISTEMAS 5º SEMESTRE TEMAS Modulo 1 Módulo 2 Primer parc. Módulo 3 Modulo 4 Segundo parcial Módulo 5 Módulo 6 Examen Final
1 2
3
4
5
6
7
SEMANAS 8 9 10 11 12 13 14 15 16
Métodos Numéricos
MODULO 1.
Rodrigo A. Ocoró
6
APROXIMACIONES Y ERRORES
1.1 INTRODUCCIÓN Dado que la mayoría de los métodos numéricos utilizan conceptos de aproximación, podemos esperar que involucren errores; por ello es indispensable definir aspectos inherentes a los errores al realizar cálculos. En muchas ocasiones se presentan varios algoritmos para producir los resultados requeridos y se necesita escoger entre ellos para lo cual es necesario tener en cuenta la rapidez y la exactitud pero la obtención de esta última evidencia la presencia de error debido a que normalmente los datos de entrada provienen de mediciones que involucran algunos errores y en ellos pueden intervenir aproximaciones. Al resolver un problema matemático por medio de una calculadora de bolsillo, estamos conscientes de que los números decimales que calculamos quizás no sean exactos. Estos números casi siempre se redondean cuando los registramos. Aún cuando los números no se redondeen de manera intencional, el número limitado de dígitos de la calculadora puede provocar errores de redondez (una calculadora de bolsillo científica puede tener 10 dígitos, mientras que una estándar puede tener 6). En una computadora electrónica los errores de redondeo aparecen por las mismas razones y afectan los resultados de los cálculos; en algunos casos causan efectos muy serios haciendo que los resultados de los cálculos carezcan por completo de sentido. De allí la importancia de aprender algunos aspectos básicos de las operaciones aritméticas en las computadoras y comprender bajo qué circunstancias pueden ocurrir severos errores de redondeo por causa de las aproximaciones. 1.2 TIPOS DE ERROR Los errores numéricos se generan con el uso de aproximaciones para representar las operaciones y cantidades matemáticas. De acuerdo con la forma en que se producen o generan los errores, tenemos: 1.2.1 Error de redondeo: Se debe a que la computadora sólo puede representar cantidades con un número limitado de dígitos, por lo tanto los valores verdaderos no son expresados exactamente ya que parte de la fracción decimal están redondeada pues ha sido acortada después del dígito final. 1.2.2 Error de truncamiento: Se debe a las aproximaciones utilizadas en la fórmula matemática del modelo. Los errores de truncamiento representan la diferencia entre una formulación matemática exacta del problema y la aproximación dada por un método numérico. Por ejemplo el resultado exacto de un problema es 7.6784231 y el método numérico arroja 7.678.
Métodos Numéricos
Rodrigo A. Ocoró
7
1.2.3 Error significativo: Ocurre cuando al ejecutar una operación en una máquina calculadora, el número de cifras válidas es menor que lo esperado. Por ejemplo la operación 7.898845 – 7.898732 = 0.000113, en una máquina que maneje notación de punto flotante para 6 dígitos genera como resultado 0.113000 x 10 –3 y los tres últimos dígitos no son significativos. 1.2.4 Error propagado: En el error que se obtiene en la salida (respuesta) generada en los pasos sucesivos debido a la ocurrencia de un error anterior. 1.3 FORMAS DE MEDIR EL ERROR 1.3.1 Error absoluto: Diferencia positiva entre los valores exacto y calculado. E ∂ = V e −V c
1.3.2 Error relativo: Es la relación aritmética entre el error absoluto y el valor exacto. E R =
V e − V c x100% V e
* Se recomienda usar el error relativo. Ejercicio 1. Supongamos que debemos evaluar la velocidad con la cual se desplaza una partícula al cabo de t=4,01 sabiendo que la función de velocidad está expresada como v(t)=10t2+14t. Si redondeamos t=4 a) Cuál es el error en t? b) Cuál es el error en v(t)? Solución: a) Error en t: Ve=4,01 Vc=4 E ∂ = 4,01 − 4 = 0,01 seg. E R =
4,01 − 4 x100% = 0,00249377 x 100% = 0,25% 4,01
b) Error en v(t) Ve V(4,01) =10 ∗ (4,01)2 +14 ∗ 4,01= 216,941 Vc V(4) =10 ∗ 42 +14 ∗ 4= 216 E ∂ = 216,941 − 216 = 0,941 E R =
216,941 − 216 ∗100% = 0,00433758 x 100% = 0,43% 216,941
Métodos Numéricos
Rodrigo A. Ocoró
8
Observe que el error absoluto que se propaga para V(t) es 94 veces; el error cometido para t (0,941 94 ∗ 0,01). Nota: Como en muchos casos no es posible tener el valor exacto (Ve) entonces en vez de establecer el error relativo (E R) podemos determinar un error de aproximación (EX) definido así: E X =
Valor actual - Valor previo x 100% Valor actual
Este es indispensable en los procesos iteractivos en los que se desea tener un valor aproximado cuando el error sea menor que un valor prefijado ( ε ). Entonces se procede hasta que se cumpla
E < ε x .
1.4 CIFRAS SIGNIFICATIVAS Cuando se emplea un número en un cálculo, debe haber seguridad que pueda usarse con confianza. El concepto de cifras o dígitos significativos se ha desarrollado para designar formalmente la confiabilidad de un valor numérico. El número de cifras significativas es el número de dígitos, más un dígito estimado que se pueda usar con confianza. Por ejemplo: Los números: 0.00001845 0.0001845 0.001845 0.01845 0.1845 Tienen 4 cifras significativas (los ceros antes de la primera cifra no nula no se consideran cifras significativas). 0.0020051 tiene 5 cifras significativas. 3.005 tiene 4 cifras significativas. El número 45300 puede tener tres cuatro o cinco cifras significativas dependiendo de los ceros que se conozcan con exactitud. La incertidumbre se puede desechar usando notación científica. 4.53 x 104 tiene 3 cifras significativas 4.530 x 104 tiene 4 cifras significativas 4.5300 x 104 tiene 5 cifras significativas
Métodos Numéricos
Rodrigo A. Ocoró
9
1.5 REGLAS DE REDONDEO Las siguientes reglas dan la pauta a seguir en el redondeo de números cuando se realizan cálculos a mano: a) El último dígito que se conserva se aumenta en uno si el primer dígito a descartar es mayor de 5 de lo contrario (si es menor a 5) se conserva igual. b) Si el primer dígito a descartar es 5 o es 5 seguido de ceros, entonces el último dígito a conservar se incrementa en 1, sólo si es impar. c) En la suma y resta, el redondeo se lleva a cabo de forma tal que el último dígito retenido en la respuesta corresponda al último dígito más significativos de las cantidades que están sumando restando (un dígito en columna de centésimas es más significativo que uno en una columna de milésimas). Ejercicios: Redondee cada uno de los números a las cifras significativas dadas: 5.6723 a 5.6723 a 10.406 a 7.3500 a 88.21650 a 7289 a
3 cifras significativas 2 cifras significativas 4 cifras significativas 2 cifras significativas 5 cifras significativas 2 cifras significativas
2. Efectúe 4.68 x 10-7 + 8.3 x 10-4 – 2.28 x 10-4 Se deben expresar los números con un mismo exponente: -4 -4 -4 0.004668 x 10 + 8.3 x 10 – 2.28 x 10 (dígito más significativo es el de 3) = 6.02468 x 10-4 = 6.0 x 10-4 3. Efectúe 0.137 x 101 + 0.269 x 10-1 Se ajustan los exponentes para poder efectuar las operaciones: 0.137 x 101 + 0.00269 x 101 = 0.13969 x 101 a) Recortado a 3 cifras significativas: 0.139 x 101 E R =
0.13969 − 0.139 x100% 0.13969 ER = 0.49 %
b)Redondeado a 3 cifras significativas: 0.140 x 101 E R =
0.13969 − 0.140 x100% 0.13969
ER = 0.22 %
Observe que el redondeo requiere de una operación extra, añadiendo con ello costo por ello es que la mayoría de las computadoras recortan (truncan) en vez de redondear.
Métodos Numéricos
Rodrigo A. Ocoró
10
1.6 NÚMEROS EN LA COMPUTADORA Por razones prácticas sólo puede manejarse una cantidad finita de bits para cada número en una computadora y esta cantidad o longitud varía de una máquina a otra. Por ejemplo cuando se realizan cálculos de ingeniería y ciencias es deseable una longitud grande pero para procesamientos administrativos es más económica y útil para los cálculos una longitud pequeña. El sistema numérico que usamos cotidianamente se llama sistema decimal cuya base es 10, pero las computadoras no usan el sistema decimal en los cálculos ni en la memoria, sino que usan el binario (base 2). Este sistema es natural para las computadoras ya que su memoria consiste en un enorme número de dispositivos de registro magnético y electrónico, en los que cada elemento solamente tiene los estados de “encendido” y “apagado” (1 ó 0). Sin embargo, en los lenguajes de máquina se usan otros sistemas numéricos como el octal (base 8) y el hexadecimal (base 16) pues las expresiones son más cortas que en binario proporcionando un uso más eficiente del espacio en memoria. 1.6.1. Números enteros: para una computadora dada, el número de bits generalmente se llama palabra y estos van desde 8 bits hasta 64 bits. Para facilitar su manejo la palabra se divide en partes más cortas denominadas bytes; por ejemplo, una palabra de 32 bits puede dividirse en 4 bytes (8 bits cada uno). Cada palabra, cualquiera sea su longitud, almacena un número aunque en ciertas circunstancias se usan varias para contener un número. Por ejemplo, consideremos una palabra de 16 bits, el primero representa el signo del número (un cero es signo más y un uno es signo menos); los 15 bits restantes pueden usarse para guardar números binarios desde 000000000000000 y 111111111111 111. Al convertirlo en decimal se obtiene: 1 x 214 + 1 x 213 + 1x 212 + ……. + 1 x 21 + 1 x 20 = 32767 (215 – 1) Por tanto cada palabra de 16 bits puede contener un número entero en el intervalo –32768 a +32767. (-215 a 215 – 1) Si la palabra es de 32 bits el primer bit representa el signo del número (0 si es positivo y 1 si es negativo) 30 29 1 0 31 1 x 2 + 1 x 2 + ……. + 1 x 2 + 1 x 2 = 2147483647 (2 –1) Por tanto cada palabra de 32 bits puede contener un número en el intervalo –2147483648 a +2147483647. (-231 a 231 –1) La menor y la mayor cantidad de un número que se puede representar en una computadora varían de acuerdo con el diseño tanto del hardware como del software (algunas pueden tener un rango de 2.9 x 10-39 hasta 1.7 x 1038).
Métodos Numéricos
Rodrigo A. Ocoró
11
n-1 En general si n es el número de2n-1 bits -2 -1 de una palabra, entonces el rango de números enteros osc ila entre
Ejemplo: Represente el número –525 en una palabra de 16 bits. –525 = -10000011012 1 0 0 0 0 0 1 0 0 0 0 0 1 1 0 1 Representa el signo – del número 1.6.2. Números reales (punto flotante): el formato para un número real en una computadora difiere según el diseño de hardware y software. Cuando se desea almacenar un número real, se emplea en su representación binaria, llamada de punto flotante, la notación en precisión simple de 32 bits: (± .abbbbbbbb …….. bbbb) x 2z El número decimal debe convertirse primero al binario más cercano en el formato normalizado(1). En la notación anterior, a es siempre 1 y cada b es un dígito binario (0 ó 1) y Z es un exponente que también se expresa en binario. Los 32 bits se distribuyen así: El primer bit: para el signo de la mantisa (0 si es +, 1 si es negativo). Los siguientes 7 bits: para el exponente z (binario) (*) incluyendo el signo del exponente (0 si es +, 1 si es negativo). Y los últimos 24 bits: para la mantisa o fracción.
1 Bit para El signo 7 bits para el exponente (característica) ♦
1
24 bits para la fracción (mantisa)
Si la palabra es de 16 bits se tiene: primer bit para el signo (0 si es +, 1 si es negativo), 7 bits para el exponente (binario) y 8 bits para la mantisa. Normalizar :
consiste en ajustar el exponente de tal manera que el primer dígito de la expresión binaria no
sea cero. * En algunas computadoras el exponente se almacena en base 8 (octal), o en base 16 (hexadecimal) en lugar de base 2.
Métodos Numéricos
♦
Rodrigo A. Ocoró
12
Algunas computadoras como las Mainframe IBM usan el formato de punto flotante normalizado en hexadecimal que se expresa como: X = (0.abbbbb) x 16k x = número en base decimal a = dígito hexadecimal diferente de cero, cada b es un dígito hexadecimal y k es un exponente binario. La mantisa tiene 6 dígitos hexadecimal.
Ejemplo: Representar el número decimal -125.32 en palabra: a) de 16 bits b)32 bits Parte entera: 125=1111101 Parte fraccional: .010100011110101 ( se toma .32 se multiplica por la base Indicada (2); se escoge el dígito que quede a la izquierda del punto y la otra parte de la derecha se multiplica nuevamente por la base y así sucesivamente. Entonces: -125.32 = -1111101.010100011110101(2) Normalizamos el número binario: desplazamos el punto flotante 7 lugares hacia la izquierda (por qué?) y lo expresamos con exponente binario (7 es 111) => -. 1111101010100011110101x2111 a)en palabra de 16 bits: 10000111 11111010 los primeros 8 campos son para signo del número, signo del exponente, y el exponente (6 campos) mientras que los otros 8 campos corresponden a la mantisa del número el cual es truncado a partir del dígito 9.
b)en palbra de 32 bits: 10000111 111110101010001111010100
los 8 primeros campos siguen la
misma estructura del a) pero la mantisa puede albergar hasta 24 dígitos.
1.6.3 Doble precisión: la doble precisión es un esfuerzo para aumentar la exactitud de los cálculos adicionando más bits a la mantisa. Esto se hace utilizando dos palabras: la primera en la forma expuesta anteriormente y los bits de la segunda para aumentar a mantisa de la primera. Así, una palabra de 16 bits al usa en doble precisión queda la mantisa de 8 + 16 = 24 bits, mientras que una de 32 bits en doble precisión permite a la mantisa 24 + 32 = 56 bits. La desventaja de la doble precisión es que se emplean más palabras con lo cual se consume más memoria para un programa. 1.6.4 Error al almacenar un número en memoria: como no es posible guardar un número binario de longitud infinita o un número de más dígitos de los que posee la mantisa de la computadora que se está empleando, se almacena sólo un número finito de estos dígitos y como consecuencia se comete automáticamente un pequeño error conocido como error de truncamiento, que al repetirse muchas veces puede llegar a se considerable. Ejercicio. Qué número en base decimal representan los siguientes dígitos almacenados en palabra de 16 bits?. a) 1100001110101001 se desnormaliza (0bserve que el exponente es negativo) = -.10101001x2-11 = -.00010101001 = -(2-4+2-6+2-8+2-11) = -0.08251958 b) 1000001110101001 Se desnormaliza: (observe que el exponente es positivo) = -.10101001x211 = -101.01001 = -(22+20+2-2+2-5)= -5,28125
Métodos Numéricos
Rodrigo A. Ocoró
13
1.7. LA SERIE DE TAYLOR. Permite epresar funciones en forma polinómica y provee una formulación para predecir el valor de la función en un punto x i+1 een términos del valor de la función y de sus derivadas en otro punto xi . Teorema de Taylor : Si la función f y sus primeras n+1 derivadas son continuas en un intervalo que contiene a y x, entonces el valor de la función de x está dada por: f ( x) = f ( a ) + f ' ( a )( x − a ) +
f ' ' ( a )
2
( x − a ) +
2!
f ' ' ' ( a )
3
( x − a) + ......
3!
f
( n)
(a )
n!
( x − a ) n + Rn
(1) f
Rn =
( n +1)
(c )
( x − a) n +1
(n + 1)! donde Rn es el residuo que se define como si evaluamos la expresión (1) en el punto x=xi+1 siendo a=xi y hacemos h=(x-a)=(xi+1 –xi) se obtiene la expresión de la serie de Taylor así: f ( xi +1 ) = f ( xi ) + f ' ( xi )h +
f ' ' ( xi )
2!
2
h +
f ' ' ' ( xi )
3!
3
h + ......
f
(n)
( xi )
n!
n
h + Rn
(2) Rn =
f
( n +1)
(c )
h
n +1
(n + 1)! donde: f(n) es la derivada de orden n y , c ε [xi+1,xi] el paso inicial o aproximación de orden cero (n=0) es: f(x i+1)=f(xi) f ( xi +1 ) = f ( xi ) + f ' ( xi )h el primer paso o aproximación de orden 1 (n=1) es:
f ( xi +1 ) = f ( xi ) + f ' ( xi )h +
el segundo paso,aproximación orden n=2: y así sucesivamente hasta el paso de aproximación de orden n:
f ( xi +1 ) = f ( xi ) + f ' ( xi ) h +
f ' ' ( xi )
2!
que podemos resumir así:
2
h +
f ' ' ' ( xi )
3!
3
h + ......
f ( xi +1 ) = f ( xi ) +
f
( n)
( xi )
n!
n
f ( i ) ( xi )h i
i =1
i!
∑
h
f ' ' ( xi )
2!
n
Ejemplo 1.Mediante la serie de Taylor hasta orden 4, encuentre el valor de f(x)= -0.1x4-0.15x3-0.5x2 -0.25x +1.2 en x=1 alrededor del punto x=0 necesitamos estimar f(1) tomando como base x=0 xi+1=1 y xi=0 por lo tanto h= xi+1 - xi =>h=1-0 =1 El valor exacto Ve es: f(1)= -0.1(1)4-0.15(1)3-0.5(1)2 -0.25(1) +1.2=0.2 La seriie de Taylor hasta n=4 es: f ( xi +1 ) = f ( xi ) + f ' ( xi ) h +
f ' ' ( xi )
2!
2
h +
f ' ' ' ( xi )
3!
3
h +
f
(4)
( xi )
4!
4
h + R4
h
2
Métodos Numéricos
Rn =
f
( 5)
Rodrigo A. Ocoró
(c)
14
h5
(5)! donde calculemos los valores de las derivadas en xi=0
f(x)= -0.1x4-0.15x3-0.5x2 -0.25x +1.2 f(0)=1.2 f’(0)=-0.25 f ‘(x)=-0.4x3-0.45x2-1.0x-0.25 2 f’’(0)=-1.0 f ‘’(x)=-1.2x -0.90x-1.0 f’’’(0)=-0.9 f ‘’’(x)= -2.4x -0.90 (4) (4) f (0)= -2.4 f (x)=-2.4 f(5)(x)=0 por lo tanto Rn=0 Veamos que sucede en cada paso de Taylor y hallemos el Error relativo en c/u: Recuerde que el ve es f(1)= 0.2 Orden de aprox n=0 n=1 n=2
n
f ( xi +1 ) = f ( xi ) +
∑
Vc
n
f ( ) ( xi )
i =1
n!
hn
f(1)=f(0)=1.2 f(1)=f(0)+f’(0)x1=1.2-0.25x1 f(1)=f(0)+f’(0)xh+(f”(0)xh 2)/2!=1.2-0.25x1(1.0x1)/2 f(1)=f(0)+f’(0)xh+(f”(0)xh 2)/2!+(f’’’(0)xh3)/3! =1.2-0.25x1-(1.0x1)/2-(0.9x1)/6 f(1)=f(0)+f’(0)xh+(f”(0)xh 2)/2!+(f’’’(0)xh3)/3! +(f(4)(0)xh4)/4!=1.2-0.25x1-(1.0x1)/2(0.9x1)/6 –(2.4x1)/24
n=3 n=4
E R =
ve − v c ve
x100%
1.2 0.95 0.45
500% 375% 125%
0.30
50%
0.20
0%
Observe que en el paso 4 se obtiene f(1)=0.2 que corresponde al valor exacto. Ejemplo 2. Usando la serie de Taylor con 5 pasos, estime el valor de f(x)=cos x en π
π
x= 3 tomando como base x= 4 π
Xi+1= 3 (60º) Veamos:
π
π π
π
Xi= 4 (45º) h= Xi+1- Xi = 3 - 4 = 12 π
f(x)=cosx
f( 4 )=0.7071068 π
f’(x)= -sen x
f’( 4 )=-0.7071068 π
f’’(x)= -cos x
f’’( 4 )=-0.7071068
(debe utilizar Radianes (RAD) en la calculadora
Métodos Numéricos
Rodrigo A. Ocoró
15
π
f’’’(x)= sen x
f’’’( 4 )=0.7071068 π
f(4)(x)=cos x
f(4)( 4 )= 0.7071068 π
f(5)(x)=-sen x f(4)( 4 )= -0.7071068 La serie de Taylor para n=5 pasos de aprox. es: f ( x i +1 ) = f ( x i ) + f ' ( x i ) h +
f ' ' ( x i )
2!
2
h +
f ' ' ' ( x i )
3!
3
h +
f
( 4)
( x i )
4!
4
h +
f ( 5) ( x)
5!
h
5
Entonces: π π
π π
π π
π π
f ' ' ( )( ) 2 f ' ' ' ( )( ) 3 f ( 4 ) ( )( ) 4 f (5)( )( ) 5 π π π π 4 12 + 4 12 + 4 12 + 4 12 f ( ) = f ( ) + f ' ( )( ) + 3 4 4 12 2! 3! 4! 5!
π
Al reemplazar y calcular estos valores se obtiene que f( 3 )=0.5000003 E R =
ve − v c ve
x100%
0.5 − 0.5000003 0.5
Y el error relativo es: = Si se agregan mas pasos pues se obtiene mas exactitud.
x100% = 0.00006%
Métodos Numéricos
Rodrigo A. Ocoró
16
ACTIVIDADES Nº 1 1. Calcule el error absoluto y relativo cometido en c/u de las siguientes situaciones: a)El número 3.0715 es redondeado a 3.07 b)Se redondea π =3.1 para evaluar f (π ) = 6π 4 − 5π 3 + 2π − 3 c)Se toma x= 2 ≅ 1.41 para evaluar la función f(x)=x 4-x2+ lnx d)Si f(x)=ex – 1, cuál es el error relativo en cada caso para valores de x={0.5; 1; -1; } al tomar el valor de f(x) aproximado a dos decimales mientras el valor exacto lo considera con 5 decimales?. 2. Se tiene para medir la longitud de un puente y de un remache, obteniéndose en cada caso 9999 y 9 cm. Si los valores verdaderos son 10000 y 10 cm respectivamente, calcule el error absoluto y relativo en cada uno y explique cómo considera la estimación de medición hecha para caso. 3. Redondéese los siguientes números a tres cifras significativas: a) 8755 b) 5555x103 c) 0.368124x102 d)4225.0002 e)0.999500 4. Encuentre el error relativo tanto “recortado” como “redondeado” a tres dígitos significativos: a) 3.28x10-2 + 2.07x103 b) (3.26x10-3 + 2.07x10-4) – 2.01x10-4 d) 3.26x10-3 + (2.07x10-4 – 2.01x10-4 ) e) (3.28x10-2 x 6.98x103) ÷ 4.82x10-8 f) 3.28x10-2 x (6.98x103 ÷ 4.82x10-8 ) g) (4.82x102 ÷ 8.81x108) ÷ 4.06x10-2 5. Repase sus apuntes de semestres anteriores para recordar las conversiones de los siguientes sistemas numéricos: a) de base decimal a binaria b) de base binaria a decimal c) de base decimal a octal d) de base octal a decimal e)de base hexadecimal a decimal f) de base decimal a hexadecimal g)de base octal a binario h)de base binaria a octal 6. Convierta c/u de los siguientes números decimales a binarios, octal y hexadecimal: a)536 b)923 c)-1536 d)5000 e)0.8 f)0.973 g)125.35 7. Convierta los siguientes números del sistema octal al binario: a)777 b)573 c)3052 d)416.42 8. Convierta los siguientes números dados en binario al sistema decimal: a)1000 b)10101 c)11001 d)-111111 e)0.0001 f)0.11111 g)1011.0101 h)-1001011.011001 9. Diseñe algoritmos para convertir números: a) del sistema decimal a binario y compruébelo aplicándolo a los ejercicios del punto 7. b) del sistema binario a decimal y compruébelo aplicándolo a los ejercicios del punto 9. c) del sistema decimal a octal y compruébelo aplicándolo a los ejercicios del punto 7. d) del sistema octal a decimal
Métodos Numéricos
Rodrigo A. Ocoró
17
e) del sistema decimal a hexadecimal f) del sistema hexadecimal a decimal 10. Considere una máquina con una palabra de 8 bits. Qué rango de números enteros puede contener dicha palabra?. Represente el número –26 en una palabra de 8 bits. 11. Qué decimales representan cada uno palabra de 16 bits? : a) 01 b) 1 0 c) 0 0 Repuestas: 0.19921875;
de los siguientes números de máquina en una 0 0 0 0 1 0 1 10 0110 0 0 0 1 0 1 1 1 0 0 1010 0 0 1 1 0 0 0 1 0 0 0111 0
-160;
9306112
12. Normalice los siguientes números y almacene en memoria con palabra de 32 bits a)723.5578 b)-15.324 c)0.003485 sugerencia: Pase los números a binarios y luego normalícelos. 13.Elabore un programa de modo que: a) 0.01 se sume 100 veces b) 0.001 se sume 1000 veces c) 0.0001 se sume 10000 veces Imprímase los valores para las sumas y compare sus resultados con los valores exactos. Esperaba que en todos los casos diera 1?. Qué puede influir para no obtener el valor de 1 en todos los casos?. 14. Usando la serie de Taylor, estime en el valor indicado y calcule el ER: a) f(x)=25x3-6x2-7x-88 hasta el tercer orden para f(3) con punto base x=2 b) f(x)=ln x hasta el paso 4 para f(4) alrededor de x=1 π
π
3
4
c) f(x)=Sen x hasta el paso 4 para f( ) con punto base x= x
d) f(x)=e -1 hasta n=5 para f(-3) alrededor de x=0
15. La forma polinomial (serie de Maclaurin) para ex se expresa así: f ( x) = 1 + x +
x
2
2!
+
x
3
3!
+ ...... +
x
n
n!
Demuestre que la expresión anterior se obtiene aplicando la serie de Taylor para f(x)=ex siendo xi=0 y h=x. 16. Obtenga la forma polinomial para f(x)=cos x (aplique la serie de Taylor con xi=0 y h=x ). 17. Obtenga la forma polinomial para f(x)=Sen x
Métodos Numéricos
Rodrigo A. Ocoró
18
MÓDULO 2. SOLUCIÓN DE ECUACIONES NO LINEALES 2.1 RAÍCES DE UNA ECUACIÓN Desde hace años se aprendió a usar la fórmula cuadrática
− b ± b 2 − 4ac x= 2a
para
resolver f(x) = ax 2 + bx + c = 0. A los valores x calculados se les llama “raíces” de la ecuación y representan los valores que hacen la ecuación igual a cero. Por lo tanto se puede definir la raíz de una ecuación al valor τ de la función f(x) si y sólo sí f( τ ) = 0. A las raíces de una ecuación también se les llama ceros de la función. Aunque la fórmula cuadrática es útil para resolver ecuaciones de 2º grado hay muchas funciones diferentes que no se pueden resolver de manera tan fácil. En estos casos los métodos numéricos que estudiaremos nos proporcionan medios eficientes para obtener respuesta. 2.2 FUNDAMENTACIÓN PRELIMINAR Cuando se requiere encontrar las raíces reales de una función algunos pasos o consideraciones previas las situarán dentro de ciertos intervalos del eje real; el objetivo será reducir considerablemente el ancho de cada intervalo. Teóricamente si estrechamos el intervalo hasta reducirlo a un punto, ese punto será una raíz de la función. En la práctica no será necesario que cada raíz se determine de manera exacta a cambio de eso se aceptarán pequeños intervalos que contengan una o más raíces. Como el intervalo no se reduce a un punto, la raíz es obtenida con error de truncamiento, el cual puede ser disminuido en la medida del aumento de las iteraciones. Pero al final, por el número de cifras decimales que puedan ocurrir, tendremos un error de redondeo. La elección de un método para disminuir la longitud del intervalo, sólo depende de qué tan pequeño deba ser éste. Teorema: si para una función f(x) continua en un intervalo cerrado [a, b] se cumple que f(a)f(b) < 0, entonces existe por lo menos un r ε [a, b] tal que f(r) = 0 f(a)
a
r1
r2
r3
b
f(b)
Figura No. 1 Raíces de la función f(x)
Métodos Numéricos
Rodrigo A. Ocoró
19
La figura 1 presenta un caso en donde f(a)f(b) < 0 con tres raíces sobre el intervalo [a, b]. El teorema anterior únicamente manifiesta que existe al menos una raíz, pero no define cuántas hay ni donde están localizadas. Los métodos para encontrar las raíces aproximadas de un polinomio se fundamentan gráficamente en este teorema para generar un algoritmo de búsqueda. 2.3 MÉTODO DE BISECCIÓN También se conoce en la literatura sore el tema como método de reducción de intervalos a la mitad o “Método Bolzano”. Se fundamenta en el conocido método de búsqueda binaria utilizado para resolver muchos problemas de la programación. El algoritmo de búsqueda consiste por lo tanto de los siguientes pasos: ♦ Paso 1: dar f(x); Localizar un intervalo [a, b] que cumpla f(a)f(b) < 0. Longitud máxima del intervalo de aceptación de la raíz = ε , por ejemplo ε = 0.001. ♦ Paso 2: I=1; Hacer X 1 = a, x2 = b ♦ Paso 3: Calcular :Xm = ½ (X 1 + X2), evalúe f(Xm) ♦ Paso 4: Si |f(Xm)| ≤ ε , tomar a Xm como una raíz aproximada de f(x). ♦ Paso 5: I=I+1 ♦ Paso 6: Si f(x1). f(xm) < 0, entonces x2 = xm, de lo contrario hacer x1 = xm, regresar al paso 3. 2.3.1 Ejemplo: si buscamos una raíz para la ecuación: f(x) = xe x -2 = 0 y evaluamos f(x) en algunos enteros encontramos que f(0) = -2 y f(1) =0.718288 Como los polinomios son funciones continuas para cualquier rango, podemos asegurar que existe al menos una raíz en algún punto del intervalo [0, 1]. De conformidad con el algoritmo propuesto, la tabla No. 1 muestra todos los cálculos iterativos realizados para encontrar una raíz aproximada en este intervalo con ε < 0.005. Tabla No. 1. Método de Bisección para xex -2 = 0 Iteración 1 2 3 4 5 6 7
X1
0 0.5 0.75 0.75 0.8125 0.84375 0.84375
X2
1 1 1 0.875 0.875 0.875 0.859375
Xm
0.5 0.75 0.875 0.8125 0.84375 0.859375 0.8515625
f(Xm) -1.175639 -0.412250 0.099016 -0.169003 -0.038222 0.029572 -0.004528
f(X1)
-2 -1.175639 -0.412250 -0.412250 -0.169003 -0.038222 -0.038222
Después de siete iteraciones sabemos que hay una raíz entre 0.84375 y 0.859375 y podemos aceptar a 0.8515625 como una raíz suficientemente aproximada. Aunque éste método de bisección es adecuado para calculadoras científicas, siempre será preferible el utilizar programas de computadora para la búsqueda de raíces.
Métodos Numéricos
Rodrigo A. Ocoró
20
2.4 MÉTODO DE REGULA FALSI Se fundamenta en los conceptos elementales de interpolación lineal. Esto es, una vez se tenga un intervalo [a, b] para el que f(a)f(b) < 0, se supone un segmento de recta que une los puntos (a,f(a)) y (b,f(b)). El cual corta el eje de las abcisas en el punto x = c. Este valor se utiliza en el algoritmo como una posible raíz. La figura No. 2 muestra gráficamente la raíz r1 y el valor c próximo a ésta para una función. Mediante propiedades de triángulos semejantes podemos escribir la siguiente proporción: b − c f(b) =− , entonces : c − a f(a) c =
c =
af(b) − bf(a) f(b) − f(a)
x 1 f(x 2 ) − x 2 f(x 1 ) f(x 2 ) − f(x 1 )
Si x1 = a, x2 = b ====> Y
f(b)
a
c
b r1
X
f(a) Figura No. 2 Gráfico para encontrar una Raíz de la función f(x) El algoritmo para este caso consiste por lo tanto de los siguientes pasos: Paso 1: Dar f(x); Localizar un intervalo [a, b] que cumpla f(a)f(b) < 0. Longitud máxima ♦ del intervalo de aceptación de la raíz = ε , por ejemplo ε = 0.001. Paso 2: I=1; Hacer x1 = a, x2 = b ♦ ♦
Paso 3: Calcular
c=
x1 f ( x 2 ) − x 2 f ( x1 ) f ( x 2 ) − f ( x1 )
, evalúe f(c)
Paso 4: Si |f(c)| ≤ ε , tomar c como una raíz aproximada de f(x). ♦ Paso 5: I=I+1 ♦ Paso 6: Si f(x1). f(c) < 0, entonces x 2 = c, de lo contrario hacer x 1 = c y regresar al paso 3. 2.4.1 Ejemplo: supongamos que necesitábamos buscar una raíz para la ecuación: f(x) = xex – 2 = 0. Al evaluar f(x) en algunos enteros encontramos que: f(0) = -2; f(1) = 0.71828183. podemos asegurar que existe al menos una raíz en algún punto de intervalo [0, 1]. ♦
Métodos Numéricos
Rodrigo A. Ocoró
21
De conformidad con el algoritmo propuesto, la tabla No. 2 muestra todos los cálculos iterativos realizados para encontrar una raíz aproximada en este intervalo con ε < = 0.001. Tabla No. 2. Método de Regula Falsi para xe x – 2 = 0 Iteración X1 X2 f(X1) f(X2) 1 0 1 -2 0.71828183 2 0.73575888 1 -046442322 0.71828183 3 0.83952077 1 -0.05629351 0.71828183 4 0.85118386 1 -0.00617135 0.71828183
C f(c) 0.73575888 -0.46442322 0.83952077 -0.05629351 0.85118386 -0.00617135 0.85245157 -0.00066888
Después de cuatro iteraciones vemos que hay una raíz entre 0.85118386 y 1; podemos aceptar a 0.85245157 como una raíz suficientemente aproximada. 2.5 MÉTODO DE LA SECANTE Es una variante del método Regular Falsi donde la convergencia hacia la raíz se hace buscando un nuevo valor llamado xk+2 a partir de los puntos anteriores xk y xk+1. Transformando la expresión que calcula el punto c para el método Regula Falsi, obtenemos; x k + 2 =
x k f(x k +1 ) − x k +1 f(x k ) f(x k +1 ) − f(x k )
x 3 =
Si x1 = xk, x2 = xk+1, x3 = xk+2 ====>
x 1f(x 2 ) − x 2 f(x 1 ) f(x 2 ) − f(x 1 )
El algoritmo para este método consistiría en los siguientes pasos: Paso 1: Dar f(x); Localizar un intervalo [a, b] que cumpla f(a)f(b) < 0. Longitud máxima ♦ del intervalo de aceptación de la raíz = ε , por ejemplo ε = 0.001. Paso 2: I=1; Hacer x1 = a, x2 = b ♦ ♦
♦
♦ ♦
Paso 3: Calcular
x3 =
x1 f ( x 2 ) − x 2 f ( x1 ) f ( x 2 ) − f ( x1 ) ≤ ε , tomar x3
, evalúe f(x3)
Paso 4: Si |f(x3) - f(x 2)| como una raíz aproximada de f(x). Termina el proceso de convergencia. De lo contrario vaya al paso siguiente. Paso 5: I=I+1 Paso 5: Hacer x1 = x2, x2 = x3, regresar al paso 3.
2.5.1 Ejemplo: supongamos que necesitábamos buscar una raíz para la ecuación: f(x) = xex – 2 = 0. Tomando como intervalo [0, 1] con ε < = 0.001. Iterac X1 1 0 2 1 3 0.7357589 4 0.8395208 5 .85383269
X2 1 0.735759 0.839521 0.853833 0.852593
F(x1)
F(x2) -2 0.71828183 0.71828183 0.46442322 -0.46442322 0.05629351 -0.05629351 0.00533812 0.00533812 -5.3928E-05
X3 0.73575888 0.83952077 0.85383269 0.85259309 0.85260549
F(x3) |f(x3)-f(x2)| -0.46442322 1.182705 -0.05629351 0.40812970 0.00533812 0.06163163 5.3928E-05 0.00539 -5.0932E-08 0.00005
Después de 5 iteraciones podemos aceptar como raíz aproximada x3= 0.85260549
Métodos Numéricos
Rodrigo A. Ocoró
22
2.7 MÉTODO DE NEWTON RAPSHON Uno de los métodos más ampliamente utilizados para resolver ecuaciones es el método de Newton. Comenzando a partir de una estimación inicial que no está lejos de la raíz, se extrapola a lo largo de la tangente hasta su intersección con el eje x y se toma esa como la siguiente aproximación. Esto se continúa hasta que los valores sucesivos de x estén suficientemente cercanos o el valor de la función sea lo suficientemente próximo a cero. El punto donde la tangente cruza al eje x, representa una aproximación de la raíz.
f(xk )
Xk+1
θ
xk
La pendiente de la recta será: tanθ = f ’(xk) = xk – xk+1 =
f ( x k ) f ' ( x k )
xk+1 = xk –
f ( x k ) x k − x k +1
f ( x k ) f ' ( xk )
Una de las ventanas de este método consiste en la convergencia rápida que se tiene cuando se está próximo a una raíz. ♦ ♦
♦ ♦
♦
♦
El algoritmo de este método consiste en los siguientes pasos: Paso 1: Dar f(x); Localizar un intervalo [a, b] que cumpla f(a)f(b) < 0. Longitud máxima del intervalo de aceptación de la raíz = ε , por ejemplo ε = 0.001. Paso 2: I=1; Iniciar con algún valor x kE [a, b], para el cual f ’(xk) ≠ 0; calcule f ’(x). Paso 3: Calcular
xk+1 = xk –
f ( x k ) f ' ( x k )
Paso 4: Si |xk+1 - xk| ≤ ε , tomar xk+1 como una raíz aproximada de f(x). Termina el proceso de convergencia. De lo contrario vaya al paso siguiente. Paso 5: I=I+1; Hacer xx= xk+1 y regresar al paso 3.
2.7.1 Ejemplo 1: Se propone buscar una raíz para la ecuación: f(x)=xex –2 = 0. Siendo f ’(x) = e x (x + 1), no se puede iniciar la evaluación para x = -1, puesto que f’(-1)= 0. Evaluamos f(x) en otros enteros: f(0) = -2 f(1) = 0.71828183. Podemos asegurar que existe al menos una raíz en algún punto del intervalo [0,1]. De conformidad con el algoritmo propuesto, tenemos la siguiente tabla:
Métodos Numéricos
Rodrigo A. Ocoró
Tabla No. 4. Método de Newton Rapshon para xex – 2 = 0 Iteración Xk f(xk) f’(xk) Xk+1 1 0 -2 1 2 2 2 12.7781122 22.1671683 1.42355686 3 1.42355686 3.9104113 10.0622731 1.03493579 4 1.03493579 0.91326712 5.7281926 0.87550206 5 0.87550206 0.10127495 4.50135492 0.85300329 6 0.85300329 0.0017292 4.34841325 0.85260562 En la sexta iteración se obtiene una raíz r ≈ 0.85260562.
23
|xk+1 – xk| 0.57644314 0.38862107 0.15943373 0.02249877 0.00039766
2.7.2.Ejemplo 2. Utilizando el método de Newton-Raphson, halle una raíz para f(x)=3x+senx –ex . NOTA: al utilizar funciones trigonométricas, el valor de x debe darse en radianes (RAD en la calculadora y no DEG). Con [1, 2] y tomando Xk =2 tenemos: f ‘(x)=3+cosx-ex Tabla No. 5. Método de Newton Rapshon para f(X)=3x+senx –ex=0 Iteración Xk f(xk) f’(xk) Xk+1 1 2 -0.47975 -4.80520 1.90016 2 1.90016 -0.04023 -4.01041 1.89013 3 1.89013 -0.00039 -3.93416 1.89003
|xk+1 – xk| 0.09984 0.01003 0.00010
X= 1.89003 es una raíz aproximada para f(X)=3x+senx –ex
2.8 MÉTODO DE PUNTO FIJO Se comienza con la ecuación general f(x) = 0 y se reorganiza en una expresión equivalente de la forma x = g(x) de manera que f(x0) = 0 siendo x 0 = g(x0). Algoritmo: se proporciona g(x) ♦ ♦ ♦
♦
Paso 1: Dar valor inicial de x 0, ε , número max de iteraciones (MAXI) Paso 2: Hacer I = 1 Paso 3: Haga mientras I < MAXI (realiza los pasos 4 a 7) ♦ Paso 4: Calcule g(x0) ♦ Paso 5: Si |g(x 0) – x0) < ε entonces x = g(x0) es una raíz aproximada de f(x) y termina; de lo contrario continúa Paso 6: I = I + 1 ♦ Paso 7: Se actualiza x 0 x0 = g(x0) ♦ Paso 8: Mensaje de falla: “El método no converge a una raíz”.
Métodos Numéricos
Rodrigo A. Ocoró
24
2.8.1 Ejemplo: consideremos f(x) = x 2 – 2x – 3 = 0 a) Reacomodamos así: x2 = 2x + 3 x = 2x + 3 g(x) = 2x + 3 determinamos un intervalo para tomar un valor inicial : [4, 2] Consideremos ε = 0.001 Se comienza con x0 = 4 y se evalúa g(x) = 2x + 3 Tabla No. 5. Método de punto fijo para f(x) = x 2 – 2x – 3 = 0, con g(x) = 2x + 3 Iteración X0 G(X0) |g(X0) – X0| 1 4 3.31662479 0.6833752 2 3.31662479 3.103747667 0.21287712 3 3.103747667 3.034385495 0.06936217 4 3.034385795 3.011440019 0.02294547 5 3.011440019 3.003810919 0.00762910 6 3.003810919 3.001270038 0.00254881 7 3.001270038 3.000423316 0.00084672 2 X = 3.000423316 ≈ 3 es una raíz aproximada de f(x) = x – 2x – 3 b) Reordemos x2 – 2x – 3 = 0 así: 2
x – 3 = 2x
x=
x 2 − 3
2
g(x) =
x 2 − 3
2
2
Tabla No. 6. Método de punto fijo para f(x) = x – 2x – 3 = 0, con g(x) = Iteración 1 2 3
X0
4 6.5 19.625
g(X0)
|g(X0) – X0| 6.5 2.5 19.625 13.125 191.0703 171.445
g(X0) no converge sino que diverge por lo tanto con g(x) =
una raíz aproximada para f(x) = x 2 – 2x – 3. Nota:
♦
2
2
tomo x0 = 4
no es posible hallar
g(x) es convergente si g’(x0) ≤ 1
Comprobemos: ♦
x 2 − 3
x 2 − 3
Para g(x) = 2x + 3 entonces converge
Para g(x) =
x 2 − 3
2
q’(x) =
q’(x) = x
1 2x + 3 g’(x0) = g’(4) = 0.3015 < 1
g’(4) = 4 > 1 por lo tanto no converge
Métodos Numéricos
Rodrigo A. Ocoró
25
c) Reordenemos x2 – 2x – 3 = 0 x(x – 2) = 3
x =
3
x − 2
g(x) =
3 x − 2
Compruebe si es convergente y establezca las iteraciones para hallar una raíz aproximada de f(x) = x2 – 2x – 3. g´( x) =
−3 ( x − 2) 2
===> g´(4) =
−3 (4 − 2)
2
= 0.75 < 1
Entonces si converge pues g´(x)<1 Tabla No. 7. Método de punto fijo para f(x) = x 2 – 2x – 3 = 0, con g(x) = Iteración 1 2 3 4 5 6 7 8 9 10 11
X0
4 1.5 -6.0 -0.375 -1.263158 -0.9193548 -1.0276243 -0.9908759 -1.0030506 -0.9989841 -1.0003387
g(X0)
1.5 -6.0 -0.375 -1.263158 -0.9193548 -1.0276243 -0.9908759 -1.0030506 -0.9989841 -1.0003387 -0.9998871
|g(X0) – X0| 2.5 7.5 5.625 0.888158 0.3438030 0.1082695 0.0367484 0.0121747 0.0004065 0.0013546 0.0004516
3 x − 2
X = -0.9998871≈ -1 es una raíz aproximada de f(x) = x2 – 2x – 3
Como puede observarse, con el mismo punto de inicio pero con g(x) diferentes, obtuvimos dos raíces para f(x) = x 2 – 2x – 3; x=3 y x= -1
Métodos Numéricos
Rodrigo A. Ocoró
26
ACTIVIDADES Nº 2 1.La ecuación x2-2=0, tiene por raíces ± 2 = ± 1.414214. Utilice el método de Bisección para evaluar la raíz positiva comenzando con el intervalo [1, 2]. 2.Con el método de bisección, encuentre una raíz aproximada para las siguientes ecuaciones. En cada caso determine un intervalo adecuado y luego calcule la raíz con una precisión del 0.5%.: a)x3-x2-2x+1=0 b)2e-x-senx=0 c)2x3+4x2-2x-5=0 d) ex-3x2=0 e)x3-x2+4x-7=0 3.Resuelva por el método de Newton –Raphson: a)ln x-x+2=0 b)x-2cosx=0 c)f(x)= x3-x- 3 e)x3-5x= -1
d)f(x)= x3-ex+3
f)tanx-x-1=0
4.Resuelva por el método de la secante: a)xlogx-10=0 b)senx-cscx+1=0 c)f(x)=x3+2x2+10x-20=0
d)ex –3x
5.Usando el método de Regula-Falsi determine las raices de las sigientes ecuaciones con e=0.5%: a)f(x)= x3+2x2+10x-20=0 b)f(x)=x4-2x-10 c)f(x)=x3-ex+3 6.Diseñe y desarrolle los respectivos software para c/u de los métodos. 7 .Para c/u de siguientes ecuaciones, encuentre todas las raices por el método que usted considere. Sugerencia: halle la gráfica para observar los puntos de intersección y resuelva utilizando los programas diseñados para computador considerando e=0.1%. a) ex +2x+4=0 b)tan x - x + 1=0 c) –x3 + x +1=0 d) 16x5 –20x3 + x2 +5x - 0.5= 0 e) x4 - 4.4x3 + 9.43x2 – 14.86x +7.15 =0 8.Utilice el método de punto fijo para hallar raíces de: a) x3 –10x –5=0 b) x2 + x –1=0 c) cos x – 3x =0 con i)x=cos x - 2x ii)x=(cos x)/3 d) Para 2x3 +4x2 –2x –5=0, establezca 3 reordenaciones y utilice en c/u x0=1 e) F(x)= x3-10x-5 =0, tiene 3 raices ubicadas en los intervalos [-3, -2], [-1,0] y [3,4] Hállelas. 9.Desarrolle un programa para correr el métdo de punto fijo y compruébelo con los ejercicios del punto 8.
Métodos Numéricos
Rodrigo A. Ocoró
27
MODULO 3. ÁLGEBRA LINEAL NUMÉRICA 3.1 INTRODUCCIÓN La solución de sistemas de ecuaciones es un tema clásico de las matemáticas y de gran utilidad en diversas ramas del conocimiento como la economía, la biología, física, ingenierías, etc. La resolución de sistemas de cualquier número de ecuaciones (10, 50, 100, 500, etc.) es una realidad hoy en día, gracias a los computadores, lo cual proporciona un atractivo especial a las técnicas de solución directa e interactivas tanto en su programación, cálculos necesarios, propagación de errores, etc. 3.2 FUNDAMENTOS PRELIMINARES 3.2.1 Notación matricial: Una matriz es un arreglo rectangular de números en el cual no sólo es importante el valor del número sino también su posición dentro del arreglo. El tamaño de una matriz se indica por el número de filas (n) y columnas (m): n x m.
a11 a12 ..... a1m a a .....a 2m A = 21 22 = [a ij ] M a a .....a nm n1 n2
i = 1,2,…,n j = 1,2,…,m
3.2.2 Suma de matrices: Para sumar dos matrices, éstas deben ser del mismo tamaño (igual número de filas y columnas) siendo la matriz resultante de igual tamaño. A + B = [aij + bij] = [cij] 4 8.5 - 3 − 1 2 - 4 3 10.5 - 7 = 2 - 1.3 7 + 5 8 3 7 6.7 10 2x3 2x3 2x3 Ejemplo 1:
Ejemplo 2:
4 8.5 - 3 − 1 2 - 4 5 6.5 1 2 - 1.3 7 − 5 8 3 = - 3 - 9.3 4
Métodos Numéricos
Rodrigo A. Ocoró
28
3.2.3 Producto de matrices por un escalar: Si ∝∈R, entonces: a 11 a 12 .....a 1n αa 11 αa 12 .....α. a a .....a αa 21 22 2n 21 αa 22 .....α. αA = α = M M a n1 a n2 .....a nm αa n1 αa n2 .....α. ⇒
2m nm
1m
α [a ij ] = [αa ij ]
4 - 6 - 2 − 2 3 1 Ejemplo: Multiplique por –2 la matriz A= 0 − 4 − 1 ===> − 2 A = 0 8 2 - 10 6 - 4 5 − 3 2
3.2.4 Multiplicación de matrices: Dada las matrices A y B, se pueden multiplicar en ese orden si el número de columnas de A es igual al número de filas de B. Anxm
Bmxp = Cnxp n
Donde Cij = ai1 b1j + ai2 b2j + aim bmj =
Ejemplo: Si − 1 2 A = 0 1 4 - 1
∑ a ik b kj
k =1
i = 1,2,…..n j = 1,2,…..p
2 4 0 B= - 1 2 3
− 1 2 - 4 0 6 2 4 0 = - 1 2 3 A ⋅ B = 0 1 • - 1 2 3 2x3 4 - 1 3x2 9 14 - 3 3x3
TAREA: Elabore un programa para sumar y multiplicar matrices.
Métodos Numéricos
Rodrigo A. Ocoró
29
3.2.5 Matrices especiales: Matriz vector: Cuando una de las dimensiones de la matriz es 1. Ejemplo:
A = [1 2 -3 5]1x4
Vector fila
x y B= z 3x1 Vector columna
Matriz cuadrada: Cuando el número de filas es igual al número de columnas. 3 7 x= 4 - 2
- 1 3 y= 5 - 7
D = [α ij ] nxn Triangular superior: Cuando todos los elementos bajo la diagonal principal son nulos A = [aij]nxn, aij = o ∀i>j a11 a12 a13 ...... a1n o a a ......a 22 23 2n M M M o ....... M ...........an-1n o o .......... .... o........ann
3 4 0 1 0 0 0 0
5 -2 5 0
5 - 3 7 - 4
Triangular inferior: Todos los elementos por encima de la diagonal principal son nulos A = [aij]nxn, aij = o ∀i
Matriz diagonal: Todos los elementos por fuera de la diagonal son todos nulos. A = [aij]nxn
aij = o
∀i≠ j
Métodos Numéricos
Rodrigo A. Ocoró
a11 o ..... o o M a ...... 22 M o .....a nn o
4 0 0 - 0
0 -1 0 0
0 0 5 0
30
0 0 0 3
Matriz identidad: Es la matriz diagonal donde los elementos de la diagonal son 1. A = [aij]nxn
aij = 1 si i=j aij = 0 si i≠ j
1 0 0 ...0 0 1 0 ...0 I = 0 0 1...0 M 0 ............1
Matriz transpuesta: Es aquella que resulta de intercambiar las filas por las columnas de una matriz dada. A = [aij]nxm
⇒ At = [a ji]mxn
1 0 -2 3 A= - 5 4 6 7 2x4
1 - 5 0 4 t ⇒A = - 2 6 3 7 4x2
3.3. SOLUCIÓN DE SISTEMAS DE ECUACIONES LINEALES Un gran número de problemas prácticos de ingeniería se reduce a resolver un sistema de ecuaciones. Un sistema de m ecuaciones lineales con n incógnitas, tiene la forma: a11x1 + a12x2 + ….. + a 1nxn = b1 a21x1 + a22x2 + ….. + a 2nxn = b2 M
am1x1 + am2x2 + …. + a mnxn = bm Con notación matricial se escribe así:
Métodos Numéricos
a11 a 21 M a m1
a 12 ..... a 1n a 22 ......a 2n a m2 .....a mn
Rodrigo A. Ocoró
x 1 b 1 x b 2 = 2 M M x m b m
31
(El número de incógnita debe ser igual al número de ecuaciones)
⇒ AX=B
Donde
A es la matriz de coeficientes del sistema X es el vector incógnita B es el vector de términos independientes.
3.3.1 Eliminación Gaussiana: Es el método más utilizado y consiste en la eliminación hacia adelante y la sustitución hacia atrás a partir de la matriz de coeficientes aumentada con el vector de términos independientes. La eliminación hacia adelante consiste en obtener una matriz triangular superior. Ejemplo 1: Resolver por eliminación de Gauss el sistema: 4x1 – 9x2 + 2x3 = 5 2x1 – 4x2 + 6x3 = 3 x1 – x2 + 3x3 = 4
⇒
4 − 9 2 x 1 5 2 − 4 6 x = 3 2 1 − 1 3 x 3 4
La matriz aumentable será: 4 − 9 2 5 2 − 4 6 3 1 − 1 3 4
i)
La matriz de coeficientes la convierto en triangular superior:
a 21 2 − Multiplico por 4 la 1ª fila y se lo sumo a la 2ª (Elemento a 11 )
Multiplico por
−
a 31 1 4 la 1ª fila y se lo sumo a la 3ª (Elemento a 11 )
Métodos Numéricos
4 − 9 2 5 5 0.5 0 0.5 0 1.25 2.5 2.75
ii)
Rodrigo A. Ocoró
akj = akj -
a ki aii
32
aij
Multiplica la 2ª fila por –2.5 y se la sumo a la 3ª.
4 − 9 2 5 0 0.5 5 0.5 0 0 - 10 1.5 Hasta aquí el proceso de eliminación hacia delante.
Ahora veamos la sustitución hacia atrás: El anterior sistema matricial se puede representar así: 4x1 – 9x2 + 2x3 = 5 0.5x2 + 5x3 = 1 / 2 -10x3 = 3 / 2 ⇒
(1) (2) (3)
1.5 Ecuación 3: x3 = − 10 = –0.15
x3 = –0.15
Ecuación 2: sustituyo x3: 0.5x2 + 5 (–0.15) = 0.5 0.5x2 = 0.5 + 0.75 = 1.25 1.25 x2 = 0.5
x2 = 2.5
Ecuación 1: se sustituye x2 y x3: 4x1 – 9 (2.5) + 2 (–0.15) = 5 4x1 = 5 + 22.5 + 0.30 = 27.8 27.8 x1 = 4
x1 = 6.95
Solución: [ 6.95 2.5 -0.15 ] NOTA: El determinante de la matriz lo podemos establecer fácilmente en la matriz triangular multiplicando los elementos de la diagonal. ⇒
Det. A = 4x 0.5x (–10) = –20 Si el det A = 0 ⇒ no es posible establecer la solución única.
Métodos Numéricos
Rodrigo A. Ocoró
3.3.2 Efecto de las cifras de redondeo en la eliminación Gaussiana: Resolver por Gauss: 3x1 –0.1x2 –0.2x3 = 7.85 0.1x1 + 7x2 –0.3x3 = –19.3 0.3x1 –0.2x2 + 10x3 = 71.4
Cifras significativas = número de cifras a partir de la 1ª cifra no nula.
a) Usando 6 cifras significativas:
3 − 0.1 - 0.2 7.85 7 - 0.3 - 19.3 0.1 0.3 - 0.2 10 71.4 ⇒
- 0.1
- 0.3
= - 0.0333333 3 = -0.1 1ª fila * ( - 0.0333333) + 2ª fila 1ª fila * (-0.1) + 3ª fila 3
3 - 0.2 7.85 − 0.1 0 7.00333 - 0.293333 - 19.56170 0 - 0.19 10.02 70.615
0.19 7.003333 =0.0271300 2ª fila * 0.0271300 + 3ª fila
3 - 0.2 7.85 − 0.1 0 7.00333 - 0.293333 - 19.56170 0 0 10.0120 70.0843
70.0843 x3 = 10.0120
⇒ x3 = 7.00003
− 19.56170 + 0.2933333( 7.00003)
x2 =
7.00333
7.85 + 0.1 (-2.50000) + 0.2 (7.00003) 3 x1 =
⇒
⇒
Evaluamos en cada ecuación y calculamos el error relativo:
(1) (2) (3)
7.85 − 7.84999 7.85 da 7.84999 ⇒ ER = =0.0001%
da –19.3000 ⇒ ER = 0% 71.4 − 71.4003 71.4 da 71.4003 ⇒ ER = = 0.0004%
x2 = - 2.50000 x1 = 3.00000
33
Métodos Numéricos
Rodrigo A. Ocoró
34
b) Usando 3 cifras significativas: 3x1 –0.1x2 –0.2x3 = 7.85 0.1x1 + 7x2 –0.3x3 = –19.3 0.3x1 –0.2x2 + 10x3 = 71.4 3 − 0.1 - 0.2 7.85 7 - 0.3 - 19.3 0.1 0.3 - 0.2 10 71.4
= - 0.033 - 0.3 3 =- 0.033 3 =- 0.100 1ª fila * ( - 0.033) + 2ª fila 1ª fila * ( - 0.100) + 3ª fila - 0.1
3 − 0.1 - 0.2 7.85 0 7.00 - 0.293 - 19.6 0 - 0.210 10.0 70.6
0.210 7.0 =0.03 2ª fila * (0.03) + 3ª fila
3 − 0.1 - 0.2 7.85 0 7.00 - 0.293 - 19.6 0 0 9.99 70.1 ⇒ x 3 =
70.1 9.99
⇒
x3 = 7.02
⇒
x2 = -2.51
x2 =
− 19.6 + 0.293(7.02) 7
x1 =
7.85 + 0.2 (7.02) + 0.1 (-2.51) 3
⇒
x1 = 3
Evaluamos en cada ecuación: 7.84 − 7.85 7.85 (1) da 7.84 ⇒ ER = =0.1% - 19.4 + 19.3 19.3 (2) da –19.4 ⇒ ER = = 0.5% 71.6 − 71.4 71.4 (3) da 71.6 ⇒ ER = = 0.3%
Luego el error de redondeo influye notoriamente en los resultados finales, pues un error en los primeros pasos tiende a propagarse.
Métodos Numéricos
Rodrigo A. Ocoró
35
3.4 MÉTODO DE GAUSS-JORDAN Es una variación de la eliminación Gaussiana y consiste en convertir la matriz de coeficiente en una matriz identidad en vez de una matriz triangular por consiguiente no es necesario emplear la sustitución hacia atrás para obtener la solución. Cada número de la diagonal debe convertirse en 1 y el resto de elementos de la columna deben ser ceros a 11 a 12 a 21 a 22 a 31 a 32
a 13 C1 a 23 C 2 a 33 C 3
1 0 0 C 1 * 0 1 0 C 2 * 0 0 1 C 3 *
⇒
x1
= C1 *
x2 ⇒
= C2
* x 3 = C3 *
Ejemplo resolver por Gauss-Jordan: 4x1 – 9x2 + 2x3 = 5 2x1 – 4x2 + 6x3 = 3 x1 – x2 + 3x3 = 4
⇒
4 - 9 2 5 2 - 4 6 3 1 - 1 3 4
Se pueden intercambiar filas de modo que la de máximo valor absoluto quede de 1ª (aquí no es necesario) se convierte al a11 en 1 y se hace operación filas para convertir a21 y a31 en 0. Luego en la fila a 22 se convierte en 1 ….. hasta obtener: 1 0 0 6.95 x1 x2 0 1 0 2.5 0 0 1 - 0.15 ⇒ x3
= = =
6.95 2.5 − 0.15
TAREA: Consultar cómo se establece la inversa de una matriz.
Métodos Numéricos
Rodrigo A. Ocoró
36
3.5 PROBLEMAS SIN SOLUCIÓN ÚNICA No siempre es posible resolver un conjunto de ecuaciones lineales en forma numérica: Por ejemplo: 1)
-x + y = 1 -2x + 2y = 2
2) -x + y = 1 -x + y = 0
3) -x + y = 1 x + 2y = -2 2x – y = 0
- Son rectas infinitas - El número de soluciones es infinito
1
Las rectas son linealmente dependientes
-1
Son paralelas No existen solución
-2 -1 -1
No hay solución única para las 3 ecuaciones y esto es debido a que hay más ecuaciones que incógnitas.
Sólo se encuentra la solución en pares de ecuaciones.
Métodos Numéricos
Rodrigo A. Ocoró
37
3.6. METODOS ITERATIVOS Al resolver un sistema de ecuaciones lineales por eliminación, la memoria de maquina requerida es proporcionar al cuadrado del det A, y el trabajo computacional es proporcional al cubo del orden de la matriz coeficiente A. Debido a esto, la solución de sistemas lineales grandes (n ≥ 50), se vuelve costoso y difícil en una computadora con los métodos de eliminación, ya que se requiere ampliar memoria; además, como él numero de operaciones que debe ejecutar es muy grande, se pueden producir errores de redondeo también muy grandes. Sin embargo, se han resuelto sistemas de orden 1000, y aun mayor, con los métodos que se estudiaran en esta sección. Estos sistemas de un numero muy grande de ecuaciones se presentan en la solución numérica de ecuaciones difereciales parciales, en la solución de los modelos resultantes en la simulación de columnas de destilación, etc. A favor de estos sistemas, puede decirse que tienen matrices con pocos elementos distintos de cero y que estas poseen ciertas propiedades (simétricas, bandeadas, diagonal dominantes, etc.), que permiten garantizar el éxito en la aplicación de los métodos de esta sección. 3.6.1. METODOS DE JACOBI Y GAUSS-SEIDEL Los métodos iterativos más sencillos y conocidos son una generación del método de punto fijo, estudiado en él modulo 2. Se puede aplicar la misma técnica a fin de elaborar métodos para la solución de A x = b, de la manera siguiente. Se parte de A Χ=b para obtener la ecuación
A x - b = 0,
(1 )
Ecuación vectorial corresponde a f (x) = 0. Se busca ahora una matriz B y un vector c, de (2 ) manera que la ecuación vectorial x = B Χ + c, sea solo un arreglo de la ecuación (1); es decir, de manera que la solución de una sea también la solución de la otra. La ecuación (2) correspondería a x=g(x). A continuación se propone un vector inicial x(0) como primera aproximación al vector solución x. Luego, se calcula con la ecuación (3) la sucesión vectorial x(1) , x(2),... de la siguiente manera donde
X (x+1) = B x (k) + c, k = 0, 1, 2, ... X (k) = X1 k X 2 k ...Xn k
(3)
T
Para que la sucesión X(0) , X (1), ..., X (n),..., converja al vector solución X es necesario que eventualmente X j m, 1 ≤ j ≤ n (los componentes del vector X (m), se aproximen tanto axJ,1 ≤ j ≤ n (los componentes correspondientes a x) que todas las diferencias x j m - x j , 1 ≤ j ≤ n sean menores que un valor pequeño previamente fijado, y que se conserven menores para todos los vectores siguientes de la iteración; es decir, Lim x jm = x j 1≤ j ≤ n m ∞
Métodos Numéricos
Rodrigo A. Ocoró
38
La forma como se llega a la ecuación (2) define el algoritmo y su convergencia. Dado el sistema A x = b, la manera más sencilla es despejar x 1 de la primera ecuación, x2 de la segunda, etc. Para ello, es necesario que todos los elementos de la diagonal principal de A, por razones obvias, sean distintos de cero. Para ver esto en detalle considérese el sistema general de tres ecuaciones (naturalmente puede extenderse a cualquier numero de ecuaciones). Sea entonces,
a11 x1 + a12 x2 + a13 x3 = b1 a21 x1 + a22 x2 + a23 x3 = b2 a31 x1 + a32 x2 + a33 x3 = b3
con a11, a22 y a33 distintos de cero. Se despeja x1 de la primera ecuación, x 2 de la segunda y x3 de la tercera con lo que se obtiene. X1 =
X2 = -
X3 = -
a12
a21
a33
x2 -
x1
a22 a31
a11
a13 a11
-
x1 -
a32 a33
a23 a11
b1
x3 +
a11 b2
x3 +
x2
a22
+
b3 a33
que en notación matricial queda X1 X2 X3
a12
0 a21 = -a22 a31 a33
X =
BX
a11 0 -
a32 a33 + C
-
a13 a11 a23 a22 0
b1
X1 X2 X3
+
a11 b2 a22 b3 a33
Métodos Numéricos
Rodrigo A. Ocoró
a12
0 -
B=
-
-
-
a11
a21
0
a22 a31
-
a32 a33
a33
a13
b1
a11
a11
a23
b2
a22
C=
39
a22 b3
0
a33
Una vez que se tiene la forma en notación matricial, se propone un vector inicial x (0) que puede ser x(0) = 0, o algún otro que sea aproximado al vector solución x. Para iterar existe dos variantes. 3.6.1.1. Iteración de Jacobi (método de desplazamiento simultáneos) Xk1 Xk) = Xk2 Xk3
Si
es el vector aproximación a la solución x después de k iteraciones, entonces se tiene para la siguiente aproximación. 1
X1k+1
(k+1)
X
=
X2k+1
a11 1
=
a22 1
X3k+1
a33
(b1 – a12 x k2 – a13 Xk3)
(b2 – a21 xk1 – a13 Xk3)
(b3 – a31 x k3 – a13 Xk3)
O bien, para un sistema de n ecuaciones con n incógnitas y usando notación más compacta y de mayor utilidad en programación, se tiene
Xik+1 =
-
1 aii
- bi + ∑ aij x jk j=1 j ≠ i
para 1 ≤ i ≤ n
Métodos Numéricos
Rodrigo A. Ocoró
40
3.6.1.2. Iteración de Gauss-Seidel (método de desplazamientos sucesivos) En este método los valores que se van calculando en la (k + 1)-ésima iteración se emplean para calcular los valores faltantes de esa misma iteración; es decir, con x (k) se calcula x(k+1) de acuerdo con: x1k+1
(k+1)
X
a11 (b1 – a12 x2 – a13 x 3)
x2k+1
=
1
=
a22 1
x3k+1
a33
(b2 – a21 x1k+1 – a23 xk3)
(b3 – a31 x1k+1 – a32 x 2k+1)
O bien, para un sistema de n ecuaciones
xik+1 =
-
1
i -1
aii
- bi + ∑ aij x j j=1
k+1
n
+ ∑ aij x jk , para 1 ≤ i ≤ n j = i +1
Sugerencia: para los siguientes ejemplos, el empleo de una calculadora programable, atenuaría considerablemente el trabajo de los cálculos. 3.6.1.3. Ejemplo : Resuelva el siguiente sistema por los métodos de Jacobi y Gauss-Seidel. Con ε = 0,5% 4x1 - x2 =1 - x1 + 4x2 - x3 =1 - x2 + 4x3 - x4 = 1 - x3 + 4x4 = 1 Solución: Despejando x1 de la primera ecuación, x 2 de la segunda, etc., se obtiene: x1 = x2 /4 x2 = x1 /4 + x3 /4 x3 = x2 /4 x4 = + x3 /4
+¼ +¼ + x4 /4 + ¼ +¼
Métodos Numéricos
Rodrigo A. Ocoró
41
Vector Inicial: Cuando no se tiene una aproximación al vector solución, se emplea generalmente como vector inicial el vector cero, esto es X(0) = [ 0 0 0 0 ] T a) Por Método de Jacobi El cálculo de x(1) en el método de Jacobi se obtiene remplazando x(0) en cada una de las ecuaciones despejado anteriormente. x1 = x2 = 0/4 x3 = x4 =
0/4
+¼ + 0/4 +¼ 0/4 + 0/4 + ¼ + 0/4 +¼
y entonces x(1) = (¼ ¼ ¼ ¼) T Para calcular x(2) se sustituye x(1) en cada una de las ecuaciones despejadas. Para simplificar la notación se han omitido los superíndices. x1 = 1/16 + ¼ = 0.3125 x2 = 1/16 + 1/16 + ¼ = 0.3750 x3 = 1/16 + 1/16 + ¼ = 0.3750 x4 = + 1/16 + ¼ = 0.3125 A continuación se presenta los resultados de subsecuentes iteraciones, en forma tabular. K 0 1 2 3 4 5 6
x1k 0.0000 0.2500 0.3125 0.3438 0.3555 0.3604 0.3623
x2k 0.0000 0.2500 0.3750 0.4219 0.4414 0.4492 0.4524
x3k 0.0000 0.2500 0.3750 0.4219 0.4414 0.4492 0.4524
x4 0.0000 0.2500 0.3125 0.3438 0.3555 0.3604 0.3623
Tabla 3.1. Solución del ejemplo 1 por método de Jacobi Solución: X = 0.3626 0.4524 0.4524 0.3623 b) Por Método de Gauss-Seidel. Para el cálculo del primer elemento del vector x (1), se sustituye x(0) en la primera ecuación despejada; para simplificar la notación se ha omitido los superíndices. x1 = 0/4 + ¼ = ¼
Métodos Numéricos
Rodrigo A. Ocoró
42
Para el cálculo de x2 de x(1), se emplea el valor de x1 ya obtenido (1/4) y los valores de x 2, x3 y x4 de x(0). Así: 1 x2 = + 0/4 + ¼ = 0.3125 4 (4) Con los valores de x1 y x2 ya obtenidos y con x3 y x4 de x(0) se evalúa x3 de x(1) X3 = 0.3125/4 + 0/4 + ¼ = 0.3281 Finalmente, con los valores de x 1, x2 y x3 calculados previamente y con x4 de x(0), se obtiene la última componente de x (1) x4 = 0.3281/4 + ¼ = 0.3320 Entonces x(1) = [ 0.25 0.3125 0.3281 0.3320 ]T Para la segunda iteración (cálculo de x (2) ) se procede de igual manera. x1 = 0.3125/4 + ¼ = 0.3281 x2 = 0.3281/4 + 0.3281/4 + ¼ = 0.4141 x3 = 0.4141/4 + 0.3320/4 + ¼ = 0.4365 x4 = 0.4365/4 + ¼ = 0.3591 Con lo que x
(2)
= [ 0.3281 0.4141 0.4365 0.3591 ]
T
En la tabla 3.2. se presentan los resultados de las iteraciones subsecuentes. K 0 1 2 3 4 5
x1k 0.0000 0.2500 0.3281 0.3335 0.3619 0.3633
x2k 0.0000 0.3125 0.4141 0.4475 0.4534 0.4544
x3k 0.0000 0.3281 0.4365 0.4517 0.4541 0.4545
x4K 0.0000 0.3320 0.3591 0.3629 0.3635 0.3636
Tabla 3.2. Solución del ejemplo 1 por el método de Gauss-Seidel. Solución: X = [ 0.3633 0.4544 0.4545 0.3636 ] Si se tomará el criterio con ε = 10 –2 y el método de Jacobi, se satisface en la sexta iteración de la tabla 3.1.; en cambio si ε = 10 –3, se necesitan de 10 iteraciones. Si se toma ε =10-2, el método de Gauss-Seidel , se requerirían sólo cinco iteraciones, mientras que con ε = 10 –3 se requieren 6 iteraciones.
Métodos Numéricos
Rodrigo A. Ocoró
43
Aunque hay ejemplos en los que Jacobi converge y Gauss-Seidel diverge y viceversa, en general puede esperarse convergencia más rápida por Gauss-Seidel, o una manifestación más rápida de divergencia. Esto se debe al hecho de ir usando los valores más recientes de x(k+1) que permitirán acercarse o alejarse más rápidamente de la solución. 3.6.1.4. Reordenaciones de Ecuaciones. Para motivar el rearreglo de ecuaciones, se propone resolver el siguiente sistema con el método de Gauss-Seidel y con ε = 0.1% aplicado a x (k+1) - x(k) donde: x(k+1) – x(k)
∑(xi k+1 – xik )2
=
< ε
- x1 + 3x2 + 5x3 + 2x4 = 10 x1 + 9x2 + 8x3 + 4x4 = 15 x2 + x4 = 2 2x1 + x2 + x3 - x4 = -3 Al resolver para x1 de la primera ecuación, para x 2 de la segunda, x3 de la cuarta y x4 de la tercera se obtiene. x1 = x2 = - x1 /9 x3 = -2x1 x4 =
3x2 + - x2 - x2
5x3 + 2 x4 - 10 - (8/9)x3 - (4/9)x4 + 15/9 + x4 - 3 + 2
Con el vector cero como vector inicial, se tiene la siguiente sucesión de vectores. Nótese que el proceso diverge. k 0 1 2 3
X1k 0.000 -10.000 67.8889 -631.1
x2k 0.000 2.7778 -18.172 170.7
x3k 0.000 14.222 -121.2 1108.0
x4K) 0.000 -0.7778 20.17 -168.71
x4(k+1) – x(k) 17.62 159.0 1439.05
Pero si reordenamos el sistema de ecuaciones tal que aii sea el mayor valor absoluto en la fila, entonces converge. Para ello, se toma la fila que tenga el pivote de mayor valor absoluto con relación a los demás coeficientes de la respectiva fila. (DIAGONAL DOMINANTE: aii ≥ aiJ El ejemplo anterior el reordenar las ecuaciones queda así: 2x1 + x2 + x3 - x4 = - 3 x1 + 9x2 + 8x3 + 4x4 = 15 - x1 + 3x2 + 5x3 + 2x4 = 10 x2 + x4 = 2
Métodos Numéricos
Rodrigo A. Ocoró
44
y la solución se presenta en el siguiente esquema: X2K
X3K
0 0.0000
0.0000
0.0000
0.0000
0.0000
1 -1.5000
1.9412
0.7073
2.0000
3.2432
2 -1.8242
0.2582
1.0000
2.0000
1.7387
3 -1.1291
-0.0000
1.0000
2.0000
0.7415
4 -1.0000
-0.0000
1.0000
2.0000
0.1291
5 -1.0000
-0.0000
1.0000
2.0000
0.0000
k
X1K
X4K
| k k+1 – X k |
La solución es X= [ -1 0 1 2] En la aplicación de estas dos variantes son válida las preguntas siguientes: 1. ¿La sucesión de vectores x(1), x(2), x(3),..., converge o se aleja del vector solución x= [ x1 x2...xn ] T ? 2. ¿Cuándo se detendrá el proceso iterativo? Las respuestas correspondientes, conocidas como criterio de convergencia, se dan a continuación. 1.
Si la sucesión converge a x, cabe esperar que los elementos de x(k) se vayan acercando a los elementos correspondientes de x, es decir, x1k, a x1, x2k a x2, etc., o que se alejen en caso contrario.
a). Los valores absolutos x1 k+1 - x1k , x2 k+1 - x2k , etc., sean todos menores de un número pequeño ε cuyo valor será dado por el programador. b). Si el número de iteraciones ha excedido un máximo predeterminado MAXIT. Si, por el contrario, esta sucesión de números diverge, entonces puede pensarse que el proceso diverge. Con esto, un criterio más es c)Detener el proceso una vez que x(k+1) – x(k) =V ∑(xi k+1 – xik )2 < ε Al elaborar un programa de cómputo para resolver sistemas de ecuaciones lineales, generalmente se utilizan los criterios (a), (b) y (c) o la combinación de (a) y (b), o la (b) y (c).
Métodos Numéricos
Rodrigo A. Ocoró
45
3.6.1.4 .Problema de aplicación: aplicación: Un granjero desea preparar una fórmula alimenticia para engordar ganado. Dispone de maíz, desperdicios, pasto y cebada; cada uno con ciertas unidades de ingredientes nutritivos de acuerdo acuerd o a la siguiente tabla: Unidades de ingredientes nutritivos por Kg. de cada alimento alimento MAIZ DESPERDICI PASTO CEBADA Requerimiento Ingr. nutritivo O Diario(unid / kg) Carbohidrato 80 15 32 60 230 s Proteínas 28 72 57 25 180 Vitaminas 20 12 20 20 80 Celulosa 50 10 20 60 160 COSTO $ 18 5 7 20 a)Determine los Kg. necesarios de cada alimento para p ara satisfacer el requerimiento diario. b)Determine el costo de la mezcla.
Métodos Numéricos
Rodrigo A. Ocoró
46
ALGORITMO Multiplicación de matrices Para multiplicar las matrices A y B, proporcionar proporcionar : DATOS:
Número de filas y columnas de A y B; N, M, N1, M1, respectivamente y sus elementos.
RESULTADOS:
La matriz producto C de dimensiones N x M1 o el mensaje "LAS MATRICES A y B NO PUEDEN MULTIPLICARSE".
PASO 1. PASO 2. PASO 3.
Si M = N1 continuar, de otro modo IMPRIMIR "LAS MATRICES A y B NO SE PUEDEN MULTIPLICAR" y TERMINAR. Hacer I = 1 Mientras ≤ N, repetir los pasos 4 a 12. PASO 4. PASO 5.
Hacer J = 1 Mientras J ≤ M1, repetir los pasos 6 a 11. PASO 6. Hacer C (I, J) = 0 PASO 7. Hacer K = 1
PASO 8.
PASO 11.
Mientras K ≤ M, repetir los pasos 9 y 10. PASO 9. Hacer C (I, J) = C (I, J) + A(I, K) * B(K, J) PASO 10. Hacer K = K + 1 Hacer J = J + 1
PASO 12. Hacer I = I + 1 PASO 13.
IMPRIMIR las matrices A, B y C y TERMINAR.
Métodos Numéricos
Rodrigo A. Ocoró
47
ALGORITMO Eliminación de Gauss Para obtener la solución de un sistema de ecuaciones lineales A x = b y el determinante de A, proporcionar DATOS:
N número de ecuaciones, A matriz coeficiente y b vector de términos independientes. RESULTADOS: El vector solución x y el determinante de A o mensaje de falla "HAY UN CERO EN LA DIAGONAL PRINCIPAL". PASO 1. Hacer DET = 1 PASO 2. Hacer I = 1 PASO 3. Mientras I ≤ N - 1, repetir los pasos 4 a 14. PASO 4. Hacer DET = DET * A (I,I) PASO 5. Si DET DET = 0 IMPRIMIR mensaje "HAY UN CERO EN LA LA DIAGONAL PRINCIPAL" y TERMINAR. De otro modo continuar. PASO 6. Hacer K = I + 1 PASO 7. Mientras K ≤ N, repetir los pasos 8 a 13. PASO 8. Hacer J = I + 1 PASO 9. Mientras J ≤ N, repetir los pasos 10 y 11. PASO 10. Hacer A(K, J) = A (K, J) -A(K, I)* A (I, J)/ (I, I) PASO 11. Hacer J = J + 1 PASO 12. Hacer b (K) = b (K) - A(K, I) * b ( I ) / A (I, I)I) PASO 13. Hacer K = K + 1 PASO 14. Hacer I = I + 1 PASO 15. Hacer DET = DET * A (N, N) PASO 16. Si DET = 0 IMPRIMIR mensaje "HAY UN CERO EN LA DIAGONAL PRINCIPAL" y TERMINAR. De otro modo continuar. PASO 17. Hacer x (N) = b (N) / A (N, N) PASO 18. Hacer I = N - 1 PASO 19. Mientras I ≥ 1, repetir los pasos 20 a 26. PASO 20. Hacer x ( I ) = b ( I ) PASO 21. Hacer J = I + 1 PASO 22. Mientras J ≤ N, repetir los pasos 23 y 24. PASO 23. Hacer X ( I ) = x ( I ) - A ( I, J) * x ( J ) PASO 24. Hacer J = J + 1 PASO 25. Hacer x ( I ) = x ( I ) / A (I, I) PASO 26. Hacer I = I - 1 PASO 27. IMPRIMIR x y DET y TERMINAR.
Métodos Numéricos
Rodrigo A. Ocoró
48
ACTIVIDADES Nº 3.
1.Elabore un programa de computación para sumar y multiplicar matrices y aplíquelo considerando las siguientes matrices: A=
1 2 -1 3 0 4 -2 1 3 -1 1 1
Efectuar: a) A+C
0 B= 1 1 2 b)C-A
1 0 2 -1 -1 3 -2 1
0 1 4 2 C= -1 -2 3 1 0 1 2 -1
c)A-2C
d)A.B
e)B.A
-1.2 X= 3.7 -2.0
f)A.C g)B.X
h)C.X
2.Calcule la inversa de c/u de las siguientes matrices y utilizando su programa de multiplicar matrices, compruebe que el producto de la matriz por su inversa da como resultado la matriz identidad: 1 -1 0 0 1 4 5 0 5 1 A= -1 2 -1 0 B= 2 1 2 C= -1 6 3 0 -1 2 -1 8 1 1 3 -9 5 0 0 -1 2 3.Resuelva por Eliminación Gaussiana: a) 2x1 + x2 – 3x3 = -1 -x1 +3x2+2x3 =12 3x1+ x2 –3x3= 0 3x2 – 13x3 = -50 b) 2x1 – 6x2 +x3 =44 4x1 - 8x3 = 4
b) 0.1x1 –0.6x2 + x3 = 0 -2x1 + 8x2 +0.3x3 =1 x1 + 6x2 + 4x3 =2 d)
2x2 + x4 = 0 2x1+ 2x2 +3x3 +2x4 = -2 4x1- 3x2 + x4 = -7 6x1+ x2 –6x3 –5x4 = 6
4.Muestre que el siguiente sistema de ecuaciones no tiene solución única: 3x1 + 2x2 – x3 – 4x4 = 10 x1 - x2 + 3x3 – x4 = -4 2x1+ x2 – 3x3 = 16 - x2 + 8x3 – 5x4 = 3 5.Utilice el programa que desarrolló para el método de Eliminación Gaussiana y resuelva: a)el ejercicio 3d. c)
b)el ejercicio 4.
0.5x1 +0.25x2 = 0.32 0.3x1 +0.8x2 +0.43x3 = 0.77 0.2x2 + x3 +0.6x4 = -0.6 x3 - 3x4 = -2
d)
4x1 + x2 = -1 -8x1 – x2 + x3 = 13 3x2 –2x3 +4x4 = -3 x3 - x4 + x5 = 2.1 2x4+ 6x5 = 3.4
Métodos Numéricos
Rodrigo A. Ocoró
49
6. Utilice cualquiera de métodos iterativos (Jacobi o Gauss-Seidel) y con ε =0.5% resuelva: a) 5x1 – x2 = 9 - x1 +5x2 – x3 = 4 -x2 +5x3 =- 6
b) 4x1 + x2 – x3 = -8 2x1 +5x2 =5 3x1 +8x2 +9x3=0
c) 2x1+ x2 – 8x3 = -15 6x1- 3x2 + x3 = 11 x1 –7x2 + x3 = 10
d) x1 –3x2+12x3 = -86 4x1 –2x2 – x3 = 39 x1 -6x2 +2x3 = -28 7. Desarrolle el algoritmo y un programa de computación para el método de Guss-Seidel y aplíquelo para resolver con ε =0.1% los siguientes sistemas: a) 5.8x1 + 3.2x2 +11.24x3 =20.24 4.3x1 + 3.4x2 +9.625x3 =17.325 2.5x1 + 5.2x2 +9.625x3 =17.325
b) 4x1+ 5x2+ 2x3 – x4 =3 5x1+ 8x2+ 7x3 +6x4 =2 3x1+ 7x2 – 4x3 –2x4=0 - x1+ 6x2 - 2x3 + 5x4=1
b) 2.156x1 + 4.102x2 – 2.3217x3 + 6x4 =18 -4.102x1 + 6 x2 +1.2x4=6.5931 - x1 - 5.7012x2+ 1.222x3 =3.4 6.532x1 + 7x2 -4x4 =0 8. Un ingeniero industrial supervisa la producción de cuatro tipos de computadoras. Se requieren 4 clases de recursos en la producción (horas/hombre, metales, plásticos y componentes electrónicos). En el cuadro se resumen las cantidades necesarias para cada uno de estos recursos en la producción de cada tipo de computadora. Si se dispone diariamente de 504 horas/hombre, 1970 kg. de metal, 970 kg de plástico y 601 componentes electrónicos, cuántas computadoras de cada tipo se pueden construír por día?. Cuadro: Recursos necesarios para producir cuatro tipos de computadoras Compudora horas/hombre metales plásticos componentes kg/computad. kg/computad. kg/computad. kg/computad. 1 3 20 10 10 2 4 25 15 8 3 7 40 20 10 4 20 50 22 15 Respuesta: [10 12 18 15] Si las ganancias correspondientes a cada tipo de computadora en dólares es:
Métodos Numéricos
Rodrigo A. Ocoró
50
500, 370, 650 y 420 respectivamente, Cuál sería la ganancia total en esta producción?.
9. Un ingeniero supervisa la producción de tres tipos de automóviles. Se requieren tres clases de materiales para la producción: metal, plástico y caucho. La cantidad necesaria para producir cada automóvil es: Automóvil 1 2 3
Metal Kg/auto 1500 1700 1900
Plástico Kg/auto 25 33 42
Caucho Kg/auto 100 120 160
Si se dispone diariamente de un total de 106 toneladas de metal, 2.17 toneladas de plástico y 8.2 toneladas de caucho. Cuántos automóviles se pueden producir diario?. 10. Un ingeniero requiere 4800 m3 de arena, 5810 m3 de grava fina y 5690 m3 de grava gruesa para la construcción de un proyecto. Existen tres bancos de materiales donde se pueden obtener estos implementos. La composición de cada banco es: Banco 1 2 3
arena % 52 20 25
grava fina % 30 50 20
grava gruesa % 18 30 55
Cuántos metros cúbicos se debe tomar de cada banco para cumplir con las necesidades del ingeniero?.
4. AJUSTE DE CURVAS
(Aproximación polinomial e interpolación)
4.1 INTRODUCCIÓN Los datos que se obtienen mediante mediciones fluctúan, esto se debe a errores aleatorios del sistema de medición aplicada al comportamiento del sistema en observación. Cualquiera que sea la razón, es frecuente que surja la necesidad de ajustar una función a los datos de una medición, por ejemplo, un economista desearía ajustar una curva a una tendencia económica para poder predecir el futuro. La interpolación es de gran importancia en el campo de la ingeniería ya que al consultar fuentes de información presentadas en forma tabular, es frecuente no encontrar el valor buscado como un punto en la tabla. Una función de interpolación es aquella que pasa a través de puntos dados como datos los cuales se muestran comúnmente por medio de una tabla de valores o se toman directamente de una función dada. La interpolación de los datos puede hacerse mediante un polinomio la cual consiste en ajustar a un polinomio los puntos dados. 4.2 REGRESIÓN LINEAL Consiste en ajustar a una línea recta un conjunto de parejas de datos observados: (X 1, Y1), (X2, Y2),………, (Xn, Yn) y para ello utilizaremos el método de mínimos cuadrados. La expresión matemática de la función lineal es g(x)=a 0 + a1x constantes a determinar.
donde a0 y a1 son
La desviación de la recta con respecto a cada dato se define como: di = yi – g(xi) y como g(xi) = a0 + a1 xi entonces di = yi – a0 – a1xi donde i = 1,2………,n siendo n número de datos El cuadrado total de las desviaciones para obtener un buen ajuste está dado por: n D=
2
∑ (di)
i=1
n
= ∑ (yi – a0 – a1xi)2 cuando este valor es mínimo. i=1
Para determinar los valores de las constantes a 0 y a 1 y hacer mínimo a D, se deriva la ecuación anterior con relación a cada uno de los coeficientes (derivados parciales) y se igualan a cero.
Métodos Numéricos
Rodrigo A. Ocoró
52
n
∂ D = −2Σ (y i − a 0 − a 1x i ) = 0 ∂ a 0 i=1 n
∂ D = −2 Σx i (y i − a 0 − a 1x i ) = 0 ∂ a1 i=1
∑yi - ∑a0 - ∑a1xi = 0
⇒
∑yixi - ∑a0xi -
∑a1xi2
n
donde ∑ = ∑ =0
i=1
n Considerando que ∑ a0 =na0, entonces las anteriores ecuaciones se pueden expresar i=1 como un conjunto de 2 ecuaciones lineales con 2 incógnitas: na0 + a1 ∑xi = ∑yi a0 ∑xi + a1 ∑xi2 = ∑xiyi
(Ecuaciones normales)
al expresar estas 2 ecuaciones en forma matricial, tenemos: n ∑xi
∑xi ∑xi2
a0 a1
∑yi ∑xiyi
=
Que se resuelven por los métodos conocidos (eliminación gaussiana o gauss-seidel) Ejemplo 1. Para el siguiente conjunto de datos observados, determine la línea de ajuste. X y
0.1 0.61
0.4 0.92
0.5 0.99
0.7 1.52
0.7 1.47
0.9 2.03
n=6 (Parejas (x,y))
Solución: Organicemos así: xi 0.1 0.4 0.5 0.7 0.7 0.9 3.3
∑
6 3.3
3.3 2.21
yi 0.61 0.92 0.99 1.52 1.47 2.03 7.54
a0 a1
=
7.54 4.844i
Xi2 0.01 0.16 0.25 0.49 0.49 0.81 2.21
xiyi 0.061 0.368 0.495 1.064 1.029 1.827 4.844
Resolviendo por eliminación Gausiana, Se tiene:
Métodos Numéricos
Rodrigo A. Ocoró
6
3.3
7.54
0
0.395
0.697
a1 =
a0 =
0.697 0.395
a 1 = 1.764557
7.54 - 3.3 (1.764557) 6
Luego la ecuación lineal de ajuste será:
53
a 0 = 0.2861603
g(x) = a0 + a1x g(x) = 0.2861603 + 1.764557x D= ∑ di 2 = ∑ ( g ( xi ) − yi ) 2
S=∑(yi - y )2
g(x)
..
2.0
.. ..
1.5 1.0
.
0.5
g(x)
0 0.5
x
1.0
i
xi
yi
1 2 3 4 5 6
0.1 0.4 0.5 0.7 0.7 0.9
0.61 0.92 0.99 1.52 1.47 2.03
g(xi)
di
2
2
(yi - y ) 0.4626 0.02172 0.4182 0.9917 0.00511 0.1134 1.1684 0.03183 0.0711 1.5214 0.00000 0.0693 1.5214 0.00264 0.0455 1.8743 0.02424 0.5979 D=0.08574 S=1.3155
Otra forma de hallar los valores de a 0 y a1 sin resolver la ecuación es: a1 =
nΣx i y i − Σx i y i nΣx 2i − ( Σx i ) 2
a 0 = y − bx
donde y =
Σy i n
x=
Σx i n
.2.1Cuantificación de error en regresión lineal: Se establece mediante el coeficiente de determinación r2 que indica el grado de correlación entre las variables, dependiendo del valor en que se ubique r2, así: (0 a 0.60) el ajuste es malo (0.60 a 0.80) el ajuste es moderado (0.80 a 0.90) el ajuste es bueno (0.90 a 0.99) el ajuste es excelente Si r2 = 1 el ajuste es perfecto 2
r =
S−D
n
S
n d on de
D = Σ (d i ) 2
(T ota l d e la s d esvia cio nes)
i=1
S = Σ (y i − y )2 (Variación total)
i=1
Para el ejemplo anterior
⇒ r 2 =
1.315533 − 0.086779 1.315533
S = 1.315533
= 0.9340
Por lo tanto el ajuste lineal realizado con g(x) = 0.2861603 + 1.764557x es excelente.
Métodos Numéricos
Rodrigo A. Ocoró
54
4.3. APLICACIÓN DE LA REGRESIÓN LINEAL La regresión lineal proporciona una técnica muy poderosa para ajustar datos a una “mejor” línea. Pero no siempre los datos se comportan en forma lineal pues en algunos casos son en forma curvilínea pero se pueden hacer las transformaciones para que sean compatibles con el ajuste.
...... . . . . .. . . .. .. . ... ... .. ...
4.3.1 El modelo exponencial g(x)=aebx donde a y b son constantes, lo podemos linealizar así: Lny = Ln(ae bx) Lny = Lna + bx Lne Lny = Lna + bx n ∑ x i
∑ x ∑ x
Ln a ∑ Ln y i = b ∑ ( x i * Ln y i )
i
2 i
4.3.2 El modelo potencial g(x)=ax b donde a y b son constantes, lo linealizamos así: Log y = Log (ax b) Log y = Log a + b Log x n ∑ Log x i
∑ Log x ∑ ( Log x i
i
)2
Log a ∑ Log y i = b Log x Log y * ∑ i i
Ejemplo 1: ajuste los siguientes datos mediante un modelo potencial: X y xi 1 2 3 4 5 ∑
yi 0.5 1.7 3.4 5.7 8.4
1 0.5 Log xi 0 0.3010 0.4771 0.6021 0.6990 2.0792
2 1.7 Log yi -0.3010 0.2304 0.5315 0.7559 0.9243 2.1411
3 3.4
4 5.7 2
(Log xi) 0 0.090601 0.227624 0.362524 0.488601 1.16930
5 8.4
Log xi . log yi 0 0.06935 0.253379 0.455127 0.644608 1.42414
n=5
G(xi) 0.5011 1.6874 3.4328 5.6819 8.3992
D
2
S
2 di (yi - y ) 0.0000 11.8336 0.0002 5.0176 0.0010 0.2916 0.0005 3.0976 0.0000 19.8976 0.0017 40.1320
Métodos Numéricos
n ∑logxi
Rodrigo A. Ocoró
∑logxi ∑(logxi)2
5 2.0792
log a b
=
2.0792 2.1411 1.16935 1.4241
De allí se obtiene: b =
0.53374 0.30474
∑log yi ∑log xi.logyi
resolviendo
= 17515 .
55
log a =
5 0
2.0792 0.30474
2.1411- 2.0792 (1.7515) 5
2.1411 0.53374
= −0.3001
La ecuación de regresión será: Logy= loga+b.log x => Log y = -0.3001 + 1.75 Log x a=antilog(-0.3001)
a=0.5010
Como g(x) = axb tomando antilog queda : g(xi) = 0.5010 X1.75 40.1320 − 0.0017 2
r = S-D S
2
r =
40.1320
= 0.99
Ajuste excelente
Ejemplo 2. Ajuste los datos del ejemplo anterior a un modelo exponencial. G(x)=aebx Ln y = ln a + bx 2 Xi yi ln yi xi2 xi.ln yi g(xi) di2 (yi - y ) 1 0.5 -0.6931 1 -0.6931 0.6806 0.0326 11.8336 2 1.7 0.5306 4 1.0612 1.3506 0.1221 5.0176 3 3.4 1.2238 9 3.6714 2.6802 0.5182 0.2916 4 5.7 1.7405 16 6.9620 5.3184 0.1456 3.0976 5 8.4 2.1282 25 10.6410 10.5537 4.6383 19.8976 _______________________________________________________________ 15 19.7 4.9300 55 21.6426 5.4568 40.1320 n ∑ x i
∑ x ∑ x i
2 i
Ln a ∑ Ln y i = b ∑ xi * Ln y i
5 15
Ln a 4.9300 = 55 b 21.6426
15
Que al resolver se obtiene: b= 0.68526 y Ln a= -1.0698 a=antiln(-1.0698) a=0.3430 G(x)= aebx g(x)= 0.3430e0.6853X r2=0.86 el ajuste exponencial es bueno
Métodos Numéricos
4.4.
Rodrigo A. Ocoró
56
REGRESIÒN LINEAL MÚLTIPLE
Una extensión útil en la regresión lineal es el caso en que y es una función lineal de dos o más variables. Por ejemplo, y pudiera ser una funciónlineal de x 1 y x2, de la forma: y = a0 + a1x1 + a2x2 Tal ecuación es útil particularmente cuando se ajustan datos experimentales en donde la variable que está analizando, a menudo es función de otras dos variables. En este caso bidimensional, la “línea” de regresión viene a ser un “plano” . Como con los casos anteriores, los “mejores” valores de los coeficientes se determinan agrupando la suma de los cuadrados de los residuos: n
D = ∑ (yi - a0 - a1x1 + a2x2)2 i=j
y derivando con respecto a cada uno de los coeficientes: ∂D1 = - 2 ∑ (yi - a0 - a1x1 - a2x2) ∂a0 ∂D1 ∂a1 ∂D1 ∂a2
= - 2 ∑ x1 (yi - a0 - a1x1 - a2x2)
(1)
= - 2 ∑ x2 (yi - a0 - a1x1 - a2x2)
Los coeficientes que generan la suma mínima de los cuadrados de los residuos se obtienen igualando cada una de las derivadas parciales a cero y expresando la ecuación (1) como un conjunto de ecuaciones lineales simultáneas, de la forma: nao + ∑ x1 a1 + ∑ x2 a2 = ∑ y1 ∑ x1ao + ∑ x12 a1 + ∑ x1 x2 a2 = ∑ x1 y1 ∑ x2ao + ∑ x1 x2 a1 + ∑ x22 a2 = ∑ x2 y1 o como una matriz: n
∑ x1 ∑ x2
∑ x1 ∑ x2 2 ∑ x1 ∑ x1 x2 ∑ x1 x2 ∑ x22
a0 a1 a2
=
∑ y1 ∑ x1 y1 ∑ x2 y1
(2)
Métodos Numéricos
Rodrigo A. Ocoró
57
4.4.1. Ejemplo: Dados los siguientes datos, determine el modelo de regresión lineal múltiple: x1
x2
y
0 2 2.5 1 4 7
0 1 2 3 6 2
5 10 9 0 3 27
Cálculos necesarios para desarrollar las ecuaciones normales del ejemplo y
x1
x12
x2
x22 x1 x2
5 0 0 0 0 10 2 1 4 1 9 2.5 2 6.25 4 0 1 3 1 9 3 4 6 16 36 27 7 2 49 4 ∑ 54 16.5 14 76.25 54
0 2 5 3 24 14 48
x1y
G(xi)
di2
4.91920 9.94012 8.98061 0.06547 3.17830 26.9109
0.00652 0.00358 0.00019 0.00428 0.03179 0.00794 0.0543
x2y
0 0 20 10 22.5 18 0 0 12 18 189 54 243.5 100
S
2
(yi - y )
16 1 0 81 36 324 458 D
Solución: Estas sumatorias necesarias para desarrollar la ecuación (1), se sustituyen en (2) y se obtiene:
6 16.5. 14 16.5. 76.25 48 14 48 54
a0 54 a1 = 243.5 a2 101
que se puede resolver usando un método como la eliminación Gaussiana para obtener: ao = 4.91922
a1 = 3.98329
a2 = -2.94568
g(x)=4.91922 +3.98329X1 -2.94568X2 −
y =
∑ yi n
=
54 6
=9
r 2 =
S − D S
=
458 − 0.0543 458
r 2 = 0.99
El ajuste es Excelente
Métodos Numéricos
Rodrigo A. Ocoró
58
La regresión lineal múltiple se puede formular en el caso más general como: y = a0 + a1x1 + a2x2 +……..+amxm en donde los coeficientes que minimizan la suma de los cuadrados de lo residuos se determinan resolviendo el sistema:
n
∑ x1,i ∑ x2,i
. . . ∑ xm,i
∑ x1,i ∑ x2,i . . . ∑xm,i 2 ∑ x 1,i ∑ x2,i x1,i . . . ∑ x1,ixm,i ∑ x2,i x1,i ∑x22,i . . . ∑x2,ixm,i
. . . ∑xm,i x1,i ∑xm,i x2,i
...
∑x2m,i,
a0
a1 a2 . . . am
∑ yi
∑ x1,i yi = ∑ x2,i yi
∑ xm,i yi
Métodos Numéricos
Rodrigo A. Ocoró
59
4.5. REGRESIÓN POLINOMIAL El método de mínimos cuadrados puede extenderse para ajustar un polinomio de cualquier valor a los datos de una medición. g(x) = a0 + a1x + a2x2 + amxm
di=yi – g(xi)
Desviación de la curva con res ecto a cada Polinomio.
La suma de los cuadrados de la desviación así: D = ∑ (yi – a0 – a1 xi – a2xi2 + ……… amxim)2 Siguiendo el mismo procedimiento hecho para la regresión lineal, tenemos: derivados parciales con respecto a los coeficientes al polinomio. ∂ D = −2Σ (y i − a 0 − a 1x i − a 2 x 2i ..........a mx m i ) ∂ a 0 ∂ D = −2 Σx i ( Σy i − a 0 − a 1x i − a 2 x 2 .......a m x 2m ) ∂ a1 ∂ D = −2Σx 2i ( Σy i − a 0 − a 1x i − a 2 x 2i .......a mx m i ) ∂ a 2
∂ D 2 m = −2Σx m i ( Σy i − a 0 − a 1x i − a 2 x i .......a m x i ) ∂ a m
Al igualar estas ecuaciones a cero y reordenando, se tiene el siguiente conjunto de ecuaciones. a0n + a1∑xi + a2∑xi2 + ……………+am ∑xim = ∑yi a0∑xi + ai∑xi2 + a2∑x23 + …………+ am∑xim+1 = ∑xiyi a0∑xi2 + a1∑xi3 + a2∑xi4 +…………+ am∑xim+2 = ∑xi2 yi a0∑xim + a1∑xim+1 + a2∑xim+2 +……+ am∑xi2m = ∑ximyi (Todos los ∑ van desde i =1 hasta n)
Métodos Numéricos
Rodrigo A. Ocoró
60
el anterior sistema lo podemos expresar así: ∑xi ∑ … xi2
∑xi ∑xi2 ∑ … xi3
∑xi2 ∑xi3 ∑…xi4
∑xim
∑xim+1
∑xim+2 …… ∑xi2m
n .
.
.
…… ∑xim …… ∑xim+1 ∑… xim+2
a0 a1 a…2
.
∑yi ∑xiyi ∑x…i2yi
=
.
.
∑ximyi
am
El sistema de ecuaciones se resuelve por los métodos vistos para resolver sistema de ecuación. Ojo: si el polinomio es de grado m => necesitamos m + 1 ecuaciones lineales con m+1 incógnitas. 4.5.1.Ejemplo: ajuste a un polinomio de 2º grado el siguiente conjunto de datos: x y
0 2.1
1 7.7
2 13.6
3 27.2
Como el polinomio de ajuste es de grado 2 (a0, a1, a2)
4 40.9
5 61.1
grafique
necesitamos 3 ecuaciones y 3 incógnitas
Organicemos la siguiente tabla: yi
xi2
xi3
xi4
Xiyi
xi2yi
2.1 7.7 13.6 27.2 40.9 61.1 152.6
0 1 4 9 16 25 55
0 1 8 27 64 125 225
0 1 16 81 256 625 979
0 7.7 27.2 81.6 163.6 305.5 585.6
0 7.7 54.4 244.8 654.4 1527.5 2488.8
xi 0 1 2 3 4 5 ∑ 15 n=6 n ∑xi ∑xi2
∑xi ∑xi2 ∑xi3
∑xi2 ∑xi3 ∑xi4
a0 a1 a2
∑yi ∑xiyi ∑xi2yi
=
Reemplazando y resolviendo por curvas: 6 15 55
15 55 225
55 225 979
152.6 585.6 2488.8
=
a0=2.47857 a1=2.35929 a2=1.86071
La ecuación cuadrática de ajuste con mínimo cuadrado es: g(x) = a0 + a1x + a2x2 == g(x) = 2.47857 + 2.35929x + 1.86071x2
Métodos Numéricos
Rodrigo A. Ocoró
61
Calculemos las desv. de cada punto; grafiquemos g(x) para la gráfica de dispersión y hallemos el coeficiente de determinación para mirar cómo es el ajuste con este método para dichos datos. di2 = [yi - g(x)]2 y=
Σy
n
=
152.6 6
xi 0 1 2 3 4 5 ∑ 15
= 25.433
yi 2.1 7.7 13.6 27.2 40.9 61.1 152.6
(yi – y)2 544.44 314.47 140.03 3.12 239.22 1272.11 2513.39
g(xi) 2.4786 6.6986 14.6406 26.3077 41.6882 60.7941 152.6038
(yi – g(xi)2 0.14332 1.00286 1.08158 0.80491 0.61951 0.09439 3.74657
S=∑ (yi – y)2 variación total D=∑ (yi – g(x))2 Total de los cuadrados de las desviaciones r 2 =
S−P S
=
2513.39 − 3.74657 2513.39
= 0.99851
r2 = 0.998 excelente (ajuste casi perfecto)
Métodos Numéricos
Rodrigo A. Ocoró
62
4.6. POLINOMIOS DE INTERPOLACIÓN La interpolación de los datos puede hacerse mediante un polinomio la cual consiste en ajustar a un polinomio los puntos dados. 4.6.1. POLIMONIO DE LAGRANGE El método de aproximación Polimonial requiere de la solución de un sistema de ecuaciones algebraicas lineales que cuando el grado del polinomio es alto, puede presentar inconvenientes. Con el método de Lagrange los cálculos se realizan directamente. Con frecuencia, se tienen que estimar valores intermedios entre valores conocidos y el método más común empleado para este propósito es la interpolación polinomial. La interpolación es una aplicación de la aproximación polinómica que permite el cálculo de algunos valores de una función y = f(x) para la que se conocen otros valores (x o, yo), (xi, yi), (x2, y2) .........., (xn,yn). Generalmente consiste en calcular y=f(x) para xE [ xk, xk+1] , suponiendo conocidos yk, yk+1 Recuerde que la fórmula general de un polinomio de n-ésimo orden es: n 2 3 n P(x) = a0 + a1x + a2x + a3x + ........ + anx = ∑ akxk k=0
El polinomio de interpolación consiste en determinar el único polinomio n-ésimo que ajusta a los n+1 puntos dados. Así por ejemplo hay solo una recta (polinomio de primer orden) que conecta dos puntos; hay solo una parábola que conecta a tres puntos (polinomio de segundo orden) y así sucesivamente. Este polinomio proporciona una fórmula para calcular valores intermedios.
(a)
(b)
(c)
Ejemplos de interpolación polinomial: a) Primer orden (lineal), conexión de dos puntos; b) conexión de tres puntos, segundo orden (cuadrática o parabólica) y c) conexión de cuatro puntos, tercer orden (cúbico) Existen gran variedad de fórmulas matemáticas mediante las cuales se puede expresar este polimonio entre ellas utilizaremos el método de Lagrange. Si dados en el plano cartesiano los puntos (x o, yo), (xi, yi), (x2, y2) .........., (xn,yn) el polimonio equivalente pasa por todos y cada uno de ellos. Observemos que para encontrar los parámetros a0, a1, a2,.....,an1 debe resolverse el siguiente sistema de ecuaciones lineales:
Métodos Numéricos
Rodrigo A. Ocoró
63
+ a1x0 + a2x02 + a3x03 + ...............+ arx0n = y0 a0 + a1x1 + a2x12 + a3x13 + ...............+ anx1n = y1 a0 + a1x2 + a2x22 + a3x23 + ...............+ anx2n = y2 ....................................................................=.... ....................................................................=.... a0 + a1xn + a2xn2 + a3xn3 + ...............+ anxnn = yn 0
con el objeto de evitar el desarrollo del sistema anterior, Lagrange demostró que este polinomio es equivalente con
p(x)
=
n
∑ yk.Lk(x),
cuando
se
encuentra
k =0
apropiadamente Lk(x) , para cada k.
Lk(x) =
(x-x0) (x-x1)........(x - xk-1) (x-xk+1).........(x-xn) (xk-x0) (xk-x1)......(xk - xk-1) (xk-xk+1).......(xk-xn)
n
(x-xJ) j=0 (xk-xJ)
= π
J ‡k
La expresión anterior dispone de n factores en el numerador que generan como resultado un polinomio Lk(x) de grano n. 4.6.1. Ejemplo 1: Encuentre el polinomio de colocación para los puntos (-1,4), (0,1) y (3,4) usando el método directo (sistema de ecuaciones lineales) y el método de Lagrange, interpole en x = 2 Como hay 3 puntos entonces, el polinomio será de 2º grado. El polinomio a buscar es P(x) = a 0 + a1x + a2x2 a)Método Directo. Reemplazando cada punto en P(x), obtenemos: (-1,4): a0 + a1 + a2 = 4 ( 0,1): a0 = 1 ( 3,4): a0 + 3a1 + 9a2 = 4 Resolviendo este Sistema por eliminació gaussiana obtenemos: a0 = 1, a1 = -2, a2 = 1 Con ello el polinomio de colocación resultante es P(x) = 1 –2x + x 2 b)Método de Lagrange: Los Lagrangianos resultan ser los siguientes: (xo, yo) = (-1,4) (x1, y1) = (-0,1) (x2, y2) = ( 3,4)
Métodos Numéricos
L0 (x) =
L1 (x) =
L1 (x) =
Rodrigo A. Ocoró
(x-x1) (x-x2)
L0(x) =
(x-0) (x-3)
=
(x0-x1) (x0-x2)
(-1-0) (-1-3)
(x-x0) (x-x2)
(x+1) (x-3)
L1(x) =
x2-3x 4
=
x2-2x - 3
(x1-x0) (x1-x2)
(0+1) (0-3)
-3
(x-x0) (x-x1)
(x+1) (x-0)
x2 +x
L2(x) =
(x2-x0) (x2-x1)
=
(3+1) (3-0)
12
El polinomio de Colocación se calcula como: 2 P(x) = ∑ ykLk (x) = y0L0(x) + y1L1(x) + y2L2(x) K=0
P(x)= 4
( x2 – 3x) 4
= x2 – 3x -
+1
( x2 – 2x – 3) -3
( x2 – 2x - 3 )
+
3
+ 4
3 3x2 – 6x +3
3
4.6.2. Ejemplo 2: i Xi F(xi)
0 0 -3
12
x2 – x
3x2 – 9x - x2 – 2x + 3 + x2 + x
p(x) = x2 – 2x + 1 p(x) = 22 – 2(2) +1 = 1
( x2 – x)
3
= x2 –2x +1
(2, 1)
Para la tabla que se presenta a continuación: 1 1 0
2 3 5
3 6 7
a) Obtenga el polinomio de Lagrange con todos los puntos. b) Interpole el valor de la función para x =1.8
64
Métodos Numéricos
Rodrigo A. Ocoró
65
4.6.3.Algoritmo para la regresión polinomial: Paso 1. Paso 2: Paso 3: Paso 4: Paso 5: Paso 6: Paso 7:
Introduzca el orden del polinomio a ser ajustado, m Introduzca el número de puntos, n y los puntos (xi , yi) si n ≤ m imprima: error el polinomio de regresión es imposible y terminar. Si n ≥ m continúe. calcule las sumas y producto de potencia (∑xi, ∑xi2...... ∑xi2m, ∑yi, ∑xiyi,...... ∑ximyi). arregle estas sumas y productos en forma de matriz aumentada. reducir la matriz suministrada por un método (llámela subrutina) imprimir resultado. ________________________________________________ 4.6.4. ALGORITMO Interpolación con polinomios de Lagrange
Para interpolar con polinomios de Lagrange de grado N, proporcionar los DATOS: RESULTADOS: PASO 1. PASO 2. PASO 3.
PASO 11.
El grado del polinomio N, las N + 1 parejas de valores (X (I ), FX (I), I = 0,1 ..., N) y el valor para el que se desea la interpolación XINT. La aproximación FXINT, el valor de la función en XINT.
Hacer FXINT = 0 Hacer I = 0 Mientras I ≤ N, repetir los pasos 4 a 10 PASO 4. Hacer L = 1 PASO 5. Hacer J = 0 PASO 6. Mientras J ≤ N, repetir los pasos 7 y 8 PASO 7. Si I ≠ J Hacer L = L*(XINT - X(J)/(X(I)-X(J)) PASO 8. Hacer J = J + 1 PASO 9. Hacer FXINT = FXINT + L*FX(I) PASO 10. Hacer I = I + 1 IMPRIMIR FXINT y TERMINAR.
Métodos Numéricos
Rodrigo A. Ocoró
ACTIVIDADES Nº 4 1.Para los siguientes datos en cada caso determine: i)La ecuación de ajuste lineal mediante el método de mínimos cuadrados ii)La gráfica de los puntos y la línea de ajuste iii)La cuantificación del error mediante el coeficiente de determinación a) Xi 1.0 1.5 2.0 2.5 3.0 Yi 2.0 3.2 4.1 4.9 5.9 b)
Xi 4 6 8 10 14 16 20 22 24 28 28 34 36 38 Yi 30 18 22 28 14 22 16 8 20 8 14 14 0.1 8
c) Xi 0.1 Yi 9.9
0.2 9.2
0.3 8.4
0.4 6.6
0.5 5.9
0.6 5.0
0.7 4.1
0.8 3.1
0.9 1.9
1.0 1.1
2.Ajuste a una ecuación de potencias: Xi 2.5 3.5 5 6 7.5 10 12.5 15 17.5 20 Yi 5 3.4 2 1.6 1.2 0.8 0.6 0.4 0.3 0.3 Grafique los datos y la ecuación de ajuste. 3.Ajuste a un modelo exponencial: Xi 0.05 0.4 0.8 1.2 1.6 2.0 2.4 Yi 550 750 1000 1400 2000 2700 3750 Grafique los datos y la ecuación de ajuste. 4.Mediante el modelo de regresión lineal múltiple, ajuste: a) X1 0 1 2 0 1 2 X2 2 2 4 4 6 6 Y 19 12 11 24 22 15 b) X1 1 X2 1 Y 18
1 2 2 3 3 4 4 2 1 2 1 2 1 2 12.8 25.7 20.6 35.0 29.8 45.5 40.3
5. Ajuste a un modelo cuadrático (parabólico) los datos del ejerccio 3. Y grafique los datos y la ecuación de ajuste. Para estos datos, cual ajuste considera mejor?.
66
Métodos Numéricos
Rodrigo A. Ocoró
67
6. Dados los siguientes datos: Xi 0.1 0.2 0.4 0.6 0.8 yi 3.25 7.78 10.68 8.37 3.97
1.0 2.5
Grafique los datos y ajuste a los siguientes modelos: a) lineal b)cuadrático c)cúbico d)exponencial e)Potencial. Cuál de ellos considera el mejor ajuste para dichos datos? (realice la gráfica de c/u en la gráfica de los datos y compare). Establezca el coefic. de determinación en c/u. 7. Dados los datos: Xi Yi
5 10 15 20 25 30 35 40 45 50 17 25 30 33 36 38 39 40 41 42
Úsese regresión con mínimos cuadrados para ajustar: a)a una línea recta b)a una parábola o cuadrático c)a una ecuación de potencia d)Grafique los datos con las curvas de ajuste y establezca cuál de ellas da mejor ajuste. 8. Desarrolle los respectivos programas para regresión lineal por mínimos cuadrados y para la regresión Polinomial y aplíquelo a los ejercicios 2, 5 y 6.
9. Al medir la velocidad (con un tubo de pitot) en una tubería circular de 20 cms de diámetro interior, se encontró la siguiente información:
V(cm/s) 600 550 450 312 240 r(cm)
0
3
5
7
8
Donde r es la distancia en cm medida a partir del centro del tubo a)Obtenga la curva v=f( r) que aproxima estos datos experimentales mediante un modelo parabólico. b)Calcule la velocidad en el punto r=4cm.
Respuesta: 503.89 cm/s
Métodos Numéricos
Rodrigo A. Ocoró
68
10. Se ha llevado a cabo un estudio para determinar la relación entre la fuerza de fricción que actúa hacia arriba y la velocidad de caída de un paracaidista obteníendose los siguientes resultados (v medida en cm/s y Fr medida en 106 dinas) v 1000 2000 Fr
5
3000 4000 5000
15.3
29.3 46.4 66.3
Grafíquese Fr contra v y mediante un modelo de regresión apropiado, determine la expresión que permita relacionar la fuerza de rozamiento y la velocidad. 11.Utilizando el método de Lagrange encuentre el polinomio de colocación para cada grado de puntos: a) (0,1) , (1, 0)
e interpole en x=0.8
b) (0,16), (-4, 0), (4,0) e interpole en x= 2.5 c) (0,0), (-2, -8), (1, 1), (2, 8) e interpole en x= -1.7 12.Dados los puntos (x k, yk)={(1,1),(2,1),(3,1) y (4,5)}, a)obtenga una estimación de y para x= -1 b)Qué efecto tendría la anterior estimación al omitir el punto (4, 5)?. 13.Dados los puntos (xk, yk)={(0,1),(1,-1),(2,1),(3,-1) y (4,1)}, obtenga una estimación para de y para a) x=0.5 b)x=-1 14.Dados los datos: x
1
2
f(x)
4.75
4.0
3
5
6
5.25 19.75 36
Calcule f(4.5) usando polinomios de interpolación de Lagrange de orden 1 hasta 4. 15.Desarrolle el software correspondiente para la interpolación de Lagrange y aplíquelo para resolver los ejercicios 11, 12 y 13. 16.Se mide la caída de voltaje v através de una resistencia para cierto número de valores de la corriente i. Los resultados obtenidos son: i
0.25
0.75
1.25
v -0.23 -0.33 0.70
1.5
2.0
1.88 6.00
Usando interpolación polinomial, calcule la caída de voltaje cuando la corriente es 0.9.
Métodos Numéricos
Rodrigo A. Ocoró
69
5. INTEGRACIÓN NUMÉRICA 5.1. INTRODUCCIÓN Una vez que se ha determinado un polinomio p(x), bien sea por criterio del ajuste exacto o por el método de mínimos cuadrados, de manera que aproxime satisfactoriamente una función dada f(x) sobre un intervalo de interés, puede esperarse que al diferenciar Pn(x) o integrarlo en forma definida también aproxime satisfactoriamente la derivada o integral definida correspondiente a f(x). Con frecuencia es necesario integrar una función cuando solo es conocida una tabulación de datos. Los métodos de integración numérica se pueden utilizar para integrar funciones dadas ya sea mediante una tabla o en forma analítica ahorrando tiempo y esfuerzo si solo se desea conocer el valor numérico de la integral. 5.2. MÉTODOS DE NEWTON – COTES las fórmulas de integración de Newton – Cotes son los esquemas más comunes dentro de la integración numérica y se basan en la estrategia de reemplazar una función complicada o un conjunto de datos tabulares con alguna función aproximada que sea más fácil de integrar: I = ∫ b∂ f(x)dx ≈ ∫ ∂ b Pn (x)dx
donde Pn(x) es un polinomio de la forma: P(x) = ao + azx + …… + an-1 xn-1 + anxn n = orden del polinomio Entre las fórmulas de Newton – Cotes tenemos: – Regla del trapecio – Reglas de Simpson – Cuadraturas de Gauss 5.3. REGLA DE TRAPECIO La regla trapezoidal corresponde al caso en donde el polinomio es de primer orden (lineal) I = ∫ b∂ f(x)dx ≈ ∫ ∂ b P1 (x)dx fx fa
a
fb
b
x
El área bajo la línea recta es una aproximación de la integral de f(x) entre los límites a y b. f(b) - f(a) (x - a) dx ⇒ I ≈ ∫ b∂ f(a) + b − a
Métodos Numéricos
Rodrigo A. Ocoró
70
que después de hacer diferentes arreglos y ajustes (ver cuadro 2), nos queda como resultado: I ≈ (b - a) ⋅
f (a ) + f (b ) 2
Geométricamente, la regla trapezoidal es equivalente a aproximar el área del trapecio bajo la línea recta que une a f(a) y f(b) en la figura pero el trapecio se encuentra sobre uno de sus lados. Recuerde que en un trapecio
A=h
( b1 + b 2 ) 2
f(a)
Donde: h= b – a (amplitud del intervalo) b1 = f(a) b2 = f(b)
f(b)
A
1
b
a
2
b
h
b
5.3.1. Error de truncamiento en la regla trapezoidal: Cuando se emplea la integral con un segmento de línea recta para aproximar la integral bajo una curva, se comete un error de truncamiento dado por: ET = −
1 12
f ' ' (c)(b − a) 3
donde c es un punto cualquiera dentro del intervalo [a,b].
5.3.2. Ejemplo 1: Utilice la regla trapezoidal para integrar: f(x) = 0.2 + 25x – 200x2 + 675x3 – 900x4 + 400x5 entre 0 y 0.8 a=0 b=0.8
f(0) = 0.2 f(0.8) = 0.232
I = (b - a) E R =
f(a) + f(b)
= (0.8 - 0) ⋅
2 1.640533 − 0.1728 1.640533
0.2 + 0.232 2
= 0.1728
= 89.468%
Para hallar el error de trucamiento, calculemos f ’’(x) F ’(x) = 25 – 400x + 2025x2 – 3600x3 + 2000x4 f ’’(x) = -400 + 4050x – 10800x2 + 8000x3 calculemos el error en 0.8 ⇒ f’’ (0.8) = 24 ET =
1 12
f ' ' (0.8)x(b − a) 3 = −
1 12
24 (0.8) 3 = −1.024
Métodos Numéricos
Rodrigo A. Ocoró
5.3.3.Cuadro 4: Derivación de la regla trapezoidal Antes de integrar, la ecuación anterior se puede expresar como: f ( x ) =
f (b) − f ( a ) b−a
x + f ( a ) −
af (b) − af (a ) b−a
Agrupando los dos últimos términos se obtiene f ( x ) =
+ f ( x ) =
f (b ) − f ( a )
x b−a bf (a ) − af ( a ) − af (b) + af ( a ) b−a f (b) − f ( a ) b−a
x + f (a ) −
bf (b) − af (b) b−a
Que se puede integrar entre x = a y x = b y obtener: f (b) − f (a) x2 bf (a) −af (b) b I ≅ + x a b−a 2 b−a
Este resultado se puede evaluar, obteniendo 2
I ≅
+
2
f (b) − f ( a ) (b − a ) b−a bf ( a ) − af (b) b−a
2 (b − a)
Ahora, considerando que b2 - a2 = (b - a) (b + b) I ≅ [ f (b) − f ( a )]
(b − a) 2
+ bf (a ) − af (b)
Multiplicando y agrupando términos se obtiene: I ≅ (b − a )
f ( a ) + f (b)
2
Que es la fórmula de la regla trapezoidal.
b
71
Métodos Numéricos
Rodrigo A. Ocoró
72
5.3.4. Ejemplo 2: a) Dada la siguiente información: Puntos 0 1 X 500 900 f(x) 9.0 13.4
2 1400 18.7
3 1800 23.0
4 2000 25.1
5 2200 27.2
Aproxime el área bajo la curva en el integral [500, 1800] A 1 = (1800 − 500)
(9 + 23)
= 20800
b = 1800 a = 500 ⇒ f(b) = 23.0 f(a) = 9.0 _____________________________________________ 5 b) Aprox. A 2 = ∫ 0 (2 + 3x)dx f(x) = 2+3x 2
a= 0 ⇒ f(0) = 2+3(0) = 2 b= 5 ⇒ f(5) = 2+3(5) = 17 2 + 17 ⇒ A 2 = (5 − 0) ⋅ = 47.5 ____________________________________________ 2
____
c) Aproxime A 3 = ∫ (1 + 2x + 3x )dx a= -2 f(-2)= 1 + 2(-2) + 3(-2)2 = 9 b= 4 f(4) = 1 + 2(4) + 3(4)2 = 57 4 −2
⇒ A 3 = (4 + 2)
2
( 9 + 57) 2
= 198
________________________________________________ π /2 d) Aproxime A 4 = ∫ 0 senx dx a= 0 ⇒ f(0) = senσ = 0 b =
π 2
⇒ f(π / 2) = sen
⇒ AX
π
=1 2 π 0 + 1 π 1 π = − 0 = ⋅ = 2 2 2 2 4
__________________________________________________ 5.4. REGLA DEL TRAPECIO USANDO SEGMENTOS MÚLTIPLES (Método Tropeziodal Compuesto) Algunas veces el intervalo de integración es tan amplio que resulta conveniente dividirlo en sub-intervalos y aproximar cada uno por medio de un polinomio con el fin de obtener mejor exactitud. Se divide [a,b] en n intervalos y se aplica a cada uno de ellos la regla trapezoidal y la suma de todos ellos de la integral sobre el intervalo completo.
Métodos Numéricos
f(x)
Rodrigo A. Ocoró
f(x2) f(x2)
x0
73
f(x)
f(x4)
f(x2)
x1 x2 x3 x4 x5 dividida en 5 segmentos
x0
x1 x2 dividida en 2 segmentos
Si tenemos n+1 puntos igualmente espaciados (x0, x1, x2 …,xn) Por lo tanto hay n segmentos de igual anchura: h=
b − a n
si a y b se igualan a x0 y a xn respectivamente, entonces la integral total será: x
x
x
I = ∫ x10 f(x) dx + ∫ x12 f(x) dx + ........+∫ x nn -1 f(x) dx
y reemplazando la regla trapezoidal en cada una de las integrales, se tiene: (x 1 ) + f(x 0 ) f(x 2 ) + f(x 1 ) f(x n ) + f(x n −1 ) I = h. f +h +.........+ h 2 2 2 que al agrupar queda: I=
h
n −1
∑ f(x
f(x 0 ) + 2
2
i =1
b - a ) f(x ) donde h = + i n n
El error de la regla trapezoidal múltiple se obtiene sumando los errores individuales de cada uno de los segmentos. ET = −
(b − a) 3 12n 3
n
∑ f' ' (c) 2=1 n
∑ f' ' (c) y como el valor promedio de la 2 a. derivada sobre el intervalo completo es f' ' =
⇒ ET = -
(b - a) 3 12n 2
i=1
n
f ' '
Observe que si el número de segmento se duplica ⇒ el error de truncamiento disminuye a un cuarto de su valor.
Métodos Numéricos
Rodrigo A. Ocoró
5.4.1.Ejemplo 3: Para f(x) = 0.2 + 25x – 200x2 + 675x3 – 900x4 + 400x5 Calcule la integral entre 0 Λ 0.8 tomando: a) 2 segmentos: n= 2 f(0) = 0.2 f(x0) I=
h=
h=0 b=0.8 f(0.4) = 2.456 f(x1)
h
b − a n
=
0.8 2
= 0.4
f(0.8)=0.232 f(x2)
0.4 f(x 0 ) + 2∑ f(x 1 ) + f(x n ) = [ 0.2 + 2(2.456) + 0.23] = 1.0688 2 2 i =1 n −1
b) 3 segmentos: n=3 x0 = a = 0 x3 = b = 0.8
h=
0.8 − 0 3
= 0.2667 ← ancho de cada intervalo
f(x0)=0.2 f(x1)= f(x0 + 0.2667)= f(0.2667) = 1.4329 f(x2)= f(x1 + 0.2667) = f(2h) = f(0.5333) = 3.4874 f(x3) = f(0.8) = 0.232 ⇒ I=
0.2667 2
[ 0.2 + 2(1.4329 + 3.4874) + 0.232] = 1.3698
c) 4 segmentos n=4 a = x0 = 0 b = x 4 = 0.8
h=
0.80 4
f(x0) = f(0) = 0.2 f(x1) = f(0.2) = 1.288 f(x2) = f(0.4) = 2.456 f(x3) = f(0.6) = 3.464 f(x4) = f(0.8) = 0.232 ⇒ I=
0.2 2
= 0.2
0
0.2
0.4
0.6
x0
x1
x2
x3
[ 0.2 + 2(1.28 + 2.456 + 3.464) + 0.232] = 1.4848
0.8 x4
74
Métodos Numéricos
Rodrigo A. Ocoró
75
Teniendo en cuenta que el valor exacto ∫ 0.8 0 f(x) dx = 1.6405
⇒ El E R =
VE − VC VE
Para n = 2
E R =
Para n = 4
E R =
x100%
1.605 − 1.0688 1.6405 1.6405 − 1.4848 1.6405
= 34.8% = 9.5%
d) 10 segmentos: h=
n=10 f(x0) = f(0)= 0.2 f(x1) = f(0.08) = f(x2) = f(0.16) = f(x3) = f(0.32) = f(x10) = f(0.32) =
0.8 − 0 10 0.08
I=
= 0.08
[ 0.2 + 2(........ ) + 0.232]
2 I = 1.615 E R =
1.640 − 1.615 1.640
x100% = 1.6%
Recomendación: use el programa de computación para obtener un cuadro como este y trate de hallar hasta que ER=0.5% n 2 3 4 5 6 7 8 9 10
h 0.4 0.2667 0.2 0.16 0.1333 0.1143 0.1 0.0889 0.08
I
1.068 1.3695 1.4848 1.5399 1.5703 1.5887 1.6008 1.609 1.6150
ER 34.9% 16.5% 9.5% 6.1% 4.3% 3.2% 2.4% 1.9% 1.6%
5.4.2. Ejemplo 4: Dada la siguiente función en forma tabular, aproxime el área bajo la curva entre -1 y 4. Como se tienen 6 puntos ⇒ n=5 intervalos x -1 F(x) 8 x0 = -1
0 10
1 10
2 20
x5 = 4 ⇒
3 4 76 238 h=
4 − ( −1) 5
h =1
Métodos Numéricos
f(x0) = 8 A=
1 2
Rodrigo A. Ocoró
76
f(x5) = 23
[ 8 + 2(10 + 10 + 20 + 76) + 238] = 239
La función analítica para este dato es: f(x)= x 4 - 2x2 + x +10 4 4 2 y el área ∫ −1 x − 2x + x + 10 = 219.17 x0 = -1
h=
x5 = 4 ⇒
⇒ E R =
219.17 − 239 219.17
4 − ( −1)
h =1
5
= 9%
5.5. REGLA DE SIMPSON Además de aplicar la regla trapezoidal con segmento cada vez más finos, otra manera de obtener una estimación más exacta de una integral es la de usar polinomios de orden superior para conectar los puntos. Así por ejemplo, si hay un punto medio entre f(a) y f(b), entonces se pueden conectar los 3 puntos con una parábola (figura a) y si hay 4 puntos igualmente espaciados esta f(a) y f(b) entonces los 4 puntos se pueden conectar a un polinomio de 3er orden (figura b). f(x)
f(x)
Figura a
x
Figura b
x
A las fórmulas de calcular la integral bajo estos polinomios, se les llama Regla de Simpson. 5.5.1. Regla de Simpson de 1/3: Resulta cuando se sustituye un polinomio de 2º grado en la ecuación: I = ∫ xx 20 f(x) dx = ∫ ba f 2 (x) dx si a y b se denominan como x 0 y x2 y f2(x) se representa mediante un polinomio de Lagranje de 2º orden. (x − x 1 )(x − x 2 ) (x − x 0 )(x − x 2 ) (x − x 0 )(x − x 1 ) ⇒ I = ∫ xx f(x 0 ) + f(x 1 ) + f(x 2 ) dx (x 2 − x 0 )(x 2 − x 1 ) (x 1 − x 0 )(x 1 − x 2 ) (x 0 − x 1 )(x 0 − x 2 ) 2 0
Después de integrar y de reordenar términos, resulta la siguiente ecuación: I=
h 3
[ f(x 0 ) +
d on de h =
4 f(x 1 ) + f(x 2 ) ]
b - a 2
(La etiqueta de “1/3”, viene de que se le divide por 3, mientras que a la trapezoidal se divide entre 2).
Métodos Numéricos
Rodrigo A. Ocoró
77
5.5.1.1. Ejemplo 5: Utilice Simpson de 1/3 simple para hallar: ∫ 0.8 donde f(x) = 0.2 + 25x - 200x 2 + 675x 3 − 900 X 4 + 400x 5 0 f(x) se toman 3 puntos (x 1 esta en la mita) x0 = a = 0
x 1 = 0.4
x 2 = b = 0.8 f(x 0 ) = 0.2 f(0.4) = 2.456 f(0.8) = 0.232 h=
0.8 − 0 2
⇒ I=
0.4 3
= 0.4
[ 0.2 + 4(2.456) + 0.232] = 1.367467
Comparando con VE = 1.6405333 ⇒ ER =
1.6405333 − 1.36747 1.6405333
= 16.6%
5.5.2. Regla de Simpson de 1/3 de segmento multiple (Método de Simpson compuesto): Al igual que la regla tropezioal, la de Simpson se puede mejorar dividiendo el intervalo de integración en segmento de igual amplitud. h=
b − a n
f(x)
a
b
x
La integral se representa como: I=
∫
x2
x0
f(x)dx +
∫
x4
x2
f(x)dx + ... +
∫
xn
x n −2
f(x)dx
al sustituir la regla de Simpson en cada uno e los integrales individuales se obtiene: I
f(x ) + 4f(x ) + f(x ) f(x ) + 4f(x ) + f(x ) f(x ) + 4f(x ) + f(xn ) − − 0 1 2 2 3 4 n 2 n 1 =h +h + ... + h
3
3
Reordenando términos queda: n −1 n −2 h + + + f(x ) 4 f(x ) 2 f(x ) f(x ) 0 ∑ i j=∑ j n 3 i =1,3,5 2,4,6 I= Nota: Se debe usar un número par de segmentos:
3
Métodos Numéricos
Rodrigo A. Ocoró
78
5.5.2.1..Ejemplo 6: Para f(x) = 0.2 + 25x – 200x2 + 675x3 – 900x4 + 400x5
∫ Utilice la regla de Simpson de 1/3 con n = 4 para calcular
0,8
0
0.8 − 0
x0 = a = 0
f(0) = 0.2
h=
4
f(x1) = f(0.2) = 1.288 f(x2) = f(0.4) = 2.456 f(x3) = f(0.6) = 3.464 f(x4) = f(0.8) = 0.232
0.2
I=
n=4
0
0.2
0.4
x0
x1
x2
0.6
x3
0.8
x4
n−1 n−2 h f(x0 ) + 4 f(xi ) + 2 f(xj ) + f(xn ) 3 i =1,3,5 j=2,4,6
∑
I=
= 0.2
f(x)dx
3
∑
[ 0.2 + 4(1.288 + 3.464) + 2(2.456) + 0.232] = 1.623466 − 162234662 1640533 . . 16405333 .
Comparando con VE = 1.6405333 ⇒ ER =
x100% => ER=1.04%
5.5.2.2. Ejemplo 7:
1 2 π
Encuentre la integral aproximada de la función (expresión de la curva normal tipificada).
e
− x 2 /2
a) Utilizando la regla trapezoidal con n = 4 b) Utilizando la regla de Simpson con n = 4 En ambos casos compare con el resultado de tablas que es 0,682 Solución: b=1 a = -1 h
a) I =
f(x 0 ) + 2
2
n −1
∑
b − a
f(x i ) + f(x n )
i =1
1 f(x0) = f(-1) = 2π
(−1)2 − 2
e
=
h=
1 2π
e−05 =
n
=
0.606 2π
1 − (−1) 4
= 0,5
entre –1 y 1
Métodos Numéricos
Rodrigo A. Ocoró
1 2π
f(x1) = f(-0.5) = 1 2π
f(x2) = f(0) =
e
1 2π
f(x3) = f(0.5) = 1
f(x4) = f(1) =
e
2π
−
(0.5)2
02
e −
e
−
2
−
1
= (0.5)
2π 2
2
12 2
2
=
= 1 2π
1
=
2π
2π
2π
2π
e − 0125 =
e − 05 =
0.882
1
e0 = 1
e −0.125 =
79
0.882 2π
0.606 2π
0.5 0.60
0.882 1 0.882 0.606 + 2 + + + = 0.672 2π 2 2π 2π 2π 2π ⇒I= 0.682 − 0672
ER =
0.682
n−1 n−2 h f(xi ) + 2 f(x j ) + f(xn ) f(x0 ) + 4 3 i =1,3,5 j=2,4,6
∑
b)
I=
0.5
I=
= 0.0147 = 1.47%
3 2π
∑
[ 0.606 + 4(0882 + 0882) + 2(1) + 0.606] = 0.683 0.682 − 0.683 0.682
ER =
= 0.0015 = 0.15%
5.5.3. Regla de Simpson de 3/8. Es de utilidad en aplicaciones de segmentos múltiples cuando el número de segmentos es impar. De manera similar a la derivación de la regla trapezoidal y a la regla de simpson de 1/3, se pueden ajustar polinomios de lagrange de tercer orden a 4 puntos e integrar: I = ∫ b∂ f(x)dx ≈ ∫ ∂ b P1 (x)dx para obtener:
3h f(x0 ) + 3 8
∑ f(x1 ) + 3∑ f(x2 ) + f(x3 ) donde h=(b-a)/3
I= La ecuación se llama regla de simpson de 3/8 por que h es múltiplo de 3/8 2 3 5 4 5.5.3.1.Ejemplo 8. para f(x) = 0.2 + 25x - 200x + 675x − 900 X + 400x , ∫ 0.8 0 f(x) Requiere 4 puntos con n=3 |____|___|____|
hallar
Métodos Numéricos
Rodrigo A. Ocoró
x0
x0 = a = 0 x3 = b = 0.8 h=
0.8 − 0 3
x1 x2
80
x3
= 0.2667 ← ancho de cada intervalo
f(x0)=f(0)=0.2 f(x1)=f(0.2667) = 1.43272428 f(x2)= f(0.5333) = 3.48717696 f(x3) = f(0.8) = 0.232 I =
3x 0.2667 8
[0.2 + 3(1.43272428 ) + 3(3.4817696) + 0.232] 1.6405333 − 1.519170
I=1.519170
1.6405333
Er=
x100% = 7.4%
5.6. CUADRATURA DE GAUSS Con los métodos trapezoidal y de Simpson la estimación de la integral se basa en puntos igualmente espaciados por lo tanto la posición de los puntos base usados en estas ecuaciones estaba predeterminado o fijo, generando con ello un error muy grande. C
A
D B
x0
x1
b Fig. a: Método trapezoidal
a
b
Fig. b: Método de Gauss
Gauus investigó y encontró que es factible disminuir el error en la integración cambiando la localización de los puntos sobre la curva de integración (figura b) obteniéndose un valor más exacto de la integral. El método de Gauss consiste en seleccionar los puntos C y D adecuados. 5.6.1. Obtención del método de integración Gaussiana usando dos puntos: f(x2) D
f(x)
f(x1) C -1
x1
x2
1
Supongamos que deseamos integrar la función f(x) entre –1 y 1 Gauss se propuso desarrollar una fórmula del tipo
Métodos Numéricos
I=
∫
1
−1
Rodrigo A. Ocoró
81
f(x) = C1 f(x 1 ) + C 2 f(x 2 )
Para simplificar el cálculo del área en donde C 1 y C2 son los coeficientes de las incógnitas y x1 y x2 no están fijos a los puntos extremos sino que son incógnitas por lo tanto se tienen 4 incógnitas que se deben evaluar y por consiguiente se requieren 4 condiciones que se pueden establecer así. f(x) = 1 f(x) = x f(x) = x2 f(x) = x3
función constante función lineal función cuadrática función cúbica
Los valores al integrar estas 4 funciones entre -1 y 1 son: I1 =
∫
1
−1 1
I2 =
∫ xdx =
I4 =
x3
2
∫
x dx =
=0
1
2
= 3 −1 3
−1
3
1
2 −1
∫ x dx = 1
=2
−1
x2
−1
1
I3 =
1
]
1dx = x
x
−1
4
1
] −1 = 0
4
⇒ tenemos al reemplazar en I = C1f (x1) + C2f(x2): I1 = C1 (1) + C2 (1) = 2 ⇒ C1 + C2 = 2 (1)
I2 = C1 x1 + C1 x2 = 0 I3 = C1 x12 + C2 x22 = 2/3 I4 = C1 x13 + C2 x23 = 0
(2) (3) (4)
Si en la primera ecuación C1 = C2 = 1 Entonces en la segunda ecuación se tendría que x 1 + x2 = 0 ⇒ x1 = –x2 que también satisface para la cuarta ecuación, al sustituir en la tercera ecuación, se obtiene: x12 + (-x1)2 = 2/3 2x12 = 2/3 ⇒ x12 = 1/3 2
Luego x 1 = − x2 =
3 1 3
Luego x1 =
= −0.57735....
= 0.57735....
Por lo tanto, la fórmula quedará:
±
1 3
=±
1 3
Métodos Numéricos
I=
1
∫
−1
Rodrigo A. Ocoró
82
f(x) dx = c 3 f(x1 ) + c 2 f(x 2 )
-1 1 ⇒ I = f + f 3 3 Si los límites de integración son de a a b y no –1 a 1 para los cuales de dedujo la fórmula entonces se debe efectuar un cambio de variables. Se reemplaza la variable dada por otra en la que esté relacionada linealmente. x* =
Si se tiene entonces
b
∫
a
b
∫
a
f(x)dx =
(b − a)x + (b + a) 2
f(x)dx =
b − a 2
∫
1
−1
b − a de manera que dx * = dx 2
(b − a)x + (b + a) dx 2
f
t= ±0.57735
b − a b − a a + b b − a a + b f ( 0.57735) f (0.57735) − + + + 2 2 2 2 2
Que es la fórmula de Gauss Legendre para dos puntos. Si se tienen más de dos puntos. I =
b−a
2
[c1 f ( x1) + c2 f ( x2) +.......... ......+cn f ( xn)
donde los valores de c (factores de peso) y los x* (argumentos de la función) se establecen en el cuadro 3 según el número de puntos considerados.
Donde:
X i =
(b − a ) 2
x * +
(b + a ) 2
Métodos Numéricos
Rodrigo A. Ocoró
83
5.6.2.Cuadro 5. Factores de peso c y argumentos x de la función usados en las fórmulas de Gauss-Legendre Puntos 2
c1 c2
Factores de peso = 1.000 000 000 = 1.000 000 000
Argumentos de la función (X*) = - 0.577 350 269 = 0.577 350 269
x1 x2
3
c1 c2 c3
= = =
0.555 555 556 0.888 888 889 0.555 555 556
X1 x2 x3
= = =
- 0.774 596 669 0.0 0.774 596 669
4
c1 c2 c3 c4
= = = =
0.347 854 845 0.652 145 155 0.652 145 155 0.347 854 845
x1 x2 x3 x4
= = = =
- 0.861 136 312 - 0.339 981 044 0.339 981 044 0.861.136.312
5
c1 c2 c3 c4 c5
= = = = =
0.236 926 885 0.478 628 670 0.568 888 889 0.478 628 670 0.236 926 885
x1 x2 x3 x4 x5
= = = = =
- 0.906 179 846 - 0.538 469 310 0.0 0.538 469 310 0.906 179 846
6
c1 c2 c3 c4 c5 c6
= = = = = =
0.171 324 492 0.360 761 573 0.467 913 935 0.467 913 935 0.360 761 573 0.171 324 492
x1 x2 x3 x4 x5 x6
= = = = = =
- 0.932 469 514 - 0.661 209 386 - 0.238 619 186 0.238 619 186 0.661 209 386 0.932 469 514
Métodos Numéricos
Rodrigo A. Ocoró
84
Ejemplo 1: Calcular ∫
5
0
e − x dx con 2 puntos b − a
b=5
a=0
⇒
2
=
: c1=1 x1*= - 0.57735
5
b + a
2
2
5
=
2
5 * 5 b − a * b + a x= x + ⇒ x= x + 2 2 2 2 5
f(x 1 ) = e x1 = e − 2
f(x 2 ) = e
⇒ I=
x 2
b - a
I=
2 5 2
( − 0.57735) +
5 2
5 5 − 2 (0.57735) + 2
=e
0.917525
f(x 1 ) = e
f(x 2 ) = e
x 2
f(x 3 ) = e
x 3
I =
I =
b − a 2 5 2
= e
= e = e
5 5 ( − 0.77460) + 2 2
−
5 5 x0.0 + 2 2
−
1
2
= e - 3.943375 = 0.019383
Si consideramos 3 puntos: ⇒
x 1
x * = −0.57735 x = * x = 0.57735 *
= e -1.056625 = 0.347627
[c 1 f(x 1 ) + c 2 f(x 2 ) ] =
[0.36701 ] =
c2=1 x2*= 0.57735
5 2
[1 x 0.347627
⇒
∫
5
0
+ 1 x 0.019383
]
e − x dx = 0.917525
x1*= -0.77460 x2*= 0.0 x3*= 0.77460
c1=0.55556 c2=0.88889 c3=0.55556
= e − 0.5635 = 0.569213
2.5 = e− = 0.082085
5 5 x0.77460 + 2 2
−
= e − 4.4365 = 0.011837
[c 1 f(x 1 ) + c 2 f(x 2 ) + c 3 f(x 3 ) ]
[0.55556
x 0.569213
El valor exacto de
∫
5
0
+ 0.88889
x 0.08208 + 0.55556
e − x dx = 0.993262
⇒
E R =
x 0.01837
0.993262
]=
0.998505
− 0.998505
0.99362
= 0.5%
Métodos Numéricos
Rodrigo A. Ocoró
85
ALGORITMO Método trapezoidal Para aproximar el área bajo la curva de una función analítica f(x) en el intervalo [a, b], proporcionar la función por integrar F(x) y los DATOS:
El número de trapecios N, el límite inferior A y límite superior B. RESULTADOS: El área aproximada AREA. PASO 1. Hacer X = A PASO 2. Hacer S = 0 PASO 3. Hacer H = (B - A)/N PASO 4. Si N = 1, ir al paso 10. De otro modo continuar. PASO 5. Hacer I = 1 PASO 6. Mientras I ≤ N - 1, repetir los pasos 7 a 9. PASO 7. Hacer X = X + H PASO 8. Hacer S = S + F(X) PASO 9. Hacer I = I + 1 PASO 10. Hacer AREA = H/2 * (F(A) + 2*S+F(B)) PASO 11. IMPRIMIR AREA y TERMINAR.
ALGORITMO Método de Simpson Para aproximar el área bajo la curva de una función analítica f(x) en el intervalo [a,b], proporcionar la función por integrar F(X) y los DATOS:
El número (par) de subintervalos N, el límite inferior A y el límite superior B. RESULTADOS: El área aproximada AREA. PASO 1. Hacer S1 = 0 PASO 2. Hacer S2 = 0 PASO 3. Hacer X = A PASO 4. Hacer H = (B-A)/N PASO 5. Si N = 2, ir al paso 13. De otro modo continuar. PASO 6. Hacer I = 1 PASO 7. Mientras I ≤ N/2-1, repetir los pasos 8 a 12. PASO 8. Hacer X = X + H PASO 9. Hacer S1 = S1 + F(X) PASO 10. Hacer X = X + H PASO 11. Hacer S2 = S2 + F(X) PASO 12. Hacer i = I + 1 PASO 13. Hacer X = X + H PASO 14. Hacer S1 = S1 + F(X) PASO 15. Hacer ÁREA = H/3 * (F(A) + 4*S1 + 2*S2 + F(B)) PASO 16. IMPRIMIR AREA y TERMINAR.
Métodos Numéricos
Rodrigo A. Ocoró
86
ACTIVIDADES Nº 5 1.Evalúe cada una de las integrales trapezoidal y regla de simpson : a) X 0 0.1 F(x) 1 7 b) X -3 -1 F(x) 1 -4
de los siguientes datos tabulares usando: regla 0.2 0.3 0.4 0.5 0.6 4 3 5 8 9 1 3 -5 2
5 4
7 8
9 6
2.Dadas las siguientes integrales: 10 (10 + 2 x − 6 x 2 + 5 x 4 )dx 0
5 (1 − x − 4 x 3 + 3 x 5 )dx −3
π
∫
a) ∫ b) c) ∫0 Evalúe cada una de ellas, utilizando: 1.1.La regla trapezoidal simple 1.2.La regla trapezoidal de segmentos múltiples con n=2, 4 y 6. 1.3.La regla simple de simpson de 1/3 1.4.La regla de simpson de 1/3 de segmentos múltiples con n=4 y 6. 1.5.La regla de simpson de 3/8. 1.6. Halle el Error relativo en cada caso
(8 + 5 sen x ) dx
3.Integre analíticamente la siguiente función y luego resuelva usando la regla de simpson de 1/3 con n=4 y 5 y determine el error en cada caso con relación al valor analítico. 6 [(4 x + 8)3]dx −4
∫
4 2 x xe dx 0
4.Dada : ∫ realice mediante: a)La regla trapezoidal con n=4 b)La regla de simpson con n=4 c)determine el valor analítico y establezca el error para a) y b). 5.Usando las fórmulas de Gauss-Legendre dos, tres y cuatro puntos, obtenga el valor de: 4 2 x xe dx 0 y calcule el Er en c/u.
∫
6.Mediante la cuadratura Gausiana de tres términos, evalúe: respuesta: 0.94608
1 sen x dx x
∫0
Métodos Numéricos
Rodrigo A. Ocoró
87
6. ECUACIONES DIFERENCIALES NUMÉRICAS
6.1 INTRODUCCIÓN Con mucha frecuencia aparecen problemas en ingeniería, física, química, ecología meteorología, sociología, etc., que exigen el manejo de ecuaciones diferenciales, muchas de las cuales no se pueden resolver pos los métodos convencionales, teniéndose que recurrir a métodos numéricos aproximados. Se llama ecuación diferencial aquella ecuación que contiene una variable dependiente y sus derivadas con respecto a una o más variables independientes. Se dividen en dos grandes grupos: Ordinarias, si contienen una sola variable independiente y Parciales, cuando contienen varias variables independientes. Estudiaremos únicamente las ecuaciones diferenciales ordinarias (EDO). Estas se clasifican y estudian según el orden de la mayor derivada que aparece en la respectiva ecuación diferencial. dy
Ejemplo:
dt
m
= − ky
d2y dt 2
d3y dx 3
es de 1er orden
= ky
−5
dy dx
+ 6y = 0
es de 2º orden es de 3er orden
Los problemas que encierran el uso de ecuaciones diferenciales ordinarias, constarán de: dy * *
Una ecuación diferencial ordinaria: Las condiciones iniciales.
dx
= y 1 = f(x, y)
El intervalo para la variable independiente x, para el cual debe determinarse la función y, si existe. *
Métodos Numéricos
Rodrigo A. Ocoró
88
6.2 MÉTODO DE EULER h=
xn − xo
n Consiste en dividir el intervalo de X 0 a Xn en n subintervalos de ancho h. obteniéndose un conjunto discreto de (n + 1) punto: x 0, x1, x2 …… xn en el intervalo [x0 , xn] generándose la sucesión de aproximaciones siguientes:
y1 = y0 + hf (x0 , y0) y2 = y1 + hf (x1 , y1) y3 = y2 + hf (x2 , y2) . . yn+1 = yn + hf(xn , yn) ⇒
yi+1 = yi + hf(xi, yi)
con i = 0 hasta n
siendo f(xi , yi) la ecuación diferencial evaluada en x i y yi. Ejemplo 1: Utilice el método de Euler para integrar numéricamente la ecuación: y’=f(x,y) = -2x3 + 12x2 – 20x + 8.5 de x=1.0 hasta x=3.0 considerando 4 intervalos y el valor inicial en x=1 es y=3 X0 = 1
h=
f(x0 , y0) = f(1 , 3)
3 −1 4
= 0.5
1.0
1.5
2.0
x0
x1
x2
2.5 x3
3.0 x4
yi+1 = yi + hf(xi , yi) i 0 1 2 3 4
xi 1 1.5 2.0 2.5 3.0
yi 3 2.25 1.625 1.875 3.000
F(xi , yi) -1.5 -1.25 0.50 2.25 2.5
∫ f ( x, y)dx = ∫ (−2 x
La solución exacta está dada por => y= -0.5x4+4x3 –10x2+8.5x+1 y(3)=4
⇒ E R =
4 − 4.25 4
x100% = 6.25%
yi+1=yi+hf(xi ,yi) 2.25 1.625 1.875 3.000 4.25
3 + 12 x 2 − 20 x + 8.5)dx
si se aumenta n entonces h es menor y da mas exacto
Métodos Numéricos
Rodrigo A. Ocoró
dy
89
= f(x, y)
Ejemplo 2. Dada la ecuación diferencial dx = x – y , resuelva por el método de Euler para el intervalo [0, 1] considerando 5 intervalos y siendo y(0)=2
Y(1)=?
h=
1− 0 5
0
= 0.2
x0
0.2
0.4
x1
x2
0.6
0.8
x3
x4
yi+1 = yi + hf(xi , yi) i 0 1 2 3 4 5
xi 0 0.2 0.4 0.6 0.8 1.0
yi 2 1.6 1.32 1.136 1.0288 0.98304
El valor exacto es 1.10364
E R =
F(xi , yi) -2 -1.4 -0.92 -0.536 -0.2288 -0.01696
1.10364 − 0.986432 1.10364
1.0 x5 yi+1=yi+hf(xi ,yi) 1.6 1.32 1.136 1.0288 0.98304 0.986432
x100% = 10.93%
3. ANÁLISIS DEL ERROR EN EL MÉTODO DE EULER MEDIANTE EL MÉTODO DE TAYLOR La solución numérica de una ecuación diferencial ordinaria incluye dos tipos de errores: 1. Errores de truncamiento causados por la naturaleza de los métodos empleados en la aproximación a los valores de y. 2. Errores de redondeo causados por el número limitado de dígitos o de cifras significativas que pueda tener la computadora. El conocimiento de la magnitud y propiedades del error de truncamiento se puede obtener derivando el método de Euler directamente de la expansión de la serie de Taylor. Recuerde que la ecuación diferencial que se está integrando tiene la forma general dy
y’=f(x,y) donde y’= dx , respectivamente.
x
y son las variables independientes y dependientes
y
La serie de Taylor se expresa como: f(xi+1 ) = f(xi ) + f ' (xi )h + R n =
f (n + 1) (c) (n + 1)!
f ' ' (xi ) 2!
2
h +
f ' ' ' (xi ) 3!
3
h + ......+
f (n)(xi )hn
h n +1
en donde Rn es el término residual
n!
+ R n
Métodos Numéricos
Rodrigo A. Ocoró
90
Si para la ecuación diferencial y’’ =f ‘(x,y) el punto inicial es (xi ,yi) siendo y’ =f(x,y) yi+1 = yi + yi ' h +
yi ' ' h2 2!
(n)
+......+
, yi hn n!
+ R n
R n =
con
f
(n + 1 )
(c )
(n + 1 )!
h
n +1
⇒ como y’=f(x,y), entonces podemos reacomodar así:
yi+1 = yi + f(xi , yi )h +
f ' (xi , yi )h2 2!
+......+
f (n−1) (x , yi )hn n!
+ R n
Este es el método de Taylor que al comparar con el método de Euler y i+1=yi+hf(xi,yi), se observa que este equivale a la serie de Taylor truncada hasta el término f(x i, yi)h que corresponde al método de Taylor de 1 er orden dejando por lo tanto una parte de la solución verdadera. Luego los términos restantes que quedan constituyen el error de truncamiento E T. ⇒ ET =
f ' (x i , y i )h 2 2!
+
f ' ' (xi , y i )h3 3!
+ ......+
f (n−1) (xi , y i )h n n!
Ejemplo: para la ecuación y’=f(x,y)= -2x3+12x2 –20x + 8.5 utilice el método de Taylor para evaluar en el primer paso (n=1), x=1.0 siendo h=0.5 valor inicial en x=1 y=3 y compare con el resultado obtenido para el mismo paso con el método de Euler. Solución:
1.0
1.5
f(xi, yi) = -2x3+12x2 –20x + 8.5 f’(xi, yi) = -6x2 +24x –20 f’’(xi, yi) = -12x + 24 f’’’(xi, yi) = - 12 de aquí en adelante las f (n)(xi, yi)=0
x0
x1
Por Taylor:
yi +1 = yi + f ( xi , yi )h +
f ' ( xi , yi ) h 2
2!
+
f " ( xi , yi )h 3
En (x0,y0) ⇒ f(1,3) = -2(1)3 + 12(1)2 – 20(1) + 8.5 = -1.5 f’(1,3) = -6(1)2 + 24(1) – 20 = -2 f’’(1,3) = -12(1) + 24 = 12 f’’’(1,3) = -12
3!
+
f ' ' ' ( xi , yi )h 4
4!
Métodos Numéricos
Rodrigo A. Ocoró
⇒ y 1 = y 0 + f(x 0 , y 0 )h +
y 1 = 3 + ( −1.5)x0.5 +
f ' (x 0 , y 0 )h 2 2!
− 2(0.5) 2
2
+
+
f ' ' (x 0 , y 0 )h 3 3!
12(0.5) 3 6
+
+
91
f ' ' ' (x 0 , y 0 )h 4 4!
− 12(0.5) 4
24
y 1 = 2.25 − 0.25 + 0.25 − 0.03125 = 2.21875
El valor exacto es x=1.5 para ∫f(x,y) dx ⇒ y= -0.5x4 + 4x3 – 10x2 + 8.5X + 1 es y= 2.21875 ⇒ E R =
2.21875 − 2.21875 = 0% 2.21875
Con el método de Euler: Y1=y0 + f(x0,y0)h =3+(-1.5)x0.5=2.25
Ea=2.21875 –2.25 = -0.03125
2.21875 − 2.25
ER=
2.21875
x100%=1.408%
f ' ( x , y )h 2 f " ( x , y )h 3 f ' ' ' ( x , y )h 4 − 2(0.5) 2 12(0.5)3 − 12(0.5) 4 0 0 0 0 0 0 + + = + + 2 ! 3 ! 4 ! 2 6 24 ET=
ET= -0.25+0.25-0.03125 = -0.03125 que corresponde a Ea ⇒ Por Euler: yi+1 = yi + f(xi,yi)h + ET
La serie de Taylor es un medio para cuantificar el error en el método de Euler. Sin embargo, presenta algunos inconvenientes como: 1) Solamente proporciona una aproximación local del error de truncamiento durante el primer paso. 2) La extensión a 4, 5 o más pasos de la serie de Taylor significa obtener métodos con mayor exactitud pero menos prácticos ya que incluirían diferenciaciones complicadas para f(x,y). 3) Las derivadas necesarias para evaluar la serie de Taylor no siempre son fáciles de obtener.
Métodos Numéricos
Rodrigo A. Ocoró
92
6.4 MÉTODO DE EULER MODIFICADO Una fuente fundamental de error en el método de Euler es que la derivada al principio del intervalo se supone que se aplica a través del intervalo entero y para obtener una exactitud razonable se utiliza un intervalo muy pequeño a cambio de un error de redondeo mayor. El método de Euler modificado trata de evitar este problema utilizando un valor promedio de la derivada tomada en los dos extremos del intervalo en lugar de la derivada tomada en un solo extremo. El método de Euler modificado consta de dos pasos básicos: 1. Se parte de (x0,y0) y se utiliza el método de Euler a fin de calcular el valor de y correspondiente a x1. Este valor de y se denotará como y 1 ya que solo es un valor transitorio para y1. Esta parte del proceso se conoce como paso predictor. 2. El segundo paso se llama corrector, pues trata de corregir la predicción. En el nuevo punto obtenido (x 1 , y 1 ) se evalúa la derivada de f (x 1 , y1 ) usando la ecuación ordinaria que se está resolviendo. Se obtiene la media aritmética de esta derivada y la derivada en el punto inicial (x0,y0). 1 2
[f(x
0
,y 0 ) + f(x 1 ,y 1 )] = derivada promedio
se usa la derivada promedio para calcular un nuevo valor de y i con la ecuación de Euler yn+1 = yi + hf(xi, yi) que deberá ser más exacto que y 1 . ⇒ y1 = y 0 +
h 2
[f(x
0
,y 0 ) + (f(x 1 ,y 1 )) ]
que se tomará como valor definitivo de y1.
xn − x0 h = n
Este procedimiento se repite hasta llegar a yn. El esquema iterativo para este método quedaría en general así: 1. Usando el paso de predicción resulta: y i + 1 = y i + hf(x i , y i ) 2. Se calcula la derivada f(x i + 1 , y i + 1 ) 3. Se establece la derivada promedio (llamémosla B) B=
1 2
[ f(x
i
, y i ) + f(x i + 1 , y i + 1 ) ]
4. Se sustituye f(xi,yi) con este valor promedio en la ecuación de iteración de Euler y h y i +1 = y i + [ f(x i , y i ) + f(x i +1 ,y i +1 )] 2 se obtiene o simplificado: yi+1= yi+ hB
Métodos Numéricos
Rodrigo A. Ocoró
93
Ejemplo 1: Resuelva Resuelva f(x,y) = x - y en el intervalo [0,1] Con y(0)=2 y(1)=? y Con 5 intervalos. x0=0 y0=2 xn=1 yn=? 0 x0
n=5
0.2
0.4
x1
x2
h =
0.6
0.8
x3
xn − x0 n
=
1− 0 5
h = 0 .2
1
x4
x5
En forma iterativa, hagamos la siguiente tabla: de i=0 hasta n=5 Recuerde que:
xi+1=xi+h
f(xi,yi) = xi - yi
yi+1 = yi + hf(xi , yi )
i
xi
Yi
0 1 2 3 4 5
0 0.2 0.4 0.6 0.8 1.0
2 1.66 1.4172 1.254104 1.156365 1.112222
Para i=5
f(xi,yi)
B=
1 2
[f(x
y i +1
-2 -1.46 -1.0172 -0.654104 -0.356365
1.6 1.368 1.21376 1.12328 1.085092
i
, y i ) + f(x i +1 , y i +1 )
Xi+1 0.2 0.4 0.6 0.8 1.0
y5=y(1)=1.112222 El valor exacto es 1.10364 ⇒ E R =
1.10364 − 1.112222 1.10364
E R = 0.78%
x100%
]
f(x i +1 , y i +1 )
-1.4 -0.968 -0.61378 -0.32328 -0.08509
B
Yi+1=yi+hB
-1.7 -1.214 -0.81548 -0.488694 -0.2207286
1.66 1.4172 1.254104 1.156365 1.112220
Métodos Numéricos
Rodrigo A. Ocoró
94
6.5 .METODOS .METODOS DE RUNGE-KUTTA (RK) Los métodos de Runge-Kutta tienen la exactitud del esquema de la serie de Taylor sin necesidad de cálculo de derivadas superiores y consisten en obtener un resultado que se obtendría al utilizar un número finito de términos térmi nos en una serie de Taylor de la forma: y i +1 = y i + f ( x i , y i ) h +
f ' ( x i , y i ) h
2
2!
+
f " ( x i , y i ) h
3!
3
+
f ' ' ' ( x i , y i ) h 4
+..........
4!
(1)
con una aproximación en la cual cual se calcula yi+1 de la fórmula: Y i +1 = yi + h(α 0 f ( xi , yi ) + α 1 f ( xi + µ 1h, yi + b1h) + ....... + α p f ( xi + µ p h, yi + b p h ) )
(2)
f ( x + µ h, y + b h)
i j i j donde los α , µ , b se determinan de modo que si expandiera con i ≤ j ≤ p en series de Taylor alrededor de (xi, yi), se observaría que los coeficientes de h, h2 , h3 , etc. Coincidirían con los coeficientes correspondientes de la ecuación (1).
Después de efectuar los cálculos convenientes a la ecuación (2) para establecer los valores de α , µ , b , tenemos las siguientes expresiones de Runge-Kutta (RK): 1.Método RK de segundo orden: yi+1 = yi +
h
[k1 +k2]
2
donde: k1=f(xi, yi) k2=f(xi+h, yi+hk1)
2.Método RK de tercer orden: yi+1 = yi +
h
6
[k1 +4k2 +k3] donde: k1=f(xi, yi) h
h
2
2
k2=f(xi+ , yi+ k1) k3= f(xi+h, yi -hk1+2hk2) 3. Método RK de cuarto orden: orden: yi+1 = yi +
h
6
[k1 +2k2 +2k3 +k4] donde: k1=f(xi, yi) h
h
2 h
2 h
2
2
k2=f(xi+ , yi+ k1) k3=f(xi+ , yi+ k2) k4=f(xi+h, yi+hk3) *El método más utilizado utilizado es el de cuarto cuarto orden ya que coincide coincide con los primeros cinco términos de la serie de Taylor lo cual significa gran exactitud sin cálculo de derivadas aunque haya que evaluar la función f(x,y) f (x,y) cuatro veces en cada subintervalo.
Métodos Numéricos
Rodrigo A. Ocoró
95
Ejemplo: Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)=x-y donde x=0 hasta hasta x=1 con 5 intervalos y siendo para x=0, y=2. 6.5.1.Ejemplo: Aplíquese el método de Runge-Kutta de cuarto orden para f(x,y)= x – y desde x=0 hasta x=1 Con 5 intervalos siendo para x=0 y=2. Solución:
1-0 (x0, y0)= (0, 2) h= ------ =0.2 5 para x=1 y=? x0=0
x1=0.2
Aplicamos: yi+1 = yi +
x2=0.4 h
6
x3=0.6
x4=0.8
x5=1.0
[k1 +2k2 +2k3 +k4]
Se calculan los valores de k1, k2, k3, k4 en cada iteración y se halla y i+1 i 0 1 2 3 4 5 para
xi 0 0.2 0.4 0.6 0.8 1.0
Yi 2 1.6562 1.410973 1.246451 1.148004 1.103656
K1 -2 -1.4562 -1.010973 -0.746451 -0.348004
K2 -1.7 -1,21058 -0.809876 -0.471806 -0.213204
K3 -1.73 -1.235142 -0.829985 -0.49927 -0.226684
x=1.0, y=1.103656
que al comparar con el valor exacto que es 1.10364, tenemos: E R =
1.10364 − 1.103656 1.10364
x100% = 0.001%
K4 -1.454 -1.009172 -0.644976 -0.346797 -0.10268
Yi+1 1.6562 1.410973 1.246451 1.148004 1.103656
Métodos Numéricos
Rodrigo A. Ocoró
96
ALGORITMO Método de Euler modificado Para obtener la aproximación YF a la solución de un PVI, proporcionar la función F(X,Y) y los DATOS: RESULTADOS: PASO 1. PASO 2. PASO 3.
PASO 8.
La condición inicial X0, Y0, el valor XF donde se desea conocer el valor de YF, y el número N de subintervalos por emplear. Aproximación a YF: Y0.
Hacer H = (XF - X0)/N Hacer I = 1 Mientras I ≤ N, repetir los pasos 4 a 7. PASO 4. Hacer Y1 = Y0 + H * F(X0, Y0) PASO 5. Hacer Y0 = Y0 + H2 * (F(X0, Y0)+F(X0+H,Y1)) PASO 6. Hacer X0 = X0 + H PASO 7. Hacer I = I + 1 IMPRIMIR Y0 y TERMINAR. ALGORITMO Método de Runge - Kutta de cuarto orden
Para obtener la aproximación YF a la solución de un PVI, proporcionar la función F(X,Y) y los DATOS: RESULTADOS: PASO 1. PASO 2. PASO 3.
PASO 11.
La condición inicial X0, Y0, el valor XF donde se desea conocer el valor de YF y el número N de subintervalos a emplear. Aproximación a YF: Y0
Hacer H = (XF - X0)/N Hacer I = 1 Mientras I ≤ N, repetir los pasos 4 a 10. PASO 4. Hacer K1 = F(X0, Y0) PASO 5. Hacer K2 = F(X0 + H/2, Y0 + H * K1/2) PASO 6. Hacer K3 = F(X0 + H/2, Y0 + H * K2/2) PASO 7. Hacer K4 = F(X0 + H, Y0 * H * K3) PASO 8. Hacer Y0 = Y0 + H/6 * (K1 + 2*K2 + 2*K3 ++ K4) PASO 9. Hacer X0 = X0 + H PASO 10. Hacer I = I + 1 IMPRIMIR Y0 y TERMINAR
