249051800-metodos-numericos.pdf

c 2014 Marlo Carranza Copyright  Publicado Publicado por Marlo Carranza

book-website.com

Material de consulta, para los estudiantes del curso de Métodos Numéricos de acuerdo al sílabus diseñado por la Facultad de Ciencias e Ingeniería de la Universida Universidadd de ciencias ciencias y Humanidades. Humanidades.

impresión, Marzo 2014

Universidad de Ciencias y Humanidades Facult acultad de Ciencias Ciencias e ingeniería ingeniería

Manual de Métodos Numéricos con implementación en Matlab

Profesor Marlo Carranza Purca

Lima - Perú Marzo -2014 i

Índice general

1. Intro ducción al Matlab

3

1.1. Oper peraciones con variables . . . . . . . . . . . . . . . . . . . . .

4

1.1.1. Formato . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.1.2. Variables . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.1. 1.1.3. 3. Algu Alguna nass func funcio ione ness mate matemá máti tica cass elem elemen enta tale less . . . . . . .

7

1.2. Números Complejos . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3. 1.3. Int Introdu roducc cció iónn al Álge Álgebr braa linea ineall numér uméric icaa . . . . . . . . . . . . . .

9

1.3. 1.3.1. 1. Div Diverso ersoss tipos tipos de matr matric ices es y pro propi pied edad ades es . . . . . . . . . 10 1.3. 1.3.2. 2. Vect ectores ores y matr matric ices es con con Mat Matlab lab . . . . . . . . . . . . . . 12 1.3. 1.3.3. 3. Vect ectores ores y matr matric ices es por por bloq bloque uess . . . . . . . . . . . . . . 18 1.3.4. Algo de sistema emas de ecua cuacion cionees . . . . . . . . . . . . . . 20 1.4. Gráficos con MatLat

. . . . . . . . . . . . . . . . . . . . . . . 21

1.4. 1.4.1. 1. Gráfi Gráfico coss en coord coorden enad adas as cart cartes esia iana nass . . . . . . . . . . . . 21 1.5. Polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.5.1. Evaluaci ación de un pol polinom nomio . . . . . . . . . . . . . . . . 28 1.5.2.

Produ oducto de pol polinomios . . . . . . . . . . . . . . . . . . 29

1.5.3.

División de pol polinomios . . . . . . . . . . . . . . . . . . 29

1.5.4. Raíces de pol polinomios . . . . . . . . . . . . . . . . . . . . 30 1.5. 1.5.5. 5. Deri Derivvada ada e inte integr gral al de un polin polinom omio io . . . . . . . . . . . . 30

ii

iii

1.6. Scripts y funciones . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.7. Matemática Simból bólica

. . . . . . . . . . . . . . . . . . . . . . 33

2. Aproximación y errores

2.1. Tipos de errores 2.1.1.

43

. . . . . . . . . . . . . . . . . . . . . . . . . . 44

Número de cifras signific nificaativas . . . . . . . . . . . . . . 47

2.2. 2.2. Repr Repres esen enta taci ción ón de númer úmeros os en la comp comput utad ador oraa . . . . . . . . . 48 2.2.1. Sistemas numéricos . . . . . . . . . . . . . . . . . . . . . 49 2.2.2. Aritmética del pun punto flotan otantte . . . . . . . . . . . . . . . 51 2.2. 2.2.3. 3. Norm Normal aliz izaci ación ón de los los númer números os de punto punto flota flotant ntee . . . . . 52 2.2. 2.2.4. 4. Oper Operac acio ione ness en pun punto flot flotante ante . . . . . . . . . . . . . . . 53 2.3. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3. Resolución de ecuaciones

57

3.1. Método odo de la Bisección

. . . . . . . . . . . . . . . . . . . . . . 59

3.1.1. Algori oritmo de la bisección ción . . . . . . . . . . . . . . . . . . 60 3.1.2. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.2. Método odo del punto fijo . . . . . . . . . . . . . . . . . . . . . . . . 64 3.2. 3.2.1. 1.

Algo Algori ritm tmoo de itera teraci ción ón de pun punto fijo fijo . . . . . . . . . . . 65

3.2.2. Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.2.3. Existencia y unicidad . . . . . . . . . . . . . . . . . . . . 66 3.2.4. Ejercicio . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.3. Méto do de Newton . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.1. Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . 72 3.3. 3.3.2. 2. Con Converge ergenc nciia del del méto método do de Newt ewton . . . . . . . . . . . 73 3.3. 3.3.3. 3. Vari ariante antess del del mét método odo de Newt ewton . . . . . . . . . . . . . 74 3.3.4. Ejercicios . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4. Sistemas de ecuaciones lineales

78

iv

4.1. Méto dos directos . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.1. 4.1.1. 1. Sist Sistem emas as diag diagon onal ales es y tria triang ngul ular ares es . . . . . . . . . . . . . 78 4.1.2. Eliminació ción Gaus aussiana ana . . . . . . . . . . . . . . . . . . . 80 4.2. Método odos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.2.1. Método odo de Jacobi . . . . . . . . . . . . . . . . . . . . . . 82 4.2.2. Método odo de Gauss Seidel . . . . . . . . . . . . . . . . . . 83 4.2.3. Implementación

. . . . . . . . . . . . . . . . . . . . . . 87

5. Resolución de sistemas no lineales

91

5.1. Méto do de Newton . . . . . . . . . . . . . . . . . . . . . . . . . 91 6. Interpolación Numérica

94

6.1. 6.1. Polin olinom omio io interpo erpola lado dorr de Lagr Lagran ange ge . . . . . . . . . . . . . . . . 94 6.1.1. Form ormulaci ación de Lagrange . . . . . . . . . . . . . . . . . . 95 6.1.2. Implementación . . . . . . . . . . . . . . . . . . . . . . . 98 6.2. Formulación de Newton . . . . . . . . . . . . . . . . . . . . . . 10 102 6.2.1. Interpol polación linea neal . . . . . . . . . . . . . . . . . . . . . 10 102 6.2. 6.2.2. 2. Int Interpo erpola laci ción ón cuad cuadrá ráti tica ca . . . . . . . . . . . . . . . . . . 104 104 6.2.3. Forma general de los polinomios de interpolación de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.2.4. Implementación . . . . . . . . . . . . . . . . . . . . . . . 108 6.3. Interpol polación por por partes . . . . . . . . . . . . . . . . . . . . . . 110 6.3.1. Trazad zadores ores linea neales . . . . . . . . . . . . . . . . . . . . . 11 110 6.3.2. Trazad zadores ores cuad uadráticos cos . . . . . . . . . . . . . . . . . . . 112 7. Integración numérica

115

7.1. Fórmula del trapec pecio . . . . . . . . . . . . . . . . . . . . . . . . 11 115 7.1. 7.1.1. 1. Fórm Fórmul ulaa del del trape trapeci cioo comp compue uest staa . . . . . . . . . . . . . . 118 118 7.1.2. Implementación . . . . . . . . . . . . . . . . . . . . . . . 120 7.2. Fórmula de Simpson . . . . . . . . . . . . . . . . . . . . . . . . 120

v

7.2.1. Fórmula de Simpson 1/3 . . . . . . . . . . . . . . . . . . 121 7.2. 7.2.2. 2. Fór Fórmula ula de Simp Simpso sonn comp compue uest staa . . . . . . . . . . . . . . 121 121 7.2.3. Implementación . . . . . . . . . . . . . . . . . . . . . . . 122 7.2. 7.2.4. 4. Ejer Ejerci cici cio, o, mecá mecáni nica ca de la frac fractu tura ra . . . . . . . . . . . . . 123 123 8. Méto dos de Runge - Kutta

126

8.1. Método odo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . 127 8.1. 8.1.1. 1. Ejer Ejerci cici cio, o, mecá mecáni nica ca de la frac fractu tura ra . . . . . . . . . . . . . 130 130

Introducción

Un método numérico es un procedimiento mediante el cual se obtiene casi siempre de manera aproximada, la solución de ciertos problemas realizando cálculos algebraicos y lógicos. Tal procedimiento consiste de una lista finita de instrucciones precisas que especifican una secuencia de operaciones algebraicas y lógicas (algoritmo), que producen una aproximación de la solución del problema (solución numérica) o un mensaje acerca del problema. En los últimos treinta años la simulación se ha convertido, por méritos propios, en uno de los tres pilares sobre los que se avanza en la ciencia (junto con la teoría y la experimentación). Esto es fruto del aumento espectacular de la capacidad de cálculo y de la disminución del coste de los ordenadores, pero también del avance en la formulación de métodos numéricos cada vez más precisos precisos Hay muchas empresas comercializando códigos de simulación e investigadores trabajando para construir modelos numéricos fiables. Por esto, los avances en las últimas décadas han sido espectaculares, no sólo en precisión, sino también en el tipo de problemas que ahora somos capaces de abordar. Pero queda mucho por hacer, la realidad es extremadamente compleja y los modelos matemáticos que se emplean han de enriquecerse. El objetivo es que la simulación sirva para comprender la realidad y para predecir su comportamiento. Los métodos numéricos, en general, tienen el atractivo de ofrecer una solución a problemas muy complejos que no pueden resolverse analíticamente, y en un entorno entorno controlado controlado de difícil difícil manipulación manipulación experimental. experimental. Su debilidad debilidad está, a nuestro criterio, en creerse cualquier solución obtenida. Todos los métodos tienen error y se necesita capacidad crítica parar valorar los resultados numéricos numéricos y comprender comprender los modelos mo delos matemáticos matemáticos subyacent subyacentes. es. Los códigos de cálculo proporcionan resultados espectaculares y muy vistosos, no necesariamente correctos. Un fallo en la selección del modelo matemático, o del método numérico puede dar lugar a resultados completamente erróneos. Entre otras de las ventajas que nos brindan los métodos numéricos, esta el ahorro de recursos, por ejemplo, ingenieros dedicados al cálculo, en 1 Mg. Marlo Carranza Purca Purca

2 instalaciones de ensayo, y en el tiempo de cálculo ha sido tremendo. En segundo lugar, los métodos numéricos y la simulación, en general, permiten estudiar situaciones inexistentes, realizar ensayos virtuales y probar prototipos, incluso en condiciones extremas sin peligro y con un costo reducido. Si bien la simulación no puede acabar con los estudios teóricos ni la experimentación, no cabe duda de que los complementa, y nos ayuda a comprender comprender mejor la naturaleza. naturaleza. El lenguaje actual de los métodos numéricos es, en general, el análisis funcional; en menor medida, y dependiendo de su aplicación, también juegan un papel muy relevante la geometría diferencial, el álgebra lineal y todo lo relacionado con la programación. El diseño y análisis de métodos numéricos es una rama de las matemáticas en la que se combinan aspectos muy diversos y cada campo de aplicación se apoya en mayor medida en una parte de las matemáticas.

Mg. Marlo Carranza Purca Purca

Manual Manual de métodos numérico numéricos s U.C.H. U.C.H.

1

Introducción al Matlab

Durante las últimas décadas, el ordenador se ha convertido en una de las herramientas más potentes y útiles de que dispone el ingeniero, científico y muchas personas. Su utilización abarca desde la fase de diseño y validación experimental en un laboratorio, hasta la fase de construcción o producción industrial, pasando por la confección de planos y la redacción de los pliegos de condiciones en los que se utilizan diferentes equipos. Paralelamente a este auge también ha aparecido la necesidad de recurrir a diferentes, y cada vez más sofisticados, métodos numéricos en varias de las anteriores fases, este conocimiento le facilitara la comprensión de los lenguajes de programación así como el análisis e interpretación tanto de los resultados obtenidos como de los posibles errores de programación. MatLab es una herramienta increíblemente poderosa, pero para usarlo con seguridad tiene que ser capaz de entender cómo funciona y ser muy preciso al introducir comandos. Cambiar la forma de introducir un comando, aunque sutilmente puede cambiar completamente su significado. El objetivo principal de esta introducción es enseñarle a conversar con MatLab y comprender sus respuestas. Es posible interactuar con MatLab utilizando un "manual de instrucciones", lo cual está bien si la respuesta es lo que usted espera. Sin embargo, es mucho mejor aprender el idioma para que pueda entender la respuesta. Además es esencial aprender la gramática ( o la sintaxis ), lo que es tal vez más importante con los lenguajes de programación programación convencionales; convencionales; MatLab utiliza utiliza un intérpret intérpretee para entender entender lo que le escribimos, que puede darnos sugerencias sobre dónde podría haber un error, a veces lo que se ha escrito tiene sentido para Matlab, pero no significa lo que usted espera, así que hay que tener cuidado. Es crucial formular las ideas con claridad ( en su cabeza o mejor aun en papel) antes de intentar traducirlas en MatLab (o cualquier otro idioma). El nombre MatLab es una abreviatura de las palabras Matrix Laboratory. MatLab es un sistema interactivo para cálculos científicos y de ingeniería basado en las matrices. Con el se pueden resolver complejos problemas numéricos sin necesidad de escribir un programa específico para ello. Además, el programa MatLab dispone, dependiendo de la versión, de diferentes módulos (Toolboxes) (Toolboxes) que permiten resolver resolver problemas problemas específicos. En esta introducción vamos a revisar las principales herramientas que posee 3 Mg. Marlo Carranza Purca Purca

1.1. Operaciones con variables

4

MatLab: manejo de variables, el uso de funciones matemáticas, las operaciones entre arrays y el diseño de gráficas, esto nos permitirá familiarizarnos de manera mas rápida de su entorno. además que posibilitará hacer un reconocimiento acerca de las reglas ortográficas y de las respectivas sintaxis de sus comandos y funciones más importantes.

1.1. 1.1. Operaci Operacione oness con variab ariables les A semejanza de C, C ++ , MatLab MatLab tiene ciertas características características 1. Hace diferencia entre las mayúsculas y las minúsculas en cuanto a los nombres de las variables, estas pueden empezar siempre con una letra y pueden constar de hasta 31 letras y números, considerándose al carácter - como como una letr letra. a. 2. Si una una expr expresi esión ón term termina ina en pun punto to y coma ; se calcu calcula la su su resul resultad tado, o, pero no se muestra en pantalla, y si esta se omite, entonces Matlab ejecutara la expresión y la mostrara en pantalla. 3. El carácter % se utiliza para insertar insertar comentarios comentarios.. Todo Todo lo que sigue (en la misma línea) es ignorado. 4. Se pueden recuperar comandos anteriores navegando con las flechas ↑ y ↓. 5. clc , borra borra el conteni contenido do de la ventan ventanaa de trabajo. 6. clear , borra todas las variables del espacio de trabajo, pero no la ventana de trabajo. 7. clear x,y, borra solo las variables x, y del espacio de trabajo, pero no la ventana de trabajo. 8. Exit , Se usa para salir de Matlab. Además hay algunas variables y constantes especiales que se utilizan por defecto: 1. ans: Es la variable que se utiliza en los resultados. En la operación siguiente se puede recuperar este resultado volviendo a escribir ans. Esta variable se modifica en cuanto haya un nuevo resultado. 2. pi: El número π . ( No hay una variable para el número e , pero se podría definir e = e = exp exp(1) (1) ) .

√ −1

3. i ó j repr repres esen enta tann la la uni unida dadd ima imagi ginar naria ia


.



5

4. eps: Es el número más pequeño que utiliza el ordenador. 5. Inf: Infinito, aparece si hacemos 1/0 . 6. NaN: Mensaje de error (Not a Number), por ejemplo, 0/ 0 /0 . Las operaciones aritméticas entre escalares son efectuadas de manera natural usando los siguientes siguientes símbolos para las respectivas respectivas operaciones: operaciones: Adición

+

Multiplicación *

Potenciación ^ División

/

Ejemplo 1.

1. >> 3+4 3+4 ans = 7

2. >> 4*8 4*8 ans = 32

3. >> 5^2 5^2 ans = 25

4. >> 1/0 1/0 Warnin Warning: g: Divide Divide by zero. zero. ans = Inf

ans = 17

8. >> pi ans = 3.1416

9. Comprueba la diferencia entre las operaciones 4/4 + 6, 6 , 4/(4 + 6) 3^5*2 , 3^(5*2)

veamos

a ) >> 4/4+6 4/4+6 ans = 7

b ) >> 4/(4+6) 4/(4+6)

5. 0/0 Warnin Warning: g: Divide Divide by zero. zero. ans = NaN

6. >> x=3+4 x=3+4

ans = 0.4000

c ) >> 3^5*2 3^5*2 ans = 486

d ) >> 3^(5*2) 3^(5*2)

x = 7

7. >> x+10 x+10


ans = 59049



1.1. 1.1.1. 1.

6

Forma ormato to

En general los resultados numéricos se muestran con cuatro cifras decimales, aunque todas las operaciones se ejecutan en doble precisión . Si se desean las salidas con toda la precisión disponible se debe insertar la instrucción instrucción format format long ; a partir de este punto punto todos los resultados resultados se mostraran con 16 cifras significativas. La instrucción, format rat , devuelve a partir de este lugar los resultados como la division de dos enteros. Además, Además, format short short , devuelve devuelve la forma estándar de los cuatro decimales. decimales.

Ejemplo 2.

Veamos como representa Matlab el número π >> form format at long long ; pi ans = 3.14159265358979

>> form format at rat; rat; pi ans = 355/113

>> form format at shor short; t; pi ans = 3.1416

>> format;pi format;pi ans = 3.1416

1.1. 1.1.2. 2.

Varia ariabl bles es

Matlab no necesita la declaración de variables como en un lenguaje tradicional. En principio todas las variables son reales, basta hacer uso de ellas para que queden declaradas.




7

Ejemplo 3.

1. >> a=1; a=1; b=2; b=2; c=3; c=3;

2. Dad Dadas as las las varia ariabl bles es y sus sus respectivos contenidos, a = 2, b = 5 y c = 10, calcule la expresión:

>> a+b a+b ans = 3 >> a*c a*c ans = 3

E = 4ab2

10c − 10c b

veamos:

>> a+b+c a+b+c ans = 6

>> a=2; a=2; >> b=5; b=5; >> c=10; c=10; >> E=4*a*b^2-10*c/b E=4*a*b^2-10*c/b E = 180

>> a*b*c a*b*c ans = 6

Matlab no hace diferencia si las variables son enteras reales o de doble precisión, como suele ocurrir en otros lenguajes de programación como Pascal, C , C ++ o Java entre los más destacados.

1.1.3.

Algunas Algunas funciones funciones matemáticas matemáticas elemen elementales tales

Matlab tiene definidos diversos tipos de funciones, todas ellas clasificadas según el objeto o tipo de problema a resolver, los detalles de estas se irán desarrollando en el transcurso de las siguientes páginas. Veamos una lista de las funciones elementales mas usadas. Función

Matlab

Función

Matlab

Función

|x|

abs(x)

ex

exp(x)

√ x

sin(x)

cos( cos(x)

cos(x) tan(x)

tan(x)

pow2(x)

loge (x)

log(x)

log10(x)

sen(x sen(x)

2z

signo(x) sign(x) senh(x)

log10 (x)

Matlab

sqrt(x)

arcos(x) acos(x) arctan(x) atan(x)

sinh(x) cosh(x cosh(x) cosh(x) tanh(x tanh(x)

tanh(x)

cos(x) acosh(x) arctanh(x) atanh(x) arcsinh(x) asinh(x) arc cos(



1.2. Números Complejos

8

Ejemplo 4.

Verificar la identidad trigonométrica: sen2 x + cos2 x = 1

Nuestro código, sería: >> x=0; x=0; >> r=sin(x)^2+cos(x)^2 r=sin(x)^2+cos(x)^2 r =1

r = 1 >> x=25; x=25; >> r=sin(x)^2+cos(x)^2 r=sin(x)^2+cos(x)^2 r =1

>> x=1; x=1; >> r=sin(x)^2+cos(x)^2 r=sin(x)^2+cos(x)^2

1.2. 1.2. Núme Número ross Comp Comple lejo joss Cuando se trata de números complejos de la forma z = a = a + bi, donde a y b son reales, Matlab pone a disposición las siguientes operaciones que reducen las operaciones: operaciones: Función

Matlab

|z|

abs(z)

Función

Matlab

Imag(z) imag(z)

Ángulo de z angle(z) Re(z) real(z) z


conj(x)


1.3. Intro ducción al Álgebra lineal numérica

9

Ejemplo 5.

Dado el número complejo z = = 3 + 4i 4i , calcule: Ángulo(z) (z),, y z z, Imag( Imag (z ), Ángulo

||

Veamos, se tiene : >> z=3+4* z=3+4*i i z = 3.0000 3.0000 + 4.0000 4.0000i i >> conjugada= conjugada=conj( conj(z) z) conjugada conjugada = 3.0000 3.0000 - 4.0000 4.0000i i >> Real=real( Real=real(z) z) Real Real =

3 >> Anglo=ang Anglo=angle(z) le(z) Angl Anglo o = 0.9273 >> modulo=ab modulo=abs(z) s(z) modulo = 5

Ejemplo 6.

Calcule el valor de la siguiente expresión: w = |z |eθi , si z = 1 + i + i además, verifique verifique si w = z = z . veamos : >> z=1+i z=1+i z = 1.0000 1.0000 + 1.0000 1.0000i i >> theta=angl theta=angle(z); e(z); >> modulo=abs modulo=abs(z); (z); >> w=modulo*exp(theta*i) w=modulo*exp(theta*i) w = 1.0000 1.0000 + 1.0000 1.0000i i

1.3. 1.3. Introd Introduc ucció ción n al Álge Álgebra bra lineal lineal numéric umérica a En la presente sección se recoge una serie de definiciones y resultados relativos a matrices, aunque una buena parte de ellos pueden ser conocidos, el hecho que esta recopilación busca fijar la notación utilizada, refrescar la memoria y tener a la mano resultados que vamos a requerir más adelante.




1.3.1. 1.3.1.

10

Divers Diversos os tipos tipos de matr matrice icess y propie propiedad dades es

En todo lo que sigue, K denota el cuerpo R de los números reales o cuerpo de los complejos, V = Kn, n ∈ N.

C el

Definición 1.

Una matriz es una colección de elementos aij siguiente forma

A =

  

a11

a12

a21

a22

.. .

...

am1 am2

··· ···

a1n

...

...

···

amn

a2n

∈ K, dispuestos de la

  

es una, matriz real si K = R y compleja cuando K = C . La matriz A tiene m filas y n columnas, es decir tiene dimension m × n en particular: 1. Un vector columna es una matriz de dimension m × 1 2. Un vector fila es una matriz de dimension 1 × n Notación 1.1.

1. Una matriz A de dimensión m × n de elementos aij ∈ K se denota m,n A = (aij )i,j=1 i,j =1 (i fila , j columna ). También denotaremos por (A)ij o A(i, j ) el elemento de de A que ocupa la fila i y la columna j 2. Escribiremos la matriz A = (a1, a2 , . . . , an) donde ai ∈ Km representa la columna i-ésima de A para i = 1, 2, . . . , n 3. El conjunto de todas las matrices de orden o dimension m × n con elementos en K con las operaciones de adición y multiplicación por un escalar lo denotaremos como Mm×n




11

Definición 2. m,n m,n Si A = (aij )i,j=1 i,j =1 y B = (bij )i,j=1 i,j =1 son matrices de dimensión m × n

1. Se define la suma de matrices como la matriz A + B = (cij )m,n + bij i,j=1 i,j =1 con c ij = a ij + b

y el producto de una matriz A = (aij )m,n i,j=1 i,j =1 por un escalar λ ∈ K como la matriz λA = λA = (λaij )m,n i,j=1 i,j =1

n,m 2. La matriz AT = (a ji ) j,i=1 j,i=1 se denomina matriz transpuesta de A.

3. Una matriz A se denomina matriz cuadrada o matriz de orden n. Diremos que una matriz es rectangular si no es cuadrada. 4. A es simétrica si A = A = A T . 5. La matriz A ∈ Mn es no singular si existe una matriz B ∈ Mn , tal que AB = BA B A = I = I , la matriz B es llamada la inversa de A y − se denota B = A 1.Las matrices no inversibl inversibles es son singulares. singulares. 6. Una matriz A ∈ Mn es no singular si det(a det(a) = 0 7. Si det(A det(A)  = 0 , la matriz inversa de A =

  b 

−x

 −y  

    a y

es A−1 =

x b

a

det(A det(A)

8. La Adjunta de una matriz A ∈ Mn es la traspuesta de la matriz que resulta de sustituir cada uno de sus elementos por su cofactor, y se denota adj adj (A) adj (A) 9. Si det(A det(A) =  0 se cumple que A−1 = adj |A| 10. A es ortogonal si A T = A −1 es decir AAT = A T A = I = I

A = diag = diag((aij ) = diag( diag (a11 , a22 , . . . , ann) =


  

a11

0

0

a22

.. . 0

...

0

··· ···

0

...

...

···

ann

0

  



12

Definición 3. m,n m,n Si A = (aij )i,j=1 i,j =1 y B = (bij )i,j=1 i,j =1 son matrices de dimensión m × n

1. La matriz A es triangular superior (respectivamente inferior si ) aij = 0 para i > j (respectivamente i < j )

 j , en este caso se 2. La matriz A es diagonal si aij = 0 cuando i = denota 3. La matriz A es de diagonal estrictamente dominante si | aii | > n Σ j=1 j =1,j ,j  =i |aij | para i = 1, 2, . . . n 4. Se denomina traza de la matriz A al número tra( tra(A) = Σni=1 aii

Teorema 7.

Sean A, B matrices cuadradas de la misma dimension, se cumple 1. tra( tra(AB) AB ) = tra tr a(BA) BA ) 2. tra( tra(A + B) B ) = tra tr a(A) + tra + tra((B ) 3. det(AB det(AB)) = det(BA det(BA)) = det(A det(A)det(B )det(B ) 4. det(λA det(λA)) = λ n det(A det(A)

1.3.2. 1.3.2.

Vectore ectoress y matr matrice icess con con Matla Matlab b

Una de las cualidades de MatLab es el manejo de los vectores y las matrices, veremos algunos ejemplos




13

Ejemplo 8.

1. Definiremos el vector fila >> a=[1 2 3 4 5 ] a = 1 2 3 4 5

para indicar la n-ésima componente de vector escribimos >> a(3) a(3) ans = 3

Obser Observe ve que que se mues muestr traa el arreglo a como un vector fila con 5 columnas, además cada elemento puede ser separado por comas o espacios como el ejemplo mostrado. Si usted desea ver los primeros m elementos del vector a escribimos >> a(1:3) a(1:3) ans = 1

2

2. A con continua inuaci ción ón tenem enemos os otros formatos para definir vectores. >> u=1:5 u=1:5 u = 1 2 3

4

5

3. Usaremos linspace para generar un vector >> a=linspace(1,5,5) a=linspace(1,5,5) a = 1 2 3 4 5

4. Pero también se pueden definir vectores columna >> v=[1; 2; 3; 4; 5 ] v = 1 2 3 4 5

3

Matrices

Para generar matrices tenemos que introducir vectores fila de la misma dimensi dimensión, ón, fila fila por fila, fila, se usa usa el punto punto y coma coma ; para para separa separarr las las filas. filas.




14

Ejemplo 9.

Generemos las matrices usando MatLab. A =

B=

  

0 2 1 4 5 0 1 0 2 0 3 7

  

  

1

2

3

4

5

6

7

8

9 10 11 12

  

tenemos

>> A=[1 2 3 4; 5 6 7 8; 9 10 11 12] A = 1 2 3 4 5 6 7 8 9 10 11 12

>> B=[0 2 1 4; 5 0 1 0 ; 2 0 3 7] B = 0 2 1 4 5 0 1 0 2 0 3 7

Operaciones

Existen en Matlat dos tipos de operaciones aritméticas, las operaciones aritméticas matriciales, que se rigen por las reglas del Álgebra lineal y las operaciones aritméticas a elemento, que se realizan elemento a elemento, veamos la tabla de que describe las operaciones. Operación

Símbolo

Operación

Símbolo

Suma de vectores y matrices

+

Resta de vectores y matrices

Producto de matrices

*

Multiplicación a elemento

.*

Potenciación de matrices

^

Potenciación a elemento

.^

-

Cociente matricial

/

División a elemento

. /

Cociente matricial

\


.\




15

Ejemplo 10.

1. Para generar matrices tenemos que que introducir vectores fila de la misma cantidad de componentes, componentes, fila por fila. se usa punto y coma para separa las filas.

     

A =

1 2 3 4 5 6 7 8 9

  

0 2 4 6

B=

1 3 5 7 0 1 2 3

  

B=[0 2 4; 1 3 5; 0 1 2] B = 0 2 4 1 3 5 0 1 2

2. Verem eremos os opera operaci cion ones es con con matrices

4 8 9

>> EA=2*A EA=2*A EA = 2 4 8 10 14 16


7 11 11

6 12 18

11 29 47

20 53 86

>> T=A’ T=A’ % transp transpues uesta ta T = 1 4 7 2 5 8 3 6 9

A=[1 2 3; 4 5 6; 7 8 9] A = 1 2 3 4 5 6 7 8 9

>> S=A+B S=A+B S = 1 5 7

>> P=A*B P=A*B P = 2 5 8

>> PM=A^3 PM=A^3 % potenc potencia ia PM = 468 576 684 1062 1062 1305 1305 1548 1548 1656 1656 2034 2034 2412 2412 >> trac trace( e(A) A) % Traz Traza a ans = 15

3. Determinante de una matriz D=[1 0 0;0 D = 1 0 0 >> det(D) det(D) ans = 6

2 0;0 0 3] 0 2 0

0 0 3

4. Matriz inversa >> inv(D) inv(D) ans = 1.0000 0 0 0 0.5000 0 0 0 0.3333



16

Ejemplo 11.

Veamos algunos ejemplos >> % oper operac acio ione nes s de suma suma y prod produc ucto to >> suma = A+B >> % transp transpues uesta ta suma suma = >> transpues transpuesta=A’ ta=A’ 1 4 4 8 transpuest transpuesta a = 10 6 8 8 1 5 9 11 10 14 19 2 6 10 3 7 11 >> % suma suma caso caso espe especi cial al 4 8 12 >> casoEspecial=10+suma casoEspecial=10+suma casoEspeci casoEspecial al = >> % producto de A por B 11 14 14 18 >> Producto= Producto= A*transpu A*transpuesta esta 20 16 18 18 Produc Producto to = 21 20 24 29 30 70 110 70 174 278 110 278 446 >> % prod prod esca escala lar r por por matr matriz iz >> escalarmat escalarmatriz=2 riz=2*A *A >> % pote potenc ncia ia E*E* E*E*E E escalarmat escalarmatriz riz = >> potencia= potencia=E^3 E^3 2 4 6 8 potenc potencia ia = 10 12 14 16 1 14 18 20 22 24 0 8

Operaciones a elementos

Ahora vamos a revisar las operaciones a elemento y vamos a ver las diferencias con las operaciones ya vistas, pero antes veamos la tabla con la descripción de las operaciones que vamos a utilizar. Operación

Símbolo

Multiplicación a elemento .*


Potenciación a elemento

.^


. /



17

Ejemplo 12.

>> E=[1 2 ; 3 4], F=[2 4 ; 8 16] E = >> % Divi Divisi sion on a elem elemen ento to 1 2 >> E./F E./F 3 4 ans = F = 0.5000 0.5000 2 4 0.3750 0.2500 8 16 >> % Prod Produc ucto to a elem elemen ento to >> E.*F E.*F ans = 2 8 24 64

>> %poten %potencia ciacio cion n a elemen elemento to >> P=E.^3 P=E.^3 P = 1 8 27 64

Ejemplo 13.

Veamos algunos ejemplos de operaciones a elemento, Dadas las matrices A y B >> A=[2 4; 8 10] A = 2 4 8 10 >> B= B=[2 3; 3; 1 B = 2 3 1 1

1]

>> %poten %potencia cia a elemen elemento to >> P=A.^2 P=A.^2 P = 4 16 64 100 >> %produ %producto cto a elemen elemento to >> A*B A*B ans = 8 10 26 34


>> %divis %division ion a elemen elemento to >> A./B A./B ans = 1.0000 1.3333 8.000 .0000 0 10.0 10.000 000 0 >> %potencia %potencia variable variable >> %a elemen elemento to >> x=[1 x=[1 2 3]; 3]; x.^x x.^x ans = 1 4 27 >> %potencia %potencia variable variable %a elemento elemento >> y=[1 2; 3 4]; y.^y ans = 1 4 27 256



1.3.3. 1.3.3.

18

Vectore ectoress y matric matrices es por bloqu bloques es

Selección de elementos de un vector

V(n) V(n)

Devu Devuel elvve en n-és n-ésim imoo elem elemen ento to del del vecto ectorr

V([n,m,p]) V([n,m,p]) Devuelv Devuelvee los elementos elementos del del vector vector situados situados en las posicion posiciones es n, m , p V(n:m) V(n:m)

Devu Devuelv elvee los elemen elementos tos del vecto vectorr situad situados os entre entre las posicio posiciones nes n y m

Ejemplo 14.

>> V= V=[ 2 4 6 8

10 ]; ];

% obtene obtenemos mos 3era 3era compon component ente e >> V(3) V(3) ans = 6 %obt %obten enem emos os la 1, 3 5 comp comp. . >> V([1 3 5])

ans = 2

6

10

% elem elemen ento tos s de 3 a 5 posi posici ción ón >> V(3:5) V(3:5) ans = 6 8 10

Selección de los elementos de una matriz

A(m,n)

Devuelve el elemento (m, (m, n) de la matriz.

A([n,m],[p A([n,m],[p,q]) ,q]) Devuelv Devuelvee la submatriz submatriz formada formada por la intersecci intersección ón de las filas filas n, m y las columnas p, q .

A(n,:)

Devuelve la fila n.

A(:,m)

Devuelve la columna m.




19

Ejemplo 15.

Sea la matriz >> A=[1 A=[1:4 :4; ; 5:8 5:8 ; 9:12 9:12 ] A = 1 2 3 4 5 6 7 8 9 10 11 12 >> % elem elemen ento to fila fila 2, col col 3 >> A(2,3) A(2,3) ans = 7 >> % extr extrae aemo mos s la fila fila 1 >> A(1,:) A(1,:) ans = 1 2 3 4 >> % obte obtene nemo mos s la col col 3 >> A(:,3) A(:,3) ans = 3 7

11 % obtene obtenemos mos la submat submatriz riz >> A([2,3],[ A([2,3],[2,3]) 2,3]) ans = 6 7 10 11 % obtene obtenemos mos la submat submatriz riz >> A(2:3,2:3 A(2:3,2:3) ) ans = 6 7 10 11 % obtene obtenemos mos la submat submatriz riz >> A(:,2:3) A(:,2:3) ans = 2 3 6 7 10 11

Ejemplo 16. Otras maneras de definir matrices . Matlat tiene otra formas de

definir matices, dado que introducirlas por el teclado no es muy práctico sobre todo cuando la matriz es muy grande. A = 1 4 7 >> B=[1 1

7 2 5 8 1; 2 2

>> col= col=[A [A B] col = 1 2 3 1 1 4 5 6 2 2


1 2

3 6 9 2; 3 3 3]

8

9

3

3

>> fil= fil=[A [A ; B] fil = 1 2 4 5 7 8 1 1 2 2 3 3

3

3 6 9 1 2 3



20

Más matrices y funciones matriciales

zeros(m,n)

Matriz nula de orden n × m

ones(m,n)

Matriz formada por unos de orden n × m

eye(n)

matriz identidad de orden n.

tril(A)

Parte triangular inferior de la matriz A

triu(A)

Parte triangular superior de la matriz A

magic(n)

Matriz mágica

vander(n)

Matriz de Vandermonde

inv(A)

Matriz inversa de la matriz A.

trace(A)

Traza de la matriz A.

[fil,col]=size(A)

Nos informa el número de filas y columnas.

length(A)

Nos da el mayor valo alor de la fila o de la columna.

max(A)

Proporciona el máximo de los elementos de A.

min(A)

proporciona el mínimo de los elementos de A.

1.3. 1.3.4. 4.

Algo Algo de de sist sistem emas as de de ecua ecuaci cion ones es

Para resolver sistemas lineales, en estos momentos se tienen muchas herramientas, pero en esta sección vamos a mostrar una forma de solución usando algunas funciones sencillas que nos ofrece MatLab.

  

2x1 + 3x 3 x2

− 4x = 3 x − x + x + x = −0,5 4x − 7x + 14 1 4x = 2 1

3

2

1

3

2

3

vamos a expresar el sistema en la forma

 −      −     −   −             2


3

4

x1

3

=

1

1

1

x2

0,5

4

7 14

x3

2

A

X

b


1.4. Gráficos con MatLat

21

Ejemplo 17.

Para resolver este sistema de ecuaciones realizamos las siguientes acciones con Matlab >> A=[2 3 -4;1 -1 1; 4 -7 14] % soluci solución ón >> x=A\b x=A\b A = x = 0.5000 2 3 -4 2.0000 1 -1 1 1.0000 4 -7 14 % Soluci Solución ón con linsol linsolve ve >> b=[3 ; -0.5 ; 2 ] >> x= linsolve( linsolve(A,b) A,b) x = b = 0.5000 3.0000 2.0000 -0.5000 1.0000 2.0000

1.4 1.4.

Gráfic ráfico os con MatLa atLatt

MatLat produce gráficos de dos dimensiones, así como contornos y gráficos de densidad. Se pueden representar los gráficos y listar los datos, permite el control de colores, sombreado y otras características de los gráficos, también soporta gráficos animados. Como podemos ver los gráficos producidos por MatLat tienen muchas muy buenas características, características, inclusive inclusive son portables a otros programas

1.4.1. 1.4.1.

Gráfico Gráficoss en coorden coordenada adass carte cartesia sianas nas

Estos gráficos se tratan como curvas que pasa por pares ordenados, pero finalmente Matlab lo que hace es trazar una poligonal lineal que pasa por estos puntos o pares ordenados.




22

Ejemplo 18.

>> x=[7 9 3 1 5 20 5]; >> plot(x plot(x) ) >>

1. plot(X) Representa los puntos (k, xk ). Si X es es una matriz hace lo mismo para cada columna de la matriz si X es un vector complejo, representa real( real(X ) frente a imag( imag(X ). 2. plot(X,Y) Representa los punto (X, Y ) Y ), si X y Y son matrices representa por filas o columnas los datos de X frente a los datos de Y , dependiendo si el otro otro vector es fila o columna. 3. plot(X,Y,S) Es la gráfica de plot(X,Y) plot(X,Y) con las opciones definidas definidas en se compone de tres caracteres entre tildas, el primero S , usualmente S se fija el color, el segundo fija la etiqueta o marca en el nodo el ultimo fija el carácter usado en el nodo. 4. plot(X 1, Y 2 , S 3, . . . , X , Y , S ) : Gráfica de las n curvas superpuestas. n

n

n

5. hold Permite Permite montar montar o sobreponer varios varios gráficos usando usando una sola ventana. hold hold on para para acti activvar hold hold hold hold off para para desa desact ctiv ivar ar hold hold 6. zoom Permite Permite ampliar o disminuir disminuir el gráfico, se activa con zoom on y se desactiva con zoom off . Los caracteres son respectivamente los siguientes.




Color y

Etiqueta

amarill lloo

m mag magenta enta c

cyan

r

23

. puntos

◦

círculos

trazo - sólido :

a puntos

x

x- marcas

-.

rojo

+

signo +

-- semisólidos

g

ve verde

*

estrellas

b

azul

s

cuadrados

w bl blanco

d

diamantes

k

p

estrella de cinco puntos

negro

guiones y puntos

Ejemplo 19.

√

Graficar f ( f (x) = x + 4x 4x3 sen(x sen(x) en el intervalo [0, [0, 5] >> x=0:0.2:5; x=0:0.2:5; >> y=sqrt(x)+4*x.^2.*sin(x); y=sqrt(x)+4*x.^2.*sin(x); >> plot(x,y) plot(x,y)




24

Ejemplo 20.

1. Graficar f ( f (x) = ex − x sen(x sen(x2) + 2 en el intervalo [−3, 3], con trazo de color rojo, etiquetas cuadradas y linea punteada >> x=-3:0.4: x=-3:0.4:3; 3; >> y=exp(x)-x.*sin(x^2)+2; y=exp(x)-x.*sin(x^2)+2; >> plot(x,y, plot(x,y,’rs:’ ’rs:’) )

2. Graficar las funciones y 1 = −x sen(x sen(x2) + 2; y2 = abs( abs(x) + 3 4 sen( sen(x x); y 3 = x + 0,3x colocando el titulo, descripción de los ejes y la leyenda. >> >> >> >> >> >> >> >> >>

x=-3:0.4: x=-3:0.4:3; 3; y1=-x.*sin(x.^2)+2; y1=-x.*sin(x.^2)+2; y2=abs(x)+4.*sin(x); y2=abs(x)+4.*sin(x); y3=x+0.3. y3=x+0.3.*x.^3 *x.^3; ; plot(x,y1,x,y2,x,y3); plot(x,y1,x,y2,x,y3); title(’Gr title(’Gráfico áfico de tres funciones funciones’); ’); xlabel(’E xlabel(’Eje je x’); ylabel ylabel(’E (’Eje je y ’); legend(’-xsen(x^2)+2’,’abs(x) legend(’-xsen(x^2)+2’,’abs(x)+4sen(x)’,’x+0. +4sen(x)’,’x+0.3x^3’); 3x^3’);




25

Ejemplo 21.

1. Graficaremos la función f ( f (x) = sen( sen(x) conjuntamente con la función h(x) = cos( cos (x) en el intervalo [−2π ; 2π] Los gráficos respectivos los mostraremos en la misma ventana. >> >> >> >> >> >> >> >>

x=linspace(-2*pi,2*pi,100); x=linspace(-2*pi,2*pi,100); y1=sin(x) y1=sin(x); ; y2=cos(x) y2=cos(x); ; plot(x,y1,’red’); plot(x,y1,’red’); hold hold on; on; plot(x,y2,’blue’); plot(x,y2,’blue’); title(’fu title(’funcion nciones: es: seno y coseno’); coseno’); legend(’ legend(’ seno(x)’, seno(x)’, ’cos(x)’); ’cos(x)’);

Gráfica de una función definida a trozos

2. Graficar la función f la cual viene definida a trozos

f ( f (x) =

   −

x2

si x < 0, 0 ,

1

si 0 ≤ x < 1 x + 2

si 1 ≤ x

Generamos una tabla de valores en el dominio en el que queramos dibujar la función >>x=linspace(-2,3,3000);

Y ahora definimos la función, multiplicando cada trozo por el índice lógico que describa el lugar en el que queremos dibujarlo, >>y=(x.^2).*(x<0)+1.*((0<=x)&( >>y=(x.^2).*(x<0)+1.*((0<=x)&(x<1))+(-x+2).* x<1))+(-x+2).*(1<=x); (1<=x);




26

Ejemplo 22.

Y ahora la dibujamos. Resulta conveniente hacerlo con puntos, asteriscos o cruces porque, de otra forma, no aparecerán las discontinuidades >>plot(x,y >>plot(x,y,’.’) ,’.’),grid ,grid on,title( on,title(’Func ’Función ión definida definida a trozos’) trozos’)

escribimos en matlab el código >>x=linspace(-2,3,3000); >>y=(x.^2).*(x<0)+1.*((0<=x)&( >>y=(x.^2).*(x<0)+1.*((0<=x)&(x<1))+(-x+2).*( x<1))+(-x+2).*(1<=x); 1<=x); >>plot(x,y >>plot(x,y,’.’) ,’.’),grid ,grid on,title( on,title(’Func ’Función ión definida definida a trozos’) trozos’)

Obtención de puntos desde el gráfico

Una vez que se ha realizado una gráfica, podemos necesitar conocer las coordenadas de algunos puntos de la misma. Por ejemplo, el lugar aproximado en el que están los máximos y mínimos, o si queremos añadir alguna recta o una poligonal al dibujo. Para conseguir esto, se puede utilizar el comando ginput . Escribiendo

[x, y ] = ginput( ginput(N )

Donde N es es el número de puntos cuyas coordenadas queremos obtener. Después de ejecutado el comando habrá que pulsar con el botón izquierdo del ratón sobre el dibujo tantas veces como puntos hayamos especificado. Las coordenadas de esos puntos quedarán almacenadas en las variables [x; y ].




27

Búsqueda de ceros de una función

La localización de ceros o raíces de una función está ligada a la al cambio de signo de dicha función, como lo muestra el Teorema Bolzano, el será visto en el Capítulo 3. Ejemplo 23.

Efectuar una localización gráfica de los ceros de la función f ( f (x) = x 2 e (x − 1) en el intervalo [−5;5] x=linspace(-5,5,100); >> y=exp(x).*(x.^2-1); y=exp(x).*(x.^2-1); >> plot(x,y,’ plot(x,y,’b’); b’);

Ahora hacemos un pequeño cambio x=linspace(-5,5,100); >> y=exp(x).*(x.^2-1); y=exp(x).*(x.^2-1); >> plot(x,y,’ plot(x,y,’b’); b’); hold on ; plot([-5,5],[0,0],’r’); zoom zoom on



1.5. Polinomios

1.5 1.5.

28

Polino linom mios ios

Definición 4 (Polinomio).

Una expresión que tiene la siguiente forma P ( P (x) = a 0 xn + a1xn−1 +

(1.1)

· · · + a − x + a + a n 1

n

es llamada polinomio.

Los polinomios son usados ampliamente en los tópicos de matemáticas y cursos de naturaleza computacional, como el curso de métodos numéricos esto se justifica por la facilidad de manipulación y de implementación en cualquier lenguaje de programación, por ejemplo es fácil de evaluar, derivar e integrar integrar analíticamente analíticamente como numéricame numéricamente, nte, claro que los polinomios polinomios tienen algunos inconvenientes como la tendencia a oscilar cuando el grado es grande, en Matlab se introduce a través de sus coeficientes pero considerando el polinomio completo y ordenado en forma decreciente Ejemplo 24.

El polinomio P ( 4x2 − 6x + 7 en Matlab se introduce así P (x) = x 3 + 4x >> p=[ 1 4 -6 7 ]

1.5.1. 1.5.1.

Evalu Evaluaci ación ón de un polinomi polinomio o

polyval(p,x) evalúa el polinomio P en x P es un vector de longitud n + 1 , cuyos elementos son los coeficientes del polinomio Ejemplo 25.

Vamos a evaluar el polinomio P ( P (x) = x 3 + 4x 4x2 − 6x + 7 en diferentes valores para x >> sumacoeficientes=polyval(p,1) sumacoeficientes=polyval(p,1) sumacoefic sumacoeficiente ientes s = 6 >> terminoindep=polyval(p,0) terminoindep=polyval(p,0) terminoind terminoindep ep =7



1.5. Polinomios

1.5. 1.5.2. 2.

29

Prod Produc ucto to de polin polinom omio ioss

Usamos la instrucción instrucción conv(p,q) que multiplica los polinomios p y q Ejemplo 26.

Cargamos Cargamos los polinom p olinomios ios p(x) = x + x + 1 y q (x) = x − 1 >> p =[ =[1 1] >> q=[1 q=[1 -1] -1] % Calculam Calculamos os p=(x+1)^2= p=(x+1)^2=x^2+2 x^2+2x+1 x+1 >> binomiocuadrado=conv(p,p) binomiocuadrado=conv(p,p) binomiocuadrado = 1 2 1 % Calcul Calculamo amos s la difere diferenci ncia a de cuadra cuadrados dos (x+1)( (x+1)(x-1 x-1) ) >> difcuadrados=conv(p,q) difcuadrados=conv(p,q) difcuadrados = 1 0 -1 % calculam calculamos os (x+2)(x+3) (x+2)(x+3) >> a=[1 a=[1 2]; 2]; >> b=[1 b=[1 3]; 3]; >> multiplicamos=conv(a,b) multiplicamos=conv(a,b) multiplicamos = 1 5

1.5. 1.5.3. 3.

6

Divi Divisi sión ón de polin polinom omio ioss

La instrucción [Q,R]=deconv(D,d) realiza realiza la división división de los polinomios polinomios D(x) entre d(x), obteniéndose el cociente Q y el residuo R Ejemplo 27.

Vamos a dividir el polinomio p(x) = x3 − 1 entre d(x) = x2 + x + x + 1 para lo cual hacemos >> p=[1 0 0 -1]; >> d=[1 1 1]; >> [q,r]=deco [q,r]=deconv(p, nv(p,d) d) q = 1 -1 r = 0 0 0


0

% cociente % residuo


1.5. Polinomios

1.5. 1.5.4. 4.

30

Raíc Raíces es de polin polinom omio ioss

roots(P)

calcula las raíces raíces del polinomio polinomio P ( P (x)

Ejemplo 28.

Vamos a calcular las raíces del polinomio p(x) = x3 + 6x2 + 11 1 1x + 6 para ello realizamos los pasos siguientes >> p=[1 6 11 6] >> raices=roo raices=roots(p) ts(p) raices raices = -3.0000 -2.0000 -1.0000

poly(vectorRaiz)

Calcula Calcula el polinomio polinomio de raíces raíces vectorRai vectorRaizz

Ejemplo 29.

Vamos a calcular el polinomio p(x) cuyas raíces vienen dadas por el vector vectorRaiz vectorRaiz vectorRaiz=[-3 =[-3 -2 -1]; >> polinomio=poly(vectorRaiz) polinomio=poly(vectorRaiz) polinomio = 1 6 11 6

1.5.5. 1.5.5.

Deriv Derivada ada e integ integral ral de de un polino polinomio mio

derivada ada de un polinomio polyder(p) calcula la deriv Ejemplo 30.

Calcular la derivada de f ( f (x) = 4x3 − x2 − 3x − 2 >> p=[4 p=[4 -1 -3 -2]; -2]; >> derivada=polyder(p) derivada=polyder(p) deriva derivada da = 12 -2 -3



1.5. Polinomios

31

integral de un polinomio. polinomio. polyint(p) calcula la integral Ejemplo 31.

Calculamos la integral de f ( f (x) = 4x3 − x2 − 3x − 2 >> p=[4 p=[4 -1 -3 -2]; -2];fo form rmat at rat rat >> integral=polyint(p) integral=polyint(p) integr integral al = 1 -1/3 -3/2 -2 0

Ejemplo 32.

Graficar el polinomio p(x) = x 3 − x p=[1 0 -1 0]; >> x=linspace(-3,3,200); x=linspace(-3,3,200); >> y=polyval( y=polyval(p,x); p,x); >> plot(x,y) plot(x,y) >> legend(’po legend(’polinom linomio io x^3-x’) x^3-x’) >> raices=roo raices=roots(p) ts(p); ; >> hold hold on >> plot(raices,zeros(1,3),’or’) plot(raices,zeros(1,3),’or’) >> legend(’po legend(’polinom linomio io x^3-x’,’r x^3-x’,’raices aices’) ’) >> plot(x,0) plot(x,0)



1.6. Scripts y funciones

32

Ejemplo 33.

Graficar el polinomio p(x) = x3 + 6x 6 x2 + 11x 11 x + 6 para ello realizamos los pasos siguientes siguientes >> >> >> >> >> >> >> >> >> >>

p=[1 6 11 6]; x=linspace(-6,2,600); x=linspace(-6,2,600); y=polyval( y=polyval(p,x); p,x); plot(x,y) plot(x,y) legend(’po legend(’polinom linomio io x^3+6x^2+ x^3+6x^2+11x+6 11x+6’) ’) raices=roo raices=roots(p) ts(p); ; hold hold on plot(raices,zeros(1,3),’or’) plot(raices,zeros(1,3),’or’) legend(’polinomio legend(’polinomio x^3+6x^2+11x+6’,’raices’) x^3+6x^2+11x+6’,’raices’) plot(x,0) plot(x,0)

1.6. 1.6. Scri Script ptss y func funcio ione ness Creación de archivos tipo m

Hasta el momento solo hemos visto pequeños ejemplos en cuanto al código, pero para códigos mas complejos donde se tenga que hacer correcciones esta tarea se tornara muy difícil, MatLab soluciona este problem a creando un archivo de de texto sin formatear (script ) en el cual se escriben todas las sentencias MatLab ejecutara todas las ordenes, especificando el archivo respectivo como si hubiesen sido escritos directamente en la ventana de comandos de MatLat. Archivos como funciones

MatLat nos proporciona un formato para crear nuestras propias funciones, estas pueden definirse en un archivo de texto (archivo con extensión M ), que a su vez lo tenemos a disposición en cualquier instante, están funciones tienen la virtud de tener argumentos de salida.



1.7. Matemática Simbólica

33

Implementación de un programa

Calcule el valor de la hipotenusa de un triángulo rectángulo a partir de sus dos catetos 1. Cree un script . Para esto vamos a crear el archivo m, desde el menu de archivo eligiendo eligiendo la file seleccionamos new y vamos a crear un nuevo archivo opción m file luego digitamos el programa

Con [ctrl]+[s] [ctrl]+[s] guardamos como hipot.m , ejecución en la ventana de comandos >> hipot hipot cate cateto to a= 3 cate cateto to b= 4 hip = 5

2. Cree una función. Para la función vamos a crear el archivo m,

Con [ctrl]+[s] [ctrl]+[s] guardamos como hipotenusa.m , ejecución en la ventana ventana de comandos comandos >> h=hipotenusa(3,4) h=hipotenusa(3,4) h = 5

1.7. 1.7.

Mate Matemá máti tica ca Sim Simbóli bólica ca

Matlab nos ofrece muchas funciones para realizar operaciones y cálculo simbólico lo cual nos permite por ejemplo factorizar, simplificar expresiones, Mg. Marlo Carranza Purca Purca



34

cálculo de raíces de polinomios, evaluar límites, derivadas integrales, transformada transformada de Laplace Laplace entre entre otros muchas muchas opciones muy importantes. importantes. Variables Simbólicas

Función

Descripción

syms yms x, y,. . . , z

Convierte Convierte las variables en simbólicas

syms yms x, y,. . . , z real real

Convie Convierte rte las variab variables les en simból simbólicas icas con valor valores es reales reales

syms yms x, y,. . . , z unreal unreal Convierte Convierte las variabl variables es en simból simbólicas icas con con valores valores no reales reales syms

Lista de variables simbólicas en el espacio de trabajo

x = syms = syms(( x )

Convierte Convierte la variable x en simbólica

pretty( pretty (x)

Convierte Convierte la expresión matemática x en escritura matemática

simplify( simplify(E )

Simplifica la expresión E

vpa( vpa (E, n)

Devuelve la expresión E con con n dígitos decimales

poly2 poly 2sym( sym (A)

Devuelve al arreglo A como un polinomio simbólico en x

Ejemplo 34.

Considere las funciones f = 2a3 − 5, g = a2 + 2 , h = a3 − 5a2 + 6a 6 a, calcular P = f + g + g + h h >> syms a >> f=2*a^3-5; f=2*a^3-5; >> g=a^2+2; g=a^2+2; >> h=a^3-5a^2 h=a^3-5a^2+6a; +6a; >> h=a^3-5*a^ h=a^3-5*a^2+6*a 2+6*a; ; >> w=f+g+ w=f+g+h h w = 3*a^3-3-4*a^2+6*a >> pretty(w) pretty(w) 3 3 a


2 - 3 - 4 a + 6 a



35

Ejemplo 35. 4 Simplificar Simplificar la expresión expresión 11−−xx2

>> >> >> >>

syms x E=simplify((1-x^4)/(1-x)); E=simplify((1-x^4)/(1-x)); E=simplify((1-x^4)/(1-x)); E=simplify((1-x^4)/(1-x)); pretty(E) pretty(E) 3 2 x + x + x + 1

Ejemplo 36.

>> P=po P=poly ly2s 2sym ym([ ([5 5 1 4 6 7]) 7]) p=5*x^4+x^3+4*x^2+6*x+7 pretty(P) 4 5 x

3 + x

2 + 4 x

+ 6 x + 7

Ejemplo 37.

Sea la matriz M = [x ed; ex d] >> syms x d >> M=[ M=[ x exp( exp(d) d); ; exp( exp(x) x) d]; d]; >> D=det(M) D=det(M) D = x*d-exp(d)*exp(x) >> pretty(D) pretty(D) x d - exp( exp(d) d) exp( exp(x) x)

Ejemplo 38.

Sea la función f ( f (x) = 13 − x3 , calcule f (2) f (2) >> f=’13-x^3’ f=’13-x^3’; ; >>p= subs(f,2) subs(f,2) p= 5




36

Función

Descripción

subs( subs(f, a)

Evalúa la función f en en el valor a

subs( subs(f,a,b) f,a,b)

Sustituye en la función f el el valor de a por el de b

compose( compose(f, g )

Función compuesta de f y g

finverse( finverse(f ) f )

Calcula la función inv inversa ersa

limit( limit(Sn,inf )

Calcula el limite de la sucesión S cuando cuando n tiende a ∞

limit( limit(f, a)

Calcula el limite de la función f cuando x tiende al valor de a

limit( limit(f, a)

Calcula el limite de la función f cuando x tiende al valor de a

diff (f, x) o diff (f ) f ) Calcula la derivada derivada de la función f respecto respecto a x diff (f, n)

Calcula la n ésima derivada de la función f respecto respecto a x

taylor( taylor(f,n,x) f,n,x)

Calcula el desarrollo de la serie de Maclaurin de orden n − 1

para la función f en en la variable x taylor( taylor(f,n,x,a) f,n,x,a)

Calcula el desarrollo de la serie de Taylor de orden n − 1

para la función f en en la variable x alrededor del punto a Ejemplo 39.

Sea la función f ( f (x) = 13 − x3 , calcule f (2) f (2) >> syms x >> f=13 f=13-x -x^3 ^3 ; >> p=subs(f,2 p=subs(f,2) ) p = 5

Ejemplo 40.

>> f=@(x)(x+6 f=@(x)(x+6) ) f = @(x)(x+6) >> f(2) f(2) ans = 8




37

Ejemplo 41.

Considere Considere las funciones f ( f (x) = x + x + 4 y g( g (x) = x 3 + 1 Calcule f ( f (g (x)) >> clear clear >> syms x >> f=x+4; f=x+4; >> g=x^3+1; g=x^3+1; >> h=co h=comp mpos ose( e(f f , g) h = x^3+5 >> h=co h=comp mpos ose( e(g g , f) h = (x+4)^3+1 >> pretty(h) pretty(h) 3 (x + 4) + 1 >> simplify(h simplify(h) ) ans = x^3+12*x^2+48*x+65

Ejemplo 42.

Considere Considere las funciones f ( f (x) =

2x3 +x2 +4 x3 +6

, calcule l´ımn→∞ f ( f (x)

syms syms x f=(2*x^3+x^2+4)/( f=(2*x^3+x^2+4)/( x^3+6); limit(f limit(f ,inf) ans = 2

Ejemplo 43. 2 Considere la función f ( f (x) = x 4 − cos(x cos(x) + e + e−x , calcular f  (x)

>> syms x >> f=x^4-cos(x)+exp(-x^2); f=x^4-cos(x)+exp(-x^2); >> df=diff(f) df=diff(f) df = 4*x^3+sin(x)-2*x*exp(-x^2)




38

Ejemplo 44. 2 Considere la función f ( f (x) = sen(xy sen(xy)) + e−x + y 3 , calcular

d2 f dx2

df dx

,

df dy

,

d2 f dxdy

,

>> syms x y >> f=sin(x*y)+exp(-x^2)+y^3; f=sin(x*y)+exp(-x^2)+y^3; >> dxf=diff(f dxf=diff(f,x) ,x) dxf = cos(x*y)*y-2*x*exp(-x^2)

>> dxyf=diff(diff(f,x),y) dxyf=diff(diff(f,x),y) dxyf dxyf = -sin(x*y)*x*y+cos(x*y)

>> dxxf=diff(diff(f,x),x) dxxf=diff(diff(f,x),x) dxxf dxxf = -sin(x*y)*y^2-2*exp(-x^2)+4*x^ -sin(x*y)*y^2-2*exp(-x^2)+4*x^2*exp(-x^2) 2*exp(-x^2)

Ejemplo 45.

Considere la función f ( f (x) = sen(x sen(x).Calcular la serie de Maclaurin de orden 5 >> syms x >> t5=taylor( t5=taylor(sin(x sin(x)) )) t5 = x-1/6*x^3+1/120*x^5

Ejemplo 46.

Considere la función f ( cos(x), calcular la serie de Taylor de orden f (x) = cos(x 3, alrededor del punto x = 2 >> syms x >> t3=taylor(cos(x),4,x,2) t3=taylor(cos(x),4,x,2) t3 = cos(2)-sin(2)*(x-2)-1/2*cos(2) cos(2)-sin(2)*(x-2)-1/2*cos(2)*(x-2)^2+1/6*si *(x-2)^2+1/6*sin(2)*(x-2)^3 n(2)*(x-2)^3




39

Función

Descripción

int( int(f, x)

Calcula la integral indefinida de de f respecto de la variable x

int( int(int( int(f, x), y )

Calcula la integral indefinida de de f respecto de x luego respecto de y

int( int(int( int(f, y ), x)

Calcula la integral indefinida de de f respecto de y luego respecto de x

int( int(

· · · int( int(int( int(f, x), y ) · · · , z )

int( int(f,x,a,b) f,x,a,b)

Calcula la integral indefinida respecto de x, luego de y ,

· · · hasta la variable z

Calcula la integral definida respecto de x de a hasta b

int( int(int( int(f, x), a , b) b), y , c , d) d)

Calcula la integral definida doble para a

≤ x ≤ b y c ≤ y ≤ d

Ejemplo 47.

Calcular la integral indefinida de la función f ( cos(x), f (x) = cos(x >> >> >> >> If

clear clear syms x f=cos(x); f=cos(x); If=int(f,x If=int(f,x) ) = sin( sin(x) x)

Ejemplo 48.

Calcular la integral indefinida de la función f ( f (x) = sen(x sen(x), >> >> >> >> Ih

clear clear syms x h=sen(x); h=sen(x); Ih=int(h,x Ih=int(h,x) ) = -cos -cos(x (x) )




40

Ejemplo 49.

Calcular la integral indefinida de la función f ( f (x) = tan(x tan(x) − x2 , >> clear clear >> syms x >> h=tan(x)-x h=tan(x)-x^3; ^3; Ih=int(h,x) Ih = -log(cos(x -log(cos(x))-1/ ))-1/4*x^4 4*x^4

Ejemplo 50.

Calcular >> >> >> >> If

 

(xy2

− tan( tan(x))dydx ))dydx,

clear clear syms x y f=x*y^3-ta f=x*y^3-tan(x); n(x); If=int(int(f,y),x) If=int(int(f,y),x) = 1/8*x^2*y^ 1/8*x^2*y^4+y*l 4+y*log(co og(cos(x)) s(x))

Función dsolve( dsolve( E  )

Descripción Resuelve la ecuación diferencial E respecto respecto de x, por defecto

dsolve( dsolve( E  , C  )

Resuelve la ecuación diferencial E con con valores iniciales C

dsolve( dsolve( E  , C 1 , C 2 , . . . , Cn )

Resuelve la ecuación diferencial E con con valores iniciales C 1 , . . . Cn

laplace( laplace(f ) f )

Calcula la transformada de Laplace de la función f

ilaplace( ilaplace(f ) f )

Calcula la transformada inv inversa ersa de Laplace de la función f




41

Ejemplo 51.

Calcular y  − by = by = 0,con valores iniciales y(2) = 1 >> clear clear >> syms y b >> y=dsolve(’Dy=b*y’,’y(2)=1’) y=dsolve(’Dy=b*y’,’y(2)=1’) y = 1/exp(2*b 1/exp(2*b)*exp )*exp(b*t) (b*t)

Ejemplo 52.

Calcular (y )2 + y 2 = 1,con valores iniciales y (0) = 0 >> clear clear >>y=dsolve(’Dy^2+y^2=1’,’y(0)=0’) y = -sin(t) sin(t)

Ejemplo 53.

Calcular y  + 5y 5 y − 6 = 0 >> clear clear >> y=dsolve(’D2y+5*Dy-6=0’) y=dsolve(’D2y+5*Dy-6=0’) y = -1/5*exp( -1/5*exp(-5*t) -5*t)*C1+6 *C1+6/5*t+ /5*t+C2 C2

Ejemplo 54.

Calcular f ( f (x) = e x >> clear clear >> syms x >> y=laplace( y=laplace(exp(x exp(x)) )) y = 1/(s 1/(s-1 -1) )




42

Ejemplo 55.

Calcular la transformada inversa de Laplace de f ( f (s) = s2 −3s−6 >> clear clear >> syms s >> y=ilaplace(3/(s^2-s-6)) y=ilaplace(3/(s^2-s-6)) y = 6/5*exp(1 6/5*exp(1/2*t) /2*t)*sinh *sinh(5/2* (5/2*t) t) >> pretty(y) pretty(y) 6/5 exp(1/2 exp(1/2 t) sinh(5 sinh(5/2 /2 t)

Ejemplo 56.

Calcular la transformada inversa de Laplace de f ( f (s) = s−1 1 >> clear clear >> syms s >> ilaplace(1 ilaplace(1/(s-1 /(s-1)) )) y = exp( exp(t) t)



2

Aproximación Aproximación y errores

Introducción P. Henrici da una definición aproximada del Análisis numérico como, la teoría de los métodos constructivos en Análisis matemático, haciendo un especial énfasis en la palabra constructivos, durante mucho tiempo, las matemáticas matemáticas fueron totalment totalmentee constructiv constructivas, as, pues su único objetiv ob jetivoo era llegar a la solución de problemas concretos. No obstante, a medida que los problemas sujetos a la investigación matemática crecían en alcance y generalidad, los matemáticos fueron interesándose, cada vez más, por cuestiones como la existencia, unicidad y propiedades cualitativas de las solución, antes que por su construcción. Una de las causas que condujeron a esta situación fue la escasa capacidad de cálculo que hacía inútil el diseño de algoritmos constructivos de la solución de problemas complejos. No obstante cuando surgieron los primeros primeros ordenadores, ordenadores, se impulso impulso nuevament nuevamentee el diseño de nuevos algoritmos numéricos.

El análisis numérico es una disciplina que contempla el desarrollo y evolución de métodos para calcular, a partir de ciertos datos numéricos, los resultados requeridos, los ingredientes del esenciales en un problema de Análisis numérico pueden resumirse en el siguiente diagrama: Como ejemplo podemos pensar en el problema del calculo del sen(x sen(x), pero con bastante bastante frecuencia frecuencia nos encontramo encontramoss con distintos distintos algoritmos algoritmos para construir la información de salida que se requiere, así volviendo al ejemplo anterior para aproximar el sen(x sen(x) podemos usar el algoritmo que se obtiene 3 5 de hacer un desarrollo de Taylor a la función sen(x sen(x) ≈ x − x3! + x5!

43 Mg. Marlo Carranza Purca Purca

2.1. Tipos de errores

44

x=linspace(-6,6,1000); y=x-x.^3/6+x.^5/120; yy=sin(x); plot(x,y,’r’); hold on ; plot(x,yy);legend(’aproximaci plot(x,yy);legend(’aproximación ón por Taylor’,’función Taylor’,’función sin’); axis([-5,5,-6,6])

Por lo tanto para escoger entre los diversos algoritmos disponibles deben estudiarse los aspectos teóricos que contribuirán a la elección del algoritmo mas adecuado a cada caso concreto. En general, los criterios fundamentales para preferir un criterio frente a otro son la rapidez y la precisión o equivalentemente el error

2.1. 2.1. Tipo Tiposs de er erro rore ress

Los errores numéricos surgen del uso de aproximaciones para representar operaciones operaciones y cantidades cantidades matemáticas exactas. 1. Error de truncamiento o discretización . Los errores de truncamiento o discretización provienen, por ejemplo, de la sustitución de una expresión continua por otra discreta (por ejemplo al aproximar la derivada de f por una expresión en diferencias), f  (x0 )

f (x + h + h)) − f ( f (x ) ≈ f ( h 0

0

Usando el desarrollo de Taylor de f : f  (x0 )h2 f  (x0)h3  f ( f (x0 + h + h)) = f ( f (x0 ) + f + f (x0 )h + + + ··· 2!


3!



45

de donde f ( f (x0 + h + h)) h

2

f (x ) f  (x )h f  (x )h − f (  = f (x ) + + + ··· 2! 3! 0

0

0

0

Así el error cometido es f ( f (x0 + h + h)) h

2

f (x ) f  (x )h f  (x )h − f (  − f (x ) = 2! + 3! + · · · = O( O (h) 0

0

0

0

2. Error de redondeo . Estos se originan debido a que la computadora emplea un numero determinado determinado de cifras significativ significativas as durante durante los cálculos. cálculos. Los numeros irracionales como π

14159265358979323846... ≈ 3, 14159265358979323846...

o e

≈ 2, 71828182845904523536028747135266249775724709369995 71828182845904523536028747135266249775724709369995...... ......

, entre muchos, no pueden representarse con un número infinito de cifras significativas, por lo tanto, no pueden ser representados de manera exacta en la computadora, además como la computadora usa la base 2 , no pueden representar exactamente algunos números en base 10. Esta discrepancia por la omisión de cifras significativas se llama error de redondeo. Para ambos tipos de errores, la relación entre el resultado exacto o verdadero y el aproximado está dado por Valor verdadero = Valor aproximado + Error de aquí se tiene (2.1)

Errort = valor = valor_verdadero

− valor_aproximado

Una desventaja en esta definición es que no se toma en consideración el orden de la magnitud del valor que se estima . Por ejemplo, un error de un centímetro es mucho mas significativo si se esta midiendo un remache o si esta una cirugía. Una manera de tomar en cuenta las magnitudes de las cantidades que se evalúan consiste en normalizar el error respecto al valor verdadero, verdadero, es decir Errorrelativofraccionalverdadero =

error_verdadero valor _verdadero

también se se puede expresar así (2.2)

εt =

error_verdadero 100% valor _verdadero

Es el error relativo relativo porcentual porcentual verdadero. verdadero. Mg. Marlo Carranza Purca Purca



46

Ejemplo 57 (Cálculo de errores).

Suponga que se tiene que medir la longitud de un puente y la de un remache, y se obtiene 9999 y 9 cm, respectivamente si lo valores verdaderos son 10000 y 10 cm. 1. Calcule el error verdadero 2. EL error relativo porcentual verdadero en cada caso. Resolución. 1. El error en la medición del puente es E t = 10000

− 9999 = 1 cm.

y en el remache es E t = 10

− 9 = 1 cm.

2. El error relativo porcentual en el puente es εt =

1 100 100 % 10000

εt =

1 100 100 % 10

= 0, 0,01 %

y para el remache es = 10 10 %

Por lo tanto, aunque ambas medidas tienen un error de 1 cm. el error relativo porcentual, es mucho mayor . se concluye entonces que se ha hecho un buen trabajo en la medición del puente, mientras en la medición del remache se hizo un mal trabajo.

El subíndice subíndice t de Errort y de εt significa que el error a sido normalizado al error verdadero; pero en situaciones reales no se conoce el valor exacto , entonces en dichos casos, una alternativa el normalizar el error usando la mejor aproximación posible al valor verdadero, es decir, para la misma aproximación (2.3)

εa =

error_aproximado 100% valor _aproximado

donde el εa significa que el error a sido normalizado a un valor aproximado.Uno de los retos que enfrentan los métodos numéricos es la estimación del error en ausencia del conocimiento de los valores verdaderos. Por ejemplo ciertos métodos numéricos usan un método iterativo para calcular los resultados, en tales casos se hace una aproximación considerando la aproximación anterior. Este proceso se efectúa varias veces o de forma iterativ iterativaa esperando cada vez una mejores aproximacion aproximaciones. es. (2.4)

εa =

aproximacion_actual aproximacion_anterior 100% aproximacion_actual

−

A menudo no nos importa el sigo del error, mas bien nos importa que su valor absoluto sea menor que una tolerancia prefijada εs , en tales casos, los Mg. Marlo Carranza Purca Purca



47

cálculos se repiten hasta que (2.5)

|ε | < ε a

s

Si se cumple lo anterior entonces se considera que el resultado obtenido está dentro dentro del nivel nivel aceptable aceptable fijado previamente previamente ε s

2.1. 2.1.1. 1.

Núme Número ro de cifr cifras as sign signifi ifica cati tiv vas

Las cifras significativas de un número son la primera no nula y todas las siguientes. Así pues, 2,350 tiene cuatro cifras significativas mientras que 0,00023 tiene sólo dos. Es conveniente relacionar los errores con el número de cifras significativas en las aproximaciones. Teorema 58 (Scarborough, 1966).

Se tendrá la seguridad de que el resultado es correcto en al menos n cifras significativas. (2.6)

εs = (0, (0,5

2 n

× 10 − ) %

Ejemplo 59 (Estimación del error con métodos iterativos).

En matemáticas con frecuencia las funciones se representan mediante series infinitas, por ejemplo, la función exponencial se calcula usando (2.7)

x

e

≈

x 2 x3 1 + x + x + + + + 2! 3!

···

x n + n!

Así cuanto más términos se le agreguen a la serie, la aproximación será cada vez más una mejor estimación del valor de ex . Empezando con el primer termino de ex = 1 y agregando término por término estime el valor de e0,5 . 1. Calcule los errores relativo, porcentual verdadero y normalizado a un valor aproximado usando las ecuaciones 2.2 y 2.4 respectivamente. Observe que el valor verdadero de e0,5 = 1,648721 . . .. 2. Agregue términos hasta que el valor absoluto del error aproximado ε a sea menor que un criterio criterio de error preestablecid preestablecidoo ε s con c on tres cifras significativas, use la ecuación 2.6

Resolución.



2.2. Representación de números en la computadora

48

La ecuación 2.6 se emplea para determinar el criterio de error que asegura un resultado sea correcto en al menos tres cifras significativas: εs = (0, (0,5

2 3

× 10 − ) % = 0,05 %

por lo tanto, se agregaran tantos términos a la serie hasta que εa sea menor que este valor. 1. La primera estimación es ex = 1 2. La segunda estimación es e x = 1 + x + x, y para e0,5 = 1 + 0, 0,5 = 1,5 De aquí se tiene un error relativo porcentual según la ecuación 2.2 εt =

1,64721 1,5 100 100 % = 9,02 % 1,648721

−

La ecuación 2.4 , se utiliza para determinar una estimación aproximada del error por 1,5 1 100 100 % = 33, 33,3 % 1,5 Com el εa no es menor que el valor requerido εs se debe seguir εa =

−

agregando términos y continuar los cálculos. tenemos la siguiente tabla. Términ Términos os Result Resultado adoss

εt ( % )

εa ( % )

1

1

39.3

2

1.5

9.02

33.3

3

1.625

1.44

7.69

4

1.645833333 0.175

1.27

5

1.648437500 0. 0.0172

0.158

6

1.648697917 0.0 0.00142 0.0 0.0158

Así después de usar seis términos, el error aproximado es menor que εs = 0,05 % y el cálculo termina, aunque esto es cierto la mayoría de las veces.

2.2. 2.2.

Repres Repr esen enta taci ción ón de númer úmeros os en la computadora

Numéricamente los errores de redondeo se relacionan de manera directa con la forma en que se guardan los números en la memoria de la computadora.C computadora.Casi asi siempre siempre la representa representación ción y la aritmética aritmética computacional computacional son satisfactori satisfactorias as y pasan inadvertidas. inadvertidas. La unidad fundamental fundamental mediante mediante la Mg. Marlo Carranza Purca Purca



49

cual se representa la información se llama palabra, esta consiste en una cadena de cadena de caracteres de dígitos binarios ( bits ), solo puede tomar dos valores 0 o 1 , encendido o apagado; positivo o negativo o cualquier otra dicotomía electrónicamente viable, es por eso que se usa la representación binaria.

2.2.1. 2.2.1.

Sistem Sistemas as numé uméric ricos os

Un sistema sistema numérico numérico es una convención convención para represent representar ar cantidades cantidades,, Debido a que tenemos diez dedos en las manos es que usamos el sistema de base 10, que utiliza 10 1 0 dígitos, ( 0, 0 , 1, 2, 3, 4, 4, 6, 7, 8, 9 ) para representar números. Ejemplo 60. 25 = 20 + 5 = 2

× 10 + 5 345 = 300 + 40 40 + 5 = 3 × 10

3

+4

× 10 + 5 86439 86439 = 80000+6000+400+30+9 80000+6000+400+30+9 = 8×10 +6 ×10 +4 ×10 +3 ×10+9 547, 547,154 = 5 × 10 + 4 × 10 + 7 × 10 + 1 × 10− + 5 × 10− + 4 × 10− 4

2

1

0

3

1

2

2

3

En el sistema binario 101, 101,101 = 1 22 +0 22 +1 20 +1 20 +1 2−1 +0 2−2 +1 2−3 = 5,625 en base 10.

×

×

×

×

×

×

×

Representación Entera

Veamos como los enteros se representan en la computadora, el método más sencillo se denomina método de magnitud con signo y emplea el primer bit de una palabra para indicar el signo, con 0 para positivo y 1 para negativo los bit sobrantes se usan para guardar el número Ejemplo 61.

La representación de un entero decimal −173 en una computadora de 16 bits usando el método de de magnitud con signo es : 1 0 0 0 0 0 0 0 1 0 1 0 1 1 0 1




50

Ejemplo 62.

determine el rango de enteros de base 10 que pueda representarse en una computadora de 16 bits. Resolución. De los 16 bits, se tiene el primer bit para el signo. Los 15 bits restantes pueden contener los números de binarios de 0 a 111111111111111 , el límite superior se convierte en un entero decimal. (1

14

×2

) + (1

13

×2

)+

3

1

· · · + (1 × 2 ) + 1 × 2

+ 1 = 32767 = 215

−1

Así en una computadora de 16 bits una palabra puede guardar en memoria un entero decimal del rango de −32767 a 32767 pero el cero esta definido como 0000000000000000, sería redundante usar el número 1000000000000000 para definirlo por lo tanto se usa para definir otro entero negativo adicional y el rango va de −32768 a 32767.

En las computadoras convencionales se prefiere usar la técnica del complement complemento o de 2 que incorpora el signo dentro de la magnitud

del número, en lugar de emplear un bit adicional para representar el más o menos.

Notación decimal en coma flotante

nos quedamos con esta ultima notación, denominada notación decimal en coma flotante normalizada, caracterizada por que la fracción es un número comprendido entre 1 y 10 y escribimos simplemente −6,2345 + 2, en general será : 1

623,45 = −62, 62,345 × 10 = −6,2345 × 10 −623,

2

10 , E ∈ ±m ± E ; 1 ≤ m < 10, ∈ N ∪ {0}

También se emplea el siguiente criterio 1 = 0,029411765 . . . se guarda en base 10 con un punto flotante que 34 únicamente únicamente puede guardar 4 lugares decimales entonces se guardaría como 1 = 0,0294 × 100 tenemos un cero inútil a la derecha del punto decimal, el 34 número número puede normalizarse normalizarse así 1 c on b : base. = 0,2941 × 10−2 en este caso 1b ≤ m < 1 con 34 Note que este proceso introduce un error de redondeo. Observación 2.1. R se

representa mediante un conjunto finito de números racionales que se representan como números en punto flotante o números máquina. Un número ±m · bc, donde b−1 < m < 1 o 1 < m < b. en punto flotante tiene la forma ± En un número de punto flotante podemos distinguir: 1. Signo




51

2. Base 3. Mantisa o parte fraccionaria Ejemplo 63.

Con el criterio anterior determine un conjunto hipotético de números con punto flotante para una máquina que usa palabras de 7 bits, el primer bit lo usa para el signo del número, los siguientes tres para el signo y la magnitud del exponente y los últimos tres para la magnitud de la mantisa. mantisa. Veamos Veamos que, el número mas pequeño posible p osible se representa representa en la figura 0 1 1 1 1 0 0

El cero inicial señala que la cantidad es positiva. El 1 en la segunda casilla indica que el exponente es tiene signo negativo. Los unos en el tercero y cuarto lugar dan un valor máximo en el exponente. 1 × 21 + 1 × 2 = 3 , por lo tanto el exponente es −3. Por lo tanto la mantisa esta especificada por el 100 10 0 en los tres últimos lugares lo cual nos da 1 × 2−1 + 0 × 2−2 + 0 × 2−3 = 0,5 En este sistema hipotético, el número mas pequeño en base 10 es +0, +0,5× 2−3 = 0,0625, esto deja en evidencia

1. El rango de cantidades que puede representarse es limitado, como el caso de los enteros, hay números grandes positivos y negativos que no pueden representarse. 2. Existe sólo un número finito de cantidades que pueden representarse dentro de un rango. Así el grado de precisión es limitado. Es evidente que los números no pueden representarse de manera exacta.

2.2.2. 2.2.2.

Aritmé Aritmétic tica a del del punto punto flotan flotante te

Por simplicidad, en lo que sigue, supongamos números máquina decimales:

±0.d d

1 2

. . . dk

n

× 10 ,

con 1 ≤ d1 ≤ 9; 0 ≤ di ≤ 9, i = 2; . . . ; k

Se puede normalizar cualquier número real positivo y para convertirlo en y = 0.d1d2 . . . dk dk+1dk+2

× 10

n

Si y está en el rango de la máquina puede ponerse en forma de punto flotante, f l(y), tomando k dígitos decimales. Hay dos maneras de elegir esos k dígitos decimales: corte y redondeo. Mg. Marlo Carranza Purca Purca



52

1. En el corte simplemente se eliminan todos los dígitos a partir del dk+1 . 2. En el redondeo redondeo se agrega agrega 5 × 10n−(k+1) a y y luego se realiza el corte. Es decir,

a ) si dk+1 ≥ 5 agregamos 1 a d k y cortamos, b ) si dk+1 < 5 < 5 simplemente cortamos. Ejemplo 64.

El número π es un número irracional, luego tiene una expresión decimal infinita de la forma π = π = 3,14159265 . . . . En forma decimal normalizada es π = 0,314159265 . . .

× 10

El número π en punto flotante con cinco dígitos y corte se representa por f l(π ) = 0,31415

Con redondeo será

× 10

1

f l(π) = (0, (0,31415 + 0, 0,00001)

= 3,1415

× 10

1

= 3,1416

El error cometido al reemplazar un número por su forma en punto flotante recibe el nombre de error de redondeo independientemente de si se ha aplicado el método de corte o de redondeo.

2.2.3. 2.2.3.

Normal Normaliza izació ción n de los númer números os de punto punto flotan flotante te

1. Para los casos donde se emplea el corte. (2.8)

|x − f l(x)| ≤ E |x|

2. Para los casos donde se emplea el redondeo. (2.9)

|x − f l(x)| ≤ E |x| 2

donde E se se denomina épsilon de la maquina, el cual se calcula como (2.10)

1 t

E = b −

donde b número base. t es el número de términos significativos en la mantisa.



2.3. Ejemplos

53

Ejemplo 65.

Épsilon la maquina Calcule el épsilon de la máquina para el ejemplo 63. El sistema de punto flotante hipotético anterior, emplea valores de base 2, y el número de bits para la mantisa t = 3 por lo tanto el épsilon de la máquina debe ser E = 21−3 = 0,25

2.2.4. 2.2.4.

Operaci Operacione oness en en pun punto to flotan flotante te

Representam Representamos os las operacion op eraciones es en la computadora computadora mediante: ⊕, , ⊗, . Supondremos Supondremos una aritmética de dígitos dígitos finitos (equivale (equivale a realizar realizar operacio op eraciones nes exactas sobre representaciones de punto flotante, y luego convertir el resultado exacto en su representación de punto flotante), viene definida por:

2.3 2.3.

Ejem jemplos los

1. Encontrar la expresión decimal de los números binarios

a ) 11, 11,11 b ) 111. 111.101 2. Encontrar la expresión binaria de los números decimales

a ) 0,1 b ) 5.3 3. Efectue 0,1557 × 101 + 0, 0,4381 × 10−1 . El número de la mantisa de menor exponente se modifica de tal forma que los exponentes sean los mismos, alineamos los puntos decimales. 0,1557

1

× 10 0,004381 × 10

1

se obiene 0,160081 × 101 y cortando queda 0,1600 × 101 4. Efectue 0,3641 × 102 − 0,2686 × 102 . 0,3641

2

× 10 −0,2686 × 10

2

se obtiene 0,0955 × 102 el cual se puede expresar así 0, 0 ,9550 × 101 Mg. Marlo Carranza Purca Purca


2.3. Ejemplos

54

5. Efectue 0,7642 × 103 − 0,7641 × 103 . 0,7642

3

× 10 −0,7641 × 10

3

se obtiene 0,0001 × 103 el cual se puede expresar así 0, 0 ,1000 × 100 . Así, en este caso, se agregan tres ceros no significativos, lo cual introduce un error sustancial de calculo debido a que las manipulaciones siguientes actúan como si los ceros fueran significativos. 6. Efectúe 0,1363 × 103 × 0,6423 × 10−1 . Aquí los exponentes se suman y las mantisas se multiplican . 3

1

2

× 10 × 0,6423 × 10− = 0,08754549 × 10 y cortando resulta 0,8754 × 10 0,1363

= 0,87545490

× 10

1

1

7. Un número grande de cálculos interdependientes . Investigar el efecto de error de redondeo en un gran número de cálculos interdependie interdependientes ntes,, desarrolle desarrolle un programa programa que sume un número número 100000 veces. Sume el número 1 con simple precisión y 0, 0 ,00001 con precisiones simple y doble . x1=1;suma1=0;suma2=0; x2=0.00001; for i=1:100000 i=1:100000 , suma1=sum suma1=suma1+x1 a1+x1; ; suma2=su suma2=suma2+x ma2+x2; 2; end suma suma1 1 = 100000 suma suma2 2 = 1.0000 suma suma2 2 = 1.0000 form format at long long e >> suma2 suma2 suma suma2 2 = 9.999999999980838e-001



2.3. Ejemplos

55

8. suma de un número grande y uno pequeño. Vamos a sumar un número pequeño 0,0010 con un número grande 4000 usando una computadora hipotética con una mantisa de 4 dígitos yun expoente de 1 dígito. Normalizamos los números, tenemos 0,4000

4

× 10 0,0000001 × 10

4

obtenemos 0,40000001 × 104 el cual se corta a 0, 0 ,4000 × 104 así resulato como i no hubieramos sumado nada. Este tipo de error uede ocurrir cuando de suman series infinitas, por ejemplo si el término inicial de la série es relativamente grande en comparación con los demás términos, luego que se han sumado unos pocos terminos, estamos en la situación de sumar una cantidad pequeña con una grande. Una manera de reducir este tipo de errores consiste en sumar la serie en sentido inverso, es decir en orden ascendente en lugar del orden descendente. 9. Cancelación por resta Resolución de la ecuación x2 + 62, 62,10x 10x + 1 = 0 , Las raíces aproximadas son

f l(x1 ) = 0, 1611


× 10−1 Manual Manual de métodos numérico numéricos s U.C.H. U.C.H.

2.3. Ejemplos

56

y el error relativo es ahora: 0,0161072 + 0, 0,1611 0,0161072

1

× 10− ≈ 2 × 10−

4

En el caso de que b fuera negativo las cosas ocurrirían exactamente al contrario. Para x1 obtendríamos una buena aproximación y para x 2 habría que considerar la segunda opción. Problemas

1. La división por un número muy pequeño (o la multiplicación por un número muy grande) da un error absoluto muy grande. 2. La sustracción sustracción de números números casi iguales da errores relativos relativos muy grandes. 3. La propagación de estos errores al resto de cálculos. Podemos evitar estos problemas: problemas: 1. Minimizando el número de operaciones, 2. Ordenando adecuadamente las operaciones, 3. Replanteando el problema en otros términos. Conclusiones

1. Los cálculos numéricos pueden ser inexactos. 2. El error final de un proceso de cálculo efectivo es fruto de la acumulación de distintos tipos de errores. Unos iniciales (de entrada, redondeo de los datos, truncadura del problema, ...) y otros generados a lo largo del proceso. 3. El aumento de la precisión (aumento de las unidades de memoria para almacenar los números), en general reduce el error final. No obstante, a veces conviene sacrificar precisión frente a economía de tiempo o de recursos. (Un equilibrio entre precisión suficiente y coste adecuado se le suele llamar eficiencia). 4. Hay operaciones o procesos de cálculo que propagan fuertemente los errores errores de redondeo, redondeo, operaciones operaciones o algoritmos algoritmos inestables (inestabilidad (inestabilidad numérica). 5. Hay problemas que tienen una naturaleza que los hace especialmente sensibles a la variación de los datos, problemas mal condicionados. Para cierto tipo de problemas se puede definir una medida de esta sensibilidad sensibilidad mediante mediante un número número de condición. condición.



3

Resolución de ecuaciones

Introducción En este capítulo, vamos a ver como podemos resolver ecuaciones con una incógnita, pero desde el punto de vista de los métodos numéricos, recordemos que desde la época escolar aprendimos resolver la ecuación = 0 f ( f (x) = ax 2 + bx + bx + c c =

(3.1) usando la fórmula

√ b ± b − 4ac − x = 2

(3.2)

2a

donde los valores calculados con la fórmula (3.2) son llamados raíces de la ecuación (3.1) , que son los valores que hacen que a la ecuación igual a cero es decir verifican la ecuación. Aunque la formula para resolver la ecuación cuadrática es útil para resolver la ecuación (3.1), existen muchas otras ecuaciones donde las raíces no se pueden calcular tan fácilmente como por ejemplo Ejemplo 66.

1. Ecuaciones de la curiosidad de los alumnos e−x

− 2x = 0,

xx = 2

2. Ecuación de Kepler x

− e sen x − b = 0

donde x: Posición de la tierra, e : Excentricidad de la orbita


3. Resolución de ecuaciones

58

Ejemplo 67.

1. Modelos de población c1 ex +

c2 x(ex

− 1) − c

3

=0

donde x:Índice de natalidad anual, c1 : Población actual, c2 : Inmigrantes en un año, c3 : Población final. 2. Raíces de polinomios x4 + 2, 2,6x3

− 3x

2

+ πx

− 0,065 = 0

Podemos demostrar que tienen raíces reales, podemos decir cuantas son, pero no existe un método para calcularlas calcularlas en forma exacta ( si se pudieran calcular en forma exacta).

Este tipo de argumentos lleva a plantearse el estudio de los métodos necesariamente iterativos, que generan una sucesión x 0 , x1 , . . . xn buscando (3.3)

l´ım xn = x,

n

→∞

tal que P ( P (x) = 0

es decir para así calcular las raíces de una ecuación P ( P (x) = 0 en el intervalo [a; b], con el grado de aproximación deseado. Definición 5 (Raíz ). f : D

f (x) = 0 en D si f ( f (r) = 0 ⊆ R → R, r ∈ D es raíz de la ecuación f (

Aspectos 1. Convergencia . Bajo que condiciones se cumple (3.3 ) 2. Orden de convergencia . Velocidad

|x − r| = k, k , →∞ |x − p|

l´ım

n

n+1 n

α

con 0 < k < 1, 1 , α > 0

Convergencia de orden α y k constante asintótica. convergenciaa es lineal. lineal. a ) Si α = 1 entonces la convergenci b ) Si α = 2entonces la convergencia es cuadrática.



3.1. Méto do de la Bisección

59

convergenciaa es superlineal. superlineal. c ) Si 1 < 1 < α < 2 entonces la convergenci 3. Análisis del error , Cotas. |k | = |xr − r| ≤ C 4. Criterios de parada

a ) Si |f ( f (xn )| ≤ tol1 . b ) Si |xn − xn−1| ≤ tol2 . c ) Si |xn −|xxnn| 1| ≤ tol3 es el número de iteraciones, d ) n < N con N es −

En lo que sigue de este documento consideraremos que las raíces de la ecuación f ( f (x) = 0 son aisladas, es decir, para cada raíz existe un intervalo que no contiene más raíces de la ecuación, teniendo en claro que no todas las ecuaciones ecuaciones tienen raíces aisladas. aisladas. Ejemplo 68.

La gráfica de la función f ( f (x) =

    x sen

1 x

0

si x = 0 si x = 0

muestra que 0 no es una raíz aislada.

3.1. 3.1. Mé Métod todo o de de la la Bis Bisec ecci ción ón Teniendo en cuenta las consideraciones hechas anteriormente el método de la bisección bisección consiste consiste en tomar un interv intervalo alo [a0 ; b0 ] ( de longitud b 0 − a0 ) donde se encuentre una raíz, luego debemos partir este intervalo en dos y ver en cual de estos dos subintervalos se encuentra pero al hacer esto el intervalo donde Mg. Marlo Carranza Purca Purca



60

se encuentra la raíz tendrá la mitad de longitud que el intervalo inicial ( b0 −2 a0 ) y si repetimos este proceso sucesivamente n veces, la raíz que buscamos estará en un intervalo de longitud b02−na0 , es decir conseguimos conseguimos asegurar asegurar que esta raíz se encuentra en un intervalo bastante pequeño y así tenemos una aproximación para la raíz buscada, un resultado de vital importancia para asegurar que en un intervalo hay una raíz es el teorema del Valor Intermedio en su forma particular llamada el teorema de Bolzano, veamos Teorema 69 (Teorema de Bolzano).

Sea f : [a, b] → R y continua sobre [a, b] y f ( f (a).f ( .f (b) < 0 entonces existe r ∈ a, b tal que f ( f (r) = 0 Notación 3.1. La función f : [a, b]

3.1. 3.1.1. 1.

→ R, continua sobre [a, b] se indicara como f ∈ ∈ C ([a, ([a, b])

Algo Algori ritm tmo o de la bis bisec ecci ción ón

Algoritmo 3.1: Bisección entrada: a, b, f, tol, n max salida : r solución aproximada a ; b0 = b = b;; f a = f = f ((a); f b = f = f ((b); 1 a0 = a; max do 2 for i = 0 : n xi = a+2bi ; 3 f i = f ( f (xi ); 4 if f i < tol then tol then 5 r = x = x i ; break 6

−

−

| |

7 8 9 10

if f ( f (ai )f i < 0 < 0 then ai+1 = a i ; bi+1 = x i ; else ai+1 = x i ; bi+1 = b i ;

Ejemplo 70.

Aplicar Aplicar el método de la bisección x x = 4 en el intervalo [1; [1; 4] y mostrar las cinco primeras iteraciones Resolución.

Sea la función f ( f (x) = x x − 4 a la cual vamos a aplicar el método de la bisección, para lo cual usamos la siguiente tabla




n an

xn =

61

an +bn 2

bn

f ( f (an )

f ( f (xn)

f ( f (bn)

252

0 1

2.5

4

-3

5.8821

1 1

1.75

2.5

-3

-1.3373 5.8821

2 1.75

2.125

2.5

-1.3373 0.9618

3 1.75

1.9375

2.125

-1.3373 -0 -0.3981 0. 0.9618

4 1.9375

2.0313

2.125

-0.3981 0.2187

5 1.9375

1.9844

2.0313 -0.3981 -0.1040 0.2187

5.8821 0.9618

Cuadro 3.1: tenemos que la solución aproximada luego de 5 iteraciones es 2,0156 Teorema 71.

Sea f : [a, b] → R y continua sobre [a, b] y f ( f (a).f ( .f (b) < 0 entonces existe r ∈ a, b tal que f ( f (r) = 0, se prueba que aplicando el algoritmo de la bisección se consigue una solución aproximada an +2 bn con un error  tal que

 −

an + b + bn 2

r

 ≤

b a < 2n+1

−

Prueba Sea f C ([a, ([a, b])tal que f ( f (a).f ( .f (b) < 0 < 0 entonces por el teorema de Bolzano, existe al menos un r Suponiendo que r es la única única a, b tal que f ( f (r ) = 0. Suponiendo

∈ ∈

∈ 

raíz de (3.4)

f ( f (x) = 0 en [a, b]

dividimos el intervalo [a, b] por la mitad. Así se puede considerar dos casos 1. f ( entonces r = f ( a+2 b ) = 0, entonces

a+b 2

2. f (  0, en este caso se elige el intervalo [a, a+2 b ] ó [ a+2 b , b] en cuyos f ( a+2 b ) = extremos la función f toma signos opuestos, para de esta forma , poder seguir aplicando el teorema de Bolzano (observe que, como r es la única raíz de f en [a, b], solo uno de estos intervalos tendrá la propiedad de que f cambie de signo en los extremos del mismo.) denotando a este intervalo por [a [ a1 , b1 ], lo dividimos por la mitad y procedemos de igual forma. Así, reiterando este procedimiento la ecuación o bien obtenemos la raíz exacta de la ecuación ( 3.4 ) o bien una sucesión de intervalos cerrados {f [ f [an , bn ]}n∞=1 encajados es decir (3.5) a

≤ a ≤ a ≤ · · · ≤ a ≤ a ≤ b ≤ b ≤ · · · ≤ b ≤ b, n ∈ N 1

2


n

n+1

n+1

n

1



62

en los cuales (3.6)

f ( f (an).f ( .f (bn )

≤ 0, n ∈ N

y

− a = b 2− a , n ∈ N la propiedad ( 3.5 ) determina que {a }∞ es una sucesión monótona, creciente y acotada y que {b }∞ es una sucesión monótona, (3.7)

bn

n

n

n n=1

n n=1

decreciente y acotada . Por lo tanto existen l 1 y l2 tales que l1 = l´ım an n

→+∞

≤ →l´ım∞ b n

+

n

= l 2

ahora bien de la ecuación 3.7 de deduce que n

l´ım (bn an ) = l´ım

→+∞

−

n

→+∞

− b

a

2n

1 = 0 entonces l 1 = l 2 = l = l n→+∞ 2n

= (b a) l´ım ım

−

por lo que de la continuidad de f y de la relación (ec:ec7 ) tenemos f ( f (l).f ( .f (l) f ( f (l)2

0 0

≤ ≤

de donde f ( = l es la única raíz de la f (l) = 0 , así de esta forma r = l ecuación (3.4) . Además

 − r

an + b + bn 2

 ≤

bn

−a ≤ b−a 2 2 ≤ b2− a , n ∈ N n

n+1

n+1

Por lo tanto, para encontrar un valor que aproxime a la raíz r con un error inferior a  > 0 basta tomar r = an +2 bn , n ∈ N tal que

 − r

es decir

n>

an + b + bn 2 ln(b ln(b

 ≤

b a < 2n+1

−

ln() − a) − ln( −1 ln 2

Observación 3.1.

El método de la bisección es útil para dar una idea rápida de la localización de las raíces, es decir, para determinar intervalos de longitud pequeña que contengan una única raíz, pero los cálculos aumentan considerablemente si queremos una buen aproximación de la raíz al ser la convergencia lenta, esto hace que este método se aplique, principalmente como un paso previo a la utilización de otros métodos iterativos de convergencia mas rápida. Mg. Marlo Carranza Purca Purca



3.1. 3.1.2. 2.

63

Ejer Ejercic cicio ioss

1. Si f tienen un único cero en [ −2, 5], ¿Cuántas ¿Cuántas iteraciones iteraciones de bisección bisección se deben hacer para aproximar este cero con error absoluto ≤ 0,5 × 10−4 ? 3(x−1) 2. Resuelva e3(x ln(x − 1)2 + 1 = 0 con al menos cinco decimales − ln( exactos.

3. Resuelva e 3x − ln( ln(x2 + 1) − 30 = 0 con al menos cinco decimales exactos. 4. (x − 1)2 = 0 tiene claramente una raíz en [0, [0, 2]. ¿Podemos usar bisección para aproximarla? 5. x = 2 es un cero del polinomio P ( P (x) = −1536 + 6272x 6272x − 11328x 11328x2 + 11872x 11872x3 − 7952x 7952x4 + 3528x 3528x5 − 1036x 1036x6 + 194x 194x7 − 21x 21x8 + x9 . Aproxime esta raíz, con un intervalo adecuado. 6. La ultima moda entre los aficionados al billar es la mesa circular. Para los principiantes, el juego consiste simplemente en golpear la bola Q con la bola P después después de un impacto impacto I en en la banda. Los parámetros del problema pueden verse en la figura, la mesa tiene radio R, la posición de las bolas P y Q queda determinada por las coordenadas cartesianas (xP ; yP ) y (x ( xQ ; yQ ), y el punto de impacto I viene definido por el ángulo θ . Mediante consideraciones geométricas sencillas (que se dejan como

Figura 3.1: Mesa de billar ejercicio) puede verse que los valores de θ que proporcionan los puntos de impacto I son son las raíces de la ecuación (3.8) f ( f (θ) =

√

xp sen θ yp cos θ

(R cos θ xp

−

−

)2 +(R +(R sen θ

−yp

)2

+

√

xQ sen θ yQ cos θ

−

(R cos θ xQ )2 +(R +(R sen θ yQ )2

−

−

=0

Nótese que R es la única incógnita de la ecuación 3.8. Los valores de R, xP , y P , xQ y yQ son datos del problema. Para resolver la ecuación 3.8 use el el método de bisección, que construye una sucesión θ k de aproximaciones a un cero de la función f .Resuelver la ecuación 3.8 tomando los valores R = R = 1, x P =.6, y P = 0, x Q =-0.6, y Q = 0. Un punto de impacto I viene viene dado entonces por 2θ , x0 = 1,5



3.2. Méto do del punto fijo

64

3.2. 3.2. Mé Métod todo o del del pun punto to fijo fijo Introducción y definiciones Definición 6 (Punto fijo).

Sea f : [a, b] → R una función. Un elemento r ∈ [a, b] es un punto fijo de f si f ( f (r ) = r

Ejemplo 72.

Jugando con una calculadora de bolsillo, uno puede verificar que aplicando repetidamente la tecla coseno al valor real x0 = 1, se consigue la siguiente sucesión de números reales x1

= cos(x cos(x0 )

= 0,54030230586814

x2

= cos(x cos(x1 )

= 0,85755321584639

.. .

...

...

x10

= cos(x cos(x9 )

= 0,74423735490056

.. .

...

...

x20 = cos(x cos(x19 ) = 0,73918439977149

que tendería al valor α = 0,73908513 . . . . Puesto que, por construcción, xk+1 = cos(x cos(xk ) para k = 0, 1, . . . con x0 = 1 el límite α satisface la ecuación cos(α cos(α) = α . Por esta razón α se llama punto fijo de la función coseno. Podemos preguntarnos cuantas de tales funciones de iteración podrán ser útiles para calcular los ceros de una función dada. En el ejemplo anterior, α no es solo un punto fijo de la función coseno sino también un cero de la función φ(x) = x − cos(x cos(x), por tanto el método propuesto puede considerarse como un método para calcular los ceros de f . Por otra parte, no toda función tiene puntos fijos. Por ejemplo, repitiendo repitiendo el experimento experimento anterior con la función exponencial y = e = e x y x0 = 1.




65

Figura 3.2: La función φ(x) = x − cos(x cos(x) admite uno y solo un punto fijo, x mientras que la función φ(x) = e no tiene ningún punto fijo

3.2. 3.2.1. 1.

Algor Algorit itmo mo de iter iterac ació ión n de de pun punto to fijo fijo

Aclaremos la idea intuitiva anterior considerando el siguiente problema. Dada una función φ : [a, b] → R , hallar α ∈ [a, b] tal que α = φ = φ((α). Si existe un tal α se llama punto fijo de φ y podría calcularse mediante el siguiente algoritmo (3.9)

φ(xk ) = x k+1

donde x0 es una conjetura inicial. Este algoritmo se llama de iteraciones de punto fijo y φ se dice que es la función de iteración de punto fijo. El ejemplo introductorio es, por tanto, un ejemplo de iteraciones de punto fijo con φ(x) = cos(x cos(x). Una interpretación geométrica de 2 se proporciona en la Figura siguiente. Se puede conjeturar que si φ es una función continua y el límite de la sucesión xk existe, entonces tal límite es un punto fijo de φ.

Figura 3.3: Representación de unas cuantas iteraciones de punto fijo para cos x




3.2. 3.2.2. 2.

66

Algo Algori ritm tmo o

Vamos encontrar una solución de f ( f ( p) p) = p , dada una aproximación inicial po Algoritmo 3.2: Punto fijo entrada: f , p0 , tol, n − max Solución aproximada aproximada o mensaje mensaje de fracaso salida : Solución 1 i = 1; while i n max do max do p = f = f (( p0 ); 3 f i = f ( f (xi ); 4 if p p p 0 tol then tol then 5 imprimir la solución p ; 6

≤ −

2

| − − |≤

7 8 9 10 11

12

el procedimient procedimientoo se realizo realizo con éxito; parar

else i = i = i + 1 ; p0 = p ;

el método fracaso después de n − max iteraciones

3.2. 3.2.3. 3.

Exis Existe tenc ncia ia y uni unici cida dad d

Definición 7 (Contracción).

Una función f : [a, b] 0 ≤ k < 1 tal que (3.10)

→ R, es una contracción en [a, b] si existe

|f ( f (x) − f ( f (y)| ≤ k|x − y |, x,y ∈ [a, b]

Ejemplo 73.

la función f : [1, [1, 2] → R dada por f ( f (x) = x2 + x1 es una contracción en [1, [1, 2] de constante k = 21 pues para todo x, y ∈ [1, [1, 2] se tiene

|f ( f (x) − f ( f (y)|

  −



 −

x 1 y 1 x y = ( + ) + ) = 2 x 2 y 2 1 1 = x y 2 xy 1 x y 2

−

 | − |

x

−y − xy



≤ | − |




67

Teorema 74 (Punto fijo de Banach).

Si f : [a, b] verifica que

[0, 1 y → R es continua y contractiva de constante k ∈ [0,

(3.11)

f ([ f ([a, a, b])

⊂ [a, b]

entonces f tiene un único punto fijo en [a, [ a, b] es decir un único r ∈ [a, b] tal que f ( f (r ) = r . Además r es el limite de una sucesión definida por

 

(3.12)

x0

∈ [a, b]

arbitrario ,

xn = f ( f (xn−1 )

n

∈ N.

y se tiene la siguiente estimación para el error kn

(3.13)

|x − r| ≤ 1 − k |x − x |, n ∈ N n

1

0

Ejemplo 75.

Vamos a aproximar las raíces de la ecuación e −x − x = 0 con un orden de error inferior a 10−3 . Para ello consideremos la función F ( F (x) = e −x

− x, x ∈ R

    ⊂   1 ,1 3

Buscamos un intervalo, como por ejemplo cumpla la hipótesis, con f ( f (x) = e −x f

Veamos se tiene f

1 ,1 3

, de tal manera que se

1 ,1 3



1 1 1 1 = 1 < 1 y f (1) f (1) = > 3 e 3 e3

entonces, por ser f estrictamente decreciente, se verifica que f

  ⊂   1 ,1 3

1 ,1 3

por otra parte al aplicar el teorema del valor medio, se tiene que x

y

n

1 3

|e− − e− | = e− |x − y| ≤ e− |x − y|, x,y Mg. Marlo Carranza Purca Purca

  ∈ 1 ,1 3


3.2. Méto do del punto fijo ya que

68

1 1 < n < 1 entonces e −n < e− 3 3

  

1 por lo tanto la función f es contractiva en 13 , 1 de constante k = e − 3 < 1 así por el teorema 74 se tiene que existe un único r ∈ 13 , 1 tal que r = e = e −r si y solo si e−r − r = 0 y además, si consideramos la sucesión

 

x0 =

1 3

xn = e −xn

1

−

n

∈ N.

se verifica l´ımn→∞ xn = r . Para aproximar r con el orden deseado, tomamos n ∈ N verificando 3.2.3, es decir n e− 3 − 13 − 1 < 10 −3 si y solo si n > 21, 21 ,6276 1 e − 1

−e

3



3



por lo tanto tomando n = 22 que x22 = 0,56714233424982 aproxima a la raíz r con un error inferior a 10−3 .Puesto que la acotación del error es simplemente eso, una acotación es probable que para valores mas pequeños también se tenga la aproximación con la precisión requerida. Vamos a ver otra versión del teorema presentado que tal vez nos ayude a aligerar algunos calculos. Teorema 76.

Si f : [a, b]

→ R continua tal que

(3.14)

f ([ f ([a, a, b])

⊂ [a, b]

entonces existe un punto punto fijo para para f ( f (x) , además supongamos que  existe f en a, b tal que (3.15)

|f (x)| ≤ k < 1para todo x ∈ a, b

entonces f tiene tiene un único punto fijo en [a, [a, b] es decir un único r ∈ [a, b] tal que f ( f (r ) = r . Además r es el limite de una sucesión definida por (3.16)

 

x0

∈ [a, b]

arbitrario ,

xn = f ( f (xn−1 )

n

∈ N.

y se tiene la siguiente estimación para el error (3.17)


kn

|x − r| ≤ 1 − k |x − x |, n ∈ N n

1

0



69

Demostración. Existencia

De la relación ( 3.14 ) se tiene f ([ f ([a, a, b]) ⊂ [a, b], es decir que: Si f ( f (a) = a ó f ( f (b) = b , se tiene la existencia del punto fijo. Supongamos que no entonces f ( f (a) > a y f ( f (b) < b , definimos h(x) = g( g (x) − x se tiene que h es continua en [a, b] y evaluando se tiene h(a) = f ( f (a) a > 0 y h(b) = f ( f (b) b < 0 se tiene que h(a) h(b) < 0 < 0 el teorema Bolzano implica que existe r a, b tal que h(r) = 0 por lo tanto h( h (r) = f ( f (r) r = 0 es decir f ( f (r) = 0, es decir existe el punto fijo de f . Unicidad Supongamos que existen dos puntos fijos p y q diferentes en [a, b], por el teorema del valor Medio existe un número α p, q [a, b] tal que

−

×

−

∈ 

−

∈ ⊂ | p p − q | = |f ( f ( p) p) − f ( f (q )| = |f  (α)|| p p − q | ≤ k | p p − q | | p p − q | ≤ k| p p − q | 0 ≤ (1 − k)| p p − q | ≤ 0 0 ≤ | p p − q | ≤ 0 p = q

Convergencia de la sucesión xn Por el teorema del valor medio y α

∈ a, b se tiene f (x − ) − f ( f (r)| ≤ |f  (α)||x − − r | ≤ k|x − − r| |x − r| = |f ( n

n 1

n 1

n 1

aplicando inductivamente este resultado y recordando que l´ımn→∞ kn = 0 se tiene

|x l´ım |x →∞ 0 ≤ l´ım |x →∞ l´ım |x →∞ n n n

n n n n

n

− r| ≤ k |x − r| ım k |x − r| − r| ≤ l´→∞ − r| ≤ 0 − r| = 0 0

n

n

0

Estimación para el error Recordemos Recordemos que xn r kn x0

|x − x | m

n

| − | ≤ | − r| −x | = |x − x − + x + x − − x − + x + x − − · · · + x ≤ |x − x − | + |x − − x − | + |x − − x − | + · · · + |x − x | ≤ k − |x − x | + k − |x − x | + k − |x − x | + · · · + k |x − x | ≤ k (1 + k + k + + k k + · · · + k − − )|x − x | ≤ l´ı→∞ m k (1 + k + k + + k k + · · · + k − − )|x − x | ≤ k (1 + k + k + + k k + k + · · · )|x − x | ≤ 1 k− k |x − x | m

m 1

m m 1

m 1

1

n

l´ım xm

m

→∞

| −x | |r − x | |r − x | n n n


m 1

0 2

n

m 1 m 2

m 2

m 2

m 2

m n

2

3

m 1 m 3

1 0 m n 1

1 m n 1

2

1

n+1

m 2

1

0

n

n+1 n 1

0

1

0

0

n

1

0


n

0


3.2. 3.2.4. 4.

70

Ejer Ejercic cicio io

La ecuación de factor de intensidad de esfuerzos para una placa de ancho w y espesor t con una grieta en el borde de largo a es √ (3.18) K = σY σ Y a Donde Y es un factor geométrico que depende del ancho de la placa y el tamaño de grieta, siendo (3.19)

Y = 1,99

−

   −  

a a 0,41 + 18, 18,7 w w

a 38, 38,48 w

2

3



a + 53, 53,85 w

4

La falla catastrófica de la placa se produce cuando el factor de intensidad de esfuerzos K iguala iguala o supera a la tenacidad a la fractura K IC IC entonces el tamaño de grieta crítica es (3.20)

af =

  K IC IC σY

2

Como el factor geométrico Y depende de a f , la ecuación ( 3.20 ) debe resolverse por el método de punto fijo. La función de iteración es entonces (3.21) ak+1 =

 

σ 1,99

K IC IC

− 0,41

 ak w

+ 18, 18,70

 −   ak 2 w

38, 38,48

ak 3 w

+ 53, 53,85

   

2

ak 4 w

Se tiene un caso de una placa sujeta a tension donde w = 2,5, σ = 2,5, K IC IC = 52 . Se elige como aproximación inicial a0 = 0,25, se muestra la solución con tres cifras decimales af = 0,620



3.3. Méto do de Newton

71

3.3. 3.3. Mé Métod todo o de Ne Newt wton on Introducción El método de Newton o de Newton Raphson es uno de los métodos numéricos más conocidos y poderosos para la resolución de un problema de búsqueda de raíces de la ecuación f ( f (x) = 0. Si el valor inicial de la raíz es x j entonces se puede trazar una tangente desde el punto (x j , f ( f (x j )) de la curva. Por lo general el punto por donde esta tangente cruza al eje x representa una mejor aproximación de la raíz .El método de Newton se deduce de esta interpretación geométrica (se puede deducir también a partir de la serie de Taylor ) De la figura 5.1 se observa que

Figura 3.4: Representación gráfica del método de Newton f  (x j ) =

f ( f (x j ) 0 x j x j+1 j +1

−

de donde se tiene x j

−x

−x

=

j+1 j +1 =

−

f ( f (x j ) f  (x j )

de aquí j+1 j +1

f (x ) −x + f f (  (x ) j

j

j

y así tenemos x j+1 j +1 = x j

f (x ) − f f (  (x ) j

j

La cuál se conoce como el método de Newton.




72

Teorema 77 (Método de Newton).

Dada una función f : [a, b] → R continua de derivada no nula en a, b, se tiene la sucesión llamada el método de Newton.

  3.3. 3.3.1. 1.

x0

∈ [a, b]

xn+1 = x n

−

f ( f (xn ) f (xn ) 

dado n



∈ N {0}

Algo Algori ritm tmo o

Nuestra intención es encontrar una solución de f ( f (x) = 0, dada una aproximación inicial p0 Algoritmo 3.3: Newton entrada: f , p0 , tol, n − max Solución aproximada aproximada o mensaje mensaje de fracaso salida : Solución 1 i = 1; while i n max do max do f ( f ( p0 ) p = p = p 0 f ( p0 ) ; 3 if p p p0 tol then tol then 4 imprimir p ; 5

≤ − − | − |≤

2



6 7 8 9 10

11

el procedimient procedimientoo se a realizado realizado satisfactori satisfactoriament amente; e; break

else i = i = i + 1 ; p0 = p ;

El método fracaso luego de n − max de iteraciones iteraciones Ejemplo 78.

Utilizar el método de Newton para calcular la raíz de f ( f (x) = e −x − x, empleando como valor inicial x0 = 0 , empleando cuatro iteraciones, indicando al final el error porcentual relativo, verifique que se obtiene la siguiente tabla de valores.




73

n xn

E r %

0 0

100

1 0.50 0.5000 0000 0000 0000 11.8 11.8 2 0.56 0.5663 6311 1100 0033 0.14 0.1477 3 0.56 0.5671 7143 4316 1655 0.00 0.0000 0022 2200 4 0.56 0.5671 7143 4329 2900 0.00 0.0000 0000 0001 01 Cuadro 3.2: tenemos que la raíz aproximada f ( f (x) = e −x − x luego de 4 iteraciones es 0,567143290 Ejemplo 79.

Utilizar el método de Newton para calcular la raíz de f ( f (x) = cos(x cos(x) − x π empleando como valor inicial x0 = 4 , empleando cuatro iteraciones, indicando al final el error porcentual relativo, verifique que se obtiene la siguiente tabla de resultados.

n

xn

0 0.78 0.7853 5398 9816 1635 35 1 0.73 0.7395 9536 3613 1337 37 2 0.73 0.7390 9085 8517 1781 81 3 0.73 0.7390 9085 8513 1332 32 4 0.73 0.7390 9085 8513 1332 32 Cuadro 3.3: tenemos que la raíz aproximada f ( f (x) = cos(x cos(x) − x luego de 4 iteraciones es 0,7390851332

3.3.2. 3.3.2.

Conv Convergenc ergencia ia del del métod método o de Newton Newton

Para garantizar la convergencia del método de Newton vamos a requerir las siguientes hipótesis sobre la función f : [a, b] → R

 ∈ ∈ C  H  f

( )

2

([a, ([a, b])

f ( f (a)f ( f (b) < 0 < 0 f  tiene signo constante en [a, b] f  tiene signo constante en [a, b]




74

Teorema 80.

Si f verifica ( H) y c ∈ [a, b] es el extremo de [a, [ a, b] tal que sign [f ( f (c)] = sign[ sign[f  (c)] entonces el método de Newton para f con x0 = c = c converge al menos cuadráticamente a una única raíz r de f en [a, [ a, b].

Ejemplo 81.

Aproximemos las raíces de la función f ( f (x) = x5 + 5x + 8, x ∈ R, mediante el método de Newton. Veremos que se obtiene la siguiente tabla de resultados para las primeras cinco iteraciones tomando como punto inicial a x0 = −2, considere la gráfica de la función f ( f (x) = x5 + 5x 5x + 8

n

xn

n

xn

0 -2.00 -2.00000 000000 000000 00000 0000 3 -1.167 -1.167693 693392 39203 03490 490 1 -1.60 -1.60000 000000 000000 00000 0000 4 -1.167 -1.167036 036664 66447 47529 529 2 -1.32 -1.32236 236390 390595 59521 2133 5 -1.167 -1.167036 036183 18370 70190 190 Cuadro 3.4: tenemos que la raíz aproximada f ( f (x) = cos(x cos(x) − x luego de 4 iteraciones es 0,7390851332

3.3.3. 3.3.3.

Varian ariantes tes del método método de Newton Newton

Como hemos visto en las condiciones de aplicación del método de Newton, éste tiene convergencia al menos cuadrática. No obstante el inconveniente que puede presentarse es que en las sucesivas iteraciones hay que evaluar la




75

derivada de f y esta puede ser difícil de evaluar, haciendo que el gasto en cada iteración sea grande. Vamos a considerar algunas variantes de este método que sustituyen la derivada de f por algo menos costoso de calcular, así de esta forma las iteraciones serán mas sencillas de calcular, aunque la convergencia sea más lenta. Teorema 82 (Método de Whittaker).

 0, En este método se toma, en lugar de f  (x) a un valor constante λ = así se tiene el método de Whittaker

 

x0

∈ [a, b]

xn = x n−1

−

dado

f (xn λ

1)

−

n

∈N

Teorema 83.

Si f verifica verifica (H) y c ∈ [a, b] es el extremo de [a, [a, b] tal que sign [f ( f (c)] =  sign[ sign[f (c)] y λ ∈ R\{0} verifica que sign[ sign[λ] = sign[ sign[f  ] y λ

| | ≥ |f (c)|

entonces el método de Whittaker para f con x0 ∈ [a, b] arbitrario converge, al menos linealmente a una única raíz r de f en [a, b].

Ejemplo 84.

Sea f : [0, [0, 1] → R dada por f ( f (x) = e x − 2cos(x 2cos(x), x ∈ [0, [0, 1], compruebe que se puede tomar los valores λ = λ = 5 y x 0 = 1 y se obtiene la siguiente tabla

n

xn

n

xn

0 1.00 1.0000 0000 0000 0000 0000 0000 00 3 0.56 0.5630 3060 6088 8817 1710 1033 33 1 0.67 0.6724 2464 6455 5566 6655 5545 45 4 0.55 0.5501 0103 0311 1180 8073 7329 29 2 0.59 0.5935 3568 6817 1715 1542 4295 95 5 0.54 0.5444 4405 0501 0143 4325 2519 19 Cuadro 3.5: tenemos que la raíz aproximada f ( f (x) = cos(x cos(x) − x luego de 4 iteraciones es 0,7390851332




76

siendo la raíz 0, 0 ,539785160080928, vea la gráfica de la función f ( f (x) = e x

2cos(x) − 2cos(x

Teorema 85 (Método de las cuerdas ).

En este método se toma, en lugar de f  (x) a la función f ( f (b) b

− f ( f (x) − x , x ∈ [a, b]

así se tiene el método de las cuerdas

 

x0

∈ [a, b]

xn = x = x n−1

−

dado

f ( f (xn 1 ) (b f ( f (b) f (xn 1 ) −

−

−

−x − ) n∈ N n 1

Teorema 86 (Método de la secante).

Este método es una version del método de Newton en el que se discretiza la derivada mediante la fórmula de derivación numérica más simple, la de dos puntos, es decir en lugar de f  (x) tomamos f ( f (xn ) xn

para n

f (x − ) − f ( −x − n 1

n 1

≥ 1, se obtiene el método de la secante

 

xn+1


∈ [a, b] = x = x −

x0 , x1

n

f (xn ) f ( f (xn ) f ( f (xn

−

1)

−

(xn

dados

−x − ) n∈N n 1



77

Teorema 87 (Método de la falsa posición).

Este método, sin dejar de ser una variante variante del método de Newton, puede interpretarse como una generalización del método de la bisección, ∈ C ([a, para aplicarlo basta que la función f ∈ ([a, b]) con f ( f (a).f ( .f (b) < 0 de forma que existe un único r ∈ a, b tal que f ( f (r) = 0. La ecuación que une a los puntos (a, f ( f (a)) y (b, f ( f (b)) es y

f ( f (b) − f ( f (a) − f ( f (a) = (x − a) b−a

y la intersecci intersección ón con el eje de las abscisas abscisas se da cuando y = 0 , tenemos x = a = a

f ( f (a) − f ( (b − a) = c f (b) − f ( f (a)

Como hacíamos con el método de la bisección, tomamos

 

a1 = a, = a, b1 = c = c

si f ( f (a).f ( .f (c) < 0 < 0

a1 = c, = c, b1 = b

si f ( f (c).f ( .f (b) < 0 < 0

la siguiente iteración se obtendría aplicando la misma estrategia al intervalo [a1 , b2 ] y así sucesivamente.

3.3. 3.3.4. 4.

Ejer Ejercic cicio ioss

1. Se desea resolver la ecuación 3x3 − 2x + 8 = 0 usando el método de Newton, con tres cifras significativas, se tiene f ( f (x) = 3x3 − 2x + 8 = 0 y f (x) = 9x2 − 2 = 0 se obtiene la sucesión 

xk + xk + 1 = xk

−

3x3 2x + 8 tomando x0 x 0 = 9x2 2

−

−

−10

muestre que la séptima iteración se obtiene −1,546 2. Se desea resolver la ecuación sen(x sen(x) + 3 cos( cos(x) = 0 en el intervalo 0 ≤ x ≤ 10 usando el método de Newton, en el intervalo 0; 3, se obtiene 1,893, verifique que otras soluciones son x = 5,034 y x = 8,176



Sistemas Sistemas de ecuaciones lineales

4 Introducción

Vamos a estudiar diversos métodos de resolución de sistemas de ecuaciones lineales.Son muchas las aplicaciones donde ( determinación de tensiones en los nodos de una red de corriente continua, calculo de estructuras reticulares definidas por vigas, vigas, modelos de económicos económicos,.,. . . ) cuya resoluci resolución ón práctica práctica involucra plantear y resolver sistemas de ecuaciones lineales: mas aun la totalidad de los problemas que surgen de en las ciencias aplicadas se resuelven mediante algoritmos que involucran, en alguna de sus etapas la resolución de un sistema lineal. Desde el punto de vista teórico los sistemas lineales no plantean ninguna dificultad, sin embargo en la práctica, los sistemas lineales involucrados en la resolución de los problemas a matrices de orden grande como 10000 o 100000 ; para estos tamaños la regla de Cramer es inviable por el elevado número de operaciones que exige, esta es la razón por la que se han diseñado otros métodos más eficientes para la resolución de sistemas sistemas de ecuaciones. ecuaciones. Dada una matriz A ∈ Mn y b ∈ V , el problema que se plantea es calcular u ∈ V tal que Au = Au = b b .

4.1. 4.1. Mé Métod todos os dire direct ctos os 4.1.1. 4.1.1.

Sistem Sistemas as diag diagona onales les y trian triangul gulare aress

n Sea la matriz A = (ai,j )i,j=1 i,j =1 ∈ Mn es triangular e inversible entonces n

det(A det(A) =



aii = a 11a22 . . . ann = 0



i=1

Lo que implica que aii  = 0 para i = 1, 2, . . . , n y por lo tanto podemos dividir por los elementos de la diagonales de la matriz A 78 Mg. Marlo Carranza Purca Purca

4.1. Méto dos directos

79

Vamos a resolver resolver el sistema lineal lineal de ecuaciones ecuaciones

Ejemplo 4.1.

  

(4.1)

2012x 2012x + 5y + 4z z = 2042 2y + 3z =

19

2z =

10

De la ultima ecuación z = 5, reemplazando en la segunda ecuación, se tiene que y = 2, ahora reemplazando en la primera ecuación tenemos x = 1, así la solución es (1, (1, 2, 5) y tenemos el cs. = cs. = {(1, (1, 2, 5)}. Es decir

  

z =

10 2

y = 21 (19 x =

1 2012

− 3z ) (2042 − (4z (4z − 5y ))

En general para resolver el sistema lineal de ecuaciones

      

a1,1 x1 + a1,2 x2 +

(4.2)

a2,2 x2 +

···+ ···+ ...

a1,n−1 xn +

a1,n xn =

b1

a1,n−1 xn +

a2,n xn =

b2

···

···

an−1,n−1 xn−1 + an−1,n xn = bn−1 an,n xn =

bn

basta considerar xn =

bn an,n

xi =

1 ai,i

 −

n j= j =i+1 ai,j x j

bi



, , i = n = n

− 1, n − 2, . . . 1

Esta técnica se conoce con el nombre de método de remonte. En el caso de que la matriz A sea triangular inferior e inversible, aplicando un remonte hacia abajo, se obtiene que la solución del sistema

(4.3)

es

    

a1,1 x1 a2,1 x1

+a2,2 x2

... an−1,1 x1 +an−1,2 x2 + an,1 n,1 x1

 

+a1,2 x2 +

x1 =

b1 a1,1

xi =

1 ai,i


··· ···

 − bi

b1

=

b2

···

···

an−1,n−1xn−1

=

bn−1

=

bn

an−1,n−1 xn−1 + an,nxn

i 1 j=1 j =1 ai,j x j

−

=



, , i = 2, 3, . . . n


4.1. Méto dos directos

4.1.2. 4.1.2.

80

Elimin Eliminaci ación ón Gaussi Gaussiana ana

Antes de dar la descripción teórica, veamos un ejemplo. Ejemplo 4.2.

Vamos a aplicar aplicar el método de Gauss para para la resolución resolución del

sistema lineal

  

x1 +

x2 +

=

x3

   

6

x1 + 2x2 + 3x3 = 13 x1 +

x2

x3

−

=

2

0x1 +

x2 +

2x3

+x4

= 1

x1 +

2x2 +

x3

+3x +3x4

= 0

x1 +

x2 +

−x

+x4

= 5

0x1 +

x2 +

8x3

3

+12x +12x4 = 2

Ax = b b, luego de efectuar este sistema se puede escribir en forma matricial Ax = 75 t los cálculos se tiene la solución (x1 , x2 , x3 , x4 ) = ( 2 , − 463 , 676 , −6)t Ejemplo 4.3.

 

2−26 x + y + y = 1 x + y + y = 2

la solución exacta del sistema viene dada por

 

x = 2

− y = 1,00000001490116 25

−

2 y = 11− −2

= 0,999999998509884

26

−

trabajando con precisión simple tenemos x = 1, y = 1. El sistema se puede escribir

 

2−26 1 1

1

        x

1

=

y

2

a) Tomando 2 −26 como primer pivot tenemos

 

26

2−

0

x

1 1

26

−2

el sistema es aproximadamente

 

2−26 0

      y

26

−2

2

26

−2

       −  x

1

1

=

 

1

=

226

y

por lo tanto de aquí la solución es (x, y ) = (0, (0, 1) Mg. Marlo Carranza Purca Purca


4.2. Méto dos iterativos

81

b) Tomando 1 como primer pivot tenemos

 

1

1

2−26 1

el sistema es aproximadamente

        x

2

=

y

1

          1 1

x

0 1

y

2

=

1

con lo cual el sistema tiene como solución a (x, ( x, y ) = (1, (1, 1) Este ejemplo muestra que los errores de redondeo con efecto desastroso provienen de la división por pivotes muy pequeños, en la practica se usa una de las estrategias siguientes al comienzo de la k − sima etapa de eliminación

1

≤ k ≤ n − 1

1. Estrategia del pivoteo parcial Se toma como pivote el elemento de mayor módulo de entre los n − k + 1 últimos elementos de la columna k − sima , es decir k k aik = máx ax a pk k p n

≤≤

| |

2. Estrategia del pivoteo total Se toma como pivote el elemento de mayor módulo de la submatriz Ak correspondiente, es decir se elige el k elemento aik , k ≤ i, j ≤ n de modo que k k = máx ax a p,q aik k p,q n

≤ ≤

| |

En general en el método de Gauss usaremos habitualmente el pivoteo parcial.

4.2. 4.2. Mé Métod todos os iter iterat ativ ivos os Introducción En esta guía vamos a estudiar dos de los métodos iterativos mas clásicos para la resolución de un sistema lineal Ax = Ax = b b , con A invertible, estos métodos comparten la misma idea básica en su construcción, el descomponer o partir la matriz del sistema en la suma de dos matrices. Es decir, expresamos la matriz A como A = M = M

− N −




82

donde la matriz M sea sea fácil de invertir, de aquí se verifica que (M

− −

Ax N )x Mx x

= = = =

b b N x + b + b − M 1 N x + M + M −1 b

      c

B

(4.4)

x = Bx + Bx + c c

De esta forma podemos considerar el método iterativo

   

(4.5)

x0

∈ V, arbitrario = Bx B x + c, k ∈ N ∪ {0}

xk+1

k

En la práctica, para calcular xk+1 se resolverá x0

(4.6)

∈ V, arbitrario

M xk+1 = N xk + b, k

en vez de trabajar directamente con 4.5

∈ N ∪ {0}

Notación 4.1. n 1. Sea una matriz A = (aij )i,j=1 = 0 para i,j =1 ∈ Mn inversible con aii 

i = 1, . . . n

2. Descomponemos la matriz A como A = D = D − E − F donde n n D = diag = diag((a1,1 , . . . , an,n ); E = (ei,j )i,j=1 i,j )i,j=1 i,j =1 , F = (f i,j i,j =1

Matrices triangulares inferir y superior respectivamente ei,j =

 − 

ai,j para i > j

f i,j i,j =

 − 

ai,j para i < j

para i ≤ j. 0 para i ≥ j. a la descomposición de la matriz A la denominaremos descomposición D − E − F por puntos de la matriz A 0

4.2. 4.2.1. 1.

Métod Método o de Jaco Jacobi bi

Consiste en tomar M = D y N = E + F + F tenemos (D

(4.7)

− E −

Ax F ) F )x Dx x

= = = =

b b (E + F + F ))x + b + b D−1 (E + F + F ))x + D + D −1b

que conduce al método iterativo de Jacobi ( por puntos ) Mg. Marlo Carranza Purca Purca



83

Definición 8 (Método de Jacobi).

 

(4.8)

x0

∈ V, arbitrario

xk+1 = D −1 (E + F + F ))xk + D −1 b, k

∈ N ∪ {0}

La matriz de este método es J = D −1 (E + F + F )) = I

1

− D− A

La iteración definida en (4.8) puede escribirse coordenada a coordenada como ai,i xki +1 = bi (4.9) = bi

k i,1 i,1 1 i 1 k j=1 j =1 i,j j

k i,1 i,1 1 xi 1 n k j= j =i+1 ai,j x j

−a x −···−a − Σ − a x − Σ

−

−

−a

k i,i+1 i,i+1 xi+1

−···−a

k i,n xn

donde xk = (xk1 , . . . , xkn ) y xk+1 = (xk1+1 , . . . , xkn+1)

4.2. 4.2.2. 2.

Métod Método o de de Gau Gauss ss Se Seid idel el

Para el calculo de la componente xki +1 se ve ve claro claro que una estrat estrategia egia adecuada para mejorar la convergencia sería utilizar las componentes ya k k calculadas {xk1+1 , . . . , xki−+1 1 } en vez de usar las antiguas {x1 , . . . , xi−1 } esta consideración nos lleva a reemplazar el sistema 4.9 por ai,i xki +1 = bi (4.10) = bi

k+1 i,1 i,1 1 i 1 k+1 j=1 j =1 i,j j

k+1 i,1 i,1 1 xi 1 n k j= j =i+1 ai,j x j

−a x −···−a −Σ− a x −Σ

−

−

k i,i+1 i,i+1 xi+1

− a

k i,n xn

−···−a

para i = 1, . . . , n matrices, estas ecuaciones se escriben Dxk+1 = Ex E xk+1 + F xk + b

es decir (D

k+1

− E )x

= F xk + b

Tenemos definido el nuevo método iterativo tomando M = D

− E,

N = F

de esta forma (D

(4.11)

−

Ax = b E )x = F x + b + b x = (D E )−1F x + (D (D

−

1

− E )− b

que conduce al método iterativo de Gauss Seidel por puntos Mg. Marlo Carranza Purca Purca



84

Definición 9 (Método de Gauss Seidel).

 

(4.12)

x0

∈ V, arbitrario = (D − E )− F x

xk+1

1

k

+ (D (D

1

− E )− b, k ∈ N ∪ {0}

La matriz de este método es G = (D

1

1

− E )− F = I − (D − E )− A

que se denomina matriz de Gauss seidel por puntos.

Contrariamente a lo que sucedía en el método de Jacobi, la n componentes del vector xk+1 deben obtenerse de manera sucesiva a partir de de las componentes ya calculadas de xk+1 y las restantes del vector xk por esta razón a este método se le denomina método de las iteraciones sucesivas, además el método de Gauss Seidel será, en principio, más rápido pues la matriz M contiene contiene más elementos de la matriz A. Ejemplo 4.4.

Consideremos la matriz

  

 −  −   − −    −  −     −  −        2

D =

     

2 0 0 0 3 0

3

α

0

2

E

F siendo

E =

A partir de la definición

J = D −1 (E + F + F )) =

Tambien

G = (D E )−1 F =

−

  


1

0

0 0 2

0

2

Donde α ∈ R, tenemos que A = D = D

  

2

0 0

2 0 0

F =

α 0 0

1 2

0 0

0

1 3

1

2 0 1

0 0

1 2

α 0 0

2 0 0 2 3 0 α 0 2

−1

0

2 0

0 0 0

0 2 0 0 0 1 0 0 0

   −  −   − −

0

=

1 2

0 2 0 0 0 1

  

=

   1 0

2 3

0

α 3

0 0

0 0

1 3

1 3

0

α 4

0

1 2

1 3

   

   0 2 0 0 0 1 0 0 0


  


85

G = (D

1

− E )− F =

  

0

1

0

− −

0

0

2 3

1 3

α 3

0

  

En la tabla se dan los valores de los radios espectrales de las matrices de los métodos de Jacobi y de Gauss-Seidel, para algunos valores de α De esta α

r(J )

r (G)

-1 0.8486 0.8486565 565391 391570 57000 0.8603 0.860379 79610 610028 02806 06 -3 0.9726 0.9726325 325833 833593 59355 1.1150 1.115069 69293 293303 30390 90 -5 1.0826 1.0826484 484563 563912 91255 1.3051 1.305158 58649 649140 14088 88

forma, a partir de estos resultados resultados se concluye concluye que 1. Para α = −1 ambos métodos son convergentes. 2. Para α = −3 el método de Jacobi converge y el método de Gauss Seidel diverge. 3. Para α = −5 ambos métodos son divergentes.

Ejemplo 4.5.

Consideremos la matriz

 −   − −  − −    −  −     2

2 0

2

Donde α ∈

1 0 R, tenemos que A = D = D E

D =

  

2 0 0 0 3 0 0 0 2

A partir de la definición J = D −1 (E + F + F )) =

También G = (D E )−1 F =

−

3

  


  

  

0

E =

α

     

2 F siendo

  −     −    −−    −    −−

0 0

2 0 0

0 2

F =

1 0 0

1 2

0 0

0

1 3

0 0

2 0 0 2 3 0 1 0 2

0

2

0

2 0

α

1 2

1 0

0

0 2 0 0 0 0

0 0 0 0

1

−1

=

1 2

0

=

α

0

0

0

α

0

  

1

0

2 3

0

1 3

−

0

0

0 0

1 3

1 3

0

1 4

0

1 2

α 3

   

   0 2

0

0 0

−

0 0

0


α 3

  


86

G = (D

1

− E )− F =

  

0 0 0

1

0

2 3

α 3

− − − 0 1 3

  

En la tabla se dan los valores de los radios espectrales de las matrices de los métodos de Jacobi y de Gauss-Seidel, para algunos valores de α De esta α

r(J )

r (G)

-1 0.8486 0.8486565 565391 391570 57000 0.4082 0.408248 48290 290463 46386 86 -4 0.0301 0.0301808 808496 496534 53411 0.8164 0.816496 96580 580927 92773 73 -7 1.1750 1.1750238 238131 131738 73833 1.0801 1.080123 23449 449734 73464 64

forma, a partir de estos resultados resultados se concluye concluye que 1. Para α = −1 ambos métodos son convergentes. 2. Para α = −4 el método de Jacobi diverge y el método de Gauss Seidel converge. 3. Para α = −7 ambos métodos son divergentes. Ejemplo 4.6.

Estudiar la convergencia de los métodos de Jacobi y Gauss Seidel por puntos para las matrices

A =

Solución

  

1 2 1 1 2 2

 −   2

y B =

1 1

  −

2

−1

1

2

2

2

−1

2

1

  

1. Para la matriz A se tiene

  − −

0

J =

−2

1

0

2

−2

  −   2

1

0

y G =

  

0

−2

0

2

0

0

  −   2

3

2

de donde se tiene que P J J (λ) = −λ3, P G (λ) = −λ(2 − λ)2 De esta forma se obtiene que sp( sp(J ) = 0 por lo que r(J ) < 1 < 1

{}

Por lo que el método de Jacobi es convergente para la matriz A. Mg. Marlo Carranza Purca Purca



87

Para sp( sp(G) = 0, 2 por lo que r(G) = 2 > 1 > 1

{ }

Por lo que el método de Gauss - Seidel para la matriz A no es convergente. 2. Para la matriz B se verifica que

  − 

0

J =

1

0,5 0

−0,5 −1

0,5 0,5

0

  

y G =

 

  

0

0,5

0

−0,5

−0,5 −0,5

0

0

0,5

 

2

  

5 2 De esta forma P J J√ y P G(λ) = −λ 12 + λ (λ) = −λ λ + 4 Tenemos r(J ) = 25 > 1 , el método de Jacobi para la matriz B no es convergente, mientras que al ser r(G) 12 < 1 , el método de Gauss-Seidel es convergen convergente. te.

4.2.3. 4.2.3.

Implem Implemen entac tación ión

1. Presentamos la implementación del método de Eliminación Gausiana %Resol %Resoluci ución ón de un sistem sistema a de ecuaci ecuacione ones s usando usando %Eliminación Gaussiana %Ingre %Ingreso so de datos datos %sol %solo o se pide pide la matr matriz iz aume aument ntad ada a sera sera llam llamad ada a a clear clear %limpia %limpia variables variables clc clc %lim %limpi pia a pant pantal alla la display(’ a=[ 1 -1 2 -1 -8; 2 -2 3 -3 -20; 1 1 1 0 -2; 1 -1 4 3 4] ’); a=inpu a=input(’ t(’Ing Ingres rese e la matriz matriz aumenta aumentada da a= ’); a [n,mm]=size(a); r=1; for k=1:n-1 k=1:n-1 w=abs(a(k,k)); %busca %buscarem remos os el mayor mayor elemen elemento to de una column columna a para para usarla usarla como como %pivot for j=k:n j=k:n if abs(a(j,k) abs(a(j,k))>w )>w w=abs(a(j w=abs(a(j,k)); ,k));r=j r=j end end %El %El mayo mayor r elem elemen ento to de la k esta esta en la fila fila r %Hacemos %Hacemos el intercambi intercambio o




88

aux=a(k,:);a(k,:)=a(r,:);a(r,:)=aux;

%Aho %Ahora ra usar usarem emos os el pivo pivot t para para hace hacer r cero ceros s deba debajo jo de el for for j=k+ j=k+1: 1:n n m=a(j,k)/a(k,k); m=a(j,k)/a(k,k); a(j,:) a(j,:)=a( =a(j,: j,:)-m )-m*a( *a(k,: k,:) ) %elimi %eliminac nación ión debajo debajo de la diagon diagonal al end end %resolver %resolver el sistema sistema triangula triangular r x(n)=a(n,mm)/a(n,n); %x(n)=b(n) %x(n)=b(n)/a(n, /a(n,n);%g n);%guarda uarda el n-esimo n-esimo valor de x suma=0; for k=n:-1:1 k=n:-1:1 suma=a(k,mm); for j=k+1:1: j=k+1:1:n n suma=suma-a(k,j)*x(j); end x(k)=suma/a(k,k); end % Impres Impresion ion de los result resultado ados s disp(’El disp(’El resultado resultado del sistema sistema es:’); es:’); for i=1:n i=1:n fprintf(’x(%i)=%8.5f\n’,i, fprintf(’x(%i)=%8.5f\n’,i, x(i)); end disp(’Pr disp(’Presion esione e una tecla tecla para continua continuar’); r’); pause

2. Presentamos la implementación del método de jacobi function function resultados resultados = jacobi(A, jacobi(A,b,x,e b,x,error) rror) %Marlo %Marlo Carranza Carranza Purca % la matriz matriz debe debe ser diagon diagonalm alment ente e domina dominante nte % sintaxix: sintaxix: % resultado= jacobi(matriz_A,vector_b,vec jacobi(matriz_A,vector_b,vector_inicial,err tor_inicial,error) or) % vpa( vpa(x’ x’,1 ,10) 0) te da el vecto vector r x’ con 10 cifr cifras as % A= A= [3 1 1 ; 1 3 1 ; 0 1 3] 3] % b= b= [5 5 4] 4]’ % x0=[1 0 0]’ % vpa(jacobi(A,b,x0,0.00001),10 vpa(jacobi(A,b,x0,0.00001),10) ) n=length(A);




89

for i=1:n, i=1:n, piv=A(i,i); b(i)=b(i)/piv; for j=1:n j=1:n A(i,j)=A(i,j)/piv; end end D=eye(n)-A; erel=1; ite=0; resultados=[ite,x’,erel]; whil while e erel erel > erro error, r, y=D*x+b; resto=y-x; erel=norm(resto,inf)/norm(y,inf); x=y; ite=ite+1; resultados=[resultados;ite,x’,erel]; end

3. Implentación del método de Gauss-Seidel function function resultados resultados = fjor(A,b, fjor(A,b,x0,nm x0,nmax,to ax,tol,ome l,omega) ga) %Marlo %Marlo Carranza Carranza Purca %JOR %JOR JOR JOR meth method od SOR SOR GAUS GAUS SEID SEIDEL EL % result resultado ados s = fjor(A fjor(A,B, ,B,X0, X0,NMA NMAX,T X,TOL, OL,OME OMEGA) GA) para para resolv resolver er el sistem sistema a % A*X= A*X=B B con con el meto metodo do JOR JOR . % TOL especi especific fica a la tolera toleranci ncia a sel metodo metodo. . % NMAX NMAX especi especific fica a el maximo maximo numero numero de iterac iteracion iones. es. % X0 especi especific fica a el valor valor inicia inicial l % OMEGA OMEGA es el parame parametro tro de relaja relajacio cion.A n.Algu lgunos nos autore autores s distin distingue guen n % OMEGA> OMEGA>1 1 SOBRE SOBRE RELAJA RELAJACIO CION N % OMEGA< OMEGA<1 1 SUB RELAJA RELAJACII CIION ON % OMEG OMEGA= A=1 1 el meto metodo do se tran transf sfor orma ma en Gaus Gauss s Seid Seidel el % ITER ITER es el nume numero ro maxi maximo mo de iter iterac acio ione nes s con con que que es calc calcul ulad ado o X % [n,m]=size(A); if (n < m) or ( n > m ) error(’O error(’Only nly square square systems’) systems’); ; end ite=0; r = b-A*x0 b-A*x0; ; r0=nor r0=norm(r m(r); ); err=no err=norm( rm(r); r); x=x0; x=x0; resultados=[ite,x’,err]; whil while e (err (err > tol) tol) & (ite (ite < nmax nmax) ) ite = ite + 1;




90

for i=1:n i=1:n s = 0; for j = 1:i-1, 1:i-1, s=s+A( s=s+A(i,j i,j)*x )*x(j) (j); ; end for j = i+1:n, i+1:n, s=s+A( s=s+A(i,j i,j)*x )*x(j) (j); ; end xnew(i,1)=omega*(b(i)-s)/A(i,i xnew(i,1)=omega*(b(i)-s)/A(i,i)+(1-omega)*x(i )+(1-omega)*x(i); ); end x=xnew; x=xnew; r=b-A*x; r=b-A*x; err=norm err=norm(r)/r (r)/r0; 0; resultados=[resultados;ite,x’,err]; end return



5

Resolución de sistemas no lineales

Introducción Vamos a estudiar los métodos que sirven para aproximar la raíces de un sistema de ecuaciones no lineales. De entre las distintas familias de métodos para abordar este problema, hemos elegido el método de Newton, la idea básica es es la misma que en el caso de una sola ecuación, linealizar el problema sustituyendo la función por su tangente, en este caso hiperplanos tangentes.

5.1. 5.1. Mé Métod todo o de Ne Newt wton on Recordemos el método de Newton se utilizo empleando la derivada de una función, para calcular su intersección con el eje x, esto es, al raíz ( vea la figura ). Dicho calculo se baso en la expansión de la serie de Taylor de primer

Figura 5.1: Representación gráfica del método de Newton orden ( xi+1 f ( f (xi+1 ) = f ( f (xi ) + (x


− x )f (x ) i

i


92

donde xi es el valor inicial de la raíz y xi+1 es el valor en el cual la recta tangente intersecta al eje x .En esta intersección , f ( f (xi+1) es por definición, igual a cero y entonces la ecuación anterior queda 0 = f ( f (xi ) + (x ( xi+1 xi )f  (xi ) f ( f (xi ) = (xi+1 xi )f  (xi ) f ( f (xi ) = xi+1 xi f  (xi ) f ( f (xi ) xi+1 = xi f  (xi )

− −

−

− −

−

que es la forma del método de Newton para una sola ecuación. En el caso de dos variables, una serie de Taylor de primer orden se escribe para cada ecuación no lineal como ( xi+1 ui+1 = ui + (x vi+1 = vi + (x ( xi+1

∂u + (y (y − y ) − x ) ∂u ∂x ∂y ∂v − x ) ∂v + (y (y − y ) ∂x ∂y i

i

i

i

i+1

i

i+1

i

i

i

de la misma manera como en la versión para una sola ecuación,la raíz aproximada aproximada corresponde a los valores x y y, donde ui+1 y vi+1 son iguales a cero 0 = ui + (x ( xi+1 0 = vi + (xi+1

∂u + (y (y − y ) − x ) ∂u ∂x ∂y ∂v + (y (y − y ) − x ) ∂v ∂x ∂y i

i

i

i

i+1

i

i+1

i

i

i

de aquí se tiene

 

(5.1)

i i xi+1 ∂u + yi+1 ∂u = ∂x ∂y

i

∂u i i ∂x

i i xi+1 ∂v + yi+1 ∂v ∂x ∂y

i

∂v i i ∂x

+ x −u + x = −v + x + x

i + yi ∂u ∂y i + yi ∂v ∂y

como se conocen los valores con subíndice i, las únicas incógnitas son xi+1 y yi+1 entonces se tiene el sistema de ecuaciones 5.1 , usando la regla de Cramer del Álgebra elemental se resuelve el sistema 5.1 y se tiene

xi+1 = x i

−

∂v i ∂y ∂u i ∂v i ∂x ∂y

ui

−vi ∂u∂yi , yi+1 = y i − − ∂u∂yi ∂v∂xi

∂u i ∂x ∂u i ∂v i ∂x ∂y

vi

−ui ∂v∂xi , − ∂u∂yi ∂v∂xi

J f ( f (x, y ) =

 

∂u i ∂x

∂u i ∂y

∂v i ∂x

∂v i ∂y

El denominador de cada una de estas ecuaciones se conoce como el Jacobiano del sistema. Donde la matriz Jacobiana de f = (ui, vi) es J f ( f (x, y ) Ejemplo 5.1.



 


93

Resolver el sistema de ecuaciones (5.2)

 

x2 + xy

= 10

y + 3xy 3 xy 2 = 57

Use como valores iniciales x = 1,5 y y = 3,5 tenga en cuenta que una solución es (2 ( 2, 3)



6

Interpolación Numérica

Introducción Son muchas y muy distintas las situaciones en las que, aparecen series de datos o resultados de mediciones experimentales de los cuales sólo se conocen como se comportan en una cierta cantidad finita de ítems ( por ejemplo, la población de un pais durante los años de una década ) y para los cuales se necesita encontrar una ley general que sirva para su tratamiento ( como conocer la población a mediados de uno de los años o la tasa de crecimiento de esta ). Esta ley general a la que se adaptan esos datos no es otra cosa que una función que tome los valores predeterminados. La función interpoladora servirá para para sustituir a la función desconocida, tanto para evaluarla en puntos donde no se conoce su valor, como para conocer su tasa e variación ( diferenciación numérica ) o su distribución acumulativa ( integración numérica ). La función interpoladora debe ser, por tanto fácil de evaluar, derivar e integrar . El método más común es la interpolación polinomial, recuerde que la fórmula general para un polinomio de grado n es P ( P (x) = a 0 + a + a1 x + a + a2 x2 +

n

· · · + a x n

dados n + 1 puntos hay un solo polinomio de interpolación de grado n que pasa por los n + 1 puntos.

6.1. 6.1.

Polino olinomio mio inter interpola polador dor de Lagran Lagrange ge

Dados los puntos

{x , x , . . . , x } 0

1

n

y los valores de una función f ( f (x)

{f ( f (x ), f ( f (x ), . . . , f ( x )} 0

1


n

6.1. Polinomio interpolador de Lagrange

95

encontrar encontrar un polinomio polinomio P n (x) de grad gradoo ≤ n que verifique P n(x0 ) = f ( f (x0 ), P n (x1 ) = f ( f (x1 ), . . . , Pn (xn) = f ( f (xn )

Dicho polinomio se conoce como polinomio interpolador de Lagrange.

6.1.1. 6.1.1.

Formula ormulació ción n de de Lagr Lagrang angee

Teorema 88 (Fórmula de interpolación de lagrange).

Sea una función f : [a, b] → R y {x0 , x1 , . . . , xn } ∈ [a, b] con xi  = x j si i =  j existe un único polinomio P n(x) que verifica P n (xi ) = f ( f (xi )

Para todo i = 0, 1, . . . , n , además este polinomio viene dado por la combinación lineal (6.1)

P n(x) = f ( f (x0 )L0 (x) + f + f ((x1 )L1 (x) +

· · · + f ( f (x )L (x) n

n

Donde para cada i ∈ {0, 1, . . . , n} se tiene n

(6.2)

Li (x) =

x xi j=0; j =0; j j  =i



−x −x

j j

además podemos agregar 1. Para que se verifique Li (x j ) = 0; j 0; j  = i n



Li (x) = ci

(x

j=0; j =0; j j =i



−x ) j

2. Para que se verifique Li (x j ) = 1; j 1; j = i n

ci =



(xi

j=0; j =0; j j =i



1

− x )− j

Definición 10.

El polinomio dado el la ecuación (2) se denomina polinomio de interpolación de Lagrange o simplemente polinomio de interpolación de f en los puntos { x0, x1, . . . , xn } y los polinomios dados en la ecuación (6.2) son los polinomios básicos de interpolación de Lagrange.




96

Ejemplo 6.1.

Encontrar el polinomio de interpolación para la siguiente función f ( f (x) = con tabla En este caso tenemos xi

f ( f (xi )

2

0.5

2.5

0.4

4

0.25

1 x

(x (2 (x L1(x) = (2, (2,5 (x L2(x) = (4

− 2,5)(x 5)(x − 4) − 2,5)(2 − 4) 2)(x − 4) − 2)(x 2)(2,5 − 4) − 2)(2, 2)(x − 2,5) − 2)(x − 2)(4 − 2,5)

L0(x) =

y el polinomio interpolador es

− 2,5)(x 5)(x − 4) (x − 2)(x 2)(x − 4) (x − 2)(x 2)(x − 2,5) + 0, 0,4 + 0, 0,25 (2, (2,5 − 2)(2, 2)(2,5 − 4) (4 − 2)(4 − 2,5) − 2,5)(2 − 4) 1 8 1 (x − 2,5)(x 5)(x − 4) + (x − 2)(x 2)(x − 4) + (x − 2)(x 2)(x − 2,5) 2 15 12

P 2 (x) = 0,5 P 2 (x) =

(x (2

Para aproximar f ( f (3) = f (3) f (3)



1 = 0.3 3 P 2 (3) = 0, 0,325

f (x) = x1 y el polinomio de interpolación Se muestra la gráfica de la función f ( P 2 (x) = 1,1167x 1167x2

− 6,8250x 8250x + 9, 9,6833

usamos el código siguiente p1=[1 -2.5],p2= -2.5],p2=[1 [1 -4],p3=[1 -4],p3=[1 -2] pl=0.5*conv(p1,p2)-8/15*conv( pl=0.5*conv(p1,p2)-8/15*conv(p2,p3)+1/12*con p2,p3)+1/12*conv(p1,p3) v(p1,p3) x=linspace(1,5,300); y=polyval(pl,x); hold hold on plot(x,0) plot(x,y,x,1./x,’:m’) legend(’p legend(’polin olinomio omio 1.1167 1.1167 x^2-6.8250x+ x^2-6.8250x+9.683 9.6833 3 ’,’funci ’,’función ón 1/x ’) X=[2 2.5 4 ]; Y=[0.5 0.4 0.25]; plot(X,Y,’or’)




97

Ejemplo 6.2.

Encontrar el polinomio de interpolación para la siguiente tabla En este caso xi

f ( f (xi )

2

1

3

2

-1

3

4

4

tenemos L0 (x) = L1 (x) = L2 (x) = L3 (x) =

(x

3)(x + 1)(x 1)(x − 4) (x − 3)(x 3)(x + 1)(x 1)(x − 4) − 3)(x = (−1)2(−4) 6 (x − 2)(x 2)(x + 1)(x 1)(x − 4) (x − 2)(x 2)(x + 1)(x 1)(x − 4) = −4 4(−1) (x − 2)(x 2)(x − 3)(x 3)(x − 4) (x − 2)(x 2)(x − 3)(x 3)(x − 4) = (−3)(−4)(−5) −60 (x − 2)(x 2)(x − 3)(x 3)(x + 1) (x − 2)(x 2)(x − 3)(x 3)(x + 1) = (2)(5)

10

y el polinomio interpolador es P 3 (x) =

(x

3)(x + 1)(x 1)(x − 4) (x − 2)(x 2)(x + 1)(x 1)(x − 4) (x − 2)(x 2)(x − 3)(x 3)(x − 4) − 3)(x − − + 6 2 20 2(x 2(x − 2)(x 2)(x − 3)(x 3)(x + 1) +

5 1 3 P 3 (x) = (x + 21x 21x2 64x 64x + 96) 60

−




98

Vamos a mostrar la gráfica del polinomio de interpolación P 3 (x) y los datos de la tabla, para lo cual damos la siguiente propuesta p=[1 p=[1/6 /60 0 21/6 21/60 0 -64/ -64/60 60 96/6 96/60] 0] x=linspace(-4,5,300); y=polyval(p,x); plot(x,y) legend(’polinomio legend(’polinomio 1/60(x^3+21x^2-64x+96)’) 1/60(x^3+21x^2-64x+96)’) roots(p) hold hold on plot(x,0) X=[2 3 -1 4 ]; Y=[1 2 3 4 ]; plot(X,Y,’or’)

6.1.2. 6.1.2.


Vamos a mostrar una función en Matlab que nos permitirá interpolar un determinado valor x , teniendo como datos las coordenadas (xi, yi ) function y=fbaselag(x,i,nodos) y=fbaselag(x,i,nodos) % funcio funcion n base base n=length(nodos)-1; y=1;




99

for j=0 : n, if i ~= j y=y.*(x-nodos(j+1))./(nodos(i y=y.*(x-nodos(j+1))./(nodos(i+1)-nodos(j+1)) +1)-nodos(j+1)); ; end end

function y=plagrange(x,xdato,ydato) y=plagrange(x,xdato,ydato) % funcio funcion n para para el polino polinomio mio de lagran lagrange ge % ej % xd=[0 0.2 0.5 0.6 0.9 % yd=[0 0.1987 0.4794 0.5646 0.7833 % y=pl y=plag agra rang nge( e(0. 0.1 1 , xd , yd ) % % n=length(xdato)-1; y=0; for i=0 : n, y=y+fbaselag(x,i,xdato).*ydato(i+1); end

1 ] 0.8415]

Ejemplo 6.3.

Calcularemos el polinomio de interpolación usando la tabla xi

0

yi 0

0.2

0.5

0.6

0.9

1

0.19 0.19887 0.479 .47944 0.56 .5646 0.78 0.78333 0.84 .8415

que se obtiene usando la función f ( f (x) = sen( sen (x) en el intervalo [−π, π] Para esto xd=[0 0.2 0.5 0.6 0.9 1 ] yd=[0 0.1987 0.4794 0.5646 0.7833 0.8415] plot(xd,yd,’or’) x=linspace(-pi,pi,100); ye=sin(x); hold hold on plot(x,ye,’r’) ya=plagrange(x,xd,yd); plot(x,ya,’b’) legend(’C legend(’Coord oordenada enadas s dato’,’ dato’,’ Función Función exacta exacta ’,... ’ polino polinomio mio de Lagran Lagrange’ ge’); ); hold hold off




10 0

Ejemplo 6.4.

Calcularemos el polinomio de interpolación usando la tabla Cuyos valores xi -1 - 1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 yi

1

0.8

0.6

0.4

0.2

0 0.2 0.4 0.6 0.8 1

corresponden a la función f ( f (x) = |x|, para lo cual emplearemos el código siguiente xd=[-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 yd=[ 1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 plot(xd,yd,’or’) x=linspace(-1,1,100); ye=abs(x); hold hold on plot(x,ye,’r’) ya=plagrange(x,xd,yd); plot(x,ya,’b’) legend(’C legend(’Coord oordenada enadas s dato’,’ dato’,’ Función Función exacta exacta ’,... ’ polino polinomio mio de Lagran Lagrange’ ge’); ); hold hold off


0.8 0.8

1 ]; 1 ];



10 1

Observación 6.1.

1. El ejemplo anterior se puede apreciar fuertes oscilaciones en los extremos del intervalo [−1, 1] conocidos como efectos de borde, esto muestra que que una función continua no puede aproximarse en general de manera arbitrariamente precisa mediante polinomios de Lagrange. 2. El cálculo del polinomio interpolador a partir de la fórmula requiere de muchas operaciones . Además una vez determinado el polinomio de {x0, x1, . . . , xn} ∈ [a, b], si añadimos interpolación de f en los puntos { un nuevo punto xn+1 disinto a los anteriores y queremos hallar el {x0, x1, . . . , xn+1} ∈ [a, b] polinomio de interpolación de f en los puntos { debemos repetir todo el proceso, dado que cambian todos los polinomios {L(x0), L(x1), . . . , L( básicos { L(xn )} por esta razón,tampoco será este el algoritmo para calcular el polinomio de interpolación de Lagrange 3. Al interpolar una función f en n + 1 puntos distintos, puede ocurrir  n, de hecho el grado de P n(x) puede ser pequeño aunque que el grad gradoo = n sea grande. Así por ejemplo P 3 (x) = 2x − 3 es el polinomio de interpolación de {−1, 0, 1, 2} Para mostrar f ( f (x) = x 4 − 2x3 − x2 + 4x 4x − 3 en los puntos {− lo dicho proponemos el siguiente código p=[2 p=[2 -3], -3], f=[1 f=[1 -2 -1 4 -3] -3] x=linspace(-2,3,300); y=polyval(f,x); yy=polyval(p,x); yy=polyval(p,x); plot(x,y) legend(’polinomio legend(’polinomio x^4-2x^3-x^2+4x-3’) x^4-2x^3-x^2+4x-3’) hold hold on plot(x,0),plot(x,yy) X=[-1 0 1 2];



6.2. Formulación de Newton

102

Y=[-5 -3 -1 -1 1]; plot(X,Y,’or’)

6.2. 6.2.

Formula ormulació ción n de Ne Newto wton n

Hasta ahora hemos visto como construir el polinomio de interpolación de Lagrange, ahora vamos a ver un método mas eficiente de construcción del polinomio de interpolación. Supongamos que una vez calculado dicho polinomio, necesitamos incluir otro punto de interpolación. Las razones pueden ser diversas se han obtenido nuevos valores experimentales para la tabla que estamos estamos interpolando interpolando,, queremos queremos obtener paulatinament paulatinamentee los polinomios de distintos grados para saber cuál nos interesa más, etc. En esta situación la fórmula de Lagrange no nos sirve pues, con sólo añadir un punto,todos punto,todos los polinomios polinomios básicos Li cambian. Necesitamos una fórmula que permita encontrar el nuevo polinomio interpolación a partir a partir del que habíamos hallado, de forma que el trabajo realizado pueda aprovecharse.

6.2.1. 6.2.1.

Inter Interpolac polación ión lineal lineal

La forma más simple de interpolación consiste en unir dos puntos por una linea recta. Dicha técnica se conoce como interpolación lineal f 1 (x) x

− f ( f (x ) −x 0

0

=

f ( f (x1 ) x1

− f ( f (x ) −x

f 1 (x) = f ( f (x0 ) +

0

0

f ( f (x1 ) x1

− f ( f (x ) − x (x − x ) 0

0

0

La notación f 1 (x) designa que éste es un polinomio de interpolación de primer grado y es una aproximación en diferencia dividida finita Mg. Marlo Carranza Purca Purca



103

Ejemplo 6.5.

calcular el polinomio que interpola a la función f ( f (x) = ln(x ln(x) en los puntos el xi

f ( f (xi )

1

0

6

1.79 .79175 1759

polinomio de interpolación de Newton es 1, 1 ,791759 0 (x 6 1 1,791759 f 1(x) = (x 1) 5 f 1(x) = 0 +

−

−

− 1)

−

y ahora en los puntos el polinomio de interpolación de Newton es xi

f ( f (xi )

1

0

4

1.38 .38629 6294

1,386294 0 (x 4 1 1,386294 g1 (x) = (x 1) 3 g1 (x) = 0 +

−

−

− 1)

−

Así usando el intervalo mas corto el error relativo se reduce,ambas interpolaciones se muestran a continuación X=[1

4


6 ];



104

Y=[0 Y=[0 1.38 1.3862 6294 94 1.79 1.7917 1759 59 ]; plot(X,Y,’or’) f=[1.79175 f=[1.791759/5 9/5 -1.791759/5], -1.791759/5], g=[1.38629 g=[1.386294/3 4/3 -1.386294 -1.386294/3] /3] x=linspace(0,7,300); y=polyval(f,x); y=polyval(f,x); yy=polyval(g,x); yy=polyval(g,x); hold hold on plot(x,y,’r’) plot(x,yy,’b’) plot(x,log(x),’g’) legend(’C legend(’Coord oordenada enadas s dato’,’p dato’,’polino olinomio mio f(x)= f(x)= 1.791759/5 1.791759/5 x - 1.791759/5 1.791759/5 ’,... ’,... ’polin ’polinomi omio o g(x)= g(x)= 1.3862 1.386294/ 94/3 3 x - 1.3862 1.386294/ 94/3 3 ’,’fun ’,’funció ción n exacta exacta ln(x)’ ln(x)’) ) plot(x,0) hold hold off

6.2.2. 6.2.2.

Inter Interpolac polación ión cuadrá cuadrátic tica a

En el ejemplo anterior el error que se comete es debido a que nuestra aproximación es hecha por una recta.En consecuencia, una estrategia para mejorar la estimación consiste en introducir una curva a la linea que une dichos puntos. Si se tienen tres puntos como datos estos se pueden ajustar por un único polinomio de grado dos, una forma particularmente conveniente de escribirlo es (6.3)

f 2 (x) = b 0 + b + b1 (x

+ b (x − x )(x )(x − x ) − x ) + b 0

2

0

1

para determinar los coeficientes evaluamos 1. x = x = x 0 para obtener b0 = f ( f (x0 ) Mg. Marlo Carranza Purca Purca



105

2. x = x = x 1 para obtener b1 = 3. x = x = x 2 para obtener b2 =

f (x1 ) f (x0 ) x1 x0

− −

f (x2 )−f (x1 ) x2 −x1

− f (xx11) x2 −x0

f (x0 ) x0

− −

Ejemplo 6.6. f (x) = ln(x ln(x) en los Calcularemos el polinomio que interpola a la función f ( puntos veamos xi

f ( f (xi )

1

0

4

1.38 .38629 6294

6

1.79 .79175 1759

x0 = 1 x1 = 4 x2 = 6

f ( f (x0 ) = 0 f ( f (x1 ) = 1,386294 f ( f (x1 ) = 1,791759

Aplicando las ecuaciones anteriores tenemos b0 = 0 1,386294 b1 = 4 1

−−

− 0 = 0,4620981

1,791759 1,386294 6 4

−

b2 =

6

− ,04620981 = −0,0518731

−1

Así el polinomio de interpolación es f 2 (x) = 0 + 0, 0,4620981(x 4620981(x

− 1) − 0,0518731(x 0518731(x − 1)(x 1)(x − 4) f (x) = −0,0519x 0519x + 0, 0,7215x 7215x − 0,6696 2

2

Mostramos el gráfico del polinomio de interpolación f 2(x) X=[1 4 6 ]; Y=[0 Y=[0 1.38 1.3862 6294 94 1.79 1.7917 1759 59 ]; plot(X,Y,’or’) f=0. f=0.46 4620 2098 981* 1*[0 [0 1 -1]-1]-0. 0.05 0518 1873 731* 1*[1 [1 -5 4 ] x=linspace(0,7,300); y=polyval(f,x); hold hold on plot(x,y,’b’) plot(x,log(x),’g’) legend(’Coordenadas legend(’Coordenadas dato’,’polinomio dato’,’polinomio f2(x)=-0.0519x^2+0.7215x-0.66 f2(x)=-0.0519x^2+0.7215x-0.6696 96 ’,... ’función ’función exacta exacta ln(x)’) ln(x)’) plot(x,0) hold hold off




6.2.3. 6.2.3.

106

Forma orma general general de los los polinomio polinomioss de interpo interpolac lación ión de Newton

El análisis anterior se puede generalizar para ajustar un polinomio de n-ésimo grado a n + 1 dato, el polinomio de n-ésimo grado es (6.4) f n (x) = b 0 + b1 (x − x0 ) + · · · + bn (x − x0 )(x )(x − x1)(x )(x − x2 ) . . . (x − xn−1 ) como se hizo antes con las interpolaciones lineales y cuadráticas, los puntos asociados con los datos se utilizan para evaluar los coeficientes b0 , b1 , . . . bn . Para un polinomio de n-ésimo grado se requieren n + 1 puntos (x0 , f ( f (x0 )), (x1 , f ( f (x1 )),. . . , (xn, f ( f (xn )) . Usamos estos datos y las siguientes ecuaciones para evaluar los coeficientes:

b0 = f ( f (x0) b1 = f [ f [x1 , x0 ] b2 = f [ f [x2 , x1 , x0 ]

.. . . = ..

(6.5)

bn = f [ f [xn, . . . , x1 , x0 ]

donde las evaluaciones de la función colocadas entre paréntesis son diferencias divididas finitas . Por ejemplo, la primera diferencia dividida finita en forma general se representa como f [ f [xi , x j ] =


f ( f (xi ) xi

− f ( f (x ) −x j

j



107

La segunda diferencia dividida finita, que representa la diferencia de las dos primeras diferencias divididas, se expresa en forma general como f [ f [xi , x j , xk ] =

f [ f [xi , x j ] xi

− f [ f [x , x ] −x j

k

k

En forma similar, la n-ésima diferencia dividida finita es f [ f [xn, xn−1 . . . , x1 , x0 ] =

f [ f [xn , . . . , x1 ] f [ f [xn−1 , . . . , x0 ] xn x0

− −

estas diferencias sirven para evaluar los coeficientes en las ecuaciones 6.5 los cuales se sustituirán en la ecuación en la ecuación 6.4 para obtener el polinomio de interpolació interpolaciónn f n (x) = f ( f (x0 ) + f + f [[x1 , x0 ](x ](x x0 ) + f + f [[x2 , x1 , x0 ](x ](x x0)(x )(x x1 ) + (6.6) + f [ f [xn, . . . x0 ](x ](x x0 )(x )(x x1)(x )(x x2 ) . . . (x xn−1 )

···

− −

xi

f ( f (xi ) P ri r imer o

−

−

−

−

−

···

que se conoce como polinomio de interpolación de Newton en diferencias divididas. Veamos la representación gráfica de la naturaleza recursiva de las diferencias divididas i

S eg undo

tercero

0 x0 f ( f (x0 )

f [ f [x1 , x0 ]

f [ f [x2 , x1 , x0 ] f [ f [x3 , x2 , x1 , x0 ]

1 x1 f ( f (x1 )

f [ f [x2 , x1 ]

f [ f [x3 , x2 , x1 ]

2 x2 f ( f (x2 )

f [ f [x3 , x2 ]

3 x3 f ( f (x3 ) Ejemplo 6.7.

Calcularemos el polinomio que interpola a la función f ( f (x) = ln(x ln(x) en los puntos veamos , utilizando la ecuación 6.4 con n = 3, el polinomio de tercer xi

1

yi 0

4

5

6

1.38 1.3862 6294 94 1.60 1.6094 9438 38 1.79 1.7917 1759 59

grado es f 3 (x) = b0 + b + b1 (x x0 ) + b + b2(x x0 )(x )(x x1 ) + b + b3 (x x0 )(x )(x x1 )(x )(x + f [[x1 , x0 ](x ](x x0 ) + f + f [[x2 , x1 , x0](x ](x x0 )(x )(x x1 ) + f 3 (x) = f ( f (x0 ) + f +f [ f [x3 , x2 , x1 , x0 ](x ](x x0 )(x )(x x1 )(x )(x x2 )

−

− −

−

−

−

−

−

−

−

−

−x ) 2

Las primera diferencias divididas del problema son 1,386294 4 1 1,791759 f [ f [x2 , x1 ] = 5 1,609438 f [ f [x3 , x2 ] = 6 f [ f [x1 , x0 ] =


− 0 = 0,4620981 − − 1,386294 = 0,2027326 −4 − 1,386294 = 0,1823216 −5 Manual Manual de métodos numérico numéricos s U.C.H. U.C.H.


108

Las segundas diferencias divididas son 0,2027326 5 0,1823216 f [ f [x3 , x2 , x1 ] = 6 f [ f [x2 , x1 , x0 ] =

La tercera diferencia dividida es f [ f [x3 , x2 , x1 , x0 ] =

− 0,4620981 = −0,05187311 −1 − 0,2027326 = −0,02041100 −4

−0,02041100 − (−0,05187311) = 0,007865529 6−1

Así el polinomio de interpolación es

f 3 (x) = 0+0, 0+0,4620981(x 4620981(x 1) 0,05187311(x 05187311(x 1)(x 1)(x 4)+0, 4)+0,007865529(x 007865529(x 1)(x 1)(x 4)(x 4)(x 6)

− −

6.2.4. 6.2.4.

−

−

−

−


Vamos a mostrar una función en Matlab que nos permitirá interpolar una función x via el polinomio de Newton basado en diferencias divididas function y=diffx(xdato,orden) y=diffx(xdato,orden) % funció función n difere diferenci ncias as en los nodos nodos n=length(xdato); for i=1:n-or i=1:n-orden, den, y(i)=xdato(orden+1)-xdato(i); end

function y=ddiv(xdato,ydato,orden,i) y=ddiv(xdato,ydato,orden,i) % función función diferenci diferencias as divididas divididas % xda xdato , ydat ydato o array rray de nodo nodos s xi , yi % orde orden n : de la dife difere renc ncia ia % i-ésim i-ésimo o elemem elememto to de las difere diferenci ncias as dividi divididas das n=lenght(xdato); for j=1:orde j=1:orden, n, ydato=diff(ydato)./diffx(xdato,j); end y=ydato(i+1);

function y=pddivnewton(x,xdato,ydato) y=pddivnewton(x,xdato,ydato) % funció función n polino polinomio mio de Newton Newton %n=length(xdato)-1; n=length(xdato)-1; y=ddiv(xdato,ydato,0,0);



−


109

prod=1; for i=1:n-1, i=1:n-1, prod=prod.*(x-xdato(i)); y=yddiv(xdato,ydato,i,0).*prod; end

Ejemplo 6.8.

Calcularemos el polinomio de interpolación de Newton, usando la tabla Cuyos xi

0

yi 0

0.2

0.5

0.6

0.9

1

0.19 0.19887 0.479 .47944 0.56 .5646 0.78 0.78333 0.84 .8415

valores corresponden a la función f ( f (x) = sen( sen(x), para lo cual emplearemos el código código siguiente xd=[ xd=[0 0 0.2 0.2 0.5 0.5 0.6 0.6 0.9 0.9 1]; 1]; yd=[0 yd=[0 0.1987 0.1987 0.4794 0.4794 0.5646 0.5646 0.7833 0.7833 0.8415] 0.8415]; ; plot(xd,yd,’or’) x=linspace(-1,1,100); ye=sin(x); hold hold on plot(x,ye,’r’) ya=pddivnewton(x,xd,yd); plot(x,ya,’b’) legend(’C legend(’Coord oordenada enadas s dato’,’ dato’,’ Función Función exacta exacta ’,... ’ polino polinomio mio de Newton Newton’); ’); hold hold off



6.3. Interpolación por partes

110

6.3. 6.3. Inte Interpo rpola laci ción ón por part partes es Introducción La palabra inglesa Spline ( trazador ) denota un instrumento flexible usado en dibujo técnico que sirve para trazar curvas suaves, se trata de una regla que puede ser adaptada flexionándola, a la forma que tome la curva que se desee dibujar.Esta técnica nació de la necesidad práctica por parte de los ingenieros navales en el diseño de los cascos de navíos, dadas las características propias de los medios acuosos, los barcos deben ofrecer la menor resistencia al momento de desplazarse, esto era conseguido de manera manual, se construían las secciones transversales mediante una herramienta denominada spline, en este método implícitamente se debe resolver un sistema de ecuaciones lineales esto hace bastante laboriosa la elaboración de los splines, pero con la ayuda de una computadora se puede hacer menos trabajosa.

6.3.1. 6.3.1.

Trazado razadores res lineal lineales es

La unión más simple entre dos puntos es una línea recta. Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto conjunto de funciones funciones lineales, lineales, f ( f (x) = f ( f (x) =

.. .

...

f ( f (x0 ) + m + m0 (x

−x ) f ( f (x ) + m + m (x − x ) 1

1

...

f ( f (x) = f ( f (xn−1 ) + m + mn−1 (x

0

,

x0 < x < x1

1

,

x1 < x < x2

−x − )

, xn−1 < x < xn

n 1

donde mi es la pendiente de la línea recta que une los puntos mi =

f ( f (xi+1 ) xi+1

f (x ) − f ( −x i

i

Estas ecuaciones se pueden usar para evaluar la función en cualquier punto entre x0 y x n localizando primero el intervalo dentro del cual está el punto. Después se usa la ecuación adecuada para determinar el valor de la función dentro del intervalo. El método es obviamente idéntico al de la interpolación lineal. Ejemplo 6.9.

Ajuste los datos de la tabla con trazadores de primer grado usando la tabla para lo cual realizamos los cálculos respectivos




111

x

3

4.5

f ( f (x)

2.5

1

f ( f (x) = 2,5 + f ( f (x) = f ( f (x) =

7

9

2.5 0.5

1,5 (x 1,5

− 3) 1 + (x − 4,5) 2,5 + − (x − 7) 1,5 2,5

2 2

, 3 < x < 4, 4 ,5 , 4,5 < x < 7 ,

7 < x < 9

Así se obtiene f ( f (x) = f ( f (x) = f ( f (x) =

2,5 + (x (x

− 3) 1 + 0,6(x 6(x − 4,5) 2,5 − (x − 7)

, 3 < x < 4, 4 ,5 , 4,5 < x < 7 ,

7 < x < 9

es decir

f ( f (x) =

 − 

5,5

−x

, 3 < x < 4, 4 ,5

1,7 + 0, 0 ,6x , 4,5 < x < 7 9,5

−x

,

7 < x < 9

Para graficar el trazador lineal, digite el siguiente código en Matlab x=linspac x=linspace(3, e(3,9,300 9,3000);xx 0);xx=[3 =[3 4.5 7 9];yy=[2.5 9];yy=[2.5 1 2.5 0.5]; 0.5]; y=(5.5-x).*((3< y=(5.5-x).*((3< x)&(x<4.5))+(-1.7+0.6.*x).*((4 x)&(x<4.5))+(-1.7+0.6.*x).*((4.5
y se obtiene el gráfico




6.3.2. 6.3.2.

112

Trazado razadores res cuadrá cuadrátic ticos os

Para asegurar que las derivadas derivadas m -ésimas sean continuas en los nodos, se debe emplear un trazador de un grado, de al menos, m + 1. En la práctica se usan con más frecuencia polinomios de tercer grado o trazadores cúbicos que aseguran aseguran primera primera y segunda derivadas derivadas continuas. continuas. Aunque las derivadas derivadas de tercer orden y mayores podrían ser discontinuas cuando se usan trazadores cúbicos, por lo común no pueden detectarse en forma visual y, en consecuencia, consecuencia, se ignoran. ignoran. El objetivo de los trazadores cuadráticos es obtener un polinomio de segundo grado para cada intervalo entre los datos. De manera general, el polinomio en cada intervalo se representa como + ci f i (x) = ax i2 + bxi + c

Para 3 datos(i = 0, 1, 2) existen, 2 intervalos y, en consecuencia, 3 constantes desconocidas ( las a, byc ) por evaluar. Por lo tanto, se requieren 3 ecuaciones o condiciones para evaluar las incógnitas. Ejemplo 6.10.

Veamos un ejemplo concreto, consideremos los siguientes datos: procedamos a xi 3

4.5

yi 2. 2.5

1

7

9

2.5 0.5

calcular la interpolación por splines de grado 2. Primero que nada, vemos que se forman tres intervalos: [3;4, [3;4,5], [4, [4,5;7], [7; [7; 9], en cada uno de estos intervalos, debemos definir una función polinomial de grado 2, como sigue: F ( F (x) =

  

ax21 + bx + bx1 + c + c1

si x ∈ [3;4, [3;4,5], 5],

ax22 + bx + bx2 + c + c2

si x ∈ [4, [4,5;7], 5;7],

si x ∈ [7; [7; 9] Hacemos que la spline pase por los puntos de la tabla de datos, es decir, se debe cumplir que: ax23 + bx + bx3 + c + c3

F (3) F (3) = 2, 2,5, F (4 F (4,,5) = 1, 1, F (7) F (7) = 2, 2,5, F (9) F (9) = 0, 0,5

Así, se forman las siguientes ecuaciones:

     

F (3) F (3) = 2, 2,5 entonces 9a1 + 3b 3 b1 + c + c1 = 2,5 F (4 F (4,,5) = 1 entonces

F (7) F (7) = 2, 2,5 entonces

   

(4, (4,5)2 a1 + 4, 4 ,5b1 + c + c1 = 1 (4, (4,5)2 a2 + 4, 4 ,5b2 + c + c2 = 1 49a 49a2 + 4, 4 ,7b2 + c + c2 = 2,5 49a 49a3 + 4, 4 ,7b3 + c + c3 = 2,5

F (9) F (9) = 0, 0,5 entonces 81a 81a3 + 9b3 + c + c3 = 0,5




113

Hasta aquí, tenemos un total de 6 ecuaciones con 9 incógnitas. El siguiente paso es manejar la existencia de las derivadas continuas. En el caso de las splines de grado 2 , necesitamos que la spline tenga derivada continua de orden k − 1 = 1, es decir, primera derivada continua. Calculamos la primera derivada:

F  (x) =

  

2a1 x + b + b1

si x ∈ [3;4, [3;4,5]

2a2 x + b + b2

[4,5;7] si x ∈ [4,

2a3 x + b + b3

si x ∈ [7; [7; 9]

Vemos que esta derivada está formada por segmentos de rectas, que pudieran presentar discontinuidad en los cambios de intervalo. Es decir, las posibles discontinuidades son x = 4,5 y x = 7. Por lo tanto para que F (x) sea continua, se debe cumplir que: 

2a1 (4, (4,5) + b + b1 = 2a2 (4, (4,5) + b + b2

→ 9a + b + b 1

1

= 9a2 + b + b2

También debe cumplirse que: 2a2 (2) + b + b2 = 2a3 (7) + b + b3

14a + b + b → 14a 2

2

= 14a 14a3 + b + b3

Así, tenemos un total de 8 ecuaciones vs. 9 incógnitas; esto nos da un grado de libertad para elegir alguna de las incógnitas. Elegimos por simple conveniencia a1 = 0. De esta forma, tenemos un total de 8 ecuaciones con 8 incógnitas. Estas son las siguientes: 3b1 + c + c1 = 2,5

49a 49a3 + 7b 7 b3 + c + c3 = 2,5

4,5b1 + c + c1 = 1

81a 81a3 + 9b 9 b3 + c + c3 = 0,5

20, 20,25a 25a2 + 4, 4 ,5b2 + c + c2 = 1 b1 49a 49a2 + 7b 7 b2 + c + c2 = 2,5

− 9a − b 2

2

=0

14a 14a2 + b + b2 = 14a 14a3 + b + b3

Este sistema de ecuaciones tiene la siguiente forma matricial:

      

3

1

0

0

0

0

0

0

4,5 1

0

0

0

0

0

0

0

0 20 20,25 4,5 1

0

0

0

0

0

49

7

1

0

0

0

0

0

0

0

0

49

7

1

0

0

0

0

0

81

9

1

1

0

−9 −1

0

0

0

0

0

0

14

0


1

−14 −1

0

           

b1 c1 a2 b2 c2 a3 b3 c3

      

=

      

2,5 1 1 2,5 2,5 0,5 0 0

      



114

Se obtiene la siguiente solución: a1 = 0

a2 = 0,64

a3 =

b1 =

b2 =

b3 = 24, 24,6

−1

−6,76

c1 = 5,5 c2 = 18, 18,46

c3 =

−1,6 −91, 91,3

Sustituyendo estos valores (junto con a1 = 0 ), obtenemos la función spline cuadrática que interpola la tabla de datos dada:

F ( F (x) =

 −   −

si x ∈ [3;4, [3;4,5], 5],

x + 5

0,64x 64x22

76x + 18, 18 ,46 − 6,76x 1,6x + 24 2 4,6x − 91, 91,3 2 3

2

3

[4,5;7], 5;7], si x ∈ [4,

si x ∈ [7; [7; 9]

La gráfica que se muestra a continuación, contiene tanto los puntos iniciales de la tabla de datos, así como la spline cuadrática. Para graficar el trazador cuadrático, digite el siguiente código en Matlab x=linspac x=linspace(3, e(3,9,300 9,3000);xx 0);xx=[3 =[3 4.5 7 9];yy=[2.5 9];yy=[2.5 1 2.5 0.5]; 0.5]; y=(5.5-x).*((3
y se obtiene el gráfico



7

Integración numérica

Introducción Uno de los problemas matemáticos más antiguos es el del cálculo de área que encierra una curva; quizás el ejemplo más significativo sea el intento de la cuadratura del circulo, que condujo a la aparición y estudio del número π . Cuando deseamos calcular el valor de la integral definida de una función f ( f (x) a veces nos damos con la sorpresa que su antiderivada no se puede expresar mediante mediante funciones funciones elementales elementales conocidas, por ejemplo ejemplo



b

2

e−x dx

a

esta integral cuando es indefinida es imposible de expresarla e términos de 2 funciones elementales, pero la integral definida existe pues sabemos que e −x es continua en todo R. La integral definida de una función positiva representa el área limitada por la curva y las rectas x = a = a y x = b = b , esto nos permitirá ampliar las aplicaciones donde se necesita calcular áreas, volúmenes, arcos, trabajo y calculo de energía entre otras aplicaciones. Las fórmulas de Newton-Cotes son los tipos de integración numérica más comunes, se basan en la estrategia de de reemplazar una función complicada o datos tabulados por un polinomio de aproximación que es fácil de integrar.



b

(7.1)

a

 ≡

b

f ( f (x)dx

f n (x)dx

a

donde f n (x) = a 0 xn + a1xn−1 + · · · + an−1 x + a + an

7.1. 7.1. Fórm Fórmul ula a del del trape trapeci cio o La regla del trapecio es es la primera fórmula de integración que vamos a ver, corresponde al caso donde el polinomio de interpolación de la ecuación 7.1 es de primer grado. 115 Mg. Marlo Carranza Purca Purca

7.1. Fórmula del trapecio

116

Teorema 89.



b

 ≡

b

f ( f (x)dx

a

f 1 (x)dx

a

(b

≡

f (a) + f + f ((b) − a) f ( 2

Demostración.

Para la demostraci demostración ón de este Teorema, Teorema, vamos vamos a reemplazar reemplazar la función f 1 (x) = f ( + f [[b, a](x ](x − a), luego vamos a efectuar f (a) + f las operaciones respectivas, así que empezamos



 ≡  ≡ {

b

b

f ( f (x)dx

f 1 (x)dx

a

a

b

f ( f (a) + f + f [[b, a](x ](x

a



b

 ≡ {  ≡ b

f ( f (x)dx

a

− f ( f (a) − a (x − a)}dx f ( f (b) − f ( f (a) f ( f (a)dx + dx + (x − a)d(x − a) b−a f ( f (b) − f ( f (a) (b − a) f (a)(b )(b − a) + ≡ f ( 2 b−a (b − a) f (a)(b )(b − a) + {f ( f (b) − f ( f (a)} ≡ f ( 2 f (a) + f + f ((b) ≡ (b − a) f ( 2 f ( f (a) +

a

b

f ( f (x)dx

a

f ( f (b) b



b

b

a



− a)}dx

a

2

Ejemplo 7.1.

Dada la función f ( f (x) = 4 − 2x2 , x ∈ [0, [0, 1] el valor aproximado de la integral de la función f que se obtiene por la fórmula del trapecio es



1

0

(4

2

− 2x )dx ≡ ≡

f (0) f (0) + f + f (1) (1) 2 3

Para ver el gráfico apunte el siguiente código, y vea la figura 7.2 X=[0 1 ]; ]; Y=[4 Y=[4 2]; plot(X,Y,’or’) f=[-2 f=[-2 0 4];g=[ 4];g=[-2 -2 4];




117

x=linspace(0,1,300); y=polyval(f,x); y=polyval(f,x); yy=polyval(g,x); yy=polyval(g,x); hold hold on plot(x,y,’r’) plot(x,yy,’b’) legend(’C legend(’Coord oordenada enadas s dato’,’p dato’,’polino olinomio mio f(x)= f(x)= -2x^2-4 -2x^2-4 ’,... ’polin ’polinomi omio o f1(x)= f1(x)= -2x+4 -2x+4 ’) plot(x,0) hold hold off

Figura 7.1: Fórmula del trapecio Ejemplo 7.2.

Dada la función f ( f (x) = 0,2 + 25x 25x

2

200x − 200x

+ 675x 675x3

4

900x − 900x

+ 400x 400x5 , x

[0;0,8] el valor ∈ [0;0,

aproximado de la integral de la función f que se obtiene por la fórmula del trapecio es



0,8

0

(0, (0,2 + 25x 25x

− 200x 200x

2

+ 675x 675x3

4

− 900x 900x

+ 400x 400x5 )dx

≡ ≡ ≡

(0, (0,8

+ f (0 (0,,8) f (0) + f − 0) f (0) 2

0,2 + 0, 0 ,232 2 0,1728 0,8

Pero el valor exacto de la integral obtenido en forma analítica es 1,640533, el 89,5 % Para ver el gráfico apunte el error porcentual de la aproximación es 89, siguiente código, y vea la figura 7.2 X=[0 0.8 ]; Y=[0 =[0.2 0.23 0.232 2]; plot(X,Y,’or’) f=[4 f=[400 00 -900 -900 675 675 -200 -200 25 25 .02] .02];g ;g=[ =[0. 0.04 04 x=linspace(0,0.9,300);


0.2] 0.2]; ;



118

y=polyval(f,x); y=polyval(f,x); yy=polyval(g,x); yy=polyval(g,x); hold hold on plot(x,y,’r’) plot(x,yy,’b’) legend(’Coordenadas legend(’Coordenadas dato’,... ’polinomio f(x)=0.2+25x-200x^{2}+675x^{3} f(x)=0.2+25x-200x^{2}+675x^{3}-900x^{4}+400x^ -900x^{4}+400x^{5}’,... {5}’,... ’pol ’polin inom omio io f1(x f1(x)= )= 0.04 0.04x x + 0.2 0.2 ’) plot(x,0) hold hold off

Figura 7.2: Fórmula del trapecio

7.1.1. 7.1.1.

Fórm Fórmula de dell trapec trapecio io comp compues uesta ta

Una forma de mejorar la precisión de la fórmula del trapecio consiste en dividir dividir el interv intervalo alo de integraci integración ón de a a b en varios segmentos y aplicar la fórmula del trapecio a cada uno de ellos, las áreas de los segmentos se suman después para obtener la integral en todo el intervalo, las ecuaciones obtenidas son llamadas formulas de integración, de aplicación multiple o compuestas.

Teorema 90.

Tomamos n + 1 puntos igualmente espaciados { x0 , . . . xn } en consecuencia existen n segmentos del mismo ancho h = b−na . Si a y b se designan como x0 y xn respectivamente la integral se designara como



xn

f ( f (x)dx = dx =

x0


h 2



n 1

f ( f (x0 ) + 2

−

 i=1



+ f ((xn ) f ( f (xi ) + f



119

Demostración.

Para realizar realizar la demostraci demostración, ón, vamos vamos a aplicar la regla regla del trapecio en cada uno de los n subintervalos que se generan al tomar los n + 1 en el dominio de la función, así tenemos





xn

x1

f ( f (x)dx =

x0

x2

f ( f (x)dx + dx +

x0

= h



xn

f ( f (x)dx + dx +

x1

f ( f (x0 )



··· +

f ( f (x)dx

xn

1

−

f (x ) f ( f (x ) − f ( f (x ) f ( f (x ) − f ( f (x ) − f ( + h + · · · + h 1

1

2

2

1

2

h f ( f (x0) + 2f 2 f ((x1 ) + + 2f 2f ((xn−1 ) + f + f ((xn) 2 n−1 h = f ( f (x0 ) + 2 f ( f (xi ) + f + f ((xn ) 2 i=1 =

{

···



(7.2)





2

2

}

Ejemplo 7.3.

Dada la función f ( f (x) = sen x, x ∈ [0, [0, π ] el valor aproximado de la integral de la función f que se obtiene por la fórmula del trapecio compuesta con dos segmentos es



 ≡  ≡

π

π/2 π/2

f ( f (x)dx

0



0

sen(x sen(x)dx + dx +

0

≡ ≡

π

f ( f (x)dx

0

  

sen(x sen(x)dx

π/2 π/2

π sen π/2 π/ 2 + sen sen 0 2 2 π π + 4 4 π = 1,5708 2

π + 2

sen π/2 π/ 2 + sen π 2



Ejemplo 7.4. f (x) = sen x, x ∈ [0, [0, π ] el valor aproximado de la integral de Dada la función f ( la función f que se obtiene por la fórmula del trapecio compuesta con cuatro segmentos es



π

 ≡  ≡

π/4 π/4

f ( f (x)dx

0


3π/4 π/4


π/4 π/4

π 4

sen π/4 π/4 + sen sen 0 π + 2 4 π sen3π/ sen3π/44 + sen sen 2π/4 π/ 4 + 4 2

π

0

2π/4 π/4

0

≡



      √ √ 

f ( f (x)dx

≡ ≡

  √

√ 2



π


2π/4 π/ 4

sen(x sen(x)dx

3π/4 π/4



sen π/4 π/ 4 + sen sen 2π/4 π/ 4 2 π sen π + sen sen 3π/4 π/ 4 + 4 2



2 2 2 π + + 1+ 1+ + 8 2 2 2 2 π (2 2 + 2) 8

√

1,896123



7.2. Fórmula de Simpson

120

Ejemplo 7.5.


2

200x − 200x

+ 675x 675x3

4

900x − 900x

+ 400x 400x5 , x

[0;0,8] el valor ∈ [0;0,

aproximado de la integral de la función f que se obtiene por la fórmula del trapecio compuesta para n=2 es



0,8

(0, (0,2 + 25x 25x

0

− 200x 200x

2

+ 675x 675x3

4

− 900x 900x

+ 400x 400x5 )dx

≡ ≡ ≡

7.1.2. 7.1.2.

f (0) f (0) + 2f 2f (0 (0,,4) + f + f (0 (0,,8) 4 0,2 + 2(2, 2(2,456) + 0, 0,232 0,8 4 1,0688

(0, (0,8)


Construiremos una función en Matlab que calcule la integral de la función f ( f (x) = sen(x sen(x) en el intervalo [0, [0, π ] usando la regla compuesta del trapecio; para ello esta función debe recibir como argumentos el intervalo y el número de nodos a usarse function y=trapecio(a,b,n) y=trapecio(a,b,n) % calc calcul ula a la inte integr gral al de una una func funció ión n seno seno %con %con la regla regla del trapec trapecio io compue compuesta sta h=(b-a h=(b-a)/( )/(n-1 n-1); ); s=0; s=0; for i=1:n-1, i=1:n-1, s=s+sin(a+i*h); end y=(h/2)*(sin(a)+2*s+sin(b));

Luego la llamamos desde Matlab como for i=1:5, i=1:5, integral=vpa(trapecio(0,pi,i*10)) end

7.2. 7.2. Fórm Fórmul ula a de Simp Simpso son n Además de aplicar la regla del trapecio con una segmentación mas fina, otra forma de obtener una estimación más exacta de una integral consiste en usar un polinomio de grado superior para unir los puntos.Por ejemplo, si hay otro punto a la mitad entre f ( f (a) y f ( f (b), los tres puntos se pueden unir con una parábola, si hay dos puntos igualmente espaciados entre f ( f (a) y f ( f (b). los cuatro puntos se pueden unir mediante un polinomio de tercer grado, las formulas que resultan de tomar las integrales bajo estos polinomios se conocen como reglas de Simpson. Mg. Marlo Carranza Purca Purca



7.2. 7.2.1. 1.

121

Fórm Fórmul ula a de de Sim Simps pson on 1/3 1/3

Teorema 91.

La regla de Simpson resulta cuando un polinomio de interpolación de segundo grado se sustituye en la ecuación (7.1), con h = b−2 a



x2

f ( f (x)dx = dx =

b

− a [f ( f (x ) + 4f 4 f ((x ) + f + f ((x )] 0

6

x0

1

2

Demostración.

Para la demostr demostración ación Vamos Vamos a usar un polinomio polinomio de interpolación de segundo grado en la ecuación ( 7.1 ) , así tenemos



x2

 ≡  ≡ { x2

f ( f (x)dx

x0

f 2(x)dx

x0 x2

f ( f (x0 ) + f + f [[x1 , x0 ](x ](x

x0

0

2

1

0

0

1

f (x ) + 4f 4 f ((x ) + f + f ((x )] ≡ h3 [f ( ≡ b −6 a [f ( f (x ) + 4f 4 f ((x ) + f + f ((x )] 0

(7.3)

− x ) + f + f [[x , x , x ](x ](x − x )(x )(x − x )} dx

1

0

2

1

2

Ejemplo 7.6.


2

− 200x 200x

+ 675x 675x3

4

− 900x 900x

+ 400x 400x5 , x

∈ [0;0, [0;0,8] vamos a

obtener el valor aproximado de la integral de la función f por la fórmula de Simpson en el intervalo [0;0, [0;0,8] . Veamos f (0) f (0) = 0, 0,2 ; f (0 f (0,,4) = 2, 2,456; f (0 f (0,,8) = 0, 0,232, por lo tanto de la ecuación 7.3 tenemos



0,8

(0, (0,2 + 25x 25x

0

200x − 200x

2

+ 675x 675x3

4

900x − 900x

+ 400x 400x5 )dx

≡ ≡ ≡

f (0) f (0) + 4f 4f (0 (0,,4) + f + f (0 (0,,8) 3 0,2 + 4(2, 4(2,456) + 0, 0,232 0,8 6 1,367467 0,8

Pero el valor exacto de la integral obtenido en forma analítica es 1,640533, el error porcentual de la aproximación es 16, 16,6 o/o

7.2.2. 7.2.2.

Fórm Fórmula de Simpso Simpson n comp compues uesta ta

Así como en la regla del trapecio, la fórmula de Simpson mejora al dividir el intervalo de integración en varios segmentos de un mismo tamaño y en un número par de ellos de longitud h = b−na Mg. Marlo Carranza Purca Purca



122

Teorema 92 (Fórmula (Fórmula de Simpson compuesta).



xn x0

f ( f (x) = (b

− a)

f (x0 )+4

n−1 f (xi )+2 i=1,3,5 f (



3n

n−2 f (xj )+f )+f ((xn ) j =2,4,6 f (



Demostración.



xn

 ≡

x2

f ( f (x)

x0



x4

f ( f (x)dx + dx +

x0

 ···

xn

f ( f (x)dx + dx +

x2

f ( f (x)dx

xn

2

−

f ( f (x0 ) + 4f 4 f ((x1 ) + f + f ((x2 ) f ( f (xn−2 ) + 4f 4 f ((xn−1 ) + f + f ((xn ) + + 2h 2h 6 6 n−1 n−2 f ( f (x0 ) + 4 i=1, f (xi ) + 2 j=2 f (x j ) + f + f ((xn ) =1,3,5 f ( j =2,,4,6 f ( (b a) 3n 2h

≡ ≡

···



−



Ejemplo 7.7.


2

200x − 200x

+ 675x 675x3

4

900x − 900x

+ 400x 400x5 , x

[0;0,8] vamos a ∈ [0;0,

obtener el valor aproximado de la integral de la función f por la fórmula de Simpson compuesta en el intervalo [0;0, con n = 4 . [0;0,8] con Veamos f (0) f (0) = 0, 0,2 ; f (0 f (0,,2) = 1, 1,288 ; f (0 f (0,,4) = 2, 2,456; f (0 f (0,,6) = 3, 3,464; f (0 f (0,,8) = 0, 0,232, por lo tanto de la ecuación 7.3 tenemos



0,8 (0, (0,2 0

+ 25x 25x

200x − 200x

≡ ≡ ≡

2

+ 675x 675x3

4

900x − 900x

+ 400x 400x5)dx

≡

f (0) f (0) + 4(f 4(f (0 (0,,2) + f + f (0 (0,,6)) + 2f 2f (0 (0,,4) + f + f (0 (0,,8) 12 0,2 + 4(1, 4(1,1288 + 3, 3,464) + 2(2, 2(2,456) + 0, 0,232 0,8 12 1,623467 0,8

Pero el valor exacto de la integral obtenido en forma analítica es 1,640533, el error porcentual de la aproximación es 1,04 o/o

7.2.3. 7.2.3.


Mostraremos un código para calcular la integral de una función usando la fórmula fórmula de Simpson compuesta function area=simpson(y,a,b) area=simpson(y,a,b) % regl regla a comp compue uest sta a de Simp Simpso son n para para un nume numero ro %impar %impar de nodos nodos unifor uniformem mement ente e distri distribui buidos dos % y : son son las las absc abscis isas as de la func funció ión n y=f( y=f(x) x) % dond donde e [a,b [a,b] ] son son los los limi limite tes s de integ integra raci ción ón n =length(y) =length(y); ;




123

m=(n-1)/2; h=(b-a)/(2*m); sigma1=0; sigma2=0; %pares for for i=1: i=1:mm-1 1 , sigma1=sigma1+y(2*i+1); end displa display(’ y(’ pares pares ’) for i=0:m-1 i=0:m-1 sigma2=sigma2+y(2*i+2); end area =(h/3)*( =(h/3)*( y(1)+2*sig y(1)+2*sigma1+4 ma1+4*sigm *sigma2+y a2+y(n) (n) );

Luego la llamamos desde Matlab como simpson(sin(linspace(0,pi,11)),0,pi)

Matlab proporciona dos dos funciones basadas en reglas de Simpson 1. quad(funcx,a,b) : Calcule la integral definida de la función funcx desde a hasta b , empleando la fórmula de Simpson (regla de tres puntos) 2. quad8(funcx,a,b) : De forma similar a quad, pero emplea 8 puntos en la fórmula de Newton Cotes. Ejemplo 7.8.



Aproxime el valor 05 sen(x sen(x)dx con estas dos funciones, luego compare con el valor exacto de la integral. Para lo cual escriba integral =[quad(’sin(x)’,0,5),... =[quad(’sin(x)’,0,5),... quad8(’sin(x)’,0,5), quad8(’sin(x)’,0,5), 1-cos(5)]

7.2.4. 7.2.4.

Ejerci Ejercicio cio,, mecán mecánica ica de la la fract fractura ura

El crecimiento de una grieta en el borde de una placa por ciclo de esfuerzos viene dado por la ecuación de Paris (7.4)

√

da = A ∆σY a dN





m

Donde N es es el número de ciclos, A y m son constantes del material,∆σ es la diferencia de esfuerzos a tensión sobre la pieza y Y viene dado por la ecuación . Cuando se observa una grieta de tamaño a0, el número de ciclos




124

restante para fractura catastrófica de la pieza se obtiene separando las variables ariables e integrando integrando la ecuación (7.4) (7.5)

dN =

da m A (∆σY (∆σY a)

√

Supóngase que se tiene la placa del ejemplo 3.2.4 con ∆σ = 17, 17,7, 9 − A = 6,6 × 10 , m = 2,25 , con una grieta inicial de a 0 = 0,25 . El número de ciclos restante para falla es N f f =



0,62 0,25

da

6,6 10

×

 9

−



17, 17,78 1,99 0,41(

−

a 2,5

  2 3 4 √ 2,25 +18,70( 2a,5 ) −38, 38,48( 2a,5 ) +53, +53,85( 2a,5 ) a )+18,

En este caso, la solución en MATLAB se implementa de la manera siguiente: function [I]=integral(I0,X,Y) [I]=integral(I0,X,Y) N=numel(X) N=numel(X); ; I(1)=I0; I(1)=I0; for n=2:N n=2:N I(n)=I(n-1)+0.5*(Y(n)+Y(n-1) I(n)=I(n-1)+0.5*(Y(n)+Y(n-1))*(X(n)-X(n-1)) )*(X(n)-X(n-1)); ; end plot(X,I,’k-’);

luego hay que escribir a=0.25:0.01:0.62; for i=1:nume i=1:numel(a) l(a) f(i)=1/(6.6e-9*(17.78*(1.99-0 f(i)=1/(6.6e-9*(17.78*(1.99-0.41*(a(i)/2.5)+ .41*(a(i)/2.5)+18.70*(a(i)/2.5 18.70*(a(i)/2.5)^2-... )^2-... 38.48*(a(i)/2.5)^3+53.85*(a(i 38.48*(a(i)/2.5)^3+53.85*(a(i)/2.5)^4)*sqrt( )/2.5)^4)*sqrt(a(i)))^2.25); a(i)))^2.25); end

debemos tomar I0=0; [N]=integral(I0,a,f);

El número de ciclos hasta la falla así calculado es N f f = 37109 ciclos. La gráfica de crecimiento de la grieta con los ciclos se muestra en la figura




125

Figura 7.3: Crecimiento de la grieta con el número de ciclos hasta la falla



8

Métodos Métodos de Runge - Kutta

Introducción Las ecuaciones diferenciales ordinarias tienen una enorme importancia en la práctica de de la ingeniería lo cual se debe a que muchas leyes físicas están expresadas en la razón de cambio de una cantidad, más que en términos de la cantidad misma.Entre los ejemplos tenemos desde los modelos de de predicción demográfica (razón de cambio de la población), hasta la aceleración de un cuerpo que cae (razón de cambio de la velocidad). Aquí nos vamos a dedicar a la solución de de ecuaciones diferenciales ordinarias de la dy forma dx = f ( f (x, y ) lo cual se puede expresar así

yi+1

dy = f ( f (t, y ) dt y = f ( f (t, y ) t yi = f ( f (ti , yi ) t yi+1 = yi + f + f ((ti , yi )

  −



t

De acuerdo con esta ecuación, la pendiente estimada f ( f (t, y ) se usa para extrapolar desde un valor anterior y i a un nuevo valor yi+1 en una distancia t , está fórmula se aplica paso a paso para calcular un valor posterior, por lo tanto para trazar la trayectoria de la solución, aquí se toma la la pendiente a l inicio del intervalo como una aproximación de la pendiente promedio sobre todo el intervalo, tal procedimiento llamado método de Euler, es lo veremos a continuación,todos los métodos de un paso que se expresan de esta forma general solo van a diferir en la manera como se estima la pendiente,dando lugar a predicciones más exactas, todas estas técnicas se conocen como métodos de Runge-Kutta


8.1. Méto do de Euler

8.1 8.1.

12 7

Métod étodo o de Eule Eulerr

La primera derivada ofrece una estimación directa de la pendiente en xi , φ = f = f ((xi , ti ), donde f ( f (xi , yi ) es la ecuación diferencial evaluada en xi y yi se tiene (8.1)

yi+1 = y i + φ + φ((xi , yi )h

Esta fórmula se conoce como Método de Euler ( o Euler Cauchy o de punto pendiente ). Se predice un nuevo valor de y usando la pendiente (igual a la primera derivada en el valor original de x ) para extrapolar linealmente sobre el tamaño de paso h. Ejemplo 8.1.

Resolver le problema de valor inicial PVI.

 



y =

y 10

y (0) = 1000, 1000, t

[0; 5] ∈ [0; t

tenga en cuenta la solución exacta y (t) = 1000e 1000e 10 Solución Vamos a usar la ecuación 8.1 con h = 1 , obtenemos los resultados y (1) = y = y(0) (0) + f + f (0 (0,, 1000) = 1000 + 100 = 1100 y (2) = y = y(1) (1) + f + f (1 (1,, 1100) = 1100 + 110 = 1210 y (3) = y = y(2) (2) + f + f (2 (2,, 1210) = 1210 + 121 = 1331 y (4) = y = y(3) (3) + f + f (3 (3,, 1331) = 1331 + 133, 133,1 = 1464, 1464,1 y (5) = y = y(4) (4) + f + f (4 (4,, 1464, 1464,1) = 1464, 1464,1 + 146, 146,41 = 1610, 1610,51

ahora vamos a ver un código que permite resolver el problema planteado, para diferentes valores de h = 1, h = 0,5, h = 0,1, h = 0,001 donde 1. yo : viene dada por la condición inicial yo = yo = y y((to) to)




12 8

2. t array de nodos ti en el intervalo [to,tf ] 3. h es el tamaño de paso entre los nodos ti 4. y representa la solución aproximada en cada nodo ti Usamos el siguiente código para resolver nuestra EDO , con el método de Euler function [t,y]=edoeuler(to,tf,yo,h) [t,y]=edoeuler(to,tf,yo,h) % méto método do de eule euler r t(1)=to; y(1)=yo; i=1; while t(i) < tf , t(i+1)=to+i*h; y(i+1)=y(i)+h*fty(t(i),y(i)); i=i+1; end

Para tener una función disponible, para llamar function f=fty(t,y) f=y/10;

Escriba el siguiente código para manejar las funciones escritas anteriormente to=0; tf=5 yo=1000; h=1; [t1,y1]=edoeuler(0,5,1000,h); plot(t1,y1,’r’) hold hold on h=0.5; [t2,y2]=edoeuler(0,5,1000,h); plot(t2,y2,’b’) h=0.1; [t3,y3]=edoeuler(0,5,1000,h); plot(t3,y3,’g’) h=0.001; [t4,y4]=edoeuler(0,5,1000,h); plot(t4,y4,’y’) % soluci solución ón exacta exacta x=linspace(to,tf,1000);




12 9

y=1000*exp(x/10); plot(x,y,’m’) legend(’h=1’,’h=0.5’,’h=0.1’, legend(’h=1’,’h=0.5’,’h=0.1’,’h=0.001’,’solu ’h=0.001’,’solución ción exacta’)

Ejemplo 8.2.

Resolver le problema de valor inicial PVI.

 



y =

200 sen(10/t sen(10/t)) t3

1000 t4

− cos(10/t cos(10/t)) y(1, (1,0610329) = 0, 0, t ∈ [0;3] tenga en cuenta la solución exacta y (t) = sin( sin( ) − −

100 t2

10 t

100 t

Solución Vamos a ver un código que permite resolver el problema planteado, para diferentes valores de h = 0,1, donde Escribimos una nueva función para ser manejada con nuestro programa edoeuler function f=fty(t,y) f=-200./t.^3.*sin(10./t)-1000 f=-200./t.^3.*sin(10./t)-1000./t^4.*cos(10./ ./t^4.*cos(10./t); t);

Para correr el programa digitamos lo siguiente to=1.0610329; tf=3; yo=0; h=0.1; %solución %solución aprox euler euler [t1,y1]=edoeuler(to,tf,yo,h);




13 0

clf plot(t1,y1,’r’) hold hold on % soluci solución ón exacta exacta t=linspace(to,tf,1000); y=(100./t.^2).*sin(10./t); plot(t,y,’g’) legend(’h= legend(’h=0.1, 0.1,.. .. sol aprox. aprox. euler’,’s euler’,’soluci olución ón exacta’) exacta’)

8.1.1. 8.1.1.

Ejerci Ejercicio cio,, mecán mecánica ica de la la fract fractura ura

El crecimiento de una grieta en el borde de una placa por ciclo de esfuerzos viene dado por la ecuación de Paris (7.4). Una forma de estimar el crecimiento de una grieta con el número de ciclos diferente a la integración numérica del ejemplo ( 7.2.4) es resolviendo la ecuación (7.4) por el método de Euler, en este caso se tiene la ecuación diferencial discretizada an+1 an ∆N

−

m

m

= A(∆ A (∆σ σ ) ak

2



1,99

− 0,41

 an 2,5

2

+ 18, 18,7

− 38, 38,48

 an 2,5

3

+ 53, 53,85

  an 2,5

Y despejando an+1 se tiene el valor del tamaño de la grieta al siguiente ciclo de carga (8.2) m

m

an+1 = a n +∆N +∆N A(∆σ (∆σ ) ak

2



1,99

− 0,41

 an 2,5

2

+ 18, 18,7

38,48 − 38,

 an 2,5

3

+ 53, 53,85

4 m

 

Una forma de implementar la solución es con un algoritmo de MATLAB que se detenga al alcanzar el valor de la grieta crítica. Con los valores del ejemplo ( 7.2.4) se tiene el algoritmo en MATLAB que implementa la ecuación (8.2 ) Mg. Marlo Carranza Purca Purca


an 2,5

4 m


13 1

% crec crecim imie ient nto o de la grie grieta ta por por el méto método do de Eule Euler r function [N,a]=CrecimientoGrieta [N,a]=CrecimientoGrieta A=6.6*10^{-9}; w=2.5; w=2.5; %[in] %[in] ds=17.78; ds=17.78; %[ksi] %[ksi] m=2.25; N(1)=0; a(1)=0.25 a(1)=0.25; ; %[in] dN=1; n=1; while while a<0.62 a<0.6200 00 %[in] %[in] Grieta Grieta Crític Crítica a a(n+1)=a(n)+dN*A*ds^m*a(n)^( a(n+1)=a(n)+dN*A*ds^m*a(n)^(m/2)*(1.99m/2)*(1.990.41*(a(n)/w)+18.7*(a(n)/w)^ 0.41*(a(n)/w)+18.7*(a(n)/w)^2-38.48*(a(n)/w 2-38.48*(a(n)/w)^3+53.85*(a(n) )^3+53.85*(a(n)/w)^4)^m; /w)^4)^m; N(n+1)=N(n)+dN; n=n+1; end

El algoritmo así implementado dice que el número de ciclos para alcanzar la grieta crítica y la falla es N = = 37102 ciclos. Compárese este valor con el valor hallado por integración numérica en el ejemplo (7.2.4 ). La gráfica del crecimiento de la grieta por el método de Euler se muestra en la figura , compárese con la figura (8) del crecimiento de la grieta hallado por integración numérica en el ejemplo (7.2.4 )

Figura 8.1: Crecimiento de la grieta con el número de ciclos hasta la falla



Bibliografía

[1]

Chains Chainsk kaia aia,, Liudm Liudmila ila y Doig, Doig, Eliza Elizabeth beth.. Elementos de Análisis Numéricos: algoritmos y aplicaciones . Fondo Editorial de la PUCP. 1999.

[2] [2]

Chap Chapra ra,, Stev Steven en.. Métodos para Ingenieros . Editorial McGraw-Hill.

[3]

Dougla Douglas, s, J. Burden Burden Richad Richad.. Métodos Numéricos . Tercera edición. THOMSON Editores. 2004. Madrid España.

[4]

Maron, Maron, Me Melvi lvin; n; Robert Robert J. Lopez. Lopez. Un enfoque Práctico . Editorial CECSA.

[5]

Nakam Nak amura ura,, Shoic Shoichir hiro. o. Métodos Numéricos aplicados con Software . Editorial Prentice Hall Hispanoamericana.

[6]

Nieve Nievess Doming Dominguez uez.. Métodos Numéricos aplicados a la ingeniería . Editorial CECSA.

[7]

Scheid, Scheid, Francis, rancis, Rosa Rosa Elena Di Const Constanzo. anzo. Métodos Numéricos . Colección Schaum. Editorial Mc Graw Hill.


249051800-metodos-numericos.pdf

Recommend Documents