Análisis Numérico de Integrales y Ecuaciones Diferenciales

Este documento es de distribución gratuita y llega gracias a

“Ciencia Matemática” www.cienciamatematica.com El mayor portal de recursos educativos a tu servicio!

ISSN 1666-3845

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

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

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

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

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

i

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

An´ alisis Num´ erico de Integrales y Ecuaciones Diferenciales

¶

Tercera edici´ on

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

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

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

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

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

ISSN 1666-3845

´Indice General Prefacio

xiii

Prefacio a la segunda edici´ on

xvii

´Indice de Tablas

1

´Indice de Figuras

3

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

1 1 3 4 7 11 12 13 16 17

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

21 21 23 33

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

35 35 35

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

37 37 38 40 41

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

47 47 48 50 53 55

Ap´endice D.

59

La integral de Riemann–Stieltjes vii

´INDICE GENERAL

viii

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

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

59 61 65 65 66

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

67 67 73 74

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

75 75 75 76 77 77

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

79 79 79 79 79 80

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

81 85

´INDICE GENERAL

Para S., N., y F.

ix

´INDICE GENERAL

xi

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

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

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

(0.1)

xiv

PREFACIO

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

(0.4)

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

n→∞

(0.5)

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

(0.6)

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

(0.7)

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

(0.8)

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

n→∞

(0.9)

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

n→∞

(0.10)

Se necesita una condici´ on adicional, la condici´on de estabilidad lim kL−1 n k≤K 0, existe una partici´on Pε de [a, b] tal que, para cada partici´ on P m´as fina que Pε y para cada elecci´on de los puntos ξk del intervalo [xk−1 , xk ], k = 1, . . . , n, se tiene |S(P, f, g) − I| < ε

(1.13)

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

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

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

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

´ 1. INTEGRACION

4

Ejemplo 1.2.3. Si f es continua y g(x) = bxc entonces la integral resulta X f (i) (1.16) i∈Z a 0, h = (b−a)/n, tj = a+jh, para j = 0, 1, . . . , n. Sea f (x) definida por la propiedad de ser lineal en cada subintervalo [tj−1 , tj ], para j = 1, . . . , n. Mostrar que este conjunto de funciones f , para todo n ≥ 1, es denso en el conjunto C[a, b] de funciones continuas en el intervalo [a, b]. Ejercicio 1.9.6. Obtener f´ormulas de integraci´on de Gauß para Z 1 n X I= xf (x) dx ' wj f (xj ) 0

(1.109)

j=1

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

1.9. EJERCICIOS

n 2 4 8 16 32 64

19

In 0.28451779686 0.28559254576 0.28570248748 0.28571317731 0.28571418363 0.28571427643

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

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

(1.110)

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

CAP´ıTULO 2

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

(2.1)

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

El campo de direcciones permite estimar gr´ aficamente las trayectorias

y 6 3 2 1

0

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

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

E Er E B E Br B B

2

3

x

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

y(x) = Ce−x

/2

(2.2)

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

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

22

2. ECUACIONES DIFERENCIALES

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

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

(2.3)

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

i = 1, . . . , n

(2.4)

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

23

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

(2.5)

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

xi−1 ≤ x ≤ xi ,

i = 1, . . . , n

(2.6)

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

xi−1 < x < xi ,

i = 1, . . . , n

(2.7)

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

(2.8)

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

δ =δ M

(2.9)

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

(2.10)

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

x 6= xi ,

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

(2.11)

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

M´ etodo de Euler

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

n≥0

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

(2.13)

24

Estabilidad del m´ etodo de Euler

2. ECUACIONES DIFERENCIALES

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

(2.14)

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

(2.16)

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

(2.18)

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

(2.19)

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

(2.20)

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

1 yn 1 − λh

(2.21)

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

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

25

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

1 yn (1 + hxn+1 )

(2.23)

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

(2.24)

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

(2.25)

y que para el paso mitad

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

(2.26)

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

(2.27)

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

(2.28)

y para medio paso

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

(2.29)

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

(2.30)

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

(2.31)

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

(2.32)

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

(2.33)

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

Primera f´ ormula de Runge

26

2. ECUACIONES DIFERENCIALES

y 6

0

10 rr Runge un (xn+1 )







r Euler











xn

xn+1

x

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

n≥0

(RK)

(2.34)

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

n≥0

(2.35)

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

(2.36)

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

n≥0

(2.37)

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

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

(2.38)

Estabilidad de la primera f´ ormula de Runge

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

27

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

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

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

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

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

⇒

ϕ(x, y, h; f ) =

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

(2.46) (2.47) (2.48)

(2.49)

(2.50)

(2.51) (2.52)

se tiene

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

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

(2.53)

(2.54)

Segunda f´ ormula de Runge

28

2. ECUACIONES DIFERENCIALES

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







r Euler

 





 







  0

xn

xn+1 x  d

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

y 6

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

(Euler)

(2.56)

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

p X

γi ki

(2.57)

i=1

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

(2.58) (2.59)

(2.60)

j=1

M´ etodo cl´ asico Runge-Kutta

de

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

(2.61)

´ 2.2. METODOS EN DIFERENCIAS DE UN PASO

29

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

y 6

0

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

r RK

r un (xn+1 )









r Euler











xn

xn+1 x

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

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

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

30

2. ECUACIONES DIFERENCIALES

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

(2.74)

(2.75)

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

n≥0

(2.76)

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

n≥0

(2.77)

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

(2.78)

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

n≥0

(2.79)

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

x0 ≤x≤b −∞ ------}} 3 6 6 (* Verifico que vale para polinomios de tercer grado *) (s3=A p3[a0] + B p3[a1] + C p3[a2])/.%15 //Simplify {0} i3=Integrate[p3[x],{x,a,b}]//Simplify 0 (* Resulta as\’{\i} que la siguiente ecuaci\’on ( s3 == i3 ) vale trivialmente *) (* Calculo del error con el polinomio de cuarto grado *) (s4=A p4[a0] + B p4[a1] + C p4[a2])/.%15 //Simplify -a + b {------} 3 i4=Integrate[p4[x],{x,a,b}]//Simplify -a + b -----5 e4=(%20 - %19)//Simplify 2 (a - b) 37

38

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

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

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

*)

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

´ B.2. EL METODO DE NEWTON

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

39

40

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

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

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

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

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

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

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

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

41

42

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

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

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

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

43

44

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

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

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

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

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

. . .

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

45

46

B. COEFICIENTES INDETERMINADOS CON EL Mathematica

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

´ APENDICE C

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

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

1

0.8

0.6

0.4

0.2

0.2

0.4

0.6

0.8

Figura 1. La funci´on mihumps

47

1

48

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

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

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

0.01 2

0.04 + (−0.9 + x)

+

Rx 0

f (s) ds son

0.01 2

0.01 + (−0.3 + x)

(C.1)

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

´ REPETIDA DE RICHARDSON C.2. EXTRAPOLACION

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

C.E.Neuman, 6-5-97

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

49

50

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

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

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

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

fin de mirichar.m

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

0.29819609270927 0.29858264933679 0.29858326161685 0.29858325443074 0.29858325398467 0.29858325395684 0.29858325395510 0

0.29860841977863 0.29858330243552 0.29858325395166 0.29858325395493 0.29858325395499 0.29858325395499 0 0

0.29858290374753 0.29858325318208 0.29858325395498 0.29858325395499 0.29858325395499 0 0 0

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

0.29858325395743 0.29858325395498 0.29858325395499 0 0 0 0 0

0.29858325395498 0.29858325395499 0 0 0 0 0 0

0.29858325395499 0 0 0 0 0 0 0

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

´ C.3. METODOS DE MONTECARLO

51

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

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

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

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

% histograma de los valores

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

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

52

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

250

200

150

100

50

0 0.48

0.485

0.49

0.495

0.5

0.505

0.51

0.515

0.52

0.525

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

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

´ DE SIMPSON ADAPTATIVA C.4. INTEGRACION

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

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

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

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

53

54

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

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

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

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

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

55

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

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

MIRK

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

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

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

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

MIEQ

56

C. CUADERNOS DEL Mathematica Y ARCHIVOS SCRIPT DE MATLAB

% % % % % % % % % %

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

%%% %%% %%%

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

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

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

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

%llamada a R-K

[tt+h; xn; xdn] ];

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

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

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

57

Integracion de dx/dt=−tx 1

0.8

0.6 2

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

x

0.2

0

−0.2

−0.4

−0.6

−0.8

0

1

2

3

4

5

6

7

t

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

´ APENDICE D

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

(D.1)

Tenemos entonces la desigualdad n−1 X

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

(D.2)

k=1

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

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

k=1

n−1 X

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

(D.3)

k=1

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

(D.4)

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

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

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

k=1

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

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

60

D. LA INTEGRAL DE RIEMANN–STIELTJES

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

(D.6)

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

(D.7)

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

(D.8)

con lo cual

y, en consecuencia,

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

(D.9)

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

n X

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

(D.10)

k=1

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

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

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

D.2. LA INTEGRAL

61

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

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

(D.13)

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

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

Nota: Este resultado se demuestra en el teorema D.2.5 en la p´agina 63 Ejemplo D.2.3. Si f es continua y g(x) = bxc entonces la integral resulta X f (i) (D.16) i∈Z a≤i 0 dado, elegimos una partici´on que ella se tenga

(1) Pε

Z |S(P, f1 , g) −

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

f1 dg| < ε a

(D.19)

(D.20)

62

D. LA INTEGRAL DE RIEMANN–STIELTJES (2)

y elegimos otra partici´ on Pε tenga

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

Z |S(P, f2 , g) −

f2 dg| < ε

(D.21)

a

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

(D.22)

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

(D.23)

a

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

a

a

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

a

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

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

n X

g(ξk )∆fk =

k=1

n X

g(ξk )f (xk ) −

k=1

n X

g(ξk )f (xk−1 )

(D.27)

k=1

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

n X k=1

f (xk )g(ξk ) −

n X k=1

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

(D.28)

D.2. LA INTEGRAL

63

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

k=1

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

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

Rb a

g(x)dg(x) existe y vale

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

a

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

(D.32)

k=1

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

n X

f (ξk )∆gk

k=1

n X

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

(D.33)

k=1

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

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

lo que demuestra el teorema.

64

D. LA INTEGRAL DE RIEMANN–STIELTJES

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

a

D.2.2.2. Integral de funciones escalonadas. Teorema D.2.6. Definimos g en el intervalo [−1, 1] como sigue   c1 si −1 ≤ x < 0 c2 si x=0 g(x) =  c3 si 0 0 existe un δ > 0 tal que, si kP k < δ, entonces |f (ξk−1 ) − f (0)| < ε (D.44) y |f (ξk ) − f (0)| < ε (D.45) de modo que se tiene la desigualdad |S(P, f, g) − f (0)(g(0+) − g(0−))| ≤ ε|g(0) − g(0−)| + ε|g(0+) − g(0)|

(D.46)

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

(D.47)

D.4. PERO, ¿HAY FUNCIONES INTEGRABLES?

65

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

(D.48)

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

k=1

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



1 f (x) d x − bxc − 2



Z + a

b



1 x − bxc − 2

 df (x) =

f (a) − f (b) 2

(D.51) (D.52)

(D.53)

por u ´ltimo Z a

b

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

(D.54)

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

66

D. LA INTEGRAL DE RIEMANN–STIELTJES

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

´ APENDICE E

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

(E.1)

k=0

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

X t tk = B k et − 1 k!

(E.2)

k=0

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

x

5 x

6 x

7 x 67

68

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

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

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

E.1. POLINOMIOS DE BERNOULLI CON EL Mathematica

0.4 0.2

0.2

0.4

0.6

0.8

-0.2 -0.4

Figura 1. El polinomio de Bernoulli B1 (x)

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

1

69

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

70

0.15 0.1 0.05

0.2

0.4

0.6

0.8

-0.05

Figura 2. El polinomio de Bernoulli B2 (x)

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

1

E.1. POLINOMIOS DE BERNOULLI CON EL Mathematica

0.02 0.01

0.2

0.4

0.6

0.8

1

-0.01 -0.02

Figura 3. El polinomio de Bernoulli B6 (x)

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

71

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

72

0.06 0.04 0.02 0.2

0.4

0.6

0.8

1

-0.02 -0.04 -0.06

Figura 4. El polinomio de Bernoulli B10 (x)

6 1649024253633120910 x

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

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

18 702064548199300 x

20 - 72937940132694 x

24 445756964055 x

26 + 27074751480 x

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

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

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

´ E.2. NUMEROS DE BERNOULLI

73

-Graphics- (*berno5.wmf*)

16

2·10

16

1·10

0.2

0.4

0.6

0.8

1

16

-1·10

16

-2·10

Figura 5. El polinomio de Bernoulli B40 (x)

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

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

(E.4)

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

74

´ E. POLINOMIOS Y NUMEROS DE BERNOULLI

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

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

(E.6)

k=0

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

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

(E.7)

(E.8) (E.9)

k=0

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

´ APENDICE F

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

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

C.E. Neuman, 30 de mayo de 1997

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

f (x) =

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

720(x − 12 )

si

x<

1 2

si

x≥

1 2

(F.1)

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

75

Hay posibles singularidades en el entorno de 0 y de 0.5

76

´ ´ ADAPTATIVA F. LABORATORIO DE MATEMATICA: INTEGRACION

400

350

300

250

200

150

100

50

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

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

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

F.5. REFERENCIAS

77

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

´ APENDICE G

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

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

79

80

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

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

´ APENDICE H

Existencia y unicidad de soluciones de ODEs Dado el PVI

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

(H.1)

y las iteraciones de Picard Z

x

yn+1 = y0 +

f (s, yn (s))ds

(H.2)

x0

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

(H.3)

(x,y)∈R

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

(H.4) (H.5)

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

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

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

82

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

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

En particular Z

x

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

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

≤ L

M (s − y0 )ds = LM x0

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

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

(H.12)

x

≤ L

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

(H.13)

x0

≤ M L2

Z

x

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

x0

(H.14)

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

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

x0 ≤ x ≤ x0 + α

(H.15)

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

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

≤ M (x − x0 ) + M L

j=1

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

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

≤ Mα + ML =

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

(H.19)

x0

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

Z yn+1 (x) = y0 +

f (s, yn (s))ds

(H.20)

x0

tomamos l´ımite: Z

x

y(x) = y0 + lim

n→∞

f (s, yn (s))ds

(H.21)

x0

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

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

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

83

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

x0

pero y(s) − yn (s) =

P∞

− yj−1 (s)) y entonces

j=n+1 (yj (s)

∞ X

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

=

Lj−1

(s − x0 )j j!

(H.25)

Lj−1

αj j!

(H.26)

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

(H.27)

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

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

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

≤ Mα

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

(H.28)

(H.29)

Resulta as´ı que Z

t

lim

n→∞

Z

t

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

f (s, y(s))ds

(H.30)

t0

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

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

(H.32)

j=N +1

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

(H.33) (H.34)

(H.35)

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

84

H. EXISTENCIA Y UNICIDAD DE SOLUCIONES DE ODES

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

Z z(x) = y0 +

f (s, z(s))ds

(H.38)

x0

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

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

x0

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

(H.42)

x0

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

Rx x0

w(s)ds entonces

dU = w(x) ≤ dx

Z

x

w(s)ds = LU (x)

(H.43)

x0

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

(x ≥ x0 )

(H.44)

o bien U (x) = 0

(H.45)

de modo que Z

x

0 ≤ w(x) ≤ L

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

(H.46)

x0

o tambien w(x) = 0

(H.47)

H.1. REFERENCIAS

85

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

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

Ek+1 = (1 + hL)Ek +

+

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

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

h2 2 (fx

Dh2 2

(H.49)

(H.50)

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

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

como kh < α Ek ≤ eαL E0 +

D αL (e − 1)h 2L

(H.51)

(H.52)

(H.53)

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

H.1. REFERENCIAS

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

87