Análisis Numérico De Integrales Y Ecuaciones Diferenciales - Carlos E. Neuman - 3ra Edición.pdf

SIGUENOS EN:

LIBROS UNIVERISTARIOS Y SOLUCIONARIOS DE MUCHOS DE ESTOS LIBROS GRATIS EN DESCARGA DIRECTA VISITANOS PARA DESARGALOS GRATIS.

ISSN 1666-3845

e-Commentarii Aplicandæ Mathematicæ Editor: Carlos E. Neuman LN: Lecture Notes

LN1: Neuman, C. E.: An´ alisis Num´erico de Integrales y Ecuaciones Diferenciales, third edition, 2001.

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

Temas de An´alisis Num´erico y Programaci´on

Vol´ umenes de la serie: I : Cardona, A., y Neuman, C.E.: Temas de Programaci´ on I, 1996, (segunda edici´ on en preparaci´on, 1997) II : Neuman, C.E. (editor): Matem´ atica asistida por computadora, 1996 III : Neuman, C.E.: An´ alisis Num´erico de Integrales y Ecuaciones Diferenciales, 1997 IV : Neuman, C.E.: Algoritmos Geom´etricos y Gr´ aficos en 3D. Laboratorios de Matem´ atica, segunda edici´on, 1999. V : Neuman, C.E., Vilchez, A.G.: Algoritmos Geom´etricos y Gr´ aficos en 3D. Laboratorios de Matem´ atica, tercera edici´on (en preparaci´on), 2001.

i

e-Commentarii Aplicandæ Mathematicæ , LN: Lecture Notes Volume 1 ISSN 1666-3845

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

¶

Tercera edici´ on

Carlos E. Neuman Departamento de Matem´atica (FIQ) Universidad Nacional del Litoral

Santa Fe 2001 ∞ ¶ Trabajo realizado con fondos provenientes de la Universidad Nacional del Litoral (UNL), a trav´ es de la programaci´ on Curso de Acci´ on para la Investigaci´ on y el Desarrollo (CAI+D), Secretar´ıa de Ciencia y T´ ecnica de la UNL

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales Carlos E. Neuman ´tica, FIQ, UNL, Santiago del Estero 2829, Departamento de Matema 3000 Santa Fe, Argentina E-mail address: [email protected] URL: http://www.cen.com.ar/eCAM/

Key words and phrases. Integraci´on num´erica, Resoluci´on num´erica de ecuaciones diferenciales, M´etodos de un paso, Integraci´on adaptativa

Trabajo realizado con fondos provenientes de la Universidad Nacional del Litoral (UNL), a trav´es de la programaci´on Curso de Acci´ on para la Investigaci´ on y el Desarrollo (CAI+D), Secretar´ıa de Ciencia y T´ecnica de la UNL. Resumen. En este texto se tratan temas de integraci´ on num´ erica de funciones y de ecuaciones diferenciales, teoremas de existencia y unicidad de soluciones, convergencia, estabilidad y consistencia de soluciones aproximadas de ecuaciones diferenciales y algunos complementos te´ oricos. Se presentan numerosos ejemplos y ejercicios y dos laboratorios de integraci´ on adaptativa.

ISSN 1666-3845

´Indice General Prefacio

xiii

Prefacio a la segunda edici´ on

xvii

´Indice de Tablas

1

´Indice de Figuras

3

Cap´ıtulo 1. Integraci´ on 1.1. Introducci´ on 1.2. La integral de Riemann–Stieltjes 1.3. La f´ ormula de suma de Euler 1.4. M´etodos b´ asicos de integraci´on num´erica 1.5. M´etodos de extrapolaci´on 1.6. Cuadratura de Gauß 1.7. M´etodos de Montecarlo 1.8. M´etodos adaptativos 1.9. Ejercicios

1 1 3 4 7 11 12 13 16 17

Cap´ıtulo 2. Ecuaciones diferenciales 2.1. Introducci´ on 2.2. M´etodos en diferencias de un paso 2.3. Ejercicios

21 21 23 33

Ap´endice A. Bibliograf´ıa y referencias A.1. Textos b´ asicos A.2. Bibliograf´ıa complementaria

35 35 35

Ap´endice B. Coeficientes indeterminados con el Mathematica B.1. El m´etodo de Simpson B.2. El m´etodo de Newton B.3. El m´etodo de cuadratura de Gauß de tres puntos B.4. Los coeficientes de la cuadratura de Gauß de cinco puntos

37 37 38 40 41

Ap´endice C. Cuadernos del Mathematica y archivos script de Matlab C.1. La funci´on mihumps C.2. Extrapolaci´ on repetida de Richardson C.3. M´etodos de Montecarlo C.4. Integraci´ on de Simpson adaptativa C.5. El m´etodo cl´ asico de Runge-Kutta

47 47 48 50 53 55

Ap´endice D.

59

La integral de Riemann–Stieltjes vii

´INDICE GENERAL

viii

D.1. D.2. D.3. D.4. D.5.

Funciones de variaci´ on acotada La integral F´ ormula de la suma de Euler Pero, ¿Hay funciones integrables? Ejercicios

59 61 65 65 66

Ap´endice E. Polinomios y n´ umeros de Bernoulli E.1. Polinomios de Bernoulli con el Mathematica E.2. N´ umeros de Bernoulli E.3. Integraci´ on por partes

67 67 73 74

Ap´endice F. Laboratorio de Matem´atica: Integraci´on Adaptativa F.1. Objetivos F.2. Antes del laboratorio F.3. Laboratorio F.4. Aplicaciones y desarrollos avanzados F.5. Referencias

75 75 75 76 77 77

Ap´endice G. Laboratorio de Matem´atica: Soluci´on adaptativa de ODEs G.1. Objetivos G.2. Antes del laboratorio G.3. Laboratorio G.4. Aplicaciones y desarrollos avanzados G.5. Referencias

79 79 79 79 79 80

Ap´endice H. Existencia y unicidad de soluciones de ODEs H.1. Referencias

81 85

´INDICE GENERAL

Para S., N., y F.

ix

´INDICE GENERAL

xi

Our ancestors were eager to understand the worlds but had not quite stumbled upon the method. They imagined a small, quaint, tidy universe in which the dominant forces were gods like Anu, Ea, and Shamash. In that universe humans played an important if not central role. We were intimately bound up with the rest of nature. The treatment of toothache with second-rate beer was tied to the deepest cosmological mysteries. Today we have discovered a powerful and elegant way to understand the universe, a method called science; it has revealed to us a universe so ancient and so vast that human affairs seem at first sight to be of little consequence. Sagan, C.: Cosmos, Ballantine, New York, 1992.

Prefacio Este texto surgi´o de la necesidad de presentar una visi´on unificada de los m´etodos elementales de integraci´on num´erica de funciones y de algunos de los que son u ´tiles para las ecuaciones diferenciales ordinarias. ♣ Los p´ arrafos que se colocan entre tr´eboles pueden ser omitidos en las primeras lecturas puesto que no son esenciales para la continuidad del texto. ♣ Seguimos en t´erminos generales las ideas planteadas en Linz∗ (1988). Para este autor el t´ermino matem´ atica computacional esta asociado con el espectro amplio de actividades relacionadas con la soluci´on aproximada de problemas cient´ıficos expresados en t´erminos de modelos matem´aticos, que, en su forma t´ıpica, est´an constituidos por ecuaciones diferenciales e integrales de las que, en general, no se conoce soluci´ on en forma cerrada. Para resolverlos las ecuaciones son convertidas, mediante discretizaci´ on, en un conjunto finito de ecuaciones m´as simples que puedan ser resueltas por m´etodos algebraicos. Se puede distinguir la metodolog´ıa num´erica, que estudia m´etodos para discretizar los operadores diferenciales e integrales y c´omo resolver los sistemas finitos resultantes, del an´ alisis num´erico que involucra el estudio riguroso de los algoritmos creados por la metodolog´ıa. La meta primaria del an´ alisis es describir la relaci´on entre la soluci´on exacta de la ecuaci´ on original, y la aproximada obtenida a partir de la versi´on discreta. El an´ alisis num´erico ha resultado muy u ´til para complementar el proceso de creaci´ on de m´etodos num´ericos generado por las aplicaciones ingenieriles. Por ejemplo las t´ecnicas de relajaci´ on o de elementos finitos fueron creadas por ingenieros basados en la intuici´ on f´ısica. El an´alisis num´erico posterior ha resultado crucial para la consolidaci´ on de estos m´etodos, la extensi´on de su aplicabilidad, y el estudio de sus propiedades de convergencia o estabilidad. En muchos casos a partir de objetivos del An´alisis Num´erico se han obtenido nuevos resultados en An´ alisis Funcional, en cuyos t´erminos suelen plantearse los problemas. Seg´ un resume Linz (1988): En art´ıculos publicados en revistas como el SIAM Journal on Numerical Analysis o Numerische Mathematik se encuentran habitualmente t´erminos como espacios de Hilbert, clausura compacta, y convergencia d´ebil. Estos conceptos le resultan muy u ´tiles al te´orico y han facilitado el establecimiento de una teor´ıa abarcativa y coherente de la soluci´on aproximada de las ecuaciones de operadores. ♣ Ilustraremos esto con el ejemplo de la ecuaci´on del operador lineal Lx = y ∗

Las referencias se enumeran en el ap´ endice A, p´ agina 35 xiii

(0.1)

xiv

PREFACIO

donde L : X → Y es un operador lineal entre los espacios normados X e Y . El miembro derecho, y, es un dato y la ecuaci´on debe ser resuelta para la inc´ognita x. Suponemos que L tiene una inversa acotada en Y . Un ejemplo muy simple de esto es la ecuaci´ on d x(t) = y(t) = f (t, x(t)) (0.2) dt con la condici´ on inicial x(t0 ) = x0 (0.3) En el proceso de discretizaci´on la ecuaci´on (0.1) es reemplazada por una sucesi´ on param´etrica de problemas Ln xn = yn

(0.4)

donde Ln es un operador en ciertos espacios n–dimensionales Xn e Yn . El s´ımbolo n se denomina par´ ametro de discretizaci´ on y mide el grado en el que el operador discreto Ln representa al operador L. Como Ln es efectivamente una matriz de n×n, es posible, en principio, resolver (0.4) en forma algebraica y obtener la soluci´on aproximada xn . El objetivo fundamental del an´alisis posterior es conocer la relaci´on entre la soluci´ on verdadera x y su aproximaci´on xn . En particular, desearemos demostrar la convergencia. Esto significa que, cuando se incrementa n, la soluci´on aproximada deber´ a acercarse m´ as y m´ as a la soluci´on verdadera, en el sentido que lim kx − xn k = 0

n→∞

(0.5)

Usualmente, el primer paso del an´alisis es definir el error de consistencia rn (x) = Lx − Ln x

(0.6)

De (0.1) y (0.4) se obtiene directamente que Ln (x − xn ) = y − yn − rn (x)

(0.7)

y obtenemos una cota para el error de la soluci´on aproximada kx − xn k ≤ kL−1 n k{ky − yn k + krn k}

(0.8)

En la mayor´ıa de los casos es relativamente elemental mostrar que lim ky − yn k = 0

n→∞

(0.9)

y que, bajo condiciones bien definidas lim krn k = 0

n→∞

(0.10)

Se necesita una condici´ on adicional, la condici´on de estabilidad lim kL−1 n k≤K <∞

n→∞

(0.11)

Si reemplazamos en la ecuaci´on (0.8), obtenemos el teorema central del an´alisis num´erico: estabilidad y consistencia implican convergencia. Se tiene un gran n´ umero de problemas no resueltos por causa de las dificultades t´ecnicas de las que el problema principal es verificar la condici´on de estabilidad (0.11). Hay otras cuestiones, por ejemplo, cu´an dificil es resolver (0.4) o cu´an velozmente converge xn a x, pero, normalmente estas son menos complejas que la estabilidad.

PREFACIO

xv

Supongamos ahora que hemos determinado que nuestro m´etodo num´erico es estable y convergente. ¿Qu´e nos dice esto? Nos expresa que, en alg´ un sentido asint´ otico, el m´etodo funciona. A menudo podemos decir m´as. En muchos casos el error de aproximaci´on puede ser acotado en forma m´ as expl´ıcita como kx − xn k ≤ n−p kL−1 n kη(x)

(0.12)

donde η es un funcional de la soluci´on desconocida x El exponente p se denomina orden de convergencia y nos da informaci´on respecto de cu´an exacto es el m´etodo. Para un dado nivel de esfuerzo computacional, un m´etodo con un orden de convergencia mayor tender´ a a dar mejores resultados que uno con un orden menor. La cota que nos da la expresi´on 0.12 tiene una utilidad limitada puesto que depende de la inc´ ognita x. Se puede sortear esa dificultad realizando un an´alisis a posteriori que utiliza la soluci´on calculada xn en lugar de la x. A partir de las ecuaciones 0.1 y 0.4 se puede obtener kx − xn k ≤ n−p kL−1 kkρn k

(0.13)

donde ρn es el residuo de la soluci´on calculada ρn (xn ) = Lxn − y

(0.14)

Es posible obtener cotas buenas para kρn k, en cambio es m´as dif´ıcil obtener cotas realistas para kL−1 k lo que se convierte en el elemento clave del que depende el an´ alisis a posteriori. ♣ Los problemas b´ asicos de la metodolog´ıa num´erica no son el establecimiento del m´etodo de discretizaci´ on, sino el manejo de los detalles involucrados en la representaci´ on de dominios, de los propios m´etodos de discretizaci´on, de los m´etodos de operaci´ on con los sistemas de ecuaciones, de las matrices ralas, y de los m´etodos para resolver los sistemas de ecuaciones con dimensiones muy altas. En el desarrollo de m´etodos para resolver estos problemas suele incluirse —cuando se trata de utilizar la metodolog´ıa en las aplicaciones— una dosis relativamente alta de lo que puede denominarse pragm´ atica num´erica la que, por tratarse de algo naturalmente impreciso, siempre puede fallar, pero que no deber´ıa ser totalmente ignorada por causa de ello. Los beneficios de la metodolog´ıa num´erica son tan considerables que se la utiliza a´ un en los casos en que no puede justificarsela rigurosamente. La idea fuerza no es rechazar los puntos de vista pragm´aticos sino desarrollar los m´as efectivos. En el primer cap´ıtulo se estudian algunos aspectos relevantes del c´alculo de integrales. Se desarrollan los m´etodos cl´asicos, los de extrapolaci´on y los de Montecarlo. Debido a que las extensiones a integrales m´ ultiples son relativamente directas, solamente se trata el problema de vol´ umenes e hipervol´ umenes utilizando el m´etodo de Montecarlo por la simple raz´on que, en muchos casos, es el u ´nico aplicable. En el segundo cap´ıtulo se introduce el problema de hallar soluciones aproximadas para problemas de valores iniciales de ecuaciones diferenciales ordinarias. Asimismo se deducen algunos m´etodos elementales de los del tipo de ‘un paso’. En los ap´endices se presentan ejemplos realizados con Matlab y Mathematica, se desarrolla una breve introducci´ on a la teor´ıa de la integral y se realizan c´alculos para obtener los polinomios de Bernoulli, que resultan necesarios para la interpretaci´on de los coeficientes en el desarrollo asint´ otico del error en los m´etodos trapezoidal y Simpson.

xvi

PREFACIO

En el u ´ltimo de los ap´endices se incluye un Laboratorio de Matem´atica en el tema de la cuadratura adaptativa. He elegido este por la simplicidad del algoritmo recursivo. Sin embargo, cualquiera de los restantes temas del texto, por su importancia matem´ atica, su riqueza y su atractivo, merece la elaboraci´on de un Laboratorio propio. El objetivo de estos laboratorios es aspirar a la recreaci´on viva de conceptos y temas que enriquecen el portafolios del pregraduado en Matem´atica. Los libros citados en la secci´on de bibliograf´ıa y el conjunto de referencias en ellos citados facilitan al lector interesado el estudio de cualquiera de los temas en que desee especializarse. Este texto, que se basa en su totalidad en ellos, y que comparte la selecci´ on de temas de inter´es que tratan, no aspira a colocarse en un plano de igualdad, sino, humildemente, enfatizar algunos aspectos y no otros. Debo agradecer a los sucesivos grupos de alumnos que han enriquecido el curso de C´ alculo Num´erico con filosas preguntas y soluciones ingeniosas a los problemas propuestos. Asimismo a mis colegas por la discusi´on de temas y el desarrollo de m´etodos aptos para la ense˜ nanza. En cambio, de los errores que pudiera contener, me hago u ´nico responsable. Carlos E. Neuman Santa Fe, 29 de abril de 1997

Prefacio a la segunda edici´ on Esta nueva edici´ on se desarrolla con varias correcciones y agregados con el fin de adaptar el texto a los cursos de C´alculo Num´erico I y II en la modalidad ‘a distancia’. Carlos E. Neuman Santa Fe, 2 de diciembre de 2001

xvii

xviii

´ PREFACIO A LA SEGUNDA EDICION

Mi padre fue un jud´ıo renegado, o un ateo, o un rebelde. Lo cierto es que siempre le hab´ıa encantado, sentado en el u ´ltimo banco, como un alumno malo, portarse mal en las bodas de la sinagoga, mirar con sorna al rabino del barrio, meterse con “la rama andaluza chupacirio de la familia” y darle vuelta la cara a los Cohen “jud´ıos ortodoxos y obtusos” de enfrente de casa. Pero mi padre le´ıa a escondidas el Antiguo Testamento, una tarde en que me encontr´o un libro de misa que me hab´ıan regalado en una fiesta de primera comuni´on me mir´o confundido y al fin — como esto pueda entenderse, como esto pueda entenderlo yo—, se sent´ıa jud´ıo. Mi padre fue un hombre y, por tanto, fue tantas cosas que no me atrevo a enumerar, porque un hombre —ahora lo s´e— contiene tantos hombres como vidas tiene un gato y tanto m´ as. En principio, mi padre fue un coleccionista. Durante a˜ nos me pas´e tardes enteras examinando una por una su colecci´ on de cajas y mu˜ necas musicales. Las mu˜ necas eran las m´ as dif´ıciles de maniobrar porque no hab´ıa que despeinarlas ni descuidar los accesorios del atuendo. La bailarina flamenca llevaba un par de casta˜ nuelas de madera —sevillanas leg´ıtimas, aclaraba mi padre—, unas miniaturas que estaban sabiamente encajadas entre las manos. La joven con mant´on de Manila sosten´ıa un abanico de encaje y piedras, la bailarina cl´asica llevaba un tut´ u que al fin acab´o manchado de caf´e con leche. Me gustaba ponerlas a todas en funcionamiento, a la vez, todas girando al son de una m´ usica extra˜ na y dislocada. Pero sobre todo lo que me importaba era desarmar. Me importaba entender el funcionamiento. Era el acto ritual ?–como ritual era la acci´ on de la anciana de la tienda cuando aparec´ıa con una caja llena de sacapuntas–?, ese acto lento, siempre el mismo, el de levantar con cuidado una tapa de terciopelo, esa operaci´on que se daba inequ´ıvoca en las cajas que realmente ten´ıan forma de caja: un peque˜ no joyero, una cigarrera en forma de dado, una caja cl´ asica dorada y de marfil. Y una vez abiertas, lentamente, como destripadas frente a mis ojos, les daba cuerda a todas a la vez para que sonase la m´ usica disonante y finalmente triste, como la de todas las mu˜ necas girando solitarias, porque lo mejor no suced´ıa al fin sino poco antes, cada vez que abr´ıa una nueva caja, cada vez que expectante, solemne, levantaba la tapa, levantaba luego el compartimento rojo, despacio, con un dedo, y otra vez la desilusi´ on. Buscaba una caja en especial, una caja que tuviera, por debajo de ese compartimento, un mecanismo distinto al del rodillo. Todas las cajas que. . . Lilian Neuman: Levantar Ciudades, Ed. Destino, Barcelona, 1999.

´Indice de Tablas 1

Extrapolaci´ on de Richardson (M´etodo de Romberg)

1

Integraci´ on de Romberg para la funci´on mihumps (ver tabla 1). La integral (con el Mathematica) resulta Q = 0.2985832539549867 50 Integraci´ on de Romberg para la funci´on mihumps (continuaci´on, ver tablas 1 y 1). La integral (con el Mathematica) resulta Q = 0.2985832539549867 50

2

1

12

´Indice de Figuras 1 2 3 4 5

Campo de direcciones de f (x, y) = −xy El rect´ angulo S est´a contenido en el R Primer m´etodo de Runge Segunda f´ ormula de Runge El m´etodo de Runge-Kutta

21 22 26 28 29

1 2

47

3

La funci´ on mihumps Histograma de distribuci´on de los valores estimados de la integral del ejemplo C.3.1 Integraci´ on por Runge-Kutta cl´asico de x0 = −tx en [0, 7]

52 57

1 2 3 4 5

El El El El El

B1 (x) B2 (x) B6 (x) B10 (x) B40 (x)

69 70 71 72 73

1

Puntos de evaluaci´on de quad para la funci´on f (x) en [0, 2]

76

polinomio polinomio polinomio polinomio polinomio

de de de de de

Bernoulli Bernoulli Bernoulli Bernoulli Bernoulli

3

CAP´ıTULO 1

Integraci´ on 1.1. Introducci´ on El objetivo del presente cap´ıtulo es estudiar m´etodos num´ericos para calcular integrales definidas de la forma Z b f (x) dx (1.1) a

Los m´etodos m´ as simples —que se originan en la propia definici´on de la integral de Riemann,— se establecen como sigue. Se elige una partici´on regular de paso h en el intervalo [a, b] de modo de obtener N subintervalos (es decir que N h = (b − a). Se toma un valor conveniente de la funci´on f en cada subintervalo y se calcula la suma N X h fi (1.2) i=1

donde fi es el valor elegido de la funci´on f en el intervalo [xi−1 , xi ]. En el caso de escoger fi = f (xi−1 ) se obtiene una “regla del rect´angulo” Z

b

f (x) dx ' h a

N X

f (xi−1 )

(1.3)

i=1

si, en cambio, se toma fi = f (xi ) se obtiene otra regla del rect´angulo Z

b

f (x) dx ' h a

N X

f (xi )

(1.4)

i=1

en realidad ‘la’ regla del rect´ angulo es m´as precisa si se elige fi = f ((xi−1 + xi )/2), en este caso se obtiene la regla del rect´angulo R(h) (o ‘midpoint’) b

Z

f (x) dx ' R(h) = h a

N X i=1

f(

xi−1 + xi ) 2

(1.5)

con fi = (f (xi−1 ) + f (xi ))/2, en cambio, se obtiene la regla del trapecio T (h) Z

b

f (x) dx ' T (h) = h a

N X f (xi−1 ) + f (xi ) i=1

2

(1.6)

es decir T (h) = h(

f (a) + 2

N −1 X i=1 1

f (xi ) +

f (b) ) 2

(1.7)

En algunos textos este m´ etodo se denomina trapezoidal compuesto, pero la distinci´ on es poco relevante

´ 1. INTEGRACION

2

y, si el n´ umero N de subintervalos es par, se pueden agrupar los subintervalos de a dos contiguos y, en cada par utilizar f2i+1 = (f (x2i ) + 4f (x2i+1 ) + f (x2i+2 ))/6, con lo que se obtiene la “regla de Simpson” S(h) Z a

La regla de Simpson es exacta para polinomios de orden cuatro

b

N

2 −1 2h X f (x) dx ' S(h) = (f (x2j ) + 4f (x2j+1 ) + f (x2j+2 )) 6 j=0

es decir 2h (f (a) + 4f (x1 ) + 2f (x2 ) + · · · + 4f (xN −1 ) + f (xN )) 6 En todos estos casos se tiene la expresi´on Z b N X f (x) dx ' A ci f (xi ) S(h) =

a

El operador ‘.*’ facilita la integraci´ on num´ erica

(1.8)

(1.9)

(1.10)

i=0

donde [Rect´ angulo]A = 2h y c = (0, 1, 0, 1, 0, 1, . . . , 0, 1, 0) [Trapezoidal]A = h/2 y c = (1, 2, 2, . . . , 2, 1) [Simpson] A = 2h/6 y c = (1, 4, 2, 4, . . . , 2, 4, 1) ´ n 1.1.1. Las expresiones anteriores sugieren un m´etodo muy simple Observacio para programar estos m´etodos elementales en Matlab

La integral de RiemannStieltjes generaliza la integral de Riemann

De la inspecci´ on de estos m´etodos surgen varias ideas interesantes (1) Se tiene una ´ıntima relaci´on entre integrales y sumas y las primeras se aproximan por las segundas. (2) Es necesario desarrollar m´etodos para estudiar los errores de aproximaci´on de estas reglas de integraci´on (es simple establecerlas pero no es tan claro c´ omo estimar los respectivos errores) (3) Se deber´ıan explorar m´etodos aplicables a particiones que no sean regulares Para la primera utilizaremos una teor´ıa de la integral que permite unificar las metodolog´ıas para tratar con integrales y sumas. Citamos textualmente a Tom Apostol en la introducci´ on al cap´ıtulo 7 (p´ag. 169) de su texto (Apostol, 1976) . . . se estudia el proceso de integraci´ on con cierto detalle. En realidad consideramos un concepto m´ as general que el de Riemann: a saber, la integral de Riemann-Stieltjes, que involucra dos funciones f y α. El s´ımbolo que utilizamos para designar tales integrales es Rb f (x)dα(x), o alguno similar, y la integral de Riemann se obtiene a como caso particular cuando α(x) = x. Cuando α tiene derivada conRb tinua, la definici´ on es tal que la integral de Stieltjes a f (x)dα(x) se Rb convierte en la intgral de Riemann a f (x)α0 (x) dx. Sin embargo, la integral de Stieltjes tiene sentido en el caso en que α no es diferenciable e incluso cuando no es continua. De hecho, es al tratar con funciones discontinuas α cuando se hace patente la importancia de la integral de Stieltjes. Eligiendo adecuadamente una funci´ on discontinua α, una suma finita o infinita puede expresarse como una integral de Stieltjes, y la sumaci´ on y la integral de Riemann ordinaria son casos especiales de este proceso m´ as general. Los problemas f´ısicos que consideran la distribuci´ on de masas que son en parte discretas

1.2. LA INTEGRAL DE RIEMANN–STIELTJES

3

y en parte continuas pueden ser abordados utilizando la integral de Stieltjes. En la teor´ıa matem´ atica de la Probabilidad esta integral es un instrumento muy u ´til que hace posible la consideraci´ on simult´ anea de variables aleatorias continuas y discretas.

En muchos casos se utilizan integrales para calcular sumas que pueden ser complicadas de obtener por otros m´etodos, en raz´on de ello las interrelaciones entre unas y otras es mucho mayor de lo que puede suponerse en una primera aproximaci´ on. En las secciones 1.4 y 1.6 estudiamos herramientas para obtener diversos m´etodos de integraci´ on y para estimar los errores de las aproximaciones conseguidas. El esquema b´ asico es aproximar la integral por una suma del tipo N X

wi f (xi )

(1.11)

i=1

sobre una partici´ on del intervalo de integraci´on, y elegir en forma conveniente los pesos wi y los puntos xi . Cuando aceptamos que las particiones no sean regulares surge naturalmente la idea de refinar la partici´ on en los subintervalos donde los errores sean muy grandes y la idea opuesta de agrupar intervalos donde los errores sean excesivamente peque˜ nos. Esto nos lleva a los denominados m´etodos adaptativos de integraci´on que estudiaremos en la secci´ on 1.8 1.2. La integral de Riemann–Stieltjes 1.2.1. Definici´ on y propiedades. ´ n 1.2.1. Sea P = {x0 , x1 , . . . , xn } una partici´on del intervalo [a, b] y Definicio sea, para cada k, ξk un punto del intervalo [xk−1 , xk ]. Una suma de la forma n X S(P, f, g) = f (ξk )∆gk (1.12) k=1

(donde ∆gk = g(xk ) − g(xk−1 )) se llama suma de Riemann-Stieltjes de f respecto de g. Diremos que f es integrable respecto de g en [a, b] si existe un n´ umero I que satisface la siguiente propiedad: para cada ε > 0, existe una partici´on Pε de [a, b] tal que, para cada partici´ on P m´as fina que Pε y para cada elecci´on de los puntos ξk del intervalo [xk−1 , xk ], k = 1, . . . , n, se tiene |S(P, f, g) − I| < ε

(1.13)

´ n 1.2.1. Si un tal n´ Observacio umero I existe, entonces es u ´nico y se representa por la expresi´ on Z b f (x)dg(x) (1.14) a

Ejemplo 1.2.1. Si g(x) = x entonces la integral se reduce a una integral de Riemann Ejemplo 1.2.2. Si g(x) es de clase C 1 en [a, b] entonces la integral resulta Z b f (x)g 0 (x) dx (1.15) a

Nota: Este resultado se demuestra en el teorema D.2.5 en la p´agina 63

La integral aproximada es un promedio pesado de los valores de la funci´ on

´ 1. INTEGRACION

4

Ejemplo 1.2.3. Si f es continua y g(x) = bxc entonces la integral resulta X f (i) (1.16) i∈Z a
y si g(x) = dxe entonces X

f (i)

(1.17)

i∈Z a≤i
Nota: Este resultado se demuestra en el teorema D.2.7 en la p´agina 65. Notar la diferencia entre los l´ımites de suma. 1.2.2. Integraci´ on por partes. Teorema 1.2.2. [F´ ormula de integraci´ on por partes] Si f es integrable respecto de g en el intervalo [a, b], en el sentido de Riemann-Stieltjes, entonces g es integrable respecto de f en el intervalo [a, b], en el sentido de Riemann-Stieltjes, y se tiene Z b Z b f (x)dg(x) + g(x)df (x) = f (b)g(b) − f (a)g(a) (1.18) a

a

´ n: Ver el teorema D.2.4 en la p´agina 62 Demostracio En el ap´endice D en p´ agina 59 se desarrollan aspectos elementales de la teor´ıa de integraci´ on de Riemann-Stieltjes siguiendo el texto de Apostol (1976) 1.3. La f´ ormula de suma de Euler Sea {x0 , x1 , . . . , xn }, donde x0 = a, x1 = a + h, . . . , xn = a + nh = b, una partici´ on regular del intervalo [a, b]. La suma trapezoidal para la integral Z b f (x) dx (1.19) a

se denota T (h) y resulta ser   1 1 T (h) = h f0 + f1 + f2 + f3 + . . . + fn−1 + fn 2 2

(1.20)

o bien n−1

T (h) = h

X 1 1 f0 + fi + fn 2 2 i=1

! (1.21)

Desear´ıamos demostrar el siguiente teorema Teorema 1.3.1. Sea T (h) la suma trapezoidal 1.21. Entonces Z b h2 T (h) = f (x) dx + (f 0 (b) − f 0 (a)) 12 a 4 h − (f (3) (b) − f (3) (a)) 720 h6 + (f (5) (b) − f (5) (a)) 30240 + . . . + c2r h2r (f (2r−1) (b) − f (2r−1) (a)) + O(h2r+2 )

(1.22) (1.23) (1.24) (1.25)

´ 1.3. LA FORMULA DE SUMA DE EULER

5

Donde los coeficientes {c2r } tienen la funci´on generatriz h eh + 1 (1.26) 2 eh − 1 (para la interpretaci´ on detallada de los coeficientes ver el ap´endice E en la p´ agina 67) 1 + c2 h2 + c4 h4 + c6 h6 + . . . =

1.3.1. El s´ımbolo ‘O’. La notaci´on ‘O grande’ es u ´til para representar y calcular en forma exacta con magnitudes aproximadas. En general la notaci´on O(f (n)) se usa —si f (n) es una funci´on definida sobre los n´ umeros naturales— con el fin de representar un n´ umero xn con el significado de que existe una constante positiva M tal que el n´ umero xn representado por O(f (n)) satisface la condici´on |xn | ≤ M |f (n)|,

∀n ≥ n0

(1.27)

No damos expl´ıcitamente —en general— las constantes M y n0 , las que suelen ser distintas en cada aparici´ on de O Ejemplo 1.3.1. Sabemos que n X 1 1 1 1 1 i2 = n(n + )(n + 1) = n3 + n2 + n 3 2 3 2 6 i=1

(1.28)

en consecuencia n X i=1 n X i=1

i2 = O(n3 ) i2 =

1 3 n + O(n2 ) 3

(1.29) (1.30)

La ecuaci´ on 1.30 es m´ as fuerte que la ecuaci´on 1.29 (Justificar la validez de estas ecuaciones y la relaci´ on entre ambas). Las siguientes son operaciones v´alidas con la notaci´on ‘O’ f (n) = O(f (n)) cO(f (n)) = O(f (n)), si c es una constante O(f (n)) + O(f (n)) = O(f (n)) O(O(f (n))) = O(f (n)) O(f (n))O(g(n)) = O(f (n)g(n)) O(f (n)g(n)) = f (n)O(g(n))

(1.31) (1.32) (1.33) (1.34) (1.35) (1.36)

Nota 1.3.1. La notaci´ on ‘ O ’ suele extenderse a la variable real x. Por ejemplo se escribe 1 1 ex = 1 + x + x2 + x3 + O(x4 ) (1.37) 2! 3! Una cierta dosis de observaci´on permite resolver la siguiente paradoja. Ejercicio 1.3.2. ¿Qu´e est´a m´al en el siguiente razonamiento? “Como n = O(n), 2n = O(n), . . . , se tiene X X kn = O(n) = 0(n2 )” (1.38) 1≤k≤n

Notar que es f´ acil demostrar que

1≤k≤n

P

1≤k≤n

kn = O(n3 ).

´ 1. INTEGRACION

6

Ejercicio 1.3.3. Demostrar que si se permite a x tomar valores arbitrariamente grandes, entonces, para toda potencia m, ex 6= O(xm ) Este tipo de desarrollos en potencias de un operador suele utilizarse en los cursos elementales de C´ alculo para la f´ ormula del polinomio de Taylor en dos variables independientes

1.3.2. lgebra de operadores. Si D representa el operador d/dx se tiene h3 h2 f (x + h) = f (x) + hf 0 (x) + f 00 (x) + f 000 (x) + . . . 2! 3!   (hD)2 (hD)3 = 1 + hD + + + . . . f (x) 2! 3! = ehD f (x)

(1.39) (1.40) (1.41)

donde el operador ehD debe entenderse definido por la serie, para la que es f´acil verificar que es normalmente convergente. En consecuencia f (x + 2h) = ehD f (x + h) = e2hD f (x)

(1.42)

y as´ı sucesivamente. Utilizamos esta representaci´on de la traslaci´on en la siguiente f´ ormula   1 1 T (h) = h f0 + f1 + f2 + . . . + fn−1 + fn (1.43) 2 2   1 1 =h + ehD + e2hD + . . . + e(n−1)hD + enhD f (x0 ) (1.44) 2 2   n−1 1 X jhD 1 nhD  = h − + e + e f (x0 ) (1.45) 2 j=0 2 luego, por la f´ormula de la suma geom´etrica,   1 enhD − 1 1 nhD T (h) = h − + hD + e f (x0 ) 2 e −1 2  
1 1 = h hD + enhD − 1 f (x0 ) e −1 2
h ehD + 1 nhD = e − 1 f (x0 ) 2 ehD − 1 pero ∞ t et + 1 X = c2j t2j 2 et − 1 j=0 para la que, con el Mathematica, se obtiene directamente: Series[(h/2)(Exp[h]+1)/(Exp[h]-1),{h,0,21}] 2 4 6 8 10 h h h h h 1 + -- - --- + ----- - ------- + -------- 12 720 30240 1209600 47900160 12 14 16 691 h h 3617 h ------------- + ----------- - ----------------- + 1307674368000 74724249600 10670622842880000 18 43867 h

20 174611 h

22

(1.46) (1.47) (1.48)

(1.49)

´ ´ ´ NUMERICA ´ 1.4. METODOS BASICOS DE INTEGRACION

7

------------------- - --------------------- + O[h] 5109094217170944000 802857662698291200000

Resulta entonces  T (h) = 1 +

∞ X

 c2j (hD)2j 

j=1

enhD − 1 f (x0 ) D

∞ X
enhD − 1 f (x0 ) + c2j h2j D2j−1 enhD − 1 f (x0 ) D j=1 Z b ∞   X = f (x) dx + c2j h2j f (2j−1) (b) − f (2j−1) (a)

=

a

(1.50)

(1.51)

(1.52)

j=1

lo que justifica el enunciado del teorema 1.3.1 1.3.3. La f´ ormula de Euler. En el ap´endice E (p´agina 67) demostramos una f´ ormula de suma de Euler que enunciamos a continuaci´on con modificaciones. Teorema 1.3.2. Supongamos que f es una funci´on suave, de clase C ∞ en el intervalo [a, b]. Sea n ≥ 1, h = (b − a)/n, xj = a + jh, para j = 0, 1, . . . , n. Sea, adem´ as, {x} = x − bxc. Entonces se tiene   Z b f (a) f (a) + f (x1 ) + . . . + f (xn−1 ) + f (x) dx = h 2 2 a   Z b x−a f 0 (x) dx (1.53) +h B1 h a La demostraci´ on est´ a presentada en el ap´endice. Se tiene la siguiente f´ ormula de integraci´on por partes (ver 74), Z b 1 Bm ({(x − a)/h})f (m) (x) dx = m! a 1 (Bm+1 (1)f (m) (b) − Bm+1 (0)f (m) (a) (m + 1)! Z b 1 − Bm+1 ({(x − a)/h})f (m+1) (x) dx (m + 1)! a

(1.54)

Utilizandola directamente se puede generalizar la f´ormula 1.53 para obtener por inducci´ on el teorema 1.3.1 de la p´agina 4 1.4. M´ etodos b´ asicos de integraci´ on num´ erica En las secciones 4.3 y 4.4 del texto de Burden y Faires (1996) se deducen y desarrollan los m´etodos b´ asicos simples y compuestos para la determinaci´on num´erica de la integral Z b f (x) dx (1.55) a

La deducci´ on de los mismos se basa en la aproximaci´on de la funci´on f en el intervalo (o, en el caso de los m´etodos compuestos, en subintervalos de) [a, b]. La lectura de las secciones citadas es instructiva puesto que suministra una metodolog´ıa alternativa para la obtenci´ on de las expresiones de los distintos m´etodos respecto

´ 1. INTEGRACION

8

de la que desarrollaremos a continuaci´on, que se denomina “m´etodo de coeficientes indeterminados”

1.4.1. El m´ etodo de coeficientes indeterminados. Para deducir los m´etodos b´ asicos de integraci´ on num´erica mediante este m´etodo reemplazamos la integral por una suma finita Z b N X f (x) dx ' wi f (xi ) (1.56) a

El m´ etodo de coeficientes indeterminados es una alternativa frente a los m´ etodos basados en la aproximaci´ on mediante polinomios

i=1

donde los N puntos a ≤ xi ≤ b se eligen en forma conveniente (lo que depender´a de la aplicaci´ on particular) y los coeficientes wi (que son pesos en un sentido general) se consideran indeterminados al principio y resultan ser las variables que debemos determinar. Para ello se reemplaza en la expresi´on una familia conveniente de polinomios de grado creciente y se utilizan las ecuaciones resultantes para determinar los wi y, eventualmente, estimar los errores de aproximaci´on del m´etodo. En una secci´ on posterior (ver 1.6 en la p´agina 12) se generalizar´a este m´etodo para el caso en que se utilizan las ecuaciones para determinar tanto los wi como los xi . No se pierde generalidad si se considera que el intervalo de integraci´on es [a, b] = [−1, 1], o bien [a, b] = [−h, h] por ello realizaremos las cuentas en alguno de estos intervalos particulares puesto que as´ı las expresiones suelen simplificarse. 1.4.2. El m´ etodo trapezoidal. Si utilizamos el producto Mathematica para obtener los coeficientes las entradas y salidas podr´ıan ser (* Trapezoidal *) a0=a;a1=b; (* Defino los puntos *) p0[x_]:=1; (* Defino los polinomios *) p1[x_]:=(2x-a-b)/(b-a); p2[x_]:=((2x-a-b)/(b-a))^2; s0=A p0[a0] + B p0[a1]//Simplify (* Integraci\’on num\’erica *) s1=A p1[a0] + B p1[a1]//Simplify; A + B i0=Integrate[p0[x],{x,a,b}]//Simplify (* Integro los polinomios *) i1=Integrate[p1[x],{x,a,b}]//Simplify; -a + b Solve[{s0==i0,s1==i1},{A,B}]//Simplify (* Resuelvo el sistema *) -a + b -a + b {{A -> ------, B -> ------}} 2 2

Se tiene, adem´ as, que el error ε1 para un paso h = b − a se puede escribir (utilizando la f´ ormula de interpolaci´on lineal o utilizando el error para el polinomio p2 (x) del cuaderno de Mathematica precedente) como ε1 (f ) = −

(b − a)3 00 f (η), 12

η ∈ [a, b]

(1.57)

´ ´ ´ NUMERICA ´ 1.4. METODOS BASICOS DE INTEGRACION

9

(Justificar), an´ alogamente el εn correspondiente al paso h = (b − a)/n resulta εn (f ) =

n X

−

j=1

h3 00 f (ηj ) 12 n

=−

h3 n 1 X 00 f (ηj ) 12 n j=1

(1.58)

y, por la suavidad de f (∈ C ∞ [a, b]) se tiene que n

∃η ∈ [a, b] : f 00 (η) =

1 X 00 f (ηj ) n j=1

(1.59)

en definitiva (b − a)h2 00 f (η) η ∈ [a, b] 12 Si utilizamos la f´ ormula del teorema 1.3.1 podemos escribir εn (f ) = −

εn (f ) = −

B 2 h2 0 (f (b) − f 0 (a)) + O(h4 ) 2!

(1.60)

(1.61)

h2 00 f (ξ)(b − a) + O(h4 ), ξ ∈ [a, b] (1.62) 12 lo que nos da un resultado del mismo orden. A partir de la expresi´ on 1.58 se puede seguir otra l´ınea deductiva: tomamos l´ımite,   n X h εn (f ) = lim − f 00 (ηj ) (1.63) lim n→∞ h2 n→∞ 12 j=1   n X 1 =− lim  f 00 (ηj )h (1.64) 12 n→∞ j=1 −

pero, en esta expresi´ on, xj−1 ≤ ηj ≤ xj , j = 1, . . . , n con lo que la f´ormula 1.64 es una suma de Riemann, en consecuencia ! Z b εn (f ) 1 1 00 lim = lim − f (x) dx = − (f 0 (b) − f 0 (a)) (1.65) n→∞ h2 n→∞ 12 a 12 es decir que h2 0 (f (b) − f 0 (a)) (1.66) 12 εf on al error asint´otico lo que podemos establecer en la sin (f ) es una aproximaci´ guiente definici´ on: εn (f ) ' εf n (f ) ≡ −

´ n 1.4.1. Sean εn (f ) y εf Definicio on del n (f ), el error exacto y una estimaci´ error respectivamente. Decimos que εf on asint´ otica del error n (f ) es una estimaci´ εn (f ) si εf n (f ) lim =1 (1.67) n→∞ εn (f )

´ 1. INTEGRACION

10

Utilizando εf etodo trapezoidal corregido ThC (f ) n (f ) es posible dar un m´ ThC (f ) = Th (f ) + εf n (f )

(1.68)

es decir ThC (f )

 =h

f0 fn + f1 + f2 + . . . + fn−1 + 2 2

 −

h2 0 (f (b) − f 0 (a)) 12

(1.69)

Ejercicio 1.4.1. Estudiar el comportamiento num´erico del m´etodo trapezoidal corregido ThC (f ) 1.4.3. El m´ etodo de Simpson. Utilizamos el producto Mathematica para determinar los coeficientes del m´etodo de Simpson (ver Ap´endice B, p´agina 37) Obtenemos     b−a b−a a+b S = f (a) + 4f ( ) + f (b) (1.70) 2 6 2 y un error de truncamiento (b − a)5 (4) a + b f ( ) 2880 2

(1.71)

lo que equivale a
b−a 5 2

a+b f (4) ( ) (1.72) 90 2 En forma an´ aloga al caso de m´etodo trapezoidal y utilizando las expresiones precedentes se tiene que
5 4 b−a f (4) (η) 2 , η ∈ [a, b] (1.73) ε2 (f ) = − 15 24  5 b−a 1 (4) =− f (η), η ∈ [a, b] (1.74) 2 90 y n
2 h5 n2 2 X f (4) (ηj ) εn (f ) = − 90 n j=1 εt =

=−

h4 (b − a) (4) f (η), 180

η ∈ [a, b]

(1.75)

y la f´ ormula asint´ otica h4 (3) (f (b) − f (3) (a)) (1.76) 180 (Esta u ´ltima tambi´en se obtiene a partir del desarrollo asint´otico de Euler) εf n (f ) = −

1.4.4. M´ etodos de Newton. En la secci´on precedente hemos obtenido la f´ ormula del m´etodo de Simpson. La misma metodolog´ıa se puede utilizar para cualquier elecci´ on de los puntos de la partici´on del intervalo de integraci´on. Por ejemplo, si desearamos calcular los que corresponden a una partici´on regular de once puntos, nuevamente utilizamos el producto Mathematica para obtener los coeficientes (ver ap´endice B en p´agina 37). Estos m´etodos llevan el nombre de m´etodos de Newton, y en algunos textos tambi´en el de f´ ormulas de Newton-Cotes.

´ ´ 1.5. METODOS DE EXTRAPOLACION

11

1.5. M´ etodos de extrapolaci´ on En esta secci´ on utilizaremos la expresi´on que calculamos para el error en el m´etodo trapezoidal para ejemplificar el uso de los m´etodos de extrapolaci´on repetida de Richardson. Supongamos que b0 = I(0) es la integral que deseamos calcular y que, para el paso h se tiene I(h) = b0 + b1 h2 + b2 h4 + b3 h6 + . . . (1.77) (f´ ormula que sabemos v´ alida para el mencionado m´etodo trapezoidal). El error de truncamiento es n−1 X I(h) − b0 = bi h2i + O(h2n ) (1.78) i=1

de modo que, con 1.77 y h h2 h4 h6 + b2 + b3 + ... I( ) = b0 + b1 2 4 16 64

(1.79)

se obtiene efectuando ‘4 por (1.79) − (1.77)’ 3 15 h 4I( ) − I(h) = 3b0 − b2 h4 − b3 h6 + . . . 2 4 16 y, dividiendo por el coeficiente de b0 I (0) ( h2 ) − I (0) (h) 1 5 h I (1) (h) = I (0) ( ) + = b 0 − b 2 h 4 − b3 h 6 + . . . 2 3 4 16

(1.80)

(1.81)

(ver la secci´ on C.2 en el ap´endice, donde se deducen estas expresiones con el Mathematica). El error de truncamiento es 1 5 I (1) (h) − b0 = − b2 h4 − b3 h6 + . . . 4 16

(1.82)

Si repetimos este c´ alculo (ver C.2) se obtiene I (0) (h) = I(h)

(1.83) I (0) ( h2 ) − I (0) (h) 22 − 1

h I (1) (h) = I (0) ( ) + 2 I (1) ( h2 ) − I (1) (h) h I (2) (h) = I (1) ( ) + 2 24 − 1 .. . h I (j+1) (h) = I (j) ( ) + 2

I (j) ( h2 ) − I (j) (h) 22j+2 − 1

(1.84) (1.85) (1.86) (1.87)

y, en general, I (j) (h) = b0 + O(h2j+2 )

(1.88)

El esquema resultante es que cada elemento de la segunda columna se obtiene con los dos inmediatos a su izquierda, cada vez que se calcula una aproximaci´on con paso m´ as fino se puede calcular una antidiagonal de la tabla, y, por u ´ltimo, el criterio de detenci´ on es que dos valores en la misma columna difieran en menos que un valor de tolerancia preestablecido.

´ 1. INTEGRACION

12

Tabla 1. Extrapolaci´on de Richardson (M´etodo de Romberg) I (0) (h)

I (1) (h)

I (2) (h)

I (3) (h)

I (0) ( h2 )

I (1) ( h2 ) I (2) ( h2 ) I (3) ( h2 )

I (0) ( h4 )

I (1) ( h4 ) I (2) ( h4 )

I (0) ( h8 )

I (1) ( h8 )

h I (0) ( 16 )

...

I (4) (h) . . . ...

...

...

.. . 1.6. Cuadratura de Gauß La expresi´ on aproximada N X

wi f (xi )

(1.89)

i=1

Rb para la integral a f (x) dx es mucho m´as u ´til si no se escogen de antemano los puntos xi . El problema de hallar los puntos y los pesos de modo que la expresi´on sea exacta para polinomios del mayor grado posible es resoluble y conduce a los m´etodos de integraci´ on de Gauß. En el ap´endice B (p´ agina 37) incluimos un cuaderno del Mathematica que calcula los coeficientes para un m´etodo de Gauß de tres puntos (se realiza en el intervalo [−1, 1], pero un simple cambio de variables permite extenderlo a cualquier intervalo [a, b]) Es conveniente, si embargo, atacar el problema de cuadratura en un contexto un poco m´ as general. Sea v una funci´ on de peso positiva en el intervalo [−1, 1], si x0 , x1 , . . . , xm se eligen como los ceros del polinomio pm+1 de grado m+1 en la familia de polinomios ortogonales asociada a v(x), entonces la f´ormula Z 1 f (x)v(x) dx ' w0 f0 + w1 f1 + . . . + wm fm (1.90) −1

La inclusi´ on de la funci´ on peso v facilita considerar distintas familias de polinomios ortogonales

es exacta para todos los polinomios de orden 2m + 2 siempre que los coeficientes satisfagan Z 1 wi = δi (x)v(x) dx (1.91) −1

donde δi (x) es el polinomio de grado m (orden m + 1) tal que  0 j 6= i δi (xj ) = j = 0, 1, 2, . . . , m 1 j=i

(1.92)

Para verificar esta afirmaci´ on tomemos un polinomio f de orden 2m + 2, existen entonces polinomios q y r de orden m + 1 tales que f = qpm+1 + r

(1.93)

´ 1.7. METODOS DE MONTECARLO

13

Por lo tanto, debido a la ortogonalidad de la familia {pi } Z 1 Z 1 Z 1 Z f (x)v(x) dx = q(x)pm+1 (x)v(x) dx + r(x)v(x) dx = −1

−1

−1

1

r(x)v(x) dx

−1

(1.94) Asimismo, dado que xi es un cero de pm+1 , i = 0, 1, . . . , m, se tiene m X i=0

wi fi =

m X

wi q(xi )pm+1 (xi ) +

i=0

m X

wi r(xi ) =

i=0

m X

wi r(xi )

(1.95)

i=0

pero, tambi´en se tiene Z

1

r(x)v(x) dx = −1

m X

wi r(xi )

(1.96)

i=0

porque los coeficientes wi fueron elegidos de modo que la f´ormula sea exacta para todos los polinomios de grado m o menor (orden m + 1) ´ n 1.6.1. Aplicando la f´ormula 1.90 al caso f (x) = (δi (x))2 se tiene Observacio (fj = 0 para j 6= i) Z 1 (δi (x))2 v(x) dx = wi (1.97) −1

por lo que los coeficientes wi en las f´ormulas de cuadratura de Gauß son positivos. ´ n 1.6.2. Para obtener los coeficientes se toman los xi como los Observacio ceros de los polinomios del orden deseado y se utiliza un calculador simb´olico como el Mathematica para obtener los coeficientes y para estimar el error de truncamiento local. En diversos manuales de f´ormulas matem´aticas figuran los puntos y coeficientes de Gauß para distintas familias de polinomios y distintos ´ordenes. El an´ alisis del error puede desarrollarse en forma an´aloga a los casos considerados previamente, por ello dejamos los detalles para el lector interesado. 1.7. M´ etodos de Montecarlo En esta secci´ on consideramos un m´etodo conceptualmente distinto, respecto de los anteriormente presentados, para obtener aproximaciones a las integrales Z b f (x) dx (1.98) a

Suponemos que se posee una definici´on de la funci´on f , por ejemplo en un m-file de Matlab function y = mifun(x) % MIFUN define mi funci\’on [mx,nx]=size(x) if (mx == 1) & (nx == 1) y = 2 * x; else disp(’error’); end % if

o, tambi´en, si la funci´ on est´ a dada por N puntos en un tabla mif de N × 2 donde mix=mif(:,1) y miy=mif(:,2), se la calcula, por ejemplo, mediante y=spline(mix,miy,x) o bien y=interp1(mix,miy,x) seg´ ´ un se desee.

´ 1. INTEGRACION

14

Para la funci´on f se determinan M f = sup f (x)

(1.99)

a≤x≤b

y mf = inf f (x) a≤x≤b

La funci´ on rand permite obtener puntos al azar en un dominio

(1.100)

(en realidad basta que se tomen cotas superior e inferior respectivamente). El m´etodo se aplica entonces a la funci´on fe(x) = f (x) − mf que resulta no negativa y con su gr´ afico contenido en el rect´angulo R = [a, b] × [0, M f − mf ], y consiste en lo siguiente. Se hallan P puntos al azar en el rect´angulo R y se determina el n´ umero Pf de los que se encuentran en el ´area que deseamos aproximar. A continuaci´on se calcula la estimaci´ on del ´ area mediante la expresi´on   Pf + mf (1.101) (b − a) (M f − mf ) P En la secci´ on C.3 incluimos ejemplos de integrales de Montecarlo. Este m´etodo no suele ser utilizado para el caso de ´areas puesto que existen buenos m´etodos de cuadratura (como los ya estudiados), sin embargo, el caso de los vol´ umenes o hipervol´ umenes es distinto. Supongamos que se desea conocer el volumen de la intersecci´ on de 5 cuerpos que se cortan en el espacio Eucl´ıdeo tridimensional, y que adem´ as, se desea distinguir la parte de la intersecci´on que pertenece a un sexto cuerpo de aquella que no est´a contenida en el mismo. El intento de obtener la integral por m´etodos cl´asicos o num´ericos puede ser de una complejidad excesiva a´ un para casos de cuerpos definidos por ecuaciones muy simples. El uso del m´etodo de Montecarlo permite superar estas dificultades. 1.7.1. Determinaci´ on de un volumen. En esta secci´on calcularemos el volumen del s´ olido que resulta de la intersecci´on, en el primer octante, de los cilindros x2 + z 2 = 1

(1.102)

y 1 y 2 + (z − )2 = 1 (1.103) 2 En primer t´ermino atacamos el problema con Matlab, a continuaci´on desarrollamos un cuaderno del Mathematica en el que obtenemos una soluci´on similar para ‘verificar’. %

MONTEK

volumen por m\’etodo de Montecarlo

%

C.E.Neuman, 11-may-97

rand(’seed’,sum(100*clock)); % inicializa el generador de n\’umeros al azar scores=[]; NN=20; N=NN^4; h=1/NN;

% en esta variable se acumulan los resultados % da el numero de cubitos por lado % da el numero total de puntos % da el paso (dimension del cubito)

for indice=1:10, X=zeros(N,3); for i=1:NN,

% construcci\’on de los puntos al azar

´ 1.7. METODOS DE MONTECARLO

im1=i-1; % en cada cubito for j=1:NN, jm1=j-1; for k=1:NN, km1=k-1; ps = rand(NN,3); ini=im1*NN*NN*NN+jm1*NN*NN+km1*NN+1; fin=ini-1+NN; X(ini:fin,:) = [ps(:,1)+i-1 ps(:,2)+j-1 ps(:,3)+k-1]*h; end % rof end % rof disp(i); end % rof x=X(:,1); y=X(:,2); % asignaci\’on de coordenadas al azar z=X(:,3); %%% %%%

la siguiente variable calcula los valores de las funciones de corte para las variables: TEST=[y.*y+4*(z-0.5).*(z-0.5) x.*x+z.*z];

%%% %%% %%%

en I se seleccionan los \’{\i}ndices correspondientes a puntos contenidos en el volumen que se desea medir (intersecci\’on de los dos cuerpos definidos por las expresiones de TEST) I=find( (TEST(:,1)<=1) & (TEST(:,2)<=1) ); [mI,nI]=size(I); scores=[scores; mI/N]; % acumulaci\’on del resultado

disp(indice); end % rof mime=mean(scores) mide=std(scores) scorord=sort(scores); [m,n]=size(scores); mimeord=mean(scores(0.1*m:0.9*m)) mideord=std(scores(0.1*m:0.9*m)) figure; hist(scores); % Salidas del programa precedente %mime = % 0.6352 %mide = % 3.8221e-004 %mimeord = % 0.6352 %mideord = % 4.0451e-004 %scores = % 0.6349 % 0.6347 % 0.6356 % 0.6356 % 0.6347 % 0.6356

15

´ 1. INTEGRACION

16

% % % %

0.6350 0.6356 0.6356 0.6352

Esta regularidad se ve confirmada por el siguiente cuaderno del Mathematica (*

M\’etodo de Montecarlo para integrales triples C\’alculo num\’erico de la integral

*)

(* Determino el l\’{\i}mite de separaci\’on *) Solve[4(1-x^2)==(1+Sqrt[1-y^2])^2,y] 2 2 {{y -> -Sqrt[-4 + 4 x - 4 Sqrt[1 - x ]]}, 2 2 {y -> Sqrt[-4 + 4 x - 4 Sqrt[1 - x ]]}, 2 {y -> -Sqrt[-4 + 4 x

2 + 4 Sqrt[1 - x ]]},

2 2 {y -> Sqrt[-4 + 4 x + 4 Sqrt[1 - x ]]}} (* Defino el l\’{\i}mite de integraci\’on *) f[x_]:=Sqrt[-4+4x x+4Sqrt[1-x x]] (* Determino el punto de separaci\’on *) Solve[f[x]==1,x] -Sqrt[3] -Sqrt[3] Sqrt[3] {{x -> --------}, {x -> --------}, {x -> -------}, 2 2 2 Sqrt[3] {x -> -------}} 2 (* Defino el punto de separaci\’on *) c=Sqrt[3]/2 Sqrt[3] ------2 (* Defino los extremos de integraci\’on *) geu[x_,y_]:=(1/2)(1+Sqrt[1-y y]) ged[x_,y_]:=(1/2)(1-Sqrt[1-y y]) fcu[x_,y_]:=Sqrt[1-x x] (* Calculo las integrales en forma num\’erica *) I1=NIntegrate[1,{x,0,c},{y,f[x],1},{z,ged[x,y],geu[x,y]}] I2=NIntegrate[1,{x,0,c},{y,0,f[x]},{z,ged[x,y],fcu[x,y]}] I3=NIntegrate[1,{x,c,1},{y,0,f[x]},{z,ged[x,y],fcu[x,y]}] I1+I2+I3 0.244751 0.358503 0.0320788 0.635332 (* valor de la integral *)

1.8. M´ etodos adaptativos En el u ´ltimo ap´endice (F, p´agina 75) se propone un laboratorio de Matem´atica orientado a estudiar algoritmos adaptativos para la integraci´on num´erica. En la secci´ on C.4, en la p´agina 53, se incluye —para ilustraci´on— el c´odigo de una funci´ on de Matlab, adaptada de la denominada quad, que implementa un algoritmo adaptativo recursivo de bajo orden basado en el m´etodo de Simpson. La

1.9. EJERCICIOS

17

funci´ on miquad es en realidad un driver que llama a la verdadera funci´on recursiva que se denomina miquadst la que realiza cada paso de la recursi´on. 1.8.1. Comparaci´ on entre NIntegrate y quad. Calculamos el ´area bajo la curva normal, que puede encontrarse en cualquier tabla estad´ıstica, utilizando las funciones NIntegrate del Mathematica y quad de Matlab. La primera da (los resultados de medici´ on de tiempos se expresan en segundos) (* Calculo con el Mathematica *) f[x_]:=(1/Sqrt[2Pi]) E^(-x^2/2);//Timing {0. Second, Null} NIntegrate[f[x],{x,0,1}]//Timing {0.11 Second, 0.341345} N[NIntegrate[f[x],{x,0,1}],15]//Timing {0.11 Second, 0.341344746068543}

Para la segunda es necesario definir en un m-archivo la funci´on que se desea integrar: function y=nintquad(x) % NINTQUAD Densidad normal %%% C.E.Neuman, 10 de mayo de 1997 y=(1/sqrt(2*pi))*exp(-x.^2/2); % fin de nintquad

La salida y la medici´ on del tiempo empleado con Matlab es Salidas de la integraci\’on adaptativa con Matlab t=cputime; Q=quad1(’nintquad’,0,1); cputime-t, Q ans = 0.06 Q = 0.34134540613909 t=cputime; Q=quad1(’nintquad’,0,1,1e-5); cputime-t, Q ans = 0.11 Q = 0.34134474647804 t=cputime; Q=quad1(’nintquad’,0,1,1e-7); cputime-t, Q ans = 0.6 Q = 0.34134474607765 t=cputime; Q=quad1(’nintquad’,0,1,1e-9); cputime-t, Q ans = 2.47 Q = 0.34134474606860

Observamos que, en esta medici´on, —la que, al margen, no es muy precisa,— parece que la primera metodolog´ıa es m´as eficiente que la segunda. 1.9. Ejercicios Ejercicio 1.9.1. Utilizar el m´etodo de Romberg para calcular la integral Z 4 f (x) dx (1.104) 0

´ 1. INTEGRACION

18

donde f (x) est´ a definida por la siguiente tabla: x f (x)

0.0 −4271

0.5 −2522

1.0 −499

1.5 1795

2.0 4358

2.5 7187

3.0 10279

3.5 13633

4.0 17247

¿Son necesarios todos los valores? Ejercicio 1.9.2. Probar que la f´ormula " r ! Z 1 1 3 f (x) dx ' 5f − + 8f (0) + 5f 9 5 −1

r !# 3 5

es exacta para polinomios de quinto grado, y aplicarlo al c´alculo de Z 1 sin x dx 1 +x 0

(1.105)

(1.106)

Ejercicio 1.9.3. (a) Deducir una f´ormula de integraci´on de Gauß de dos puntos para integrales de la forma Z 1 f (x)(1 + x2 ) dx (1.107) −1

que sea exacta cuando f (x) es un polinomio de grado 3. (b) Aplicar la f´ ormula a f (x) = x4 . Usar el resultado para obtener una aproximaci´ on del resto. Ejercicio 1.9.4. Supongamos que una integral numericamente aproximada se representa mediante el s´ımbolo In donde n es el n´ umero de subdivisiones del intervalo de integraci´ on. Si se calculan los valores de In para alg´ un m´etodo num´erico, analizar en forma emp´ırica la velocidad de convergencia de In a I (el valor de la integral calculada exactamente) mediante los cocientes Rn =

I2n − In I4n − I2n

(1.108)

Dar ejemplos y explicar los resultados. Ejercicio 1.9.5. [Optativo] Considerar el conjunto S1 de funciones definidas en el intervalo [a, b] de la siguiente forma. Sea n > 0, h = (b−a)/n, tj = a+jh, para j = 0, 1, . . . , n. Sea f (x) definida por la propiedad de ser lineal en cada subintervalo [tj−1 , tj ], para j = 1, . . . , n. Mostrar que este conjunto de funciones f , para todo n ≥ 1, es denso en el conjunto C[a, b] de funciones continuas en el intervalo [a, b]. Ejercicio 1.9.6. Obtener f´ormulas de integraci´on de Gauß para Z 1 n X I= xf (x) dx ' wj f (xj ) 0

(1.109)

j=1

con funci´ on peso w(x) = x Ejercicio 1.9.7. Considerar la siguiente tabla de integrales aproximadas In obtenidas mediante la regla de Simpson. Predecir el orden de convergencia de la sucesi´ on de In a la integral I:

1.9. EJERCICIOS

n 2 4 8 16 32 64

19

In 0.28451779686 0.28559254576 0.28570248748 0.28571317731 0.28571418363 0.28571427643

Es decir que, si I − In ' c/np , entonces, ¿cu´anto vale p? El resultado ¿resulta ser una forma v´ alida para el error de estos datos? Predecir un valor de c y del error en I64 . ¿Qu´e valor de n es necesario alcanzar para que el error en In sea menor que 10−11 ? Ejercicio 1.9.8. Denotamos InT a la regla trapezoidal para aproximar la inteRb gral a f (x) dx, y, analogamente, InM para la f´ormula de ‘midpoint’. Los respectivos errores asint´oticos —en el caso en que f sea suficientemente regular en [a, b],— son I − InT = − y

h2 0 [f (b) − f 0 (a)] + O(h4 ) 12

(1.110)

h2 0 [f (b) − f 0 (a)] + O(h4 ) (1.111) 24 Utilizar estos resultados para obtener un nuevo m´etodo num´erico, Ien , con orden de convergencia mayor, combinando InT y InM . ¿Cu´ales son los pesos de la nueva f´ ormula Ien ? I − InM =

CAP´ıTULO 2

Ecuaciones diferenciales 2.1. Introducci´ on Decimos que la ecuaci´ on diferencial ordinaria y 0 = f (x, y)

(2.1)

en una funci´ on inc´ ognita y, con f (x, y) funci´on continua en un dominio D del plano, tiene una soluci´ on y(x) en un intervalo x0 ≤ x ≤ x1 , si y(x) es una fuci´on diferenciable, (x, y(x)) est´ a en el dominio D para cada x del intervalo [x0 , x1 ] y, adem´ as, y 0 (x) = f (x, y(x)) en [x0 , x1 ]. Desde el punto de vista geom´etrico podemos considerar que la ecuaci´on y 0 = f (x, y) define un campo continuo de direcciones sobre D. Y que las funciones soluci´ on son tangentes a esas direcciones (en cada punto del dominio).

El campo de direcciones permite estimar gr´ aficamente las trayectorias

y 6 3 2 1

0

B Br B B A Ar A A @r @ 1

E Er E E C Cr C C A Ar A A

E Er E B E Br B B

2

3

x

Figura 1. Campo de direcciones de f (x, y) = −xy Ejemplo 2.1.1. Sea la ecuaci´on diferencial y 0 = −xy, entonces y 0 /y = −x y es posible integrar ambos miembros obteniendo log y = −x2 /2 + C 0 , expresi´on de la que tomando la exponencial de ambos miembros se deduce la soluci´on general 2

y(x) = Ce−x

/2

(2.2)

de la ecuaci´ on diferencial. En la figura 1 esquematizamos el campo de direcciones definido por f (x, y) = −xy. (En cada punto (x, y) dibujamos un peque˜ no segmento en la direcci´ on de la recta por el punto con pendiente −xy.) 21

Una soluci´ on ε-aproximada difiere casi uniformemente de la soluci´ on en menos de ε

22

2. ECUACIONES DIFERENCIALES

´ n 2.1.1. Sea f (x, y) una funci´on continua definida en el dominio D. Definicio Una funci´ on y(x), definida en el intervalo [x1 , x2 ], es una soluci´on (aproximada) de y 0 = f (x, y) con error menor que ε si (i) (x, y(x)) ∈ D, x1 ≤ x ≤ x2 (ii) y(x) es continua en [x1 , x2 ] (iii) y(x) tiene derivada continua a trozos en [x1 , x2 ] que puede no estar definida en un n´ umero finito de puntos del intervalo [x1 , x2 ], llam´emoslos ξ1 , ξ2 , . . ., ξn (iv) |y 0 (x) − f (x, y(x))| ≤ ε, x1 ≤ x ≤ x2 , x 6= ξi , i = 1, 2, . . . , n. Una vez definidas las soluciones aproximadas de una ecuaci´on diferencial deseamos demostrar un resultado de existencia de las mismas que, adem´as, resulta constructivo, puesto que nos da un m´etodo para obtener la aproximaci´on Teorema 2.1.1. Sea (x0 , y0 ) un punto del dominio D (de la definici´on 2.1.1) y supongamos que el rect´ angulo R = {|x − x0 | ≤ a, |y − y0 | ≤ b} est´a contenido en D. Sea |f (x, y)| ≤ M , (x, y) ∈ R. Entonces, si h = min(a, b/M ), puede construirse una soluci´ on aproximada y(x) de y 0 = f (x, y) en el intervalo |x − x0 | ≤ h, tal que y(x0 ) = y0 , donde el error ε puede ser un n´ umero positivo arbitrariamente peque˜ no. Notar que h es independiente de ε.

P # # # # # # " # " " #" s # !(x (x0 + h, y0 ) ! O# 1 , y1 ) ! (x0 , yc 0) α 6 c c c c Mh c c c c S c ?  R h Q Figura 2. El rect´angulo S est´a contenido en el R ´ n: El rect´ Demostracio angulo S = {|x − x0 | ≤ h, |y − y0 | ≤ M h} est´a contenido en R por la definici´ on de h (ver la figura 2). Supongamos dado el ε del teorema. Como f (x, y) es continua en S, resulta uniformemente continua en S; es decir que, dado ε > 0 (que lo tomamos como el del teorema) existe δ > 0 tal que si |˜ x − x| ≤ δ y |˜ y − y| ≤ δ con (˜ x, y˜) ∈ S y (x, y) ∈ S, entonces |f (˜ x, y˜) − f (x, y)| ≤ ε

(2.3)

Sea x1 , . . . , xn−1 un conjunto de puntos tales que: x0 < x1 < x2 < · · · < xn−1 < xn = x0 + h, y xi − xi−1 ≤ min(δ, δ/M ),

i = 1, . . . , n

(2.4)

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

23

Construiremos la soluci´ on aproximada en el intervalo x0 ≤ x ≤ x0 + h; un proceso similar permite definirla en el intervalo x0 − h ≤ x ≤ x0 . La soluci´ on aproximada ser´a una poligonal construida de la siguiente manera: desde (x0 , y0 ) se dibuja un segmento hacia la derecha con pendiente f (x0 , y0 ), este corta la recta x = x1 en un punto (x1 , y1 ). Desde (x1 , y1 ) se dibuja un segmento hacia la derecha con pendiente f (x1 , y1 ), que corta la recta x = x2 en y2 ; etc. El punto (x1 , y1 ) debe pertenecer al tri´angulo OQP de la figura 2 porque el ´angulo α se toma de modo que tan α = M , y |f (x0 , y0 )| ≤ M . Por motivos an´alogos (x2 , y2 ) tambi´en est´ a en OQP ; etc. En consecuencia el proceso puede continuarse hasta xn = x0 + h, porque la u ´nica raz´on por la que podr´ıa detenerse ser´ıa que f (xk , yk ) no estuviese definida, caso en el que se tendr´ıa |yk − y0 | > M h lo que ser´ıa contrario a la construcci´ on. Podemos definir y(x) por las f´ormulas recursivas y(x) = yi−1 + (x − xi−1 )f (xi−1 , y(xi−1 ))

(2.5)

donde yi−1 = y(xi−1 ),

xi−1 ≤ x ≤ xi ,

i = 1, . . . , n

(2.6)

Por su definici´ on y(x) es admisible, continua, y tiene derivada continua a trozos y 0 (x) = f (xi−1 , y(xi−1 ),

xi−1 < x < xi ,

i = 1, . . . , n

(2.7)

que no est´ a definida solamente en los puntos xi , i = 1, . . . , n − 1. Adem´as, si xi−1 < x < xi |y 0 (x) − f (x, y(x))| = |f (xi−1 , yi−1 ) − f (x, y(x))|

(2.8)

Pero, por 2.4 |x − xi−1 | < min(δ, δ/M ), y, por 2.5 |y − yi−1 | ≤ M |x − xi−1 | ≤ M

δ =δ M

(2.9)

En consecuencia, por 2.3 |f (xi−1 , yi−1 ) − f (x, y(x))| ≤ ε

(2.10)

y |y 0 (x) − f (x, y(x))| ≤ ε,

x 6= xi ,

i = 1, 2, . . . , n − 1

(2.11)

Resulta as´ı que y(x) satisface todas las condiciones de la definici´on previa y la construcci´ on requerida por el teorema ha sido completada.(ver ap´endice H) Este m´etodo de construir una soluci´on aproximada se conoce como m´etodo de Euler. Es innecesario mejorar el valor del h del teorema porque, en general, se puede demostrar que y(x) est´a definida en un intervalo mayor que |x − x0 | ≤ h. 2.2. M´ etodos en diferencias de un paso 2.2.1. M´ etodos de Euler y retroEuler. El m´etodo de Euler para la resoluci´ on aproximada de la ecuaci´on diferencial  y 0 = f (x, y) (2.12) y(x0 ) = y0 se puede formular as´ı:

M´ etodo de Euler

yn+1 = yn + hf (xn , yn ),

n≥0

donde h es de ahora en adelante el paso del m´etodo en estudio.

(2.13)

24

Estabilidad del m´ etodo de Euler

2. ECUACIONES DIFERENCIALES

En cada paso se cambia de miembro de la familia de soluciones de la ecuaci´on diferencial, de manera que la precisi´on del m´etodo deber´a depender de la estabilidad de las ecuaciones. Si las ecuaciones son estables, los errores en los primeros pasos tendr´ an peque˜ no efecto posterior. Definimos provisoriamente la noci´on de estabilidad de un m´etodo num´erico en funci´on de su comportamiento para resolver la ecuaci´ on diferencial y 0 = λy con λ un n´ umero complejo. La regi´on de estabilidad absoluta es entonces el subconjunto de puntos hλ del plano complejo con h ≥ 0 para los que una perturbaci´ on en un u ´nico valor yn produce una sucesi´on de cambios en los siguientes valores que no crecen de paso a paso. Ejemplo 2.2.1. Para el m´etodo de Euler se tiene yn+1 = yn + λhyn = (1 + λh)yn

(2.14)

de modo que este m´etodo es absolutamente estable en la regi´on |1 + λh| ≤ 1 que es el c´ırculo unitario del plano complejo centrado en el punto (−1, 0). Ejemplo 2.2.2. El m´etodo de Euler aplicado a la ecuaci´on del ejemplo 2.1.1 resulta ser yn+1 = yn − hxn yn = (1 − hxn )yn (2.15) El m´etodo de Euler puede obtenerse tambi´en del siguiente razonamiento. Se reemplaza la derivada y 0 en la ecuaci´on y 0 (x) = f (x, y)

(2.16)

por un cociente de incrementos en (xn , yn ) para obtener yn+1 − yn = f (xn , yn ) (2.17) h resultando el m´etodo al despejar la ecuaci´on para yn+1 . Si, en cambio, se utiliza el mismo cociente, pero en (xn+1 , yn+1 ), se obtiene yn+1 − yn = f (xn+1 , yn+1 ) h M´ etodo de retroEuler

(2.18)

lo que conduce al m´etodo denominado “retroEuler” (backward Euler) yn+1 = yn + hf (xn+1 , yn+1 )

(2.19)

que tiene la caracter´ıstica especial de ser impl´ıcito y requiere la soluci´on de una ecuaci´ on en cada paso. Ejemplo 2.2.3. Para el m´etodo retroEuler se tiene yn+1 = yn + λhyn+1

(2.20)

es decir que luego de despejar yn+1 yn+1 =

1 yn 1 − λh

(2.21)

de modo que este m´etodo es absolutamente estable en la regi´on |1/(1 − λh)| ≤ 1 v´ alida en el semiplano de partes reales no positivas del plano complejo. Estabilidad del m´ etodo retroEuler

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

25

Ejemplo 2.2.4. El m´etodo retroEuler aplicado a la ecuaci´on del ejemplo 2.1.1 resulta ser yn+1 = yn − hxn+1 yn+1 (2.22) es decir yn+1 =

1 yn (1 + hxn+1 )

(2.23)

2.2.2. M´ etodos de Runge–Kutta. Nuestro siguiente objetivo es generalizar el m´etodo de Euler utilizando los desarrollos en polinomios de Taylor pero limitando lo m´ as posible el c´ alculo de derivadas de las funciones involucradas. Recordemos que en este m´etodo se puede escribir yh (x) = y(x) + hD(x) + O(h2 )

(2.24)

yh/2 (x) = y(x) + (h/2)D(x) + O(h2 )

(2.25)

y que para el paso mitad

(¿porqu´e?) y que aplicando la extrapolaci´on (de Richardson) para eliminar D(x) y(x) = 2yh/2 (x) − yh (x) + O(h2 )

(2.26)

Apliqu´emosla a un paso de longitud h, se tiene yh (xn + h) = yn + hf (xn , yn )

(2.27)

yh/2 (xn + (h/2)) = yn + (h/2)f (xn , yn )

(2.28)

y para medio paso

en consecuencia yh/2 (xn + h) = yn + (h/2)f (xn , yn ) +(h/2)f (xn + (h/2), yn + (h/2)f (xn , yn ))

(2.29)

luego, aplicando la extrapolaci´on, y(xn+1 ) = 2yn + hf (xn , yn ) + hf (xn + (h/2), yn + (h/2)f (xn , yn )) −yn − hf (xn , yn ) + O(h2 )

(2.30)

y simplificando yn+1 = yn + hf (xn + (h/2), yn + (h/2)f (xn , yn ))

(2.31)

En definitiva, podemos escribir yn+1 = yn + hϕ(xn , yn , h; f )

(2.32)

ϕ(x, y, h; f ) = f (x + (h/2), y + (h/2)f (x, y))

(2.33)

con ¿En qu´e consiste esta primera f´ormula de Runge en forma gr´afica? En la figura 3 representamos el efecto de aplicar la primera f´ormula de Runge en el paso de xn a xn+1 . Notemos que para el m´etodo de Euler se puede definir en forma an´aloga ϕ(x, y, h; f ) = f (x, y). Ejercicio 2.2.5. ¿Es posible expresar en esta notaci´on el m´etodo de retroEuler?

Primera f´ ormula de Runge

26

2. ECUACIONES DIFERENCIALES

y 6

0

10 rr Runge un (xn+1 )







r Euler











xn

xn+1

x

Figura 3. Primer m´etodo de Runge Para estudiar la estabilidad aplicamos la f´ormula a la ecuaci´on y 0 = λy. Se tiene yn+1 = yn + hλ(yn + (h/2)λyn ) = yn (1 + hλ + 12 h2 λ2 ) raz´on por la cual debe tenerse que |1 + hλ + 12 h2 λ2 | ≤ 1. Este primer m´etodo de Runge es generalizaci´on del de Euler. Aunque existen otras posibles generalizaciones, por ejemplo utilizando series de potencias, el problema es que hay que derivar la funci´on f . La idea de Runge fue dar f´ormulas que se originan en tomar combinaciones de valores de f (x, y) en puntos adecuadamente elegidos para obtener, sin derivaciones, expresiones cuyos desarrollos coincidan en sus primeros t´erminos con la serie de potencias desarrollada en el entorno del punto. Esas f´ ormulas fueron mejoradas por Heun y por Kutta. Llamaremos m´etodo de Runge-Kutta (RK) a yn+1 = yn + hϕ(xn , yn , h; f ),

n≥0

(RK)

(2.34)

y en lo que sigue determinaremos diversas elecciones para ϕ. Es natural esperar que ϕ(x, y(x), h; f ) ' y 0 (x) = f (x, y(x)) si h es peque˜ no. El error de truncamiento local es en este m´etodo εn+1 (y) = y(xn+1 ) − y(xn ) − hϕ(xn , y(xn ), h; f ),

n≥0

(2.35)

♣ Ilustremos la deducci´ on de distintas ϕ en el caso de segundo orden de aproximaci´ on. La propuesta es ϕ(x, y, h; f ) = γ1 f (x, y) + γ2 f (x + αh, y + βhf (x, y))

(2.36)

donde deben determinarse las constantes γ1 , γ2 , α, y β. Para ello desarrollemos el error de truncamiento local en ´este m´etodo εn+1 (y) = Y (xn+1 ) − Y (xn ) − hϕ(xn , Y (xn ), h; f ),

n≥0

(2.37)

en potencias de h εn+1 (y) = hYn0 +

h2 00 h3 000 Y + Yn + O(h4 ) − hϕ(xn , Yn , h; f ) 2 n 6

(2.38)

Estabilidad de la primera f´ ormula de Runge

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

27

donde (la notaci´ on fx significa ∂f /∂x y analogamente para las derivadas parciales de orden superior) y0 = f y 00 = fx + fy y 0 = fx + fy f y 000 = fxx + fxy f + (fyx + fyy f )f + fy (fx + fy f ) = fxx + 2fxy f + fyy f 2 + fy fx + fy2 f

(2.39) (2.40) (2.41) (2.42)

y ϕ(x, y, h; f ) = γ1 f (x, y) + γ2 (f (x, y) + h(αfx + βf fy ) 1 1 +h2 ( α2 fxx + αβfxy f + β 2 fyy f 2 )) + O(h3 ) (2.43) 2 2 Sustituyendo y agrupando t´erminos en potencias de h se obtiene 1 1 εn+1 = h(1 − γ1 − γ2 )f + h2 (( − γ2 α)fx + ( − γ2 β)f fy ) 2 2 1 1 1 1 3 1 2 +h (( − γ2 α )fxx + ( − γ2 αβ)f fxy + ( − γ2 β 2 )f 2 fyy 6 2 3 6 2 1 1 2 3 + fy fx + fy f ) + O(h ) (2.44) 6 6 todo evaluado en (xn , Yn ). Deseamos que εn+1 (y) converja a cero tan r´apido como se pueda. Aceptando f arbitrarias no podremos, en general, anular el coeficiente de h3 (¿porqu´e?). Anulando los de h y h2 obtenemos 1 1 εn+1 (y) = O(h3 ), γ1 + γ2 = 1, γ2 α = , γ 2 β = (2.45) 2 2 cuya soluci´ on general es γ2 arbitrario γ1 = 1 − γ2 1 α=β= 2γ2 Consideremos algunos casos particulares (a) γ2 = 1, lo que implica γ1 = 0, α = β = 12 h h , y + f (x, y)) 2 2 (b) γ2 = 12 , lo que implica γ1 = 21 , α = β = 1 ⇒

ϕ(x, y, h; f ) = f (x +

1 1 f (x, y) + f (x + h, y + hf (x, y)) 2 2 Notemos que si definimos

⇒

ϕ(x, y, h; f ) =

k1 = hf (x, y) k2 = hf (x + h, y + k1 )

(2.46) (2.47) (2.48)

(2.49)

(2.50)

(2.51) (2.52)

se tiene

1 (k1 + k2 ) 2 Runge not´ o que este m´etodo puede mejorarse definiendo hϕ =

k3 = hf (x + h, y + k2 )

(2.53)

(2.54)

Segunda f´ ormula de Runge

28

2. ECUACIONES DIFERENCIALES

a r Runge 2 r u (x n n+1 )







r Euler

 





 







  0

xn

xn+1 x  d

  Figura 4. Segunda  f´ormula de Runge    y tomando   1 hϕ = (k1 + k3 ) (2.55) 2 d que es la segunda f´ ormula de Runge (ver figura 4). (c) (primer orden) γ1 = 1, lo que implica γ2 = 0, se obtiene en consecuencia

y 6

ϕ(x, y, h; f ) = f (x, y)

(Euler)

(2.56)

♣ Las f´ ormulas de mayor orden involucran un ´algebra m´as complicada, para ello se suele proponer hϕ(x, y, h; f ) =

p X

γi ki

(2.57)

i=1

con k1 = hf (x, y) k2 = hf (x + α2 h, y + β21 k1 ) .. . i−1 X ki = hf (x + αi h, y + βij kj ) i = 1, . . . , p

(2.58) (2.59)

(2.60)

j=1

M´ etodo cl´ asico Runge-Kutta

de

Los coeficientes se pueden escoger de modo que los primeros t´erminos en el error de truncamiento local sean nulos. En consecuencia, por un procedimiento an´alogo al realizado en el caso de segundo orden se obtienen diversos m´etodos de orden superior de los que el m´ as utilizado es el llamado M´etodo de Runge-Kutta 1 yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) 6

(2.61)

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

29

con k1 = hf (xn , yn ) h k1 k2 = hf (xn + , yn + ) 2 2 h k2 k3 = hf (xn + , yn + ) 2 2 k4 = hf (xn + h, yn + k3 )

y 6

0

(2.62) (2.63) (2.64) (2.65)

r RK

r un (xn+1 )









r Euler











xn

xn+1 x

Figura 5. El m´etodo de Runge-Kutta Se puede demostrar que ´esta es una f´ormula de cuarto orden: εn+1 = O(h5 ). ´ n 2.2.1. El m´etodo RK cl´asico es una generalizaci´on de la regla Observacio de Simpson en el caso que f no sea independiente de y Z xn+1 h h f (x)dx ' (f (xn ) + 4f (xn + + f (xn+1 )) (2.66) 6 2 xn Ejemplo 2.2.6. Apliquemos el m´etodo cl´asico de RK a la ecuaci´on de prueba y 0 = λy k1 = hλyn 1 k2 = hλ(yn + hλyn ) 2 1 = hλyn + h2 λ2 yn 2 1 1 k3 = hλ(yn + (hλyn + h2 λ2 yn )) 2 2 1 2 2 1 = hλyn + h λ yn + h3 λ3 yn 2 4 1 2 2 1 k4 = hλ(yn + hλyn + h λ yn + 1 h3 λ3 yn ) 2 4 1 3 3 1 2 2 = hλyn + h λ yn + h λ yn + h4 λ4 yn ) 2 4

(2.67) (2.68) (2.69) (2.70) (2.71) (2.72) (2.73)

El m´ etodo cl´ asico de Runge-Kutta generaliza el m´ etodo de Simpson

30

2. ECUACIONES DIFERENCIALES

de manera que 1 yn + (k1 + 2k2 + 2k3 + k4 ) = yn 6 1 1 + (6hλyn + 3h2 λ2 yn + h3 λ3 yn + h4 λ4 yn ) 6 4 lo que permite escribir en definitiva 1 1 1 yn+1 = yn (1 + hλ + h2 λ2 + h3 λ3 + h4 λ4 ) 2 6 24

(2.74)

(2.75)

vemos as´ı que se aproxima a ehλ hasta el cuarto orden. Ejercicio 2.2.7. [Optativo] Con la expresi´on 2.75 se puede dibujar la regi´on de estabilidad absoluta. ♣ Para estudiar la convergencia del esquema general RK yn+1 = yn + hϕ(xn , yn , h; f ),

n≥0

(2.76)

recordemos que el error local es εn+1 (Y ) = Y (xn+1 ) − Y (xn ) − hϕ(xn , Y (xn ), h; f ),

n≥0

(2.77)

conviene definir τn+1 (Y ) por la expresi´on εn+1 (y) = hτn+1 (Y )

(2.78)

y con esta τ se tiene Y (xn+1 ) = Y (xn ) + hϕ(xn , Y (xn ), h; f ) + hτn+1 (Y ),

n≥0

(2.79)

con lo que resulta Y (xn+1 ) − Y (xn ) − ϕ(xn , Y (xn ), h; f ) (2.80) h de modo que debe tenerse τn+1 (Y ) → 0 cuando h → 0, equivalentemente debe requerirse que ϕ(x, Y (x), h; f ) → Y 0 (x) = f (x, Y (x)) cuando h → 0. Definiendo δ(h) = max |f (x, y) − ϕ(x, y, h; f )| (2.81) τn+1 (Y ) =

x0 ≤x≤b −∞
se tiene que limh→0 δ(h) = 0 se denomina condici´on de consistencia para el m´etodo 2.76 y es suficiente para asegurar la convergencia si se tienen ciertas hip´otesis para ϕ, a saber: ϕ continua y |ϕ(x, y, h; f ) − ϕ(x, y˜, h; f )| ≤ L|y − y˜|,

x0 ≤ x ≤ b, −∞ < y, y˜ < ∞

(2.82)

(funci´ on de Lipschitz en la variable y) Teorema 2.2.2. Supongamos que el m´etodo 2.76 dado por yn+1 = yn + hϕ(xn , yn , h; f )

(2.83)

satisface la condici´ on |ϕ(x, y, h; f ) − ϕ(x, y˜, h; f )| ≤ L|y − y˜| con x0 ≤ x ≤ b − ∞ < y, y˜ < ∞, L constante. Entonces la soluci´on aproximada {yn } del problema de valores iniciales y 0 = f (x, y), y(x0 ) = y0 satisface max |Y (xn ) − yn | ≤ e(b−x0 )L |Y0 − y0 | +

x0 ≤xn ≤b

e(b−x0 )L − 1 τ (h) L

(2.84)

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

31

donde τ (h) = maxx0 ≤xn ≤b |τn+1 (Y )|. Si se satisface la condici´on de consistencia δ(h) =

max

x0 ≤x≤b −∞
|f (x, y) − ϕ(x, y, h; f )| −→ 0

para

h→0

(2.85)

entonces la soluci´ on num´erica {yn } converge a Y (x) ´ n: Restamos (2.79)−(2.76) Demostracio Y (xn+1 ) − yn+1 = Y (xn ) − yn + h(ϕ(xn , Y (xn ), h; f ) −ϕ(xn , yn , h; f )) + hτn+1 (Y )

(2.86)

Luego |Y (xn+1 ) − yn+1 | ≤ |Y (xn ) − yn | + h|ϕ(xn , Y (xn ), h; f ) − ϕ(xn , yn , h; f )| +h|τn+1 (Y )| ≤ |Y (xn ) − yn |(1 + hL) + hτ (h), x0 ≤ xn ≤ b (2.87) Recursivamente se tiene |Y (xn+1 ) − yn+1 | ≤ (1 + hL)n |Y (x0 ) − y0 | + (1 + (1 + hL) + . . . + (1 + hL)n−1 )hτ (h) (1 + hL)n − 1 ≤ (1 + hL)n |Y (x0 ) − y0 | + τ (h) L

(2.88) 2

pero (1 + hL)n ≤ enhL = e(xn −x0 )L ≤ e(b−x0 )L pues ex = 1 + x + x2 eξ con 0 ≤ ξ ≤ x y x > −1 ⇒ (1 + x)m ≤ emx . Luego |Y (xn ) − yn | ≤ e(b−x0 )L |Y (x0 ) − y0 | +

e(b−x0 )L − 1 τ (h), L

x0 ≤ xn ≤ b

(2.89)

En muchos casos se puede demostrar que τ (h) → 0 cuando h → 0 por c´alculo directo obteniendose as´ı la convergencia buscada. Veamos que basta saber que limh→0 δ(h) = 0, hτn+1 (Y ) = Y (xn+1 ) − Y (xn ) − hϕ(xn , Y (xn ), h; f ) 1 = hY 0 (n) + h2 Y 00 (ξn ) − hϕ(xn , Y (xn ), h; f ), 2 con xn < ξn < xn+1

(2.90)

resulta as´ı h|τn+1 (Y )| ≤ hδ(h) +

h2 00 ||Y ||∞ 2

(2.91)

obteniendose en definitiva 1 τ (h) ≤ δ(h) + h||Y 00 ||∞ 2

(2.92)

lo que completa la demostraci´on. Corolario 2.2.3. Si el m´etodo RK dado por yn+1 = yn +hϕ(xn , yn , h; f ) tiene un error de truncamiento εn+1 (y) = O(hm+1 ) entonces la velocidad de convergencia de {yn } a y(x) es O(hm ).

32

2. ECUACIONES DIFERENCIALES

Teorema 2.2.4. Si ϕ(x, y, h; f ) es continua en (x, y, h) en x0 ≤ x ≤ b, 0 ≤ h ≤ h0 , para todo y y es Lipschitz en y, entonces la convergencia es equivalente a tener ϕ(x, Y (x), 0; f ) = f (x, Y (x)) o, equivalentemente ϕ(x, y, 0; f ) = f (x, y)

Ejercicio 2.2.8. Demostrar el teorema 2.2.4 La velocidad de convergencia de yh (x) a Y (x) es O(h4 ) en el m´etodo cl´asico RK pues εn+1 (y) = O(h5 ). Se puede demostrar, en consecuencia, que Y (x) − yh (x) = D(x)h4 + O(h5 )

(2.93)

donde D(x) satisface un problema an´alogo de valores iniciales. Asimismo, como Y (x) − y2h (x) = 16D(x)h4 + O(h5 )

(2.94)

se deduce que 1 (yh (x) − y2h (x)) + O(h5 ) 15 donde el t´ermino yh (x) − y2h (x) permite estimar el error. En general vale el siguiente

(2.95)

Y (x) = yh (x) +

Teorema 2.2.5. Si el error de truncamiento εn+1 (Y ) = hϕ(x, Y (x), h; f ) − (Y (x + h) − Y (x))

(2.96)

se puede expresar εn+1 (Y ) = hr+1 ψ(x, Y ) + O(hr+2 )

(2.97)

y ϕ tiene derivadas segundas continuas, entonces el error satisface en = hr D(xn ) + O(hr+1 )

(2.98)

donde D(x) es soluci´ on de D0 (x) =

∂f (x, Y (x))D(x) + ψ(x, Y (x)) ∂y

(2.99)

´ n: Por el orden del m´etodo buscamos para en una expresi´on de la Demostracio forma del enunciado. Si reemplazamos la expresi´on en = hr δn en en+1 = en + h(ϕ(xn , yn , h; f ) − ϕ(xn , Y (xn ), h; f )) + εn+1 (Y )

(2.100)

obtenemos δn+1 = δn + h1−r (ϕ(xn , Y (xn ) + hr δn , h; f ) − ϕ(xn , Y (xn ), h; f )) +hψ(xn , Y (xn )) + O(h2 )

(2.101)

Por la suavidad supuesta para ϕ podemos escribir ϕ(xn , Y (xn ) + hr δn , h; f ) = ϕ(xn , Y (xn ), h; f ) + ϕy (xn , Y (xn ), h; f )hr δn 1 + ϕyy (xn , Y (xn ), h; f )h2r δn2 (2.102) 2

2.3. EJERCICIOS

33

en raz´ on de la convergencia este u ´ltimo sumando se puede escribir en la forma k1 h2r con ||k1 || acotada. El segundo sumando se puede reescribir as´ı ϕy (xn , Y (xn ), h; f )hr δn = ϕy (xn , Y (xn ), 0; f )hr δn +ϕyh (xn , Y (xn ), ξ; f )hr+1 δn

(2.103)

cuyo u ´ltimo sumando se puede escribir en la forma k2 hr+1 con ||k2 || acotada. En definitiva δn+1 = δn + h(fy (xn , Y (xn ))δn + hψ(xn , Y (xn ) + hk2 + hr k1 )

(2.104)

que es un m´etodo num´erico para D0 (x) = fy (x, Y (x))D(x) + ϕ(x, Y (x))

(2.105)

D(0) + e0 h−r

(2.106)

lo que permite completar la demostraci´on del teorema.

♣

2.3. Ejercicios Ejercicio 2.3.1. Dada y 0 = 1 + x2 y 2 , con y(0) = 0, calcular y(0.5) mediante el m´etodo de Euler con extrapolaci´on de repetida de Richardson. Usar aritm´etica de redondeo de cinco decimales. Ejercicio 2.3.2. Aplicando el m´etodo de Euler, se obtuvieron los resultados: 1.22726 (con paso h = 0.05), 1.22595 (h = 0.1), 1.22345 (h = 0.2). Calcule un mejor valor por extrapolaci´ on. Ejercicio 2.3.3. Utilizar el m´etodo de Euler con longitud de paso h sobre el problema test y 0 = −y, y(0) = 1 (2.107) (a) Determinar una expresi´on expl´ıcita para yn (b) ¿Para qu´e valores de h la sucesi´on {yn }∞ 0 es acotada? (c) Calcular el limh→0 (y(x, h) − e−x )/h Ejercicio 2.3.4. Determinar el desarrollo de Taylor para la soluci´on de la ecuaci´ on y 0 = y 2 , con y(0) = 1, en el entorno de x = 0. Usar esta aproximaci´on para calcular y(0.2) e y(1.2) con cuatro decimales. Comparar con la soluci´on exacta y explicar porqu´e el segundo caso (x = 1.2) no es exitoso. Ejercicio 2.3.5. La funci´on y(x) se define mediante el problema y 0 = x2 − y 2 ,

y(0) = 1

(2.108)

Calcular y(0.2) utilizando los siguientes m´etodos: (a) (b) (c) (d)

Euler-Richardson: h = 0.1 y h = 0.2; M´etodo Runge-Kutta: h = 0.1; Desarrollo en serie de Taylor hasta el cuarto t´ermino; M´etodo del trapecio: h = 0.2.

Ejercicio 2.3.6. [Optativo]Se desea realizar una tabla de la funci´on Z ∞ −u2 e y(x) = du (u + x) 0

(2.109)

34

2. ECUACIONES DIFERENCIALES

para varios valores de x. Se procede de la siguiente manera: y(x) se calcula para x = 1 usando alg´ un m´etodo para integraci´on num´erica; se calcula y(1) = 0.6051. Mostrar que y satisface la ecuaci´on diferencial dy 1 √ + 2xy = − + π (2.110) dx x Se resuelve la ecuaci´ on diferencial num´ericamente con el valor inicial y(1) = 0.6051, y as´ı se pueden calcular m´ as valores de la tabla. Determinar y(x) para x = 1.2 y x = 1.4 mediante el m´etodo de Runge-Kutta y estimar el error en y(1.4).

´ APENDICE A

Bibliograf´ıa y referencias A.1. Textos b´ asicos (1) Atkinson, K.E.: An Introduction to Numerical Analysis, 2nd edition, J. Wiley & Sons., New York, 1989. (2) Burden, R.L. y Faires, J.D.: An´ alisis Num´erico, 2da edici´on, Grupo Editorial Iberoam´erica, M´exico, 1993. (3) Dahlquist, G., and Bjrck, ˚ A.: Numerical Methods, Prentice-Hall, E. Cliffs, NJ, 1974.

A.2. Bibliograf´ıa complementaria (1) Aguilera, N.E.: Introducci´ on a la computaci´ on en Matem´ atica usando Mathematica, Red Ol´ımpica, OMA, Buenos Aires, 1994. (2) Apostol, T.M.: An´ alisis Matem´ atico, 2da edici´on, Revert´e, Barcelona, 1976. (3) Braun, M.: Differential Equations and Their Applications. An Introduction to Applied Mathematics, 4th edition, Springer, N.York, 1993. (4) Davis, Ph.J. y Rabinowitz, Ph.: Methods of Numerical Integration, 2nd edition, Academic Press, N.York, 1984. (5) Gear, C.W.: Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, E. Cliffs, 1971 (6) Golub, G.H. (Editor): Studies in Numerical Analysis, MAA Studies in Mathematics 24, Washington, 1984. (7) Graham, R.L., Knuth, D.E., and Patashnik, O.: Concrete Mathematics, Addison-Wesley, Reading, MA, 1989. (8) Hurewicz, W.: Lectures on Ordinary Differential Equations, The M.I.T. Press, Cambridge, MA, 1958. (9) Knuth, D.E.: The Art of Computer Programming Vol.1: Fundamental Algorithms Vol.2: Seminumerical algorithms Vol.3: Sorting and searching , Addison-Wesley, Reading, 1973. Versi´ on espa˜ nola: V.1, Revert´e, Barcelona, 1980. (10) Linz, P.: A Critique of Numerical Analysis, Bulletin American Mathematical Society, 19(2), 407–416, 1988 (11) Maeder, R.: Programming in Mathematica, 2nd edition, Addison-Wesley, 1991. (12) Marshall, G.: Soluci´ on Num´erica de Ecuaciones Diferenciales, Tomos 1 y 2, Revert´e, Buenos Aires, 1985 y 1986 35

36

A. BIBLIOGRAF´IA Y REFERENCIAS

(13) Press, W.H., Teukolsky, S.A., Vetterlin, W.T., y Flannery, B.P.: Numerical Recipes in C. The Art of Scientific Computing, 2nd edition, Cambridge University Press, Cambridge, 1992. (14) Rice, J.R.: Numerical Methods, Software, and Analysis, 2nd edition, Academic Press, Boston, 1993. (15) Stoer, J. y Bulirsch, R.: Introduction to Numerical Analysis, Springer, N.York, 1980. (16) Strang, G.: Introduction to Applied Mathematics, Wellesley-Cambridge, Wellesley, 1986 (17) Wolfram Research: Guide to Standard Mathematica Packages, Version 2.2, Technical Report, Wolfram Research, 1993. (18) Wolfram Research: MathLink Reference Guide, Version 2.2, Technical Report, Wolfram Research, 1993. (19) Wolfram, S.: The Mathematica Book , 3rd edition, Wolfram MediaCambridge, Cambridge, UK, 1996.

´ APENDICE B

Coeficientes indeterminados con el Mathematica B.1. El m´ etodo de Simpson (* Simpson *) a0=a;(* Defino los puntos *) a1=a/2+b/2; a2=b;(*a=-1;b=1;*) p0[x_]:=1; (* Defino los polinomios *) p1[x_]:=(2x-a-b)/(b-a); p2[x_]:=((2x-a-b)/(b-a))^2; p3[x_]:=((2x-a-b)/(b-a))^3; p4[x_]:=((2x-a-b)/(b-a))^4; (* Integraci\’on num\’erica *) s0=A p0[a0] + B p0[a1] + C p0[a2] //Simplify s1=A p1[a0] + B p1[a1] + C p1[a2] //Simplify; s2=A p2[a0] + B p2[a1] + C p2[a2] //Simplify; A + B + C (* Integro los polinomios *) i0=Integrate[p0[x],{x,a,b}]//Simplify i1=Integrate[p1[x],{x,a,b}]//Simplify; i2=Integrate[p2[x],{x,a,b}]//Simplify; -a + b (* Resuelvo el sistema *) Solve[{s0==i0,s1==i1,s2==i2},{A,B,C}]//Simplify 2 (-a + b) -a + b -a + b {{B -> ----------, A -> ------, C -> ------}} 3 6 6 (* Verifico que vale para polinomios de tercer grado *) (s3=A p3[a0] + B p3[a1] + C p3[a2])/.%15 //Simplify {0} i3=Integrate[p3[x],{x,a,b}]//Simplify 0 (* Resulta as\’{\i} que la siguiente ecuaci\’on ( s3 == i3 ) vale trivialmente *) (* Calculo del error con el polinomio de cuarto grado *) (s4=A p4[a0] + B p4[a1] + C p4[a2])/.%15 //Simplify -a + b {------} 3 i4=Integrate[p4[x],{x,a,b}]//Simplify -a + b -----5 e4=(%20 - %19)//Simplify 2 (a - b) 37

38

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

{---------} 15 (* Calculo del polinomio de Taylor *) Series[ f[x], {x, (a+b)/2, 4}] a + b -(a + b) 2 f’’[-----] (-------- + x) a + b a + b -(a + b) 2 2 f[-----] + f’[-----] (-------- + x) + -------------------------- + 2 2 2 2 (3) a + b -(a + b) 3 (4) a + b -(a + b) 4 [-----] (-------- + x) f [-----] (-------- + x) 2 2 2 2 --------------------------- + --------------------------- + 6 24 f

-(a + b) 5 O[-------- + x] 2 (* Resulta asi que el error de truncamiento es 5 (4) a + b ( b - a ) f [-----] 2 -----------------------2880

*)

B.2. El m´ etodo de Newton (* Newton de once puntos *) a0=a; (* Defino los puntos *) a1=9a/10+b/10; a2=8a/10+2b/10; a3=7a/10+3b/10; a4=6a/10+4b/10; a5=5a/10+5b/10; a6=4a/10+6b/10; a7=3a/10+7b/10; a8=2a/10+8b/10; a9=a/10+9b/10; a10=b; p0[x_]:=1; (* Defino los polinomios *) p1[x_]:=(2x-a-b)/(b-a); p2[x_]:=((2x-a-b)/(b-a))^2; p3[x_]:=((2x-a-b)/(b-a))^3; p4[x_]:=((2x-a-b)/(b-a))^4; p5[x_]:=((2x-a-b)/(b-a))^5; p6[x_]:=((2x-a-b)/(b-a))^6; p7[x_]:=((2x-a-b)/(b-a))^7; p8[x_]:=((2x-a-b)/(b-a))^8; p9[x_]:=((2x-a-b)/(b-a))^9; p10[x_]:=((2x-a-b)/(b-a))^10; p11[x_]:=((2x-a-b)/(b-a))^11; p12[x_]:=((2x-a-b)/(b-a))^12;

´ B.2. EL METODO DE NEWTON

(* Integraci\’on num\’erica *) s0=A0 p0[a0] + A1 p0[a1] + A2 p0[a2] + A3 p0[a3]+ A4 p0[a4] + A5 p0[a5] + A6 p0[a6] + A7 p0[a7]+ A8 p0[a8] + A9 p0[a9] + A10 p0[a10]//Simplify; s1=A0 p1[a0] + A1 p1[a1] + A2 p1[a2] + A3 p1[a3]+ A4 p1[a4] + A5 p1[a5] + A6 p1[a6] + A7 p1[a7]+ A8 p1[a8] + A9 p1[a9] + A10 p1[a10]//Simplify; s2=A0 p2[a0] + A1 p2[a1] + A2 p2[a2] + A3 p2[a3]+ A4 p2[a4] + A5 p2[a5] + A6 p2[a6] + A7 p2[a7]+ A8 p2[a8] + A9 p2[a9] + A10 p2[a10]//Simplify; s3=A0 p3[a0] + A1 p3[a1] + A2 p3[a2] + A3 p3[a3]+ A4 p3[a4] + A5 p3[a5] + A6 p3[a6] + A7 p3[a7]+ A8 p3[a8] + A9 p3[a9] + A10 p3[a10]//Simplify; s4=A0 p4[a0] + A1 p4[a1] + A2 p4[a2] + A3 p4[a3]+ A4 p4[a4] + A5 p4[a5] + A6 p4[a6] + A7 p4[a7]+ A8 p4[a8] + A9 p4[a9] + A10 p4[a10]//Simplify; s5=A0 p5[a0] + A1 p5[a1] + A2 p5[a2] + A3 p5[a3]+ A4 p5[a4] + A5 p5[a5] + A6 p5[a6] + A7 p5[a7]+ A8 p5[a8] + A9 p5[a9] + A10 p5[a10]//Simplify; s6=A0 p6[a0] + A1 p6[a1] + A2 p6[a2] + A3 p6[a3]+ A4 p6[a4] + A5 p6[a5] + A6 p6[a6] + A7 p6[a7]+ A8 p6[a8] + A9 p6[a9] + A10 p6[a10]//Simplify; s7=A0 p7[a0] + A1 p7[a1] + A2 p7[a2] + A3 p7[a3]+ A4 p7[a4] + A5 p7[a5] + A6 p7[a6] + A7 p7[a7]+ A8 p7[a8] + A9 p7[a9] + A10 p7[a10]//Simplify; s8=A0 p8[a0] + A1 p8[a1] + A2 p8[a2] + A3 p8[a3]+ A4 p8[a4] + A5 p8[a5] + A6 p8[a6] + A7 p8[a7]+ A8 p8[a8] + A9 p8[a9] + A10 p8[a10]//Simplify; s9=A0 p9[a0] + A1 p9[a1] + A2 p9[a2] + A3 p9[a3]+ A4 p9[a4] + A5 p9[a5] + A6 p9[a6] + A7 p9[a7]+ A8 p9[a8] + A9 p9[a9] + A10 p9[a10]//Simplify; s10=A0 p10[a0] + A1 p10[a1] + A2 p10[a2] + A3 p10[a3]+ A4 p10[a4] + A5 p10[a5] + A6 p10[a6] + A7 p10[a7]+ A8 p10[a8] + A9 p10[a9] + A10 p10[a10]//Simplify; (* Integro los polinomios *) i0=Integrate[p0[x],{x,a,b}]//Simplify i1=Integrate[p1[x],{x,a,b}]//Simplify; i2=Integrate[p2[x],{x,a,b}]//Simplify; i3=Integrate[p3[x],{x,a,b}]//Simplify; i4=Integrate[p4[x],{x,a,b}]//Simplify; i5=Integrate[p5[x],{x,a,b}]//Simplify; i6=Integrate[p6[x],{x,a,b}]//Simplify; i7=Integrate[p7[x],{x,a,b}]//Simplify; i8=Integrate[p8[x],{x,a,b}]//Simplify; i9=Integrate[p9[x],{x,a,b}]//Simplify; i10=Integrate[p10[x],{x,a,b}]//Simplify; -a + b (* Resuelvo el sistema *) Solve[{s0==i0,s1==i1,s2==i2 ,s3==i3,s4==i4 ,s5==i5,s6==i6 ,s7==i7,s8==i8 ,s9==i9,s10==i10} ,{A0,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10}]//Simplify 17807 (-a + b) -16067 (a - b) {{A5 -> --------------, A0 -> --------------,

39

40

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

24948 598752 -26575 (a - b) 16175 (a - b) A1 -> --------------, A2 -> -------------, 149688 199584 -5675 (a - b) 4825 (a - b) A3 -> -------------, A4 -> ------------, 12474 11088 4825 (a - b) -5675 (a - b) A6 -> ------------, A7 -> -------------, 11088 12474 16175 (a - b) -26575 (a - b) A8 -> -------------, A9 -> --------------, 199584 149688 -16067 (a - b) A10 -> --------------}} 598752 (* Verifico que los pesos suman uno *) 598752/24948 24 598752/149688 4 598752/199584 3 598752/12474 48 598752/11088 54 2 16067 + 2 26575 4 - 2 16175 3 + 2 5675 48 - 2 4825 54 + 17807 24 598752

B.3. El m´ etodo de cuadratura de Gauß de tres puntos (* Gauss de tres puntos *) a=-1; (* Se trabaja en el [-1,1] *) b=1; p0[x_]:=1; (* Definici\’on de los polinomios *) p1[x_]:=(2x-a-b)/(b-a); p2[x_]:=((2x-a-b)/(b-a))^2; p3[x_]:=((2x-a-b)/(b-a))^3; p4[x_]:=((2x-a-b)/(b-a))^4; p5[x_]:=((2x-a-b)/(b-a))^5; p6[x_]:=((2x-a-b)/(b-a))^6; p7[x_]:=((2x-a-b)/(b-a))^7; (* Las expresiones de Gauss *) s0=A p0[X] + B p0[Y] + C p0[Z] s1=A p1[X] + B p1[Y] + C p1[Z] s2=A p2[X] + B p2[Y] + C p2[Z] s3=A p3[X] + B p3[Y] + C p3[Z] s4=A p4[X] + B p4[Y] + C p4[Z] s5=A p5[X] + B p5[Y] + C p5[Z] s6=A p6[X] + B p6[Y] + C p6[Z] A + B + C

//Simplify //Simplify; //Simplify; //Simplify; //Simplify; //Simplify; //Simplify;

(* Las integrales *) i0=Integrate[p0[x],{x,a,b}]//Simplify i1=Integrate[p1[x],{x,a,b}]//Simplify; i2=Integrate[p2[x],{x,a,b}]//Simplify;

B.4. LOS COEFICIENTES DE LA CUADRATURA DE GAUSS DE CINCO PUNTOS

i3=Integrate[p3[x],{x,a,b}]//Simplify; i4=Integrate[p4[x],{x,a,b}]//Simplify; i5=Integrate[p5[x],{x,a,b}]//Simplify; 2 (* El sistema de ecuaciones no lineal *) Solve[{s0==i0, s1==i1, s2==i2, s3==i3, s4==i4, s5==i5},{A,B,C,X,Y,Z}]//Simplify 5 5 8 3 3 {{A -> -, B -> -, C -> -, Z -> 0, Y -> -Sqrt[-], X -> Sqrt[-]}, 9 9 9 5 5 5 5 8 3 3 {A -> -, B -> -, C -> -, Z -> 0, Y -> Sqrt[-], X -> -Sqrt[-]}, 9 9 9 5 5 5 8 5 3 3 {A -> -, B -> -, C -> -, Z -> -Sqrt[-], Y -> 0, X -> Sqrt[-]}, 9 9 9 5 5 5 8 5 3 3 {A -> -, B -> -, C -> -, Z -> Sqrt[-], Y -> 0, X -> -Sqrt[-]}, 9 9 9 5 5 8 5 5 3 3 {A -> -, B -> -, C -> -, Z -> -Sqrt[-], Y -> Sqrt[-], X -> 0}, 9 9 9 5 5 8 5 5 3 3 {A -> -, B -> -, C -> -, Z -> Sqrt[-], Y -> -Sqrt[-], X -> 0}} 9 9 9 5 5 %//N (* La versi\’on num\’erica *) {{A -> 0.555556, B -> 0.555556, C -> 0.888889, Z -> 0, Y -> -0.774597, X -> 0.774597}, {A -> 0.555556, B -> 0.555556, C -> 0.888889, Z -> 0, Y -> 0.774597, X -> -0.774597}, {A -> 0.555556, B -> 0.888889, C -> 0.555556, Z -> -0.774597, Y -> 0, X -> 0.774597}, {A -> 0.555556, B -> 0.888889, C -> 0.555556, Z -> 0.774597, Y -> 0, X -> -0.774597}, {A -> 0.888889, B -> 0.555556, C -> 0.555556, Z -> -0.774597, Y -> 0.774597, X -> 0}, {A -> 0.888889, B -> 0.555556, C -> 0.555556, Z -> 0.774597, Y -> -0.774597, X -> 0}}

B.4. Los coeficientes de la cuadratura de Gauß de cinco puntos (* Cuadratura de Gauss de cinco puntos *)

41

42

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

(* Determino las coordenadas de los puntos *) Roots[LegendreP[5,x]==0,x] N[%,30] (* Defino las variables que representan los puntos *) a0=x/.x->%%[[5]][[2]] a1=x/.x->%%%[[3]][[2]] a2=x/.x->%%%%[[1]][[2]] a3=x/.x->%%%%%[[2]][[2]] a4=x/.x->%%%%%%[[4]][[2]] 10 10 Sqrt[5 - 2 Sqrt[--]] -Sqrt[5 - 2 Sqrt[--]] 7 7 x == 0 || x == -------------------- || x == --------------------- || 3 3 10 10 Sqrt[5 + 2 Sqrt[--]] -Sqrt[5 + 2 Sqrt[--]] 7 7 x == -------------------- || x == --------------------3 3 x == 0 || x == 0.5384693101056830910363144207 || x == -0.5384693101056830910363144207 || x == 0.906179845938663992797626878299 || x == -0.906179845938663992797626878299 -0.906179845938663992797626878299 -0.5384693101056830910363144207 0 0.5384693101056830910363144207 0.906179845938663992797626878299 (* Fijo los extremos del intervalo *) a=-1; b=1; (* Defino los polinomios *) p0[x_]:=1; p1[x_]:=(2x-a-b)/(b-a); p2[x_]:=((2x-a-b)/(b-a))^2; p3[x_]:=((2x-a-b)/(b-a))^3; p4[x_]:=((2x-a-b)/(b-a))^4; p5[x_]:=((2x-a-b)/(b-a))^5; p6[x_]:=((2x-a-b)/(b-a))^6; p7[x_]:=((2x-a-b)/(b-a))^7; p8[x_]:=((2x-a-b)/(b-a))^8; p9[x_]:=((2x-a-b)/(b-a))^9; p10[x_]:=((2x-a-b)/(b-a))^10; p11[x_]:=((2x-a-b)/(b-a))^11; p12[x_]:=((2x-a-b)/(b-a))^12; (* Defino las expresiones para los pesos A0,A1,A2,A3,A4 *) s0=A0 p0[a0]+A1 p0[a1]+A2 p0[a2]+A3 p0[a3]+A4 p0[a4]//Simplify s1=A0 p1[a0]+A1 p1[a1]+A2 p1[a2]+A3 p1[a3]+A4 p1[a4]//Simplify

B.4. LOS COEFICIENTES DE LA CUADRATURA DE GAUSS DE CINCO PUNTOS

s2=A0 p2[a0]+A1 p2[a1]+A2 p2[a2]+A3 p2[a3]+A4 p2[a4]//Simplify; s3=A0 p3[a0]+A1 p3[a1]+A2 p3[a2]+A3 p3[a3]+A4 p3[a4]//Simplify; s4=A0 p4[a0]+A1 p4[a1]+A2 p4[a2]+A3 p4[a3]+A4 p4[a4]//Simplify; s5=A0 p5[a0]+A1 p5[a1]+A2 p5[a2]+A3 p5[a3]+A4 p5[a4]//Simplify; s6=A0 p6[a0]+A1 p6[a1]+A2 p6[a2]+A3 p6[a3]+A4 p6[a4]//Simplify; s7=A0 p7[a0]+A1 p7[a1]+A2 p7[a2]+A3 p7[a3]+A4 p7[a4]//Simplify; s8=A0 p8[a0]+A1 p8[a1]+A2 p8[a2]+A3 p8[a3]+A4 p8[a4]//Simplify; s9=A0 p9[a0]+A1 p9[a1]+A2 p9[a2]+A3 p9[a3]+A4 p9[a4]//Simplify; s10=A0 p10[a0]+A1 p10[a1]+A2 p10[a2]+ A3 p10[a3]+A4 p10[a4]//Simplify; s11=A0 p11[a0]+A1 p11[a1]+A2 p11[a2]+ A3 p11[a3]+A4 p11[a4]//Simplify; A0 + A1 + A2 + A3 + A4 -0.906179845938663992797626878299 A0 - 0.5384693101056830910363144207 A1 + 0.5384693101056830910363144207 A3 + 0.906179845938663992797626878299 A4 (* Calculo las integrales *) i0=Integrate[p0[x],{x,a,b}]//Simplify i1=Integrate[p1[x],{x,a,b}]//Simplify; i2=Integrate[p2[x],{x,a,b}]//Simplify; i3=Integrate[p3[x],{x,a,b}]//Simplify; i4=Integrate[p4[x],{x,a,b}]//Simplify; i5=Integrate[p5[x],{x,a,b}]//Simplify; i6=Integrate[p6[x],{x,a,b}]//Simplify; i7=Integrate[p7[x],{x,a,b}]//Simplify; i8=Integrate[p8[x],{x,a,b}]//Simplify; i9=Integrate[p9[x],{x,a,b}]//Simplify; i10=Integrate[p10[x],{x,a,b}]//Simplify; i11=Integrate[p11[x],{x,a,b}]//Simplify; 2 (* Resoluci\’on del sistema de ecuaciones *) Solve[{s0==i0 ,s1==i1 ,s2==i2 ,s3==i3 ,s4==i4},{A0,A1,A2,A3,A4}]//Simplify {{A2 -> 0.56888888888888888888888889 , A0 -> 0.236926885056189087514264041 , A1 -> 0.478628670499366468041291515 , A3 -> 0.478628670499366468041291515 , A4 -> 0.2369268850561890875142640407}} (* Calculo de los residuos para los siguientes polinomios *) s5 - i5 /. %56 -16 {-2.22045 10 } s6 - i6 /. %56 -16 {4.16334 10 } s7 - i7 /. %56 -16 {-1.66533 10 } s8 - i8 /. %56 -16 {3.46945 10 } s9 - i9 /. %56

43

44

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

-16 {-1.38778 10 } s10 - i10 /. %56 {-0.00293181} s11 - i11 /. %56 -16 {-1.249 10 } s12=A0 p12[a0]+A1 p12[a1]+A2 p12[a2]+ A3 p12[a3]+A4 p12[a4]//Simplify; i12=Integrate[p11[x],{x,a,b}]//Simplify; s12 - i12 /. %56 {0.145853}

En algunos casos se trabaja con mayor precisi´on (* Cuadratura de Gauss de cinco puntos (operaciones exactas) *) (* Determino las coordenadas de los puntos *) Roots[LegendreP[5,x]==0,x] N[%,50] (* Defino las variables que representan los puntos *) a0=x/.x->%%%[[5]][[2]] a1=x/.x->%%%%[[3]][[2]] a2=x/.x->%%%%%[[1]][[2]] a3=x/.x->%%%%%%[[2]][[2]] a4=x/.x->%%%%%%%[[4]][[2]] 10 10 Sqrt[5 - 2 Sqrt[--]] -Sqrt[5 - 2 Sqrt[--]] 7 7 x == 0 || x == -------------------- || x == --------------------- || 3 3 10 10 Sqrt[5 + 2 Sqrt[--]] -Sqrt[5 + 2 Sqrt[--]] 7 7 x == -------------------- || x == --------------------3 3 x == 0 || x == 0.53846931010568309103631442070020880496728660690556 || x == -0.53846931010568309103631442070020880496728660690556 || x == 0.90617984593866399279762687829939296512565191076253 || x == -0.90617984593866399279762687829939296512565191076253 10 -Sqrt[5 + 2 Sqrt[--]] 7 --------------------3 10 -Sqrt[5 - 2 Sqrt[--]] 7 --------------------3 0 10 Sqrt[5 - 2 Sqrt[--]]

B.4. LOS COEFICIENTES DE LA CUADRATURA DE GAUSS DE CINCO PUNTOS

7 -------------------3 10 Sqrt[5 + 2 Sqrt[--]] 7 -------------------3

. . .

(* Resoluci\’on del sistema de ecuaciones *) Solve[{s0==i0 ,s1==i1 ,s2==i2 ,s3==i3 ,s4==i4},{A0,A1,A2,A3,A4}]//Simplify; N[%,50] {{A2 -> 0.56888888888888888888888888888888888888888888888889, A0 -> 0.23692688505618908751426404071991736264326000221241, A1 -> 0.47862867049936646804129151483563819291229555334314, A3 -> 0.4786286704993664680412915148356381929122955533431, A4 -> 0.23692688505618908751426404071991736264326000221241}} s5 - i5 /. %61; N[%,35] -43 {0. 10 } s6 - i6 /. %61; N[%,35] -42 {0. 10 } s7 - i7 /. %61; N[%,35] -43 {0. 10 } s8 - i8 /. %61; N[%,35] -43 {0. 10 } s9 - i9 /. %61; N[%,35] -43 {0. 10 } s10 - i10 /. %61; N[%,35] {-0.00293181245562197943150324102705055} s11 - i11 /. %61; N[%,35] -43 {0. 10 }

45

46

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

s12 - i12 /. %61; N[%,35] {0.14585257971501357744743988130231516} s13 - i13 /. %61; N[%,35] -44 {0. 10 } s14 - i14 /. %61; N[%,35] {-0.01386690413313408190371321302706202}

´ APENDICE C

Cuadernos del Mathematica y archivos script de Matlab C.1. La funci´ on mihumps (*MIHUMPS*) f[x_]:= 0.01/((x-0.3)^2+.01) + 0.01/((x-0.9)^2+.04) - 0.06 Plot[f[x],{x,0,1},PlotRange->{{0,1},{0,1}} ,GridLines->Automatic ,AspectRatio->1 ,PlotStyle->RGBColor[1.000,0.000,0.000] ,Frame->True] -Graphics-

La figura representada por el Plot precedente es la 1, la que en el cuaderno queda insertada a continuaci´ on.

1

0.8

0.6

0.4

0.2

0.2

0.4

0.6

0.8

Figura 1. La funci´on mihumps

47

1

48

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

Integrate[f[x],x] -0.06 x + 0.05 ArcTan[5. (-0.9 + x)] + 0.1 ArcTan[10. (-0.3 + x)] N[Integrate[f[x],{x,0,1}],30] 0.2985832539549867 f[x]//TeXForm -0.06 + {{0.01}\over {0.04 + {{\left( -0.9 + x \right) }^2}}} + {{0.01}\over {0.01 + {{\left( -0.3 + x \right) }^2}}} %3//TeXForm -0.06\,x + 0.05\,\arctan (5.\,\left( -0.9 + x \right) ) + 0.1\,\arctan (10.\,\left( -0.3 + x \right) )

Las funciones f (x) —mihumps— y su integral F (x) = f (x) = −0.06 +

0.01 2

0.04 + (−0.9 + x)

+

Rx 0

f (s) ds son

0.01 2

0.01 + (−0.3 + x)

(C.1)

y F (x) = −0.06 x + 0.05 arctan(5. (−0.9 + x)) + 0.1 arctan(10. (−0.3 + x)) (C.2) C.2. Extrapolaci´ on repetida de Richardson En el siguiente cuaderno del Mathematica calculamos los errores de truncamiento de los distintos pasos de extrapolaci´on. Extrapolaci\’on repetida de Richardson Table[2^(2j)-1,{j,9}] {3, 15, 63, 255, 1023, 4095, 16383, 65535, 262143} I0[h_]=b0+b1 h^2+b2 h^4+b3 h^6+b4 h^8+b5 h^10+b6 h^12+ b7 h^14+b8 h^16 2 4 6 8 10 b0 + b1 h + b2 h + b3 h + b4 h + b5 h + 12 14 16 b6 h + b7 h + b8 h I1[h_]=I0[h/2]+(I0[h/2]-I0[h])/3//Simplify; Collect[%,h] 4 6 8 10 b2 h 5 b3 h 21 b4 h 85 b5 h b0 - ----- - ------- - -------- - --------- 4 16 64 256 12 14 16 341 b6 h 1365 b7 h 5461 b8 h ---------- - ----------- - ----------1024 4096 16384 I2[h_]=I1[h/2]+(I1[h/2]-I1[h])/15//Simplify; Collect[%,h] 6 8 10 12 b3 h 21 b4 h 357 b5 h 5797 b6 h b0 + ----- + -------- + ---------- + ----------- + 64 1024 16384 262144

´ REPETIDA DE RICHARDSON C.2. EXTRAPOLACION

14 16 93093 b7 h 1490853 b8 h ------------ + -------------4194304 67108864 I3[h_]=I2[h/2]+(I2[h/2]-I2[h])/63//Simplify; Collect[%,h] 8 10 12 b4 h 85 b5 h 5797 b6 h b0 - ----- - --------- - ----------- 4096 262144 16777216 14 16 376805 b7 h 24208613 b8 h ------------- - --------------1073741824 68719476736 I4[h_]=I3[h/2]+(I3[h/2]-I3[h])/255//Simplify; Collect[%,h] 10 12 14 b5 h 341 b6 h 93093 b7 h b0 + ------- + ---------- + ------------ + 1048576 268435456 68719476736 16 24208613 b8 h --------------17592186044416 I5[h_]=I4[h/2]+(I4[h/2]-I4[h])/1023//Simplify; Collect[%,h] 12 14 16 b6 h 1365 b7 h 1490853 b8 h b0 - ---------- - ------------- - ---------------1073741824 1099511627776 1125899906842624 % MIRICHAR aplica la extrapolaci\’on repetida de Richardson % al m\’etodo trapezoidal (Integraci\’on de Romberg) %

C.E.Neuman, 6-5-97

% inicializaci\’on i=1; a=0; b=1; h=(b-a)/10; % lazo while 1 x=a:h:b; % partici\’on l=length(x); pesos=[1 2*ones(1,l-2) 1]; % pesos de Trapezoidal II(i,1)=(h/2)*sum(pesos.*mihumps(x)); % m\’etodo Trapezoidal %%% %%% En este algoritmo se recalcula la funci\’on de nuevo, se %%% puede mejorar aprovechando las evaluaciones previas. %%% del(i)=2^(2*i)-1; % divisor para la extrapolaci\’on

49

50

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

for j=(i-1):(-1):1, % l\’{\i}nea de extrapolaciones k=i-j+1; II(j,k)=II(j+1,k-1)+(II(j+1,k-1)-II(j,k-1))/del(k-1) ; end % rof if i > 7,

% terminaci\’on ( se puede sustituir por una % que compare los valores obtenidos de I )

break; end % fi h=h/2; i=i+1; end % elihw %%%

fin de mirichar.m

La salida de mirichar (valores de II) se consigna en las tablas 1 y 2. Tabla 1. Integraci´on de Romberg para la funci´on mihumps (ver tabla 1). La integral (con el Mathematica) resulta Q = 0.2985832539549867 0.29851740902665 0.29827642178861 0.29850609244975 0.29856396932507 0.29857843315432 0.29858204877708 0.29858295266190 0.29858317863180

0.29819609270927 0.29858264933679 0.29858326161685 0.29858325443074 0.29858325398467 0.29858325395684 0.29858325395510 0

0.29860841977863 0.29858330243552 0.29858325395166 0.29858325395493 0.29858325395499 0.29858325395499 0 0

0.29858290374753 0.29858325318208 0.29858325395498 0.29858325395499 0.29858325395499 0 0 0

Tabla 2. Integraci´ on de Romberg para la funci´on mihumps (continuaci´ on, ver tablas 1 y 1). La integral (con el Mathematica) resulta Q = 0.2985832539549867 0.29858325455241 0.29858325395801 0.29858325395499 0.29858325395499 0 0 0 0

0.29858325395743 0.29858325395498 0.29858325395499 0 0 0 0 0

0.29858325395498 0.29858325395499 0 0 0 0 0 0

0.29858325395499 0 0 0 0 0 0 0

C.3. M´ etodos de Montecarlo En lo que sigue veremos dos ejemplos en el cuadrado unitario [0, 1]x[0, 1]. Ejemplo C.3.1. Calculamos mediante el m´etodo de Montecarlo el ´area bajo la funci´on f (x) = 0.5 en el cuadrado unitario. Realizamos mil veces el c´alculo del ´area estimada mediante diez mil puntos seleccionados al azar.

´ C.3. METODOS DE MONTECARLO

51

% AZAR3.M script de prueba para m\’etodo de Montecarlo % C.E.N., 05-may-97 rand(’seed’,sum(100*clock)); %inicializa los n\’umeros aleatorios scores=[]; for j=1:1000, % pru=rand(1,1) NN=10000;%0*pru; puntos=rand(NN,2); I=find(puntos(:,2) < 0.5); % determina los puntos bajo la curva [m,n]=size(I); scores=[scores; m/NN]; end mime=mean(scores) mide=std(scores)

% calcula la media aritm\’etica % calcula la desviaci\’on t\’{\i}pica

scorord=sort(scores); [m,n]=size(scores); mimeord=mean(scores(0.1*m:0.9*m)) mideord=std(scores(0.1*m:0.9*m)) figure; hist(scores);

% media podada 20% % desv\’{\i}o podado 20%

% histograma de los valores

save azar3; % Resultados de la corrida de 1000 % azar3 %mime = % 0.4997 %mide = % 0.0050 %mimeord = % 0.4999 %mideord = % 0.0050 % axis([0.475 0.525 0 300]) % grid

Vemos que el resultado —utilizando los estimadores robustos— es de 0.500 ± 0.005, intervalo de confianza que contiene el valor exacto. En la figura 2 dibujamos el histograma que representa la distribuci´on de puntajes obtenidos por la rutina de Montecarlo. Ejemplo C.3.2. Calulamos la integral de la funci´on ‘humps’ mediante el m´etodo de Montecarlo rand(’seed’,sum(100*clock)) scores=[]; for j=1:1000, % pru=rand(1,1) NN=10000;%0*pru; puntos=rand(NN,2); I=find( puntos(:,2) < 0.01*humps(puntos(:,1)) ); [m,n]=size(I);

52

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

250

200

150

100

50

0 0.48

0.485

0.49

0.495

0.5

0.505

0.51

0.515

0.52

0.525

Figura 2. Histograma de distribuci´on de los valores estimados de la integral del ejemplo C.3.1

scores=[scores; m/NN]; disp(j); end mime=mean(scores) mide=std(scores) scorord=sort(scores); [m,n]=size(scores); mimeord=mean(scores(0.1*m:0.9*m)) mideord=std(scores(0.1*m:0.9*m)) figure; hist(scores); save azar4; % Resultados de la corrida de 100 % azar4 %mime = % 2.9845e-001 %mide = % 1.3248e-002 %mimeord = % 2.9938e-001 %mideord = % 1.3119e-002

´ DE SIMPSON ADAPTATIVA C.4. INTEGRACION

% Resultados de quad % format short e % quad(’humps’,0,1,1e-5,1)/100 %ans = % 2.9858e-001 % %

C.4. Integraci´ on de Simpson adaptativa function [Q,cnt] = miquad(funfcn,a,b,tol) %MIQUAD Numerical evaluation of an integral, low order method. % Q = MIQUAD(’F’,A,B) approximates the integral of F(X) from A to B % to within a relative error of 1e-3. ’F’ is a string % containing the name of the function. Function F must return a % vector of output values if given a vector of input values. % Q = MIQUAD(F,A,B,TOL) integrates to a relative error of TOL. % % MIQUAD uses an adaptive recursive Simpson’s rule. % % See also QUAD, QUAD8. % % % %

Basada en QUAD.M C.B. Moler, 3-22-87. Copyright (c) 1984-94 by The MathWorks, Inc. C.E. Neuman, 4-29-97

% [Q,cnt] = miquad(F,a,b,tol) also returns a % function evaluation count. if nargin < 4, tol = 1.e-3; end c = (a + b)/2; % Top level initialization x = [a b c a:(b-a)/10:b]; % set up function call args = ’(x)’; y = eval([funfcn,args]); fa = y(1); fb = y(2); fc = y(3); lev = 1; % Adaptive, recursive Simpson’s quadrature if any(imag(y)) Q0 = 1e30; else Q0 = inf; end [Q,cnt] = eval([’miquadst(funfcn,a,b,tol,lev,fa,fc,fb,Q0)’]); cnt = cnt + 3;

53

54

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

%fin de la funci\’on miquad function [Q,cnt] = miquadst(FunFcn,a,b,tol,lev,fa,fc,fb,Q0) %MIQUADST Recursive function used by MIQUAD. % [Q,cnt] = miquadst(F,a,b,tol,lev,fa,fc,fb,Q0) tries to % approximate the integral of f(x) from a to b to within a % relative error of tol. F is a string containing the name % of f. The remaining arguments are generated by quad or % by the recursion. lev is the recursion level. % fa = f(a). fc = f((a+b)/2). fb = f(b). % Q0 is an approximate value of the integral. % See also QUAD and QUAD8. % % % %

Basada en QUADSTP.M C.B. Moler, 3-22-87. Copyright (c) 1984-94 by The MathWorks, Inc. C.E. Neuman, 4-29-97.

LEVMAX = 11; if lev > LEVMAX disp(’Limite de orden de recursion alcanzado en miquad.’) disp(’Probable singularidad.’) Q = Q0; cnt = 0; c = (a + b)/2; else % Evaluate function at midpoints of left % and right half intervals. h = b - a; c = (a + b)/2; x = [a+h/4, b-h/4]; f = eval([FunFcn,’(x)’]); cnt = 2; % Simpson’s rule for half intervals. Q1 = h*(fa + 4*f(1) + fc)/12; Q2 = h*(fc + 4*f(2) + fb)/12; Q = Q1 + Q2; % Recursively refine approximations. if abs(Q - Q0) > tol*abs(Q) ts2=tol/2; lm1=lev+1; [Q1,cnt1]=eval([’miquadst(FunFcn,a,c,ts2,lm1,fa,f(1),fc,Q1)’]); [Q2,cnt2]=eval([’miquadst(FunFcn,c,b,ts2,lm1,fc,f(2),fb,Q2)’]); Q = Q1 + Q2; cnt = cnt + cnt1 + cnt2; end end % fin de la funci\’on miquadst function [y]=mifun(x) % MIFUN es una funci\’on de prueba para MIQUAD % debe ser una funci\’on vectorizada % C.E. Neuman, 29-4-97

´ ´ C.5. EL METODO CLASICO DE RUNGE-KUTTA

55

y=x.^(1/2); % fin de la funci\’on mifun

C.5. El m´ etodo cl´ asico de Runge-Kutta En esta secci´ on incluimos el c´odigo Matlab que implementa el m´etodo clasico de Runge-Kutta. Se trata de los programas mirk.m, mieq.m, donde se define la funci´ on para integrar, y milo.m que es el ‘driver’ del m´etodo. function [xn,xdn]=mirk(t,x,h) % % % % MIRK Es la implementacion de un Runge-Kutta % de cuatro pasos para avanzar de t a t+h % inputs t:tiempo % x:estado % outputs xn:nuevo estado % xdn:nueva derivada del estado % La rutina mieq calcula el vector xdot % de derivadas de los estados. %%% %%% %%%

MIRK

C.E.Neuman, 10-11-1992 Revisado: 5-5-1997 Ultima revisi\’on: 7-5-1997

tv = t; xv = x; xdot=mieq(tv,xv); c1 = h*xdot; tn = tv + h/2; xn = xv + c1/2; xdot=mieq(tn,xn); c2 = h*xdot; xn = xv + c2/2; xdot=mieq(tn,xn); c3 = h*xdot; tn = tv + h; xn = xv + c3; xdot=mieq(tn,xn); c4 = h*xdot; xn = xv + (c1 + 2*c2 + 2*c3 + c4)/6; xdn=mieq(tn,xn); return %%% fin de mieq.m function xdot=mieq(t,x) % % % % MIEQ es la funcion que permite determinar la % derivada de los estados % %%% %%% %%%

C.E.Neuman, 10-11-1992 Revisado: 5-5-1997 Ultima revisi\’on: 7-5-1997

xdot=-t*x; % xdot = mihumps(x); return %%% fin de mieq.m

MIEQ

56

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

% % % % % % % % % %

MILO MILO Este es el loop central de administracion de MIRK y su driver para la integraci\’on. En MIRK.M se program\’o un paso de Runge-Kutta cl\’asico. En MIEQ.M se debe incluir la funci\’on que se desea integrar. Con modificaciones puede utilizarse en una rutina de paso variable.

%%% %%% %%%

C.E.Neuman, 10-11-1992 Revisado: 5-5-1997 Ultima revisi\’on: 7-5-1997

%%% Tiempos inicial y final y cantidad de pasos de integraci\’on a = 0; %comienzo de integraci\’on b = 7; %fin de la integraci\’on numpaso = 300; %n\’umero de pasos %%% Condici\’on inicial x0 = 1; %%% Paso de integraci\’on h = (b - a)/numpaso; %%% Definici\’on de tini y tmax tini = a; tmax = a + numpaso * h; if tmax ~= b, disp(’error en el paso’); end % fi

%%% Lista de salida (inicializaci\’on) xdn0 = mieq(tini,x0); xout = [tini; x0; xdn0];

%%% Lazo principal xn = x0; for tt=tini:h:(tmax-h), [xn,xdn] = mirk(tt,xn,h); xout = [ xout

%llamada a R-K

[tt+h; xn; xdn] ];

end % rof plot(xout(1,:),xout(2:3,:),’r’)

En la figura 3 representamos la salida de las rutinas precedentes cuando se integra la ecuaci´ on x0 = −tx (C.3)

´ ´ C.5. EL METODO CLASICO DE RUNGE-KUTTA

57

Integracion de dx/dt=−tx 1

0.8

0.6 2

x(t)=exp(−t /2) 0.4

x

0.2

0

−0.2

−0.4

−0.6

−0.8

0

1

2

3

4

5

6

7

t

Figura 3. Integraci´on por Runge-Kutta cl´asico de x0 = −tx en [0, 7]

´ APENDICE D

La integral de Riemann–Stieltjes Comencemos estudiando las propiedades de las funciones de variaci´on acotada. D.1. Funciones de variaci´ on acotada Teorema D.1.1. [Apostol (1976), p´ag. 153] Sea f una funci´on creciente definida en [a, b] y sean x0 , x1 , . . . , xn n + 1 puntos tales que a = x0 < x1 < x2 < . . . < xn = b

(D.1)

Tenemos entonces la desigualdad n−1 X

[f (xk +) − f (xk −)] ≤ f (b) − f (a)

(D.2)

k=1

´ n: Sea ξk ∈ (xk , xk+1 ). Se tiene entonces que f (xk +) − f (xk −) ≤ Demostracio f (ξk ) − f (ξk−1 ) por la monoton´ıa de la funci´on f . Resulta as´ı n−1 X

[f (xk +) − f (xk −)] ≤

k=1

n−1 X

[f (ξk ) − f (ξk−1 )]

(D.3)

k=1

= f (ξn−1 ) − f (ξ0 ) ≤ f (b) − f (a)

(D.4)

lo que completa la demostraci´on. Teorema D.1.2. [Apostol (1976), p´ag. 153] Si f es mon´otona en [a, b], entonces el conjunto de discontinuidades de f es numerable Ejercicio D.1.1. Demostrar el teorema D.1.2.

Una funci´ on mon´ otona tiene, a lo sumo, numerables saltos

´ n D.1.1. Sea f definida en [a, b]. Sea P = {x0 , x1 , . . . , xn } una Definicio partici´ on de [a, b]. Escribiremos, para k = 1, 2, . . . , n, ∆fk = f (xk ) − f (xk−1 ). Si existe un n´ umero positivo M tal que n n X X |∆fk | = |f (xk ) − f (xk−1 )| ≤ M (D.5) k=1

k=1

para toda partici´ on P de [a, b], entonces diremos que f es de variaci´ on acotada. Teorema D.1.3. [Apostol (1976), p´ag. 154] Si f es mon´otona en [a, b], entonces f es de variaci´ on acotada en [a, b]. Teorema D.1.4. [Apostol (1976), p´ag. 155] Si f es de clase C 1 en [a, b], entonces f es de variaci´ on acotada en [a, b]. Ejercicio D.1.2. Demostrar los teoremas D.1.3 y D.1.4. 59

Una funci´ on de clase C 1 es de variaci´ on acotada

60

D. LA INTEGRAL DE RIEMANN–STIELTJES

Teorema D.1.5. P [Apostol (1976), p´ag. 155] Si f es de variaci´on acotada en [a, b], es decir que |∆fk | ≤ M para toda partici´on de [a, b], entonces f est´a acotada en [a, b]. ´ n: Sea x ∈ (a, b) y la partici´on P = {a, x, b}. Se tiene Demostracio |f (x) − f (a)| + |f (b) − f (x)| ≤ M

(D.6)

|f (x) − f (a)| ≤ M

(D.7)

|f (x)| ≤ |f (a)| + M

(D.8)

con lo cual

y, en consecuencia,

Cuando x = a vale D.8 trivialmente y cuando x = b |f (b)| ≤ |f (a)| + M

(D.9)

por la hip´ otesis. De modo que |f (a)| + M es cota de |f (x)| en todo el intervalo. ´ n D.1.2. Sea f una funci´on de variaci´on acotada en [a, b] y sea P = Definicio {x0 , x1 , . . . , xn } una partici´ on de [a, b]. El n´ umero V (a, b) = sup{

n X

|∆fk | : P es partici´on de [a,b]}

(D.10)

k=1

se llama variaci´ on total de f en el intervalo [a, b]. Teorema D.1.6. Sea f una funci´on de variaci´on acotada en [a, b]. Sea V definida en [a, b] por V (a) = 0 y, para a < x ≤ b, V (x) = V (a, x), la variaci´on total de f entre a y x. Entonces (i) V es una funci´ on creciente en [a, b] (ii) W = V − f es una funci´on creciente en [a, b]

´ n: Sean x < y en [a, b]. Entonces V (a, y) = V (a, x) + V (x, y), lo que Demostracio implica que V (y) − V (x) = V (x, y) ≥ 0 lo que prueba (i). Asimismo W (y) − W (x) = V (y) − V (x) − (f (y) − f (x)) = V (x, y) − (f (y) − f (x)) ≥ 0 (D.11) lo que prueba (ii). Teorema D.1.7. Sea f definida sobre [a, b]. Entonces f es de variaci´on acotada en [a, b] sii f puede expresarse como diferencia de dos funciones crecientes. Una funci´ on es de variaci´ on acotada sii es diferencia de dos funciones crecientes

´ n: Si f es de variaci´on acotada en [a, b], podemos escribir f = V −W , Demostracio donde V es la funci´ on del teorema D.1.6 y W = V − f . Las dos son funciones crecientes. La otra implicaci´ on se obtiene de la monoton´ıa de las dos funciones.

D.2. LA INTEGRAL

61

D.2. La integral ´ n D.2.1. Sea P = {x0 , x1 , . . . , xn } una partici´on del intervalo [a, b] Definicio y sea, para cada k, ξk un punto del intervalo [xk−1 , xk ]. Una suma de la forma n X S(P, f, g) = f (ξk )∆gk (D.12) k=1

(donde ∆gk = g(xk ) − g(xk−1 )) se llama suma de Riemann-Stieltjes de f respecto de g. Diremos que f es integrable respecto de g en [a, b] si existe un n´ umero I que satisface la siguiente propiedad: para cada ε > 0, existe una partici´on Pε de [a, b] tal que, para cada partici´ on P m´as fina que Pε y para cada elecci´on de los puntos ξk del intervalo [xk−1 , xk ], k = 1, . . . , n, se tiene |S(P, f, g) − I| < ε

(D.13)

´ n D.2.1. Si un tal n´ Observacio umero I existe, entonces es u ´nico y se representa por la expresi´ on Z b f (x)dg(x) (D.14) a

Ejemplo D.2.1. Si g(x) = x entonces la integral se reduce a una integral de Riemann Ejemplo D.2.2. Si g(x) es de clase C 1 en [a, b] entonces la integral resulta Z b f (x)g 0 (x)dx (D.15) a

Nota: Este resultado se demuestra en el teorema D.2.5 en la p´agina 63 Ejemplo D.2.3. Si f es continua y g(x) = bxc entonces la integral resulta X f (i) (D.16) i∈Z a≤i
Nota: Este resultado se demuestra en el teorema D.2.7 en la p´agina 65 Teorema D.2.2. Si f1 y f2 son integrables en el sentido de Riemann-Stieltjes respecto de g en el intervalo [a, b], entonces, para todo par de constantes c1 y c2 , se tiene que c1 f1 + c2 f2 resulta es integrable respecto de g en el mismo intervalo y vale Z Z Z b

b

(c1 f1 + c2 f2 )dg = c1

b

f1 dg + c2

a

a

f2 dg

(D.17)

a

´ n: Sea f = c1 f1 + c2 f2 . Para una partici´on P del intervalo vale Demostracio n n n X X X S(P, f, g) = f (ξk )∆gk = c1 f1 (ξk )∆gk + c2 f2 (ξk )∆gk (D.18) k=1

k=1

k=1

es decir S(P, f, g) = c1 S(P, f1 , g) + c2S(P, f2 , g) Sea ε > 0 dado, elegimos una partici´on que ella se tenga

(1) Pε

Z |S(P, f1 , g) −

tal que para toda partici´on P m´as fina b

f1 dg| < ε a

(D.19)

(D.20)

62

D. LA INTEGRAL DE RIEMANN–STIELTJES (2)

y elegimos otra partici´ on Pε tenga

tal que para toda partici´on P m´as fina que ella se b

Z |S(P, f2 , g) −

f2 dg| < ε

(D.21)

a

Unimos las particiones para obtener Pε = Pε(1) ∪ Pε(2)

(D.22)

de modo que, para cada partici´on P m´as fina que Pε se obtiene Z b Z b f1 dg − c2 f2 dg| ≤ |c1 |ε + |c2 |ε |S(P, f, g) − c1 a

(D.23)

a

lo que (con una reelecci´ on del ε) prueba el teorema. Teorema D.2.3. Si f es integrable respecto de g1 en el sentido de RiemannStieltjes en el intervalo [a, b], y tambi´en es integrable respecto de g2 entonces, para todo par de constantes c1 y c2 , se tiene que f tambi´en es integrable respecto de c1 g1 + c2 g2 en el mismo intervalo y vale Z b Z b Z b f d(c1 g1 + c2 g2 ) = c1 f dg1 + c2 f dg2 (D.24) a

a

a

Ejercicio D.2.4. Demostrar el teorema D.2.3 en forma an´aloga al D.2.2. Demostrar adem´ as la aditividad respecto del intervalo de integraci´on de la integral de Riemann-Stieltjes D.2.1. Integraci´ on por partes. Teorema D.2.4. [F´ ormula de integraci´ on por partes] Si f es integrable respecto de g, en el sentido de Riemann-Stieltjes, en el intervalo [a, b] entonces g es integrable respecto de f en el sentido de Riemann-Stieltjes en el intervalo [a, b], y se tiene Z b Z b f (x)dg(x) + g(x)df (x) = f (b)g(b) − f (a)g(a) (D.25) a

a

´ n: Fijamos un n´ Demostracio umero ε > 0, por la existencia de la primera integral existe una partici´ on Pε del intervalo [a, b] tal que para toda partici´on P˜ m´as fina que Pε tenemos Z b |S(P˜ , f, g) − f (x)dg(x)| < ε (D.26) a

Para una partici´ on P m´ as fina que Pε planteamos una suma de Riemann-Stieltjes Rb cualquiera para la integral a g(x)df (x): S(P, g, f ) =

n X

g(ξk )∆fk =

k=1

n X

g(ξk )f (xk ) −

k=1

n X

g(ξk )f (xk−1 )

(D.27)

k=1

como se tiene adem´ as que f (b)g(b) − f (a)g(a) =

n X k=1

f (xk )g(ξk ) −

n X k=1

f (xk−1 )g(ξk−1 )

(D.28)

D.2. LA INTEGRAL

63

es posible restar las expresiones D.27 y D.28 para obtener f (b)g(b) − f (a)g(a) − S(P, g, f ) n n X X f (xk )[g(xk ) − g(ξk )] − f (xk−1 )[g(xk ) − g(ξk−1 )] (D.29) = k=1

k=1

Si P˜ es una partici´ on formada por los puntos xk y los ξk , entonces es m´as fina que la partici´ on P y, por ello, que la Pε entonces si llamamos S(P˜ , f, g) al segundo miembro de la ecuaci´ on D.29 se tiene que la desigualdad D.26 es v´alida y se tiene Z b g(x)df (x)| < ε (D.30) |f (b)g(b) − f (a)g(a) − S(P, g, f ) − a

para toda partici´ on P m´ as fina que Pε entonces Rb f (b)g(b) − f (a)g(a) − a f (x)dg(x)

Rb a

g(x)dg(x) existe y vale

D.2.2. Propiedades de la integral. Teorema D.2.5. Sea f una funci´on integrable en el sentido de RiemannStieltjes respecto de una funci´on g de clase C 1 en el intervalo [a, b], entonces existe la integral de Riemann del segundo miembro y vale Z b Z b f (x)dg(x) = f (x)g 0 (x)dx (D.31) a

a

´ n: Sea la suma de Riemann Demostracio n X S(P, f g 0 ) = f (ξk )g 0 (ξk )∆xk

(D.32)

k=1

y (para la misma partici´ on) la suma de Riemann-Stieltjes S(P, f, g) =

n X

f (ξk )∆gk

k=1

n X

f (ξk )g 0 (ζk )∆xk

(D.33)

k=1

donde aplicamos el teorema del valor intermedio y ζk ∈ (xk−1 , xk ), restando se obtiene n X S(P, f, g) − S(P, f g 0 ) = f (ξk )[g 0 (ζk ) − g 0 (ξk )]∆xk (D.34) k=1

Utilizando la hip´ otesis para el integrador g se puede probar que ε |S(P, f, g) − S(P, f g 0 )| < (D.35) 2 para cada partici´ on P m´ as fina que una dada Pε . Por la hip´otesis de integrabilidad de f respecto de g se tiene que, para toda partici´on P m´as fina que una dada P˜ε , vale Z b ε |S(P, f, g) − f (x)dg(x)| < (D.36) 2 a Entonces cuando P es m´ as fina que la uni´on Pε ∪ P˜ε se tiene que Z b |S(P, f g 0 ) − f (x)dg(x)| < ε (D.37) a

lo que demuestra el teorema.

64

D. LA INTEGRAL DE RIEMANN–STIELTJES

D.2.2.1. Cambio de variables en la integral de Riemann-Stieltjes. Si la funci´on h es un cambio de variables biyectivo de clase C ∞ entre el intervalo [a, b] y el intervalo [c, d] entonces vale la siguiente f´ormula de cambio de variables (que damos sin demostraci´ on) Z h(d) Z b f (t)dg(t) = f (h(x))dg(h(x)) (D.38) h(c)

a

D.2.2.2. Integral de funciones escalonadas. Teorema D.2.6. Definimos g en el intervalo [−1, 1] como sigue   c1 si −1 ≤ x < 0 c2 si x=0 g(x) =  c3 si 0<x≤1

(D.39)

(donde c1 , c2 , y c3 son arbitrarios. Sea f una funci´on definida en [−1, 1] de modo que por lo menos una de las funciones f y g sea continua a la izquierda del origen y por lo menos una de las dos lo sea a la derecha de 0. Entonces f es integrable en el sentido de Riemann-Stieltjes en el intervalo [−1, 1] y la integral vale Z 1 f (x)dg(x) = f (0)(g(0+) − g(0−)) (D.40) −1

´ n: Sea P una partici´on que contiene al origen, entonces si xk = 0 Demostracio S(P, f, g) = f (ξk−1 )(g(0) − g(0−)) + f (ξk )(g(0+) − g(0))

(D.41)

donde ξk−1 ≤ 0 ≤ ξk . Sumando y restando f (0)g(0) se puede escribir S(P, f, g) − f (0)(g(0+) − g(0−)) = (f (ξk−1 ) − f (0))(g(0) − g(0−)) +(f (ξk ) − f (0))(g(0+) − g(0))

(D.42)

|S(P, f, g) − f (0)(g(0+) − g(0−))| ≤ |f (ξk−1 ) − f (0)||g(0) − g(0−)| +|f (ξk ) − f (0)||g(0+) − g(0)|

(D.43)

por lo tanto se tiene

Si f es continua en el origen, para cada ε > 0 existe un δ > 0 tal que, si kP k < δ, entonces |f (ξk−1 ) − f (0)| < ε (D.44) y |f (ξk ) − f (0)| < ε (D.45) de modo que se tiene la desigualdad |S(P, f, g) − f (0)(g(0+) − g(0−))| ≤ ε|g(0) − g(0−)| + ε|g(0+) − g(0)|

(D.46)

En los otros casos la desigualdad sigue siendo v´alida (1) f discontinua a izquierda y derecha del origen: entonces g debe ser continua y vale S(P, f, g) − f (0)(g(0+) − g(0−)) = 0 (2) f es discontinua a la derecha y continua a la izquierda: en este caso debe valer g(0+) − g(0) = 0 y vale |S(P, f, g) − f (0)(g(0+) − g(0−))| ≤ ε|g(0) − g(0−)|

(D.47)

D.4. PERO, ¿HAY FUNCIONES INTEGRABLES?

65

(3) f es continua a la derecha y discontinua a la izquierda: en este caso debe valer g(0) − g(0−) = 0 y vale |S(P, f, g) − f (0)(g(0+) − g(0−))| ≤ ε|g(0+) − g(0)|

(D.48)

En todos los casos vale la acotaci´on por lo que el teorema queda demostrado. Este u ´ltimo resultado permite demostrar que las integrales respecto de funciones escalonadas se reducen a sumas. Teorema D.2.7. Sea g una funci´on escalonada definida en el intervalo [a, b] con salto gk en el punto xk de una partici´on {x1 , x2 , . . . , xn } del intervalo. Sea f definida en [a, b] con la propiedad que f y g no sean las dos discontinuas a la izquierda o a la derecha en cada punto xk . Entonces existe la integral de f respecto de g y vale Z b n X f (x)dg(x) = f (xk )gk (D.49) a

k=1

´ n: Basta descomponer el intervalo [a, b] en subintervalos y aplicar el Demostracio teorema D.2.6. D.3. F´ ormula de la suma de Euler Teorema D.3.1. Supongamos que f es de clase C ∞ en [a, b]. Si a y b son enteros se tiene  Z b Z b b X 1 f (a) + f (b) f (n) = f (x)dx + x − bxc − f 0 (x)dx + (D.50) 2 2 a a n=a ´ n: Por el teorema D.2.4 se tiene que Demostracio   Z b  Z b 1 1 f (x)d x − bxc − + x − bxc − df (x) 2 2 a   a   1 1 = f (b) b − bbc − − f (a) a − bac − 2 2 es decir Z b a



1 f (x) d x − bxc − 2



Z + a

b



1 x − bxc − 2

 df (x) =

f (a) − f (b) 2

(D.51) (D.52)

(D.53)

por u ´ltimo Z a

b

  Z b b X 1 f (x)d x − bxc − = f (x) dx − f (n) 2 a n=a+1

(D.54)

lo que, al combinarla con la anterior, completa la demostraci´on. D.4. Pero, ¿Hay funciones integrables? Se puede demostrar que si f es continua y g es de variaci´on acotada en el intervalo [a, b] entonces f es integrable en el sentido de Riemann-Stieltjes respecto del integrador g. La demostraci´on es muy simple y se basa en la definici´on de las denominadas sumas superior e inferior de Riemann-Stieltjes y en ver que, cuando se refina la partici´ on, ambas sumas se aproximan tanto como se desee.

66

D. LA INTEGRAL DE RIEMANN–STIELTJES

D.5. Ejercicios Ejercicio D.5.1. Una funci´on f , definida en [a, b], verifica una condici´on uniforme de Lipschitz de orden α > 0 en [a, b] si existe una constante M > 0 tal que |f (x) − f (y)| < M |x − y|α para todo x e y de [a, b] (a) Si f es una tal funci´on, probar que α > 1 implica que f es constante en [a, b], mientras que α = 1 implica que f es de variaci´on acotada en [a, b] (b) Dar un ejemplo de una funci´on f que satisfaga una condici´on uniforme de Lipschitz de orden α < 1 en [a, b] tal que f no sea de variaci´on acotada en [a, b] (c) Dar un ejemplo de una funci´on f que sea de variaci´on acotada y que, sin embargo, no satisfaga ninguna condici´on uniforme de Lipschitz en [a, b]. Ejercicio D.5.2. Probar que una funci´on polin´omica f es de variaci´on acotada en todo intervalo cerrado y acotado [a, b]. Describir un m´etodo que permita calcular la variaci´ on total de f en [a, b] conociendo los ceros de la derivada f 0 . Ejercicio D.5.3. La siguiente definici´on de la integral de Riemann-Stieltjes es bastante usual en textos matem´aticos: Se dice que f es integrable respecto de g si existe un n´ umero real A con la siguiente propiedad: para cada ε > 0 existe un δ > 0 tal que para cada partici´on P de [a, b] con norma kP k < δ y cada elecci´on de ξk en [xk−1 , xk ], tenemos |S(P, f, g) − A| < ε Rb (a) Demostrar que si a f dg existe seg´ un esta definici´on entonces tambi´en existe de acuerdo con la definici´on D.2.1 (p´agina 61) y ambas son iguales. (b) Sean f (x) = g(x) = 0 para −1 ≤ x < 0, f (x) = g(x) = 1 para 0 < x ≤ 1, Rb f (0) = 0, y g(0) = 1. Demostrar que a f dg existe de acuerdo a la definici´ on D.2.1 pero no existe seg´ un la definici´on de este ejercicio.

´ APENDICE E

Polinomios y n´ umeros de Bernoulli Los polinomios de Bernoulli, Bn (x), resultan de la serie generadora ∞ X text tk = Bk (x) t e −1 k!

(E.1)

k=0

y los n´ umeros de Bernoulli, Bn , de ∞

X t tk = B k et − 1 k!

(E.2)

k=0

(es decir, de hacer x = 0 en la anterior, raz´on por la cual Bn (0) = Bn ). A continuaci´ on desarrollamos un cuaderno del Mathematica con algunos ejemplos de estos n´ umeros y polinomios. E.1. Polinomios de Bernoulli con el Mathematica (* N\’umeros y polinomios de Bernoulli *) (* Calculamos la serie generadora de los polinomios la correspondiente serie para los n\’umeros se obtiene, simplemente, fijando x = 0 en la primera *) Series[t Exp[t x]/(Exp[t]-1),{t,0,4}] 2 2 3 1 1 x x 2 x x x 3 1 + (-(-) + x) t + (-- - - + --) t + (-- - -- + --) t + 2 12 2 2 12 4 6 2 3 4 1 x x x 4 5 (-(---) + -- - -- + --) t + O[t] 720 24 12 24 Series[t Exp[t x]/(Exp[t]-1),{t,0,12}]; CoefficientList[%,t] 2 2 3 2 3 4 1 1 x x x x x 1 x x x {1, -(-) + x, -- - - + --, -- - -- + --, -(---) + -- - -- + --, 2 12 2 2 12 4 6 720 24 12 24 3 4 5 2 4 5 6 -x x x x 1 x x x x --- + -- - -- + ---, ----- - ---- + --- - --- + ---, 720 72 48 120 30240 1440 288 240 720 3 x

x

5 x

6 x

7 x 67

68

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

----- - ---- + ---- - ---- + ----, 30240 4320 1440 1440 5040 2 4 6 7 8 1 x x x x x -(-------) + ----- - ----- + ---- - ----- + -----, 1209600 60480 17280 8640 10080 40320 3 5 7 8 9 -x x x x x x ------- + ------ - ----- + ----- - ----- + ------, 1209600 181440 86400 60480 80640 362880 2 4 6 8 9 10 1 x x x x x x -------- - ------- + ------ - ------ + ------ - ------ + -------, 47900160 2419200 725760 518400 483840 725760 3628800 3 5 7 9 10 11 x x x x x x x -------- - ------- + ------- - ------- + ------- - ------- + --------, 47900160 7257600 3628800 3628800 4354560 7257600 39916800 2 4 6 8 10 691 x x x x x -(-------------) + -------- - -------- + -------- - -------- + -------- 1307674368000 95800320 29030400 21772800 29030400 43545600 11 12 x x -------- + ---------} 79833600 479001600 c=2;cm1=c-1; BernoulliB[cm1] BernoulliB[cm1,x] cm1! %3[[c]]//Simplify D[%,x]/cm1 Integrate[%%,{x,0,1}] Plot[%%%,{x,0,1}] 1 -(-) 2 1 -(-) + x 2 1 -(-) + x 2 1 0 -Graphics- (*berno1.wmf*)

En la figura 1 y siguientes representamos algunos de los polinomios de Bernoulli. c=3;cm1=c-1; BernoulliB[cm1] BernoulliB[cm1,x]==cm1! %3[[c]]//Simplify cm1! %3[[c]]//Simplify

E.1. POLINOMIOS DE BERNOULLI CON EL Mathematica

0.4 0.2

0.2

0.4

0.6

0.8

-0.2 -0.4

Figura 1. El polinomio de Bernoulli B1 (x)

D[%,x]/cm1 //Simplify Integrate[%%,{x,0,1}] Plot[%%%,{x,0,1}] 1 6 True 1 2 - - x + x 6 1 -(-) + x 2 0 -Graphics- (*berno2.wmf*) c=7;cm1=c-1; BernoulliB[cm1] BernoulliB[cm1,x]==cm1! %3[[c]]//Simplify cm1! %3[[c]]//Simplify D[%,x]/cm1 //Simplify Integrate[%%,{x,0,1}] Plot[%%%,{x,0,1}] 1 -42 True

1

69

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

70

0.15 0.1 0.05

0.2

0.4

0.6

0.8

-0.05

Figura 2. El polinomio de Bernoulli B2 (x)

2 4 1 x 5 x 5 6 -- - -- + ---- - 3 x + x 42 2 2 3 4 -x 5 x 5 x 5 -- + ---- - ---- + x 6 3 2 0 -Graphics- (*berno3.wmf*) c=11;cm1=c-1; BernoulliB[cm1] BernoulliB[cm1,x]==cm1! %3[[c]]//Simplify cm1! %3[[c]]//Simplify D[%,x]/cm1 //Simplify Integrate[%%,{x,0,1}] Plot[%%%,{x,0,1}] 5 -66 True 2 8 5 3 x 4 6 15 x 9 10 -- - ---- + 5 x - 7 x + ----- - 5 x + x 66 2 2 2 4 6 7 8

1

E.1. POLINOMIOS DE BERNOULLI CON EL Mathematica

0.02 0.01

0.2

0.4

0.6

0.8

1

-0.01 -0.02

Figura 3. El polinomio de Bernoulli B6 (x)

x (-3 + 20 x - 42 x + 60 x - 45 x + 10 x ) ---------------------------------------------10 0 -Graphics- (*berno4.wmf*) (* Resulta que las integrales de los coeficientes son nulas, (todas, salvo la primera) porque vale la siguiente integral: *) Integrate[t Exp[t x]/(Exp[t]-1),x] t x E ------t -1 + E Integrate[t Exp[t x]/(Exp[t]-1),{x,0,1}]//Simplify 1 (* Calculamos ahora el coeficiente de t^40 en el desarrollo en serie *) n=40; n! Coefficient[Series[t Exp[t x]/(Exp[t]-1),{t,0,n}],t^n]//Simplify %==BernoulliB[40,x] Integrate[%%,{x,0,1}] Plot[%%%,{x,0,1}] 4 261082718496449122051 2 26315271553053477373 x -(---------------------) + 380899208799402670 x - ----------------------- + 13530 21

71

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

72

0.06 0.04 0.02 0.2

0.4

0.6

0.8

1

-0.02 -0.04 -0.06

Figura 4. El polinomio de Bernoulli B10 (x)

6 1649024253633120910 x

8 10 2325031004857511379 x 10708663585311718520 x - ---------------------- + ------------------------ 2 21

12 16 457533651717217348 x 14 38092256087030385 x ---------------------- + 33081876872548920 x - --------------------- + 3 7

18 702064548199300 x

20 - 72937940132694 x

24 445756964055 x

26 + 27074751480 x

22 43628525827900 x + ------------------ 7

28 30 29696275036 x 192650120 x - --------------- + ------------- 21 3

32 36 5126979 x 34 9139 x 38 39 40 ----------- + 91390 x - -------- + 130 x - 20 x + x 2 3 True 0

´ E.2. NUMEROS DE BERNOULLI

73

-Graphics- (*berno5.wmf*)

16

2·10

16

1·10

0.2

0.4

0.6

0.8

1

16

-1·10

16

-2·10

Figura 5. El polinomio de Bernoulli B40 (x)

(* Observamos que los valores son muy grandes por eso calculamos el valor del polinomio y luego los partimos por el correspondiente factorial *) BernoulliB[40,0.5] %/40! 16 1.92966 10 -32 2.36502 10

E.2. N´ umeros de Bernoulli Para demostrar que los n´ umeros de Bernoulli con ´ındice impar ≥ 3 son nulos observemos que t t t t et + 1 − B1 t = t + = (E.3) t e −1 e −1 2 2 et − 1 pero esta u ´ltima es una funci´ on par puesto que t et + 1 t e−t + 1 = − −t t 2e −1 2e −1

(E.4)

con lo cual todos los coeficientes de ´ındice impar de su desarrollo deben anularse.

74

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

Ejercicio E.2.1. Multiplicar ambos miembros de E.2 por et − 1, e igualar coeficientes para obtener n   X n Bk = Bn + δ1n (E.5) k k=0

Ejercicio E.2.2. Demostrar que n   X m Bm (x) = Bk xm−k k

(E.6)

k=0

En algunos textos la expresi´on E.6 se toma como definici´on de los polinomios de Bernoulli. Cuando m = 1, se tiene B1 (x) = B0 x + B1 = x − 12 . Si m > 1 se tiene Bm (1) = Bm = Bm (0) (ver el ejercicio E.2.1). Resulta, en consecuencia que Bm (x − bxc) no tiene saltos en las coordenadas enteras. E.3. Integraci´ on por partes Teorema E.3.1. Vale la identidad 0 Bm (x) = mBm−1 (x)

´ n: Por la expresi´on E.6 se tiene Demostracio n   X m 0 Bm (x) = (m − k)Bk xm−k−1 k k=0  n  X m−1 =m Bk xm−1−k = mBm−1 (x) k

(E.7)

(E.8) (E.9)

k=0

Por u ´ltimo se tiene la siguiente f´ormula de integraci´on por partes, Z 1 1 Bm (x)f (m) (x)dx = m! 0 1 (Bm+1 (1)f (m) (1) − Bm+1 (0)f (m) (0) (m + 1)! Z 1 1 Bm+1 (x)f (m+1) (x)dx (E.10) (m + 1)! 0 Este u ´ltimo resultado, convenientemente extendido, lo utilizamos para demostrar la f´ ormula de Euler para el error del m´etodo trapezoidal.

´ APENDICE F

Laboratorio de Matem´ atica: Integraci´ on Adaptativa F.1. Objetivos ∗

Las metas del laboratorio de Cuadratura Adaptativa son • Introducir el concepto de m´etodo adaptativo para la resoluci´on num´erica de distintos problemas de Matem´atica, y, en particular, para el caso de la integraci´ on de funciones • Comparar la eficiencia de distintos m´etodos de integraci´on num´erica • Estudiar los errores de aproximaci´on asociados a distintos m´etodos de integraci´ on num´erica F.2. Antes del laboratorio La idea b´ asica de un m´etodo adaptativo es refinar los c´alculos solamente en los dominios en los que es necesario y realizar el m´ınimo de cuentas en los restantes para lograr a la postre una distribuci´on uniforme de los errores de aproximaci´on function y=mifun(x) % MIFUN funci\’on para probar la adaptatividad de QUAD. % Tiene que ser una funci\’on vectorizada. %

C.E. Neuman, 30 de mayo de 1997

I=find(x < 0.5); J=find(x >= 0.5); y=zeros(size(x)); % mayores o iguales que 1/2 y(I)=-2.5e8*((x(I)-0.5).*x(I)).^5; % menores que 1/2 y(J)=720*(x(J)-0.5); % fin de la funci\’on mifun

f (x) =


5  − 52 108 x(x − 12 ) 

720(x − 12 )

si

x<

1 2

si

x≥

1 2

(F.1)

∗ Trabajo realizado con fondos provenientes de la Universidad Nacional del Litoral (UNL), a trav´ es de la programaci´ on Curso de Acci´ on para la Investigaci´ on y el Desarrollo (CAI+D), Secretar´ıa de Ciencia y T´ ecnica de la UNL, subsidio 1042-1042-34-182 y 94-16-4-20. Proyectos Uso del producto Mathematica en la ense nanza de la Matem´ atica y Problemas te´ oricos y num´ ericos asociados a la obtenci´ on de soluciones de ecuaciones el´ıpticas y parab´ olicas

75

Hay posibles singularidades en el entorno de 0 y de 0.5

76

´ ´ ADAPTATIVA F. LABORATORIO DE MATEMATICA: INTEGRACION

400

350

300

250

200

150

100

50

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Figura 1. Puntos de evaluaci´on de quad para la funci´on f (x) en [0, 2]

N´ otese en la figura 1 que en el entorno de 0 y de 0.5 se acumulan los puntos y se alcanza en nivel de posible singularidad en el algoritmo quad. Asimismo, en el intervalo donde la funci´ on es un polinomio de orden 4 el m´etodo de Simpson es exacto y la distribuci´ on de puntos es uniforme y mucho menos densa. Ejercicio F.2.1. Dada una funci´on f definida en un intervalo [a, b] dar un m´etodo adaptativo de interpolaci´on de f mediante polinomios a trozos de orden 2 (splines lineales). Ejercicio F.2.2. Dada una funci´on f definida en un intervalo [a, b] dar un m´etodo adaptativo de interpolaci´on de f mediante polinomios a trozos de orden 4 (splines c´ ubicos) F.3. Laboratorio Ejercicio F.3.1. Desarrollar un algoritmo —y correspondiente programa en Matlab— que implemente un m´etodo de integraci´on num´erica adaptativa de una funci´ on dada en un intervalo cerrado y acotado del eje real. Ejercicio F.3.2. Realizar experimentos num´ericos destinados a determinar las caracter´ısticas del comportamiento del programa desarrollado en F.3.1 para distintos integrandos y distintos dominios de integraci´on. Comparar con otros m´etodos de integraci´ on num´erica.

F.5. REFERENCIAS

77

F.4. Aplicaciones y desarrollos avanzados Ejercicio F.4.1. Analizar en detalle las funciones quad y quad8 de Matlab y, a partir de ellas, y del programa desarrollado en F.3.1, generar un integrador autom´ atico m´ as robusto. F.5. Referencias (1) Atkinson, K.E.: An Introduction to Numerical Analysis, 2nd edition, J. Wiley & Sons., New York, 1989. (2) Burden, R.L. y Faires, J.D.: An´ alisis Num´erico, 2da edici´on, Grupo Editorial Iberoam´erica, M´exico, 1993.

´ APENDICE G

Laboratorio de Matem´ atica: Soluci´ on adaptativa de ODEs G.1. Objetivos ∗

Las metas del laboratorio de Soluci´on adaptativa de ODEs son • Desarrollar el concepto de m´etodo adaptativo introducido en otro laboratorio de este texto para la resoluci´on num´erica de distintos problemas de Matem´ atica, y, en particular, para el caso de la soluci´on de ecuaciones diferenciales • comparar la eficiencia de distintos m´etodos de soluci´on num´erica de ODEs • estudiar los errores de aproximaci´on asociados a distintos m´etodos de soluci´ on num´erica de ODEs. G.2. Antes del laboratorio La idea b´ asica de un m´etodo adaptativo es refinar los c´alculos solamente en los dominios en los que es necesario y realizar el m´ınimo de cuentas en los restantes para lograr a la postre una distribuci´on uniforme de los errores de aproximaci´on. Analizar en detalle los resultados del laboratorio de integraci´on adaptativa. G.3. Laboratorio Ejercicio G.3.1. Desarrollar un algoritmo —y correspondiente programa en Matlab— que implemente un m´etodo de integraci´on num´erica adaptativa de una ecuaci´ on diferencial dada en un intervalo del eje real. Ejercicio G.3.2. Realizar experimentos num´ericos destinados a determinar las caracter´ısticas del comportamiento del programa desarrollado en G.3.1 para distintas funciones de carga f (x, y) y distintos dominios de soluci´on. Comparar con otros m´etodos de soluci´ on num´erica de ODEs. G.4. Aplicaciones y desarrollos avanzados Ejercicio G.4.1. Analizar en detalle las funciones ode23 y ode45 de Matlab y, a partir de ellas, y del programa desarrollado en G.3.1, generar un integrador autom´ atico m´ as robusto. ∗ Trabajo realizado con fondos provenientes de la Universidad Nacional del Litoral (UNL), a trav´ es de la programaci´ on Curso de Acci´ on para la Investigaci´ on y el Desarrollo (CAI+D), Secretar´ıa de Ciencia y T´ ecnica de la UNL, subsidio 1042-1042-34-182 y 94-16-4-20. Proyectos Uso del producto Mathematica en la ense nanza de la Matem´ atica y Problemas te´ oricos y num´ ericos asociados a la obtenci´ on de soluciones de ecuaciones el´ıpticas y parab´ olicas y sus continuaciones

79

80

´ ´ ADAPTATIVA DE ODES G. LABORATORIO DE MATEMATICA: SOLUCION

G.5. Referencias (1) Atkinson, K.E.: An Introduction to Numerical Analysis, 2nd edition, J. Wiley & Sons., New York, 1989. (2) Burden, R.L. y Faires, J.D.: An´ alisis Num´erico, 2da edici´on, Grupo Editorial Iberoam´erica, M´exico, 1993.

´ APENDICE H

Existencia y unicidad de soluciones de ODEs Dado el PVI

( y 0 = f (x, y) x ∈ I, y(x0 ) = y0 x0 ∈ I.

(H.1)

y las iteraciones de Picard Z

x

yn+1 = y0 +

f (s, yn (s))ds

(H.2)

x0

definimos el rect´ angulo R : |x − x0 | ≤ a, |y − y0 | ≤ b y M = max |f (x, y)|

(H.3)

(x,y)∈R

∂f (x, y) L = max ∂y (x,y)∈R ∂f (x, y) ∂f (x, y) D = max + f (x, y) ∂x ∂y (x,y)∈R

(H.4) (H.5)

Lema H.0.1. Si α = min(a, b/M ) entonces |yn (x) − y0 | ≤ M (x − x0 ) para x0 ≤ x ≤ x0 + α ´ n: Por inducci´on en n. Si n = 0 se toma y0 (x) = y0 . Veamos ahora Demostracio que si vale para n = k entonces se puede demostrar el caso n = Rk + 1: |yk (x) − y0 | ≤ Rx x M (x − x0 ) implica que |yk+1 (x) − y0 | = | x0 f (s, yk (s))ds| ≤ x0 |f (s, yk (s))|ds ≤ M (x − x0 ) si x0 ≤ x ≤ x0 + α lo que demuestra el lema por inducci´on. Se tiene que (yn (x)) converge sii y0 (x) + (y1 (x) − y0 (x)) + (y2 (x) − y1 (x)) + . . . + (yn (x) − yn−1 (x)) + . . . converge como serie. En consecuencia bastar´a ver que la serie es absolutamente convergente: ∞ X |yn (x) − yn−1 (x)| < ∞ (H.6) n=1

´ n H.0.2. La serie es absolutamente convergente Observacio ´ n: Demostracio Z x |yn (x) − yn−1 (x)| = | (f (s, yn−1 (s)) − f (s, yn−2 (s)))ds| x0 Z x ≤ |f (s, yn−1 (s)) − f (s, yn−2 (s))|ds x0 Z x ∂f (s, ξ(s)) |yn−1 (s) − yn−2 (s)|ds = ∂y x0 81

(H.7) (H.8) (H.9)

82

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

con ξ(s) ∈ int(yn−1 (s), yn−2 (s)) de modo que (s, ξ(s)) ∈ R, ∀s ∈ [x0 , x0 + α] Z x ≤L |yn−1 (s) − yn−2 (s)|ds, x0 ≤ x ≤ x0 + α (H.10) x0

En particular Z

x

|y2 (x) − y1 (x)| ≤ L

|y1 (s) − y0 |ds x0 Z x

≤ L

M (s − y0 )ds = LM x0

Z |y3 (x) − y2 (x)|

(H.11) (x − x0 )2 2!

(H.12)

x

≤ L

|y2 (s) − y1 (s)|ds

(H.13)

x0

≤ M L2

Z

x

(s − x0 )2 (x − x0 )3 ds = M L2 2 3!

x0

(H.14)

y por inducci´ on |yn (x) − yn−1 (x)| ≤

M Ln−1 (x − x0 )n , n!

x0 ≤ x ≤ x0 + α

(H.15)

Resulta as´ı que, para x0 ≤ x ≤ x0 + α ∞ X

|yj (x) − yj−1 (x)|

≤ M (x − x0 ) + M L

j=1

(x − x0 )3 (x − x0 )2 + M L2 +(H.16) ... 2! 3!

α2 α3 + M L2 + ... (H.17) 2! 3! M (αL)2 (αL)3 M αL (αL + + + . . .) = (e − 1)(H.18) L 2 3 L

≤ Mα + ML =

Luego yn (x) converge para cada x en x0 ≤ x ≤ x0 + α Un argumento similar muestra que yn (x) converge para cada x en x0 − α ≤ x ≤ x0 El l´ımite de (yn (x)) se denotar´a y(x) ´ n H.0.3. La funci´on y(x) satisface la ecuaci´on integral Observacio Z x y(x) = y0 + f (s, y(s))ds

(H.19)

x0

y se tiene que y(x) es continua Como x

Z yn+1 (x) = y0 +

f (s, yn (s))ds

(H.20)

x0

tomamos l´ımite: Z

x

y(x) = y0 + lim

n→∞

f (s, yn (s))ds

(H.21)

x0

y tenemos que calcular este l´ımite. Notemos que y(x) pertenece a R para x0 ≤ x ≤ x0 + α porque es el l´ımite de funciones con gr´aficos en R.

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

H.0.1. C´ alculo del l´ımite. Z x Z x f (s, y(s))ds − f (s, yn (s))ds| ≤ | x0 x0 Z x ≤ |f (s, y(s)) − f (s, yn (s))|ds x0 Z x ≤ L |y(s) − yn (s)|ds

83

(H.22) (H.23) (H.24)

x0

pero y(s) − yn (s) =

P∞

− yj−1 (s)) y entonces

j=n+1 (yj (s)

∞ X

|y(s) − yn (s)| ≤ M ≤ M M L

=

Lj−1

(s − x0 )j j!

(H.25)

Lj−1

αj j!

(H.26)

(αL)j M αL (e − 1) j! L j=n+1

(H.27)

j=n+1 ∞ X j=n+1 ∞ X

de modo que la integral es (H.24) ≤ M

Z ∞ X (Lα)j x ds j! x0 j=n+1

≤ Mα

∞ X (Lα)j → 0 (n → ∞) j! j=n+1

(H.28)

(H.29)

Resulta as´ı que Z

t

lim

n→∞

Z

t

f (s, yn (s))ds = t0

f (s, y(s))ds

(H.30)

t0

H.0.2. Continuidad de y(x). Escribimos y(x+h)−y(x) = (y(x+h)−yN (x+h))+(yN (x+h)−yN (x))+(yN (x)−y(x)) (H.31) Sea N tal que M L

∞ X (αL)j ε < j! 3

(H.32)

j=N +1

de modo que si h es suficientemente peque˜ no, h < α, se tiene ε |y(x + h) − yN (x + h)| < 3 y ε |y(x) − yN (x)| < 3 y, por la continuidad de yN (integral repetida) ε ∃δ > 0 : |yN (x + h) − yN (x)| < si |h| < δ 3 entonces

(H.33) (H.34)

(H.35)

|y(x+h)−y(x)| ≤ |y(x+h)−yN (x+h)|+|yN (x+h)−yN (x)|+|yN (x)−y(x)| ≤ ε si |h| < δ (H.36)

84

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

Hemos demostrado as´ı que el problema de valores iniciales tiene una soluci´on continua en x0 ≤ x ≤ x0 + α Teorema H.0.4. La soluci´on es u ´nica ´ n: Supongamos que adem´as de y(x) se tiene que z(x) tambi´en es Demostracio soluci´ on: Z x y(x) = y0 + f (s, y(s))ds (H.37) x0 x

Z z(x) = y0 +

f (s, z(s))ds

(H.38)

x0

entonces Z x |y(x) − z(x)| = | (f (s, y(s)) − f (s, z(s)))ds| x0 Z x ≤ |f (s, y(s)) − f (s, z(s))|ds x0 Z x ≤ L |y(s) − z(s)|ds

(H.39) (H.40) (H.41)

x0

lo que prueba el teorema una vez demostrado el siguiente lema Lema H.0.5. Sea w(x) una funci´on no negativa con Z x w(x) ≤ w(s)ds

(H.42)

x0

Entonces w(x) es id´enticamente cero. ´ n: Sea U (x) = Demostracio

Rx x0

w(s)ds entonces

dU = w(x) ≤ dx

Z

x

w(s)ds = LU (x)

(H.43)

x0

lo que implica e−L(x−x0 ) U (x) ≤ U (x0 )

(x ≥ x0 )

(H.44)

o bien U (x) = 0

(H.45)

de modo que Z

x

0 ≤ w(x) ≤ L

w(s)ds = LU (x) = 0

(H.46)

x0

o tambien w(x) = 0

(H.47)

H.1. REFERENCIAS

85

H.0.3. Errores en el m´ etodo de Euler. El m´etodo de Euler tiene la expresi´ on yk+1 = yk + hf (xk , yk ). Si desarrollamos en serie de Taylor: Y (xk+1 = 2 ∂f Y (xk ) + hf (xk , Y (xk )) + h2 ( ∂f ∂x + f ∂y )(ξk , Y (ξk )) con ξk ∈ int(xk , xk+1 ) Si restamos: Yk+1 − yk+1 = Yk − yk + h(f (xk , Yk ) − f (xk , yk )) + f fy )(ξk , Y (ξk )) y se tiene |Yk+1 − yk+1 |

≤ |Yk − yk | + h|fy (xk , ηk )||Yk − yk | + |Yk − yk | + hL|Yk − yk | +

Ek+1 = (1 + hL)Ek +

+

h2 |(fx + f fy )(ξk , Y (H.48) (ξk ))| 2

Dh2 2 y si llamamos Ek = |Yk − yk | al error (k > 0) se tiene ≤

h2 2 (fx

Dh2 2

(H.49)

(H.50)

Se prueba que Dh2 ((1 + hL)k − 1) 2 hL se tiene tambien que Dh khL Ek ≤ ekhL E0 + (e − 1) 2L

Ek ≤ (1 + hL)k E0 + y como (1 + hL) ≤ ehL

como kh < α Ek ≤ eαL E0 +

D αL (e − 1)h 2L

(H.51)

(H.52)

(H.53)

H.1. Referencias (1) Atkinson, K.E.: An Introduction to Numerical Analysis, 2nd edition, J. Wiley & Sons., New York, 1989. (2) Burden, R.L. y Faires, J.D.: An´ alisis Num´erico, 2da edici´on, Grupo Editorial Iberoam´erica, M´exico, 1993.

H.1. REFERENCIAS

Terminado de imprimir por CEN el 30 de noviembre de dos mil uno en la ciudad de Santa Fe, Rep´ ublica Argentina.

87