Analisis Numerico de Integrales y Ecuaciones Diferenciales

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

Temas de Análisis Numérico y Programación

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

ii

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

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

¶

Segunda 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.intec.unl.edu.ar/CN/CNBook3.pdf

Key words and phrases. Integración numérica, Resolución numérica de ecuaciones diferenciales, Métodos de un paso, Integración adaptativa

Trabajo realizado con fondos provenientes de la Universidad Nacional del Litoral (UNL), a través de la programación Curso de Acci´ on para la Investigaci´ on y el Desarrollo (CAI+D), Secretar´ıa de Ciencia y Técnica 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.

Índice General Prefacio

xiii

Prefacio a la segunda edici´ on

xvii

Índice de Tablas

1

Índice 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étodos b´ asicos de integración numérica 1.5. Métodos de extrapolación 1.6. Cuadratura de Gauß 1.7. Métodos de Montecarlo 1.8. Métodos 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étodos en diferencias de un paso 2.3. Ejercicios

21 21 23 33

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

35 35 35

Apéndice B. Coeficientes indeterminados con el Mathematica B.1. El método de Simpson B.2. El método de Newton B.3. El método 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éndice C. Cuadernos del Mathematica y archivos script de Matlab C.1. La función mihumps C.2. Extrapolaci´ on repetida de Richardson C.3. Métodos de Montecarlo C.4. Integraci´ on de Simpson adaptativa C.5. El método cl´ asico de Runge-Kutta

47 47 48 50 53 55

Apéndice D.

59

La integral de Riemann–Stieltjes vii

ÍNDICE 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éndice 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éndice F. Laboratorio de Matemática: Integración 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éndice G. Laboratorio de Matemática: Solución 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éndice H. Existencia y unicidad de soluciones de ODEs H.1. Referencias

81 85

ÍNDICE GENERAL

Para S., N., y F.

ix

ÍNDICE 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ó de la necesidad de presentar una visión unificada de los métodos elementales de integración numérica de funciones y de algunos de los que son u ´tiles para las ecuaciones diferenciales ordinarias. ♣ Los p´ arrafos que se colocan entre tréboles pueden ser omitidos en las primeras lecturas puesto que no son esenciales para la continuidad del texto. ♣ Seguimos en términos generales las ideas planteadas en Linz∗ (1988). Para este autor el término matem´ atica computacional esta asociado con el espectro amplio de actividades relacionadas con la solución aproximada de problemas cient´ıficos expresados en términos de modelos matemáticos, que, en su forma t´ıpica, están 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ás simples que puedan ser resueltas por métodos algebraicos. Se puede distinguir la metodolog´ıa numérica, que estudia métodos para discretizar los operadores diferenciales e integrales y cómo resolver los sistemas finitos resultantes, del an´ alisis numérico que involucra el estudio riguroso de los algoritmos creados por la metodolog´ıa. La meta primaria del an´ alisis es describir la relación entre la solución exacta de la ecuaci´ on original, y la aproximada obtenida a partir de la versión discreta. El an´ alisis numérico ha resultado muy u ´til para complementar el proceso de creaci´ on de métodos numéricos generado por las aplicaciones ingenieriles. Por ejemplo las técnicas de relajaci´ on o de elementos finitos fueron creadas por ingenieros basados en la intuici´ on f´ısica. El análisis numérico posterior ha resultado crucial para la consolidaci´ on de estos métodos, la extensión de su aplicabilidad, y el estudio de sus propiedades de convergencia o estabilidad. En muchos casos a partir de objetivos del Análisis Numérico se han obtenido nuevos resultados en An´ alisis Funcional, en cuyos términos 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érminos como espacios de Hilbert, clausura compacta, y convergencia débil. Estos conceptos le resultan muy u ´tiles al teórico y han facilitado el establecimiento de una teor´ıa abarcativa y coherente de la solución aproximada de las ecuaciones de operadores. ♣ Ilustraremos esto con el ejemplo de la ecuación 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ón debe ser resuelta para la incógnita 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ón la ecuación (0.1) es reemplazada por una sucesi´ on paramétrica 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ón aproximada xn . El objetivo fundamental del análisis posterior es conocer la relación entre la soluci´ on verdadera x y su aproximación xn . En particular, desearemos demostrar la convergencia. Esto significa que, cuando se incrementa n, la solución aproximada deber´ a acercarse m´ as y m´ as a la solución verdadera, en el sentido que lim kx − xn k = 0

n→∞

(0.5)

Usualmente, el primer paso del análisis 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ón 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ón de estabilidad lim kL−1 n k≤K <∞

n→∞

(0.11)

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

PREFACIO

xv

Supongamos ahora que hemos determinado que nuestro método numérico es estable y convergente. ¿Qué nos dice esto? Nos expresa que, en alg´ un sentido asint´ otico, el método funciona. A menudo podemos decir más. En muchos casos el error de aproximación 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ón desconocida x El exponente p se denomina orden de convergencia y nos da información respecto de cuán exacto es el método. Para un dado nivel de esfuerzo computacional, un método 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ón 0.12 tiene una utilidad limitada puesto que depende de la inc´ ognita x. Se puede sortear esa dificultad realizando un análisis a posteriori que utiliza la solución 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ón calculada ρn (xn ) = Lxn − y

(0.14)

Es posible obtener cotas buenas para kρn k, en cambio es más 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érica no son el establecimiento del método de discretizaci´ on, sino el manejo de los detalles involucrados en la representaci´ on de dominios, de los propios métodos de discretización, de los métodos de operaci´ on con los sistemas de ecuaciones, de las matrices ralas, y de los métodos para resolver los sistemas de ecuaciones con dimensiones muy altas. En el desarrollo de métodos 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érica 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érica 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áticos sino desarrollar los más efectivos. En el primer cap´ıtulo se estudian algunos aspectos relevantes del cálculo de integrales. Se desarrollan los métodos clásicos, los de extrapolación 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étodo de Montecarlo por la simple razón que, en muchos casos, es el u ńico 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étodos elementales de los del tipo de ‘un paso’. En los apéndices 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álculos para obtener los polinomios de Bernoulli, que resultan necesarios para la interpretación de los coeficientes en el desarrollo asint´ otico del error en los métodos trapezoidal y Simpson.

xvi

PREFACIO

En el u ´ltimo de los apéndices se incluye un Laboratorio de Matemática 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ón de un Laboratorio propio. El objetivo de estos laboratorios es aspirar a la recreación viva de conceptos y temas que enriquecen el portafolios del pregraduado en Matemática. Los libros citados en la sección 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és 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érico con filosas preguntas y soluciones ingeniosas a los problemas propuestos. Asimismo a mis colegas por la discusión de temas y el desarrollo de métodos aptos para la ense˜ nanza. En cambio, de los errores que pudiera contener, me hago u ńico 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álculo Numérico I y II en la modalidad ‘a distancia’. Carlos E. Neuman Santa Fe, 12 de noviembre de 2001

xvii

Índice de Tablas 1

Extrapolaci´ on de Richardson (Método de Romberg)

1

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

2

1

12

Índice de Figuras 1 2 3 4 5

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

21 22 26 28 29

1 2

47

3

La funci´ on mihumps Histograma de distribución de los valores estimados de la integral del ejemplo C.3.1 Integraci´ on por Runge-Kutta clásico 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ón de quad para la función 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étodos numéricos para calcular integrales definidas de la forma Z b f (x) dx (1.1) a

Los métodos m´ as simples —que se originan en la propia definición de la integral de Riemann,— se establecen como sigue. Se elige una partición 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ón 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ón f en el intervalo [xi−1 , xi ]. En el caso de escoger fi = f (xi−1 ) se obtiene una “regla del rectángulo” 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ángulo Z

b

f (x) dx ' h a

N X

f (xi )

(1.4)

i=1

en realidad ‘la’ regla del rect´ angulo es más precisa si se elige fi = f ((xi−1 + xi )/2), en este caso se obtiene la regla del rectángulo 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ón 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étodo muy simple Observacio para programar estos métodos elementales en Matlab

La integral de RiemannStieltjes generaliza la integral de Riemann

De la inspecci´ on de estos métodos surgen varias ideas interesantes (1) Se tiene una ´ıntima relación entre integrales y sumas y las primeras se aproximan por las segundas. (2) Es necesario desarrollar métodos para estudiar los errores de aproximación de estas reglas de integración (es simple establecerlas pero no es tan claro c´ omo estimar los respectivos errores) (3) Se deber´ıan explorar métodos 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ág. 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étodos, en razón 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étodos 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ón, 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étodos adaptativos de integración 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ón 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ón Pε de [a, b] tal que, para cada partici´ on P más fina que Pε y para cada elección 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 ńico 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ágina 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ágina 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ágina 62 Demostracio En el apéndice 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ón generatriz h eh + 1 (1.26) 2 eh − 1 (para la interpretaci´ on detallada de los coeficientes ver el apéndice E en la p´ agina 67) 1 + c2 h2 + c4 h4 + c6 h6 + . . . =

1.3.1. El s´ımbolo ‘O’. La notación ‘O grande’ es u ´til para representar y calcular en forma exacta con magnitudes aproximadas. En general la notación O(f (n)) se usa —si f (n) es una función 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ón |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ón 1.29 (Justificar la validez de estas ecuaciones y la relaci´ on entre ambas). Las siguientes son operaciones válidas con la notación ‘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ón permite resolver la siguiente paradoja. Ejercicio 1.3.2. ¿Qué está mál 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ácil 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ón de la traslación 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órmula de la suma geométrica, 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éndice E (página 67) demostramos una f´ ormula de suma de Euler que enunciamos a continuación con modificaciones. Teorema 1.3.2. Supongamos que f es una función 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éndice. Se tiene la siguiente f´ ormula de integración 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órmula 1.53 para obtener por inducci´ on el teorema 1.3.1 de la página 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étodos b´ asicos simples y compuestos para la determinación numérica de la integral Z b f (x) dx (1.55) a

La deducci´ on de los mismos se basa en la aproximación de la función f en el intervalo (o, en el caso de los métodos 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étodos respecto

´ 1. INTEGRACION

8

de la que desarrollaremos a continuación, que se denomina “método de coeficientes indeterminados”

1.4.1. El m´ etodo de coeficientes indeterminados. Para deducir los métodos b´ asicos de integraci´ on numérica mediante este método 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á 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ón 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ón del método. En una secci´ on posterior (ver 1.6 en la página 12) se generalizará este método 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ón 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ón 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órmula 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ótico 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érico del método trapezoidal corregido ThC (f ) 1.4.3. El m´ etodo de Simpson. Utilizamos el producto Mathematica para determinar los coeficientes del método de Simpson (ver Apéndice B, página 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étodo 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én se obtiene a partir del desarrollo asintótico de Euler) εf n (f ) = −

1.4.4. M´ etodos de Newton. En la sección precedente hemos obtenido la f´ ormula del método de Simpson. La misma metodolog´ıa se puede utilizar para cualquier elecci´ on de los puntos de la partición del intervalo de integración. Por ejemplo, si desearamos calcular los que corresponden a una partición regular de once puntos, nuevamente utilizamos el producto Mathematica para obtener los coeficientes (ver apéndice B en página 37). Estos métodos llevan el nombre de métodos de Newton, y en algunos textos también 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ón que calculamos para el error en el método trapezoidal para ejemplificar el uso de los métodos de extrapolación 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étodo 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éndice, 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ón 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ón de Richardson (Método 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ás 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ón sea exacta para polinomios del mayor grado posible es resoluble y conduce a los métodos de integraci´ on de Gauß. En el apéndice B (p´ agina 37) incluimos un cuaderno del Mathematica que calcula los coeficientes para un método 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órmula 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én 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órmula sea exacta para todos los polinomios de grado m o menor (orden m + 1) ´ n 1.6.1. Aplicando la fórmula 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órmulas 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ólico como el Mathematica para obtener los coeficientes y para estimar el error de truncamiento local. En diversos manuales de fórmulas matemáticas figuran los puntos y coeficientes de Gauß para distintas familias de polinomios y distintos órdenes. El an´ alisis del error puede desarrollarse en forma análoga 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étodo 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ón de la función 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én, 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ón 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étodo se aplica entonces a la función fe(x) = f (x) − mf que resulta no negativa y con su gr´ afico contenido en el rectángulo R = [a, b] × [0, M f − mf ], y consiste en lo siguiente. Se hallan P puntos al azar en el rectángulo R y se determina el n´ umero Pf de los que se encuentran en el área que deseamos aproximar. A continuación se calcula la estimaci´ on del ´ area mediante la expresión Pf + mf (1.101) (b − a) (M f − mf ) P En la secci´ on C.3 incluimos ejemplos de integrales de Montecarlo. Este método no suele ser utilizado para el caso de áreas puesto que existen buenos métodos 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ón que pertenece a un sexto cuerpo de aquella que no está contenida en el mismo. El intento de obtener la integral por métodos clásicos o numéricos puede ser de una complejidad excesiva a´ un para casos de cuerpos definidos por ecuaciones muy simples. El uso del método de Montecarlo permite superar estas dificultades. 1.7.1. Determinaci´ on de un volumen. En esta sección calcularemos el volumen del s´ olido que resulta de la intersección, 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érmino atacamos el problema con Matlab, a continuación desarrollamos un cuaderno del Mathematica en el que obtenemos una solución 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éndice (F, página 75) se propone un laboratorio de Matemática orientado a estudiar algoritmos adaptativos para la integración numérica. En la secci´ on C.4, en la página 53, se incluye —para ilustración— el código de una funci´ on de Matlab, adaptada de la denominada quad, que implementa un algoritmo adaptativo recursivo de bajo orden basado en el método de Simpson. La

1.9. EJERCICIOS

17

funci´ on miquad es en realidad un driver que llama a la verdadera función recursiva que se denomina miquadst la que realiza cada paso de la recursión. 1.8.1. Comparaci´ on entre NIntegrate y quad. Calculamos el área 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ón 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ón, —la que, al margen, no es muy precisa,— parece que la primera metodolog´ıa es más eficiente que la segunda. 1.9. Ejercicios Ejercicio 1.9.1. Utilizar el método 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órmula " 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álculo de Z 1 sin x dx 1 +x 0

(1.105)

(1.106)

Ejercicio 1.9.3. (a) Deducir una fórmula de integración 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étodo numérico, 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órmulas de integración 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ánto 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é 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órmula de ‘midpoint’. Los respectivos errores asintóticos —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étodo numérico, Ien , con orden de convergencia mayor, combinando InT y InM . ¿Cuáles 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ón 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ón 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étrico podemos considerar que la ecuación 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ón diferencial y 0 = −xy, entonces y 0 /y = −x y es posible integrar ambos miembros obteniendo log y = −x2 /2 + C 0 , expresión de la que tomando la exponencial de ambos miembros se deduce la solución 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ón continua definida en el dominio D. Definicio Una funci´ on y(x), definida en el intervalo [x1 , x2 ], es una solución (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émoslos ξ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ón diferencial deseamos demostrar un resultado de existencia de las mismas que, además, resulta constructivo, puesto que nos da un método para obtener la aproximación Teorema 2.1.1. Sea (x0 , y0 ) un punto del dominio D (de la definición 2.1.1) y supongamos que el rect´ angulo R = {|x − x0 | ≤ a, |y − y0 | ≤ b} está 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ángulo S está contenido en el R ´ n: El rect´ Demostracio angulo S = {|x − x0 | ≤ h, |y − y0 | ≤ M h} está 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á 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ángulo OQP de la figura 2 porque el ángulo α se toma de modo que tan α = M , y |f (x0 , y0 )| ≤ M . Por motivos análogos (x2 , y2 ) también est´ a en OQP ; etc. En consecuencia el proceso puede continuarse hasta xn = x0 + h, porque la u ńica razón 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órmulas 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ás, 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ón previa y la construcci´ on requerida por el teorema ha sido completada.(ver apéndice H) Este método de construir una solución aproximada se conoce como método de Euler. Es innecesario mejorar el valor del h del teorema porque, en general, se puede demostrar que y(x) está 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étodo de Euler para la resoluci´ on aproximada de la ecuación 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étodo en estudio.

(2.13)

24

Estabilidad del m´ etodo de Euler


En cada paso se cambia de miembro de la familia de soluciones de la ecuación diferencial, de manera que la precisión del método deberá 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ón de estabilidad de un método numérico en función de su comportamiento para resolver la ecuaci´ on diferencial y 0 = λy con λ un n´ umero complejo. La región 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 ńico valor yn produce una sucesión de cambios en los siguientes valores que no crecen de paso a paso. Ejemplo 2.2.1. Para el método de Euler se tiene yn+1 = yn + λhyn = (1 + λh)yn

(2.14)

de modo que este método es absolutamente estable en la región |1 + λh| ≤ 1 que es el c´ırculo unitario del plano complejo centrado en el punto (−1, 0). Ejemplo 2.2.2. El método de Euler aplicado a la ecuación del ejemplo 2.1.1 resulta ser yn+1 = yn − hxn yn = (1 − hxn )yn (2.15) El método de Euler puede obtenerse también del siguiente razonamiento. Se reemplaza la derivada y 0 en la ecuación 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étodo al despejar la ecuación 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étodo 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ón de una ecuaci´ on en cada paso. Ejemplo 2.2.3. Para el método 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étodo es absolutamente estable en la región |1/(1 − λh)| ≤ 1 v´ alida en el semiplano de partes reales no positivas del plano complejo. Estabilidad del m´ etodo retroEuler


25

Ejemplo 2.2.4. El método retroEuler aplicado a la ecuación 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étodo 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étodo 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é?) y que aplicando la extrapolación (de Richardson) para eliminar D(x) y(x) = 2yh/2 (x) − yh (x) + O(h2 )

(2.26)

Apliquémosla 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ón, 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é consiste esta primera fórmula de Runge en forma gráfica? En la figura 3 representamos el efecto de aplicar la primera fórmula de Runge en el paso de xn a xn+1 . Notemos que para el método de Euler se puede definir en forma análoga ϕ(x, y, h; f ) = f (x, y). Ejercicio 2.2.5. ¿Es posible expresar en esta notación el método de retroEuler?

Primera f´ ormula de Runge

26


y 6

0

10 rr Runge un (xn+1 )

r Euler

xn

xn+1

x

Figura 3. Primer método de Runge Para estudiar la estabilidad aplicamos la fórmula a la ecuación y 0 = λy. Se tiene yn+1 = yn + hλ(yn + (h/2)λyn ) = yn (1 + hλ + 12 h2 λ2 ) razón por la cual debe tenerse que |1 + hλ + 12 h2 λ2 | ≤ 1. Este primer método de Runge es generalización del de Euler. Aunque existen otras posibles generalizaciones, por ejemplo utilizando series de potencias, el problema es que hay que derivar la función f . La idea de Runge fue dar fórmulas 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érminos con la serie de potencias desarrollada en el entorno del punto. Esas f´ ormulas fueron mejoradas por Heun y por Kutta. Llamaremos método 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étodo ε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 éste método ε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


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érminos 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ápido como se pueda. Aceptando f arbitrarias no podremos, en general, anular el coeficiente de h3 (¿porqué?). 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étodo puede mejorarse definiendo hϕ =

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

(2.53)

(2.54)

Segunda f´ ormula de Runge

28


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

r Euler

0

xn

xn+1 x d

Figura 4. Segunda fórmula 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 álgebra más 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érminos en el error de truncamiento local sean nulos. En consecuencia, por un procedimiento análogo al realizado en el caso de segundo orden se obtienen diversos métodos de orden superior de los que el m´ as utilizado es el llamado Método de Runge-Kutta 1 yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) 6

(2.61)


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étodo de Runge-Kutta Se puede demostrar que ésta es una fórmula de cuarto orden: εn+1 = O(h5 ). ´ n 2.2.1. El método RK clásico es una generalización 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étodo clásico de RK a la ecuación 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


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ón 2.75 se puede dibujar la región 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ón ε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ón de consistencia para el método 2.76 y es suficiente para asegurar la convergencia si se tienen ciertas hipótesis 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étodo 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ón 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)


31

donde τ (h) = maxx0 ≤xn ≤b |τn+1 (Y )|. Si se satisface la condición 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érica {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álculo 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ón. Corolario 2.2.3. Si el método 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


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étodo clásico 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álogo 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érmino 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étodo buscamos para en una expresión de la Demostracio forma del enunciado. Si reemplazamos la expresión 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étodo numérico 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ón 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étodo de Euler con extrapolación de repetida de Richardson. Usar aritmética de redondeo de cinco decimales. Ejercicio 2.3.2. Aplicando el método 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étodo de Euler con longitud de paso h sobre el problema test y 0 = −y, y(0) = 1 (2.107) (a) Determinar una expresión expl´ıcita para yn (b) ¿Para qué valores de h la sucesión {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ón de la ecuaci´ on y 0 = y 2 , con y(0) = 1, en el entorno de x = 0. Usar esta aproximación para calcular y(0.2) e y(1.2) con cuatro decimales. Comparar con la solución exacta y explicar porqué el segundo caso (x = 1.2) no es exitoso. Ejercicio 2.3.5. La función 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étodos: (a) (b) (c) (d)

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

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

(2.109)

34


para varios valores de x. Se procede de la siguiente manera: y(x) se calcula para x = 1 usando alg´ un método para integración numérica; se calcula y(1) = 0.6051. Mostrar que y satisface la ecuación diferencial dy 1 √ + 2xy = − + π (2.110) dx x Se resuelve la ecuaci´ on diferencial numéricamente 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étodo 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érico, 2da edición, Grupo Editorial Iberoamérica, México, 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ón, Reverté, 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é, 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érica de Ecuaciones Diferenciales, Tomos 1 y 2, Reverté, Buenos Aires, 1985 y 1986 35

36

A. BIBLIOGRAFÍA 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


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


(* 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


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


-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ón (* 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[--]]


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


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ón 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ón. 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


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ón de Romberg para la función 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ón 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étodo de Montecarlo el área bajo la función f (x) = 0.5 en el cuadrado unitario. Realizamos mil veces el cálculo del área 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ón de puntajes obtenidos por la rutina de Montecarlo. Ejemplo C.3.2. Calulamos la integral de la función ‘humps’ mediante el método 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


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ón 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


%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ódigo Matlab que implementa el método 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étodo. 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 % %%% %%% %%%


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

MIEQ

56


% % % % % % % % % %

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.

%%% %%% %%%


%%% 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ón por Runge-Kutta clásico de x0 = −tx en [0, 7]

´ APENDICE D

La integral de Riemann–Stieltjes Comencemos estudiando las propiedades de las funciones de variación acotada. D.1. Funciones de variaci´ on acotada Teorema D.1.1. [Apostol (1976), pág. 153] Sea f una función 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ón 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ón. Teorema D.1.2. [Apostol (1976), pág. 153] Si f es monótona 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ág. 154] Si f es monótona en [a, b], entonces f es de variaci´ on acotada en [a, b]. Teorema D.1.4. [Apostol (1976), pág. 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ág. 155] Si f es de variación acotada en [a, b], es decir que |∆fk | ≤ M para toda partición de [a, b], entonces f está acotada en [a, b]. ´ n: Sea x ∈ (a, b) y la partición 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ón de variación 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ón 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ón de variación 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ón 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ón 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ón 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ón 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ón 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ón Pε de [a, b] tal que, para cada partici´ on P más fina que Pε y para cada elección 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 ńico 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ágina 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ágina 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ón 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ón que ella se tenga

(1) Pε

Z |S(P, f1 , g) −

tal que para toda partición P más 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ón P más 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ón P más 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én es integrable respecto de g2 entonces, para todo par de constantes c1 y c2 , se tiene que f también 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áloga al D.2.2. Demostrar adem´ as la aditividad respecto del intervalo de integración 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ón P˜ más 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ás 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álida 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ón integrable en el sentido de RiemannStieltjes respecto de una función 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ótesis de integrabilidad de f respecto de g se tiene que, para toda partición P más 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ón 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.2.2.1. Cambio de variables en la integral de Riemann-Stieltjes. Si la función 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órmula 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
(D.39)

(donde c1 , c2 , y c3 son arbitrarios. Sea f una función 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ón 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álida (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ón 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ón escalonada definida en el intervalo [a, b] con salto gk en el punto xk de una partición {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ón. D.4. Pero, ¿Hay funciones integrables? Se puede demostrar que si f es continua y g es de variación acotada en el intervalo [a, b] entonces f es integrable en el sentido de Riemann-Stieltjes respecto del integrador g. La demostración es muy simple y se basa en la definición 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.5. Ejercicios Ejercicio D.5.1. Una función f , definida en [a, b], verifica una condición 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ón, probar que α > 1 implica que f es constante en [a, b], mientras que α = 1 implica que f es de variación acotada en [a, b] (b) Dar un ejemplo de una función f que satisfaga una condición uniforme de Lipschitz de orden α < 1 en [a, b] tal que f no sea de variación acotada en [a, b] (c) Dar un ejemplo de una función f que sea de variación acotada y que, sin embargo, no satisfaga ninguna condición uniforme de Lipschitz en [a, b]. Ejercicio D.5.2. Probar que una función polinómica f es de variación acotada en todo intervalo cerrado y acotado [a, b]. Describir un método 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ón de la integral de Riemann-Stieltjes es bastante usual en textos matemáticos: 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ón P de [a, b] con norma kP k < δ y cada elección 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ón entonces también existe de acuerdo con la definición D.2.1 (página 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ón 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ón 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


70

0.15 0.1 0.05

0.2

0.4

0.6

0.8

-0.05


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


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


72

0.06 0.04 0.02 0.2

0.4

0.6

0.8

1

-0.02 -0.04 -0.06


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


(* 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


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ón E.6 se toma como definición 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ón 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órmula de integración 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étodo 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étodo adaptativo para la resolución numérica de distintos problemas de Matemática, y, en particular, para el caso de la integraci´ on de funciones • Comparar la eficiencia de distintos métodos de integración numérica • Estudiar los errores de aproximación asociados a distintos métodos de integraci´ on numérica F.2. Antes del laboratorio La idea b´ asica de un método adaptativo es refinar los cálculos 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ón uniforme de los errores de aproximación 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ón de quad para la función 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étodo de Simpson es exacto y la distribuci´ on de puntos es uniforme y mucho menos densa. Ejercicio F.2.1. Dada una función f definida en un intervalo [a, b] dar un método adaptativo de interpolación de f mediante polinomios a trozos de orden 2 (splines lineales). Ejercicio F.2.2. Dada una función f definida en un intervalo [a, b] dar un método adaptativo de interpolación 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étodo de integración numérica adaptativa de una funci´ on dada en un intervalo cerrado y acotado del eje real. Ejercicio F.3.2. Realizar experimentos numéricos destinados a determinar las caracter´ısticas del comportamiento del programa desarrollado en F.3.1 para distintos integrandos y distintos dominios de integración. Comparar con otros métodos de integraci´ on numérica.

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érico, 2da edición, Grupo Editorial Iberoamérica, México, 1993.

´ APENDICE G

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

Las metas del laboratorio de Solución adaptativa de ODEs son • Desarrollar el concepto de método adaptativo introducido en otro laboratorio de este texto para la resolución numérica de distintos problemas de Matem´ atica, y, en particular, para el caso de la solución de ecuaciones diferenciales • comparar la eficiencia de distintos métodos de solución numérica de ODEs • estudiar los errores de aproximación asociados a distintos métodos de soluci´ on numérica de ODEs. G.2. Antes del laboratorio La idea b´ asica de un método adaptativo es refinar los cálculos 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ón uniforme de los errores de aproximación. Analizar en detalle los resultados del laboratorio de integración adaptativa. G.3. Laboratorio Ejercicio G.3.1. Desarrollar un algoritmo —y correspondiente programa en Matlab— que implemente un método de integración numérica adaptativa de una ecuaci´ on diferencial dada en un intervalo del eje real. Ejercicio G.3.2. Realizar experimentos numéricos 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ón. Comparar con otros métodos de soluci´ on numérica 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érico, 2da edición, Grupo Editorial Iberoamérica, México, 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ón 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ón. 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á 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á y(x) ´ n H.0.3. La función y(x) satisface la ecuación 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áficos en R.


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


Hemos demostrado as´ı que el problema de valores iniciales tiene una solución continua en x0 ≤ x ≤ x0 + α Teorema H.0.4. La solución es u ńica ´ n: Supongamos que además de y(x) se tiene que z(x) también 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ón no negativa con Z x w(x) ≤ w(s)ds

(H.42)

x0

Entonces w(x) es idénticamente 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étodo 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érico, 2da edición, Grupo Editorial Iberoamérica, México, 1993.

H.1. REFERENCIAS

Terminado de imprimir el doce de noviembre de dos mil uno en el Departamento de Matemática (FIQ), Universidad Nacional del Litoral, en la ciudad de Santa Fe, Rep´ ublica Argentina.

87

Analisis Numerico de Integrales y Ecuaciones Diferenciales

Recommend Documents