Cálculo fraccionario y dinámica newtoniana

Trabajo de Fin de Máster Cálculo fraccionario y dinámica newtoniana REALIZADO POR: Antón Lombardero Ozores DIRIGIDO POR: Juan José Nieto Roig DEPARTAMENTO: Análisis Matemático Para los estudios de: Máster Universitario en Matemáticas Curso 2012 - 2013 Santiago de Compostela, septiembre de 2013 Trabajo de Fin de Máster Cálculo fraccionario y dinámica newtoniana REALIZADO POR: Antón Lombardero Ozores DIRIGIDO POR: Juan José Nieto Roig DEPARTAMENTO: Análisis Matemático Para los estudios de: Máster Universitario en Matemáticas Curso 2012 - 2013 Santiago de Compostela, septiembre de 2013 Juan José Nieto Roig, Catedrático del Departamento de Análisis Matemático de la Universidad de Santiago de Compostela, ha dirigido el presente Trabajo de Fin de Máster de título Cálculo Fraccionario y dinámica newtoniana , proyecto destinado a los estudios del Máster Universitario en Matemáticas de la USC (curso 2012-2013). El Director del TFM: Prof. Juan José Nieto Roig El alumno: Antón Lombardero Ozores Índice general 1. Introducción 9 2. Cálculo Fraccionario 11 2.1. Historia .................................. 11 2.2. Operadores fraccionarios . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1. Integral fraccionaria de Riemann-Liouville . . . . . . . . . . . 13 2.2.2. Derivada fraccionaria de Riemann-Liouville . . . . . . . . . . . 15 2.2.3. Derivada fraccionaria de Caputo . . . . . . . . . . . . . . . . . 20 2.2.4. Derivada fraccionaria de Grünwald-Letnikov . . . . . . . . . . 22 2.2.5. Otros operadores fraccionarios . . . . . . . . . . . . . . . . . . 24 2.3. Ecuaciones diferenciales fraccionarias . . . . . . . . . . . . . . . . . . 24 2.3.1. Deniciones y conceptos . . . . . . . . . . . . . . . . . . . . . 25 2.3.2. Teoremas de existencia y unicidad . . . . . . . . . . . . . . . . 30 2.3.3. Soluciones de casos particulares . . . . . . . . . . . . . . . . . 31 3. Aplicaciones dinámicas 37 3.1. Péndulosimple .............................. 37 3.1.1. Ecuación diferencial fraccionaria . . . . . . . . . . . . . . . . . 38 3.1.2. Solución.............................. 39 3.2. Proyectil.................................. 41 3.2.1. Ecuaciones diferenciales fraccionarias . . . . . . . . . . . . . . 42 3.2.2. Soluciones............................. 43 3.2.3. Características . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3. Resorteomuelle ............................. 47 3.3.1. Ecuación diferencial fraccionaria . . . . . . . . . . . . . . . . . 47 3.3.2. Soluciones............................. 48 Bibliografía 51 7 8 Capítulo 1 Introducción Desde el punto de vista matemático, el estudio del movimiento de un cuerpo a partir de las fuerzas que actúan sobre él se reduce a aplicar la ecuación fundamental de Newton F=ma , teniendo en cuenta que F es la fuerza resultante o suma vectorial de las fuerzas que intervienen sobre el cuerpo. Se trata, en suma, de resolver la ecuación diferencial de segundo orden F=mx00, (1.0.1) donde x es la función que determina el movimiento. Puede ser interesante, como experimento físico-matemático, seguir este proceso desde una perspectiva más amplia, que utilice como punto de partida una versión generalizada de la segunda Ley de Newton. Esto es lo que se pretende con este trabajo, y para ello haremos uso del Cálculo Fraccionario. El Cálculo Fraccionario es una activa rama del Análisis Matemático que nace de una idea muy básica: de la misma forma que una función se puede derivar o integrar un número entero de veces, ¾podemos hallar su derivada de orden 1 /2 ? ¾Y su integral π - ésima? Surgen así operadores diferenciales e integrales de orden racional, real e incluso complejo. Con estos cimientos se puede construir un nuevo edicio del Cálculo, que, sin perder de vista al Cálculo tradicional (al que no sustituye, sino que extiende), posee sus propios conceptos fraccionarios . Un ejemplo, entre otros muchos, es el de ecuación diferencial fraccionaria. Como en cualquier área matemática viva, el interés del Cálculo Fraccionario no se restringe al plano teórico, siendo este notable. Los campos de aplicación práctica se incrementan año tras año. En particular, los modelos de derivadas fraccionarias proporcionan mejores representaciones que los modelos tradicionales en aquellos materiales con amortiguación interna. 9 Riemann-Liouville de orden α∈C(<(α)≥0) de f se dene, si existe 1 , por (Dα af) (x):=d dxnIn−α af(x) =1 Γ (n−α)d dxnZx a f(t) (x−t)α−n+1 dt (x>a) (2.2.4) donde n= [<(α)] + 1 y d dxn es la n -ésima derivada usual. Observación 13 . Al contrario que la integral fraccionaria, la derivada fraccionaria de Riemann-Liouville sí permite órdenes imaginarios puros. Cuando <(α) = 0 , la derivada de orden α=iθ se expresa de la siguiente manera: Diθ af(x) = 1 Γ (1 −iθ) d dx Zx a f(t) (x−t)iθ dt (θ∈R\{0};x>a) (2.2.5) Observación 14 . Cuando α es un número natural, α∈N , se tiene que n=α+ 1 y la derivada fraccionaria coincide con la derivada usual: (Dα af) (x) = d dxα+1 I(α+1)−α af(x) = d dxα+1 (Iaf) (x) = d dxα f(x). En particular D0 a será el operador identidad, i.e., D0 af(x) = f(x) . Observación 15 . Algo llama poderosamente la atención en la denición (2.2.4): se especica un límite inferior de integración a . La derivada fraccionaria de RiemannLiouville, al igual que otras deniciones de derivada fraccionaria, es un operador no local . Contrariamente al caso entero, esta derivada queda denida por medio de una integral que depende de los valores que la función tome a lo largo de un intervalo . Solamente cuando α es natural la derivada fraccionaria se convierte en un operador local. Las derivadas de orden fraccionario contienen parcial o totalmente la historia temporal o el comportamiento espacial de la función, promediadas de una cierta forma. Esto convierte a las ecuaciones diferenciales fraccionarias en candidatas idóneas para la modelización de fenómenos con memoria , aquellos en los que lo que ocurre en un punto del espacio o en un instante de tiempo depende de un intervalo (espacial o temporal) que contiene al punto o al instante. Al igual que en el caso entero, la derivada fraccionaria de una función f(x) existe bajo unas condiciones más restrictivas que la correspondiente integral. No basta, como 1 En el Teorema 20 se dan condiciones sucientes de existencia. 16 ocurría con Iα af , con que f(x) sea de L1 . Antes de enunciar el resultado principal sobre existencia, necesitamos la noción de continuidad absoluta. Denición 16. Una función f: [a, b]−→ C se dice que es absolutamente continua en [a, b], f ∈AC [a, b] , si para cualquier ε > 0 existe un δ > 0 tal que para toda familia nita de intervalos disjuntos [ak, bk]⊂[a, b], k = 1,2, ..., n , que verique Pn k=1 (bk−ak)< δ se cumple que Pn k=1 |f(bk)−f(ak)|<  . Observación 17 . El espacio AC [a, b] coincide con el espacio de primitivas de funciones Lebesgue integrables: f∈AC [a, b]⇔f(x) = k+Zx a ϕ(t)dt (ϕ∈L(a, b)) y en consecuencia una función f absolutamente continua tiene derivada f0(x) = ϕ(x) en casi todo punto. Además, se deduce trivialmente que k=f(a) , con lo que f(x) = f(a) + Zx a f0(t)dt Observación 18 . Ser absolutamente continua es una propiedad bastante fuerte: implica continuidad uniforme (y por tanto continuidad). Sin embargo es menos restrictiva que ser Lipchiziana: toda función de Lipschitz es absolutamente continua. Observación 19 . Denotamos por ACn[a, b] (n∈N) al espacio de funciones f con derivadas continuas hasta el orden n−1 en [a, b] tales que f(n−1) ∈AC [a, b] . Es claro que AC1[a, b] = AC [a, b] y que el espacio ACn[a, b] consiste en todas las funciones representables por una integral de Lebesgue n -múltiple sumada a un polinomio de orden n−1 . Teorema 20. Sea f∈ACn[a, b], n ∈N y α∈C vericando <(α)≥0 y [<(α)]+1 ≤ n . Entonces Dα af existe en casi todo punto de [a, b] . Demostración. Ver en Samko, Kilbas y Marichev [1]. Corolario 21. Si f∈AC [a, b] , entonces Diθ af existe en casi todo punto de [a, b] para cualquier θ∈R\{0} y se representa en la forma (2.2.5) . Enumeramos ahora las propiedades más relevantes de la derivada fraccionaria de Riemann-Liouville (en todas ellas suponemos que existen las derivadas fraccionarias implicadas). La derivada fraccionaria de Riemann-Liouville carece de algunas cualidades importantes, como la de semigrupo (ver Proposición 9), que sí cumple la integral. 17 La relación Dα aDβ af=Dα+β af sólo es válida en casos muy especícos. Otra importante diferencia con el caso entero es que, en general, la derivada de una constante no nula no vale cero. Las demostraciones pueden encontrarse en Samko, Kilbas y Marichev [1] y Kilbas, Srivastava y Trujillo [2]. Proposición 22 ( linealidad ) . Sean f, g ∈L1(a, b) (−∞ <a<b<∞) y α∈C (<(α)≥0) . Entonces Dα a[λf +µg] = λDα af+µDα ag∀λ, µ ∈R Proposición 23. Sea f∈L1(a, b) (−∞ <a<b<∞) . Entonces, tomando valores de α∈C tales que <(α)≥0 , l´ım α→0+(Dα af) (x) = f(x) en casi todo punto de [a,b]. El siguiente resultado reeja que, análogamente a lo que ocurre en el caso entero, la derivada fraccionaria de Riemann-Liouville es el operador inverso izquierdo de la integral fraccionaria de Riemann-Liouville. Proposición 24. Sea f∈L1(a, b) (−∞ <a<b<∞) y α, β ∈C(<(α),<(β)>0) . Entonces las igualdades (Dα aIα af) (x) = f(x) (2.2.6) Dα aIβ af(x) = Iβ−α af(x), si α ≤β (2.2.7) Dα aIβ af(x) = Dα−β af(x), si α ≥β (2.2.8) son ciertas en casi todo punto de [a,b]. Sin embargo, también como en el caso entero, no es cierto el recíproco: la integral no es la inversa por la izquierda de la derivada. En consecuencia no se puede hablar estrictamente de operadores inversos. Proposición 25. Sea f∈L1(a, b) (−∞ <a<b<∞) y α∈C(<(α)>0) . Entonces si fn−α=In−α af∈ACn[a, b] (Iα aDα af) (x) = f(x)− n X k=1 (x−a)α−k Γ (α−k+ 1)f(n−k) n−α(a), (2.2.9) donde n= [<(α)] + 1 , es cierto en casi todo punto de [a,b]. 18 Corolario 26. En las condiciones de la Proposición 25, se verica que (Iα aDα af) (x) = f(x) + c0xα−1+c1xα−2+···+cn−1xα−n (2.2.10) con ci∈C∀i= 1,2, ..., n −1 . Ya adelantábamos el hecho sorprendente de que la derivada fraccionaria no es nula para las funciones constantes. Esto se concluye del valor de la derivada de una potencia: Proposición 27. Sea [a, b]⊂R(−∞ <a<b<∞) , α∈C(<(α)>0) y β > 0 . Se verica que Dα a(x−a)β=Γ (β+ 1) Γ (β+ 1 −α)(x−a)β−α (2.2.11) Figura 2.2.1: Evolución de las derivadas fraccionarias Dα 0x2 , con α variando entre 0 (parábola) y 2 (función constante) 19 Observación 28 . Como caso particular de (2.2.11), se tiene que cuando β=α−1, α− 2, ..., n (n= [<(α)] + 1) la derivada de orden α -ésimo se anula, por lo que Dα a anula las potencias (x−a)α−1,(x−a)α−2, ..., (x−a)α−n . Observación 29 . Para el cálculo de la derivada de una función constante f(x) = k hacemos β= 0 e introducimos la constante multiplicativa k : Dα ak=kΓ (1) Γ (1 −α)(x−a)−α=k Γ (1 −α) (x−a)α 2.2.3. Derivada fraccionaria de Caputo La derivada fraccionaria de Riemann-Liouville jugó un papel determinante en el desarrollo del cuerpo teórico del Cálculo Fraccionario, y se utilizó con éxito en aplicaciones estrictamente matemáticas. Pero al tratar de realizar modelizaciones matemáticas de fenómenos físicos reales por medio de ecuaciones diferenciales fraccionarias, surgió el problema de que las condiciones iniciales eran también de orden fraccional. Este tipo de condiciones no son físicamente interpretables y presentan un obstáculo considerable a la hora de hacer uso práctico del Cálculo Fraccionario. El operador diferencial de Caputo, en contraste con el de Riemann-Liouville, emplea como condiciones iniciales derivadas de orden entero, es decir, valores iniciales que son físicamente interpretables a la manera tradicional. La denición que sigue representó pues un notable avance práctico en el estudio de fenómenos físicos como los de tipo viscoelástico y otros. Denición 30. Sea α∈C(<(α)≥0) y f∈ACn[a, b] (−∞ <a<b<∞) , con n= [<(α)] + 1 . La derivada fraccionaria de Caputo de orden α de f se dene por CDα af(x) : = In−α ad dxn f(x) =1 Γ (n−α)Zx a f(n)(t) (x−t)α−n+1 dt (x>a) (2.2.12) donde f(n) es la n -ésima derivada usual de f . Observación 31 . Al contrario que en la derivada fraccionaria de Riemann-Liouville, en la que primero se integra y luego se deriva, en la derivada de Caputo primero derivamos (n veces) y seguidamente integramos. En consecuencia, se trata de una denición más restrictiva, ya que requiere la integrabilidad de f(n) . A pesar de ello, la hipótesis del Teorema 20 (que f pertenezca a ACn[a, b] ) sigue siendo una condición suciente para garantizar la existencia de CDα af con cualquier orden α tal que [<(α)] + 1 ≤n . 20 Las derivadas fraccionarias de Caputo y de Riemann-Liouville son buenas generalizaciones de la derivada ordinaria, en el sentido de que respetan los valores de las derivadas enteras usuales, concordando así entre ellas. Pero en el caso no entero no coinciden, como pone de maniesto el siguiente resultado, cuya demostración puede encontrarse en Kilbas, Srivastava y Trujillo [2]. Proposición 32. Sea α /∈N(<(α)≥0) , n = [<(α)] + 1 y f∈L1[a, b] una función para la que existen las derivadas fraccionarias de Caputo CDα af y Riemann-Liouville (Dα af) . Entonces se verica la siguiente relación: CDα af(x) = (Dα af) (x)− n−1 X k=0 f(k)(a) Γ (k+ 1 −α)(x−a)k−α(x > a) (2.2.13) Observación 33 . Por lo tanto, para órdenes de derivación no enteros, las derivadas de Caputo y Riemann-Liouville coincidirán cuando se cumpla f(a) = f0(a) = ··· =f(n−1) (a) = 0 (2.2.14) Observación 34 . Como ejemplo de la no coincidencia de estos dos tipos de derivada fraccionaria, observemos las grácas de D0,5 0cos(x) y CD0,5 0cos(x) : 21 La divergencia era de esperar ya que no se cumple la condición (2.2.14): cos (0) 6= 0. La integral de Riemann-Liouville tampoco es la inversa por la izquierda de la derivada de Caputo. Pero en el polinomio resto obtenido al menos aparecen derivadas enteras, a diferencia de (2.2.9). Esta característica nos será de utilidad a la hora de resolver analíticamente algunas Ecuaciones Diferenciales Fraccionarias: Proposición 35. Sea f∈L1(a, b) (−∞ <a<b<∞) y α∈C(<(α)>0) . Entonces si f∈ACn[a, b] o f∈Cn[a, b] IαC aDα af(x) = f(x)− n−1 X k=0 (x−a)k k!f(k)(a), (2.2.15) donde n= [<(α)] + 1 para α /∈N y n=α para α∈N . Corolario 36. En las condiciones de la Proposición 35, se verica que IαC aDα af(x) = f(x) + c0+c1x+c2x2+···+cn−1xn−1 (2.2.16) con ci∈C∀i= 1,2, ..., n −1 . La derivada de Caputo puede considerarse por tanto como una regularización de la derivada de Riemann-Liouville efectuada mediante la sustracción de un polinomio de Taylor, a través de la cual se obtienen ciertas ventajas importantes, como las citadas condiciones iniciales enteras, unos requisitos menos restrictivos para el cumplimiento de la propiedad de semigrupo CDαC aDβ af=CDα+β af , y la no menos signicativa característica que enunciamos a continuación. Proposición 37. Sea α∈C(<(α)>0) y sea n= [<(α)] + 1 . Se verica que CDα a(x−a)k= 0 ∀k= 0,1,··· , n −1 . Esta propiedad, demostrada en Kilbas, Srivastava y Trujillo [2], evidencia que el valor de la derivada fraccionaria de Caputo operada en una constante es nulo. 2.2.4. Derivada fraccionaria de Grünwald-Letnikov Al contrario que los operadores de Riemann-Liouville o Caputo, que obtienen sus deniciones de una integral repetida, el enfoque de Grünwald-Letnikov toma como 22 punto de partida la derivada, en concreto de la siguiente fórmula para la derivada n -ésima, fácilmente deducible: d dxn f(x) = l´ım h→0 1 hn n X k=0 (−1)kn kf(x−kh). (2.2.17) Esta expresión puede ser generalizada para valores no enteros de n , teniendo en cuenta que para α > 0 podemos interpretar α k como α k=Γ (α+ 1) k!Γ (α+ 1 −k). (2.2.18) Entonces, si consideramos una serie innita se obtiene d dxα f(x)≈l´ım h→0+ 1 hα ∞ X k=0 (−1)kα kf(x−kh). Por otra parte, dado un h determinado, no tienen sentido términos de la serie para valores de k superiores a x−a h , ya que en ese caso f(x) tomaría valores fuera del intervalo de denición [a, b] . Se llega entonces a la denición que sigue. Denición 38. Sea f una función acotada en [a, b] (−∞ <a<b<∞) . La derivada fraccionaria de Grünwald-Letnikov de orden α∈R(α > 0) de f se dene, si existe, por fα) a(x) := l´ım h→0+ 1 hα [x−a h] X k=0 (−1)kα kf(x−kh) (x > a) (2.2.19) donde α k signica (2.2.18). Observación 39 . La derivada Grünwald-Letnikov sólo está denida para órdenes de derivación reales. Observación 40 . En Podlubny [3] se demuestra que los operadores diferenciales de Riemann-Liouville y Grünwald-Letnikov dan los mismos resultados para idénticos límites inferiores a y órdenes de integración reales α > 0 . Cabe preguntarse cuál utilizar ante un problema concreto. En la literatura relacionada con la resolución de ecuaciones diferenciales fraccionarias se utiliza preferentemente la denición de Riemann-Liouville para la formulación de los problemas, y luego, a la hora de obtener la solución numérica, se pasa a la denición de Grünwald-Letnikov, muy apropiada para ser implementada en cálculos numéricos de tipo iterativo. 23 2.2.5. Otros operadores fraccionarios Todos los operadores fraccionarios que hemos tratado se denían sobre un intervalo acotado de la recta real [a, b] . Las formas integrales de Riemann-Liouville y diferenciales de Caputo, con esta característica, serán las que utilizaremos en los cálculos a lo largo de este trabajo. Pero existen otros tipos de operadores, denidos sobre intervalos no acotados de la recta real, o incluso sobre el plano complejo. Otros están diseñados para funciones de varias variables. El operador integral de Weil (también llamado de Liouville) es una extensión de la integral fraccionaria de Riemann-Liouville al semieje no acotado (−∞, x] : Iα +f(x) := 1 Γ (α)Zx −∞ f(t) (x−t)1−αdt (−∞ <x<∞). Esta forma, y su correspondiente forma diferencial, son especialmente adecuadas para realizar integración y derivación fraccionaria de funciones periódicas. Pero hay una gran variedad de operadores fraccionarios. Otros registrados en la literatura especializada son el de Riesz, para derivadas e integrales fraccionarias de funciones de varias variables, Erdélyi-Kober, Hadamard, Bessel, Chen o Dzherbashyan. 2.3. Ecuaciones diferenciales fraccionarias De la misma forma que en el Cálculo ordinario, las ideas de derivación e integración fraccionarias conducen al concepto más avanzado de ecuación diferencial. Una relación involucrando uno o más operadores fraccionarios aplicados a una función desconocida f se conoce como ecuación diferencial fraccionaria 2 (EDF, de ahora en adelante). Gran parte de la teoría de ecuaciones diferenciales ordinarias tiene una correspondencia fraccionaria más o menos natural. Por ejemplo, es conocido que una EDF de orden α necesita [<(α)] + 1 condiciones iniciales para ser resuelta de manera única. Sin embargo surgen también notables diferencias, la primera de las cuales es el hecho de que una misma EDF tendrá diferentes signicados según los operadores fraccionarios implicados sean los de Riemann-Liouville, Caputo u otros. 2 En la obra pionera de Oldham y Spanier [5] aparecen con el curioso nombre de ecuaciones diferenciales extraordinarias , en oposición a las EDOs. 24 2.3.1. Deniciones y conceptos Análogamente a la teoría clásica de ecuaciones diferenciales, las EDF se dividen en ecuaciones lineales, homogéneas, no homogéneas, con coecientes constantes o variables. Empezamos dando la denición más general para luego ir concretando casos más especícos. Denición 41. Una ecuación diferencial fraccionaria de orden α∈C(<(α)>0) es una relación del tipo (Dα ay) (x) = f[x, y (x),(Dα1 ay) (x),(Dα2 ay) (x), ..., (Dαm−1 ay) (x)] (2.3.1) donde y(x) es una función compleja desconocida de dominio real, f[x, y, y1, y2, ..., ym−1] es una función conocida y Dαk a(k= 1,2, ..., m −1) son operadores diferenciales fraccionarios vericando 0<<(α1)<<(α2)<··· <<(αm−1)<<(α) y m≥2 . Observación 42 . Establecemos varias condiciones sobre los operadores Dαk a de la de- nición anterior. Han de ser operadores de tipo diferencial (en otro caso hablaríamos de ecuaciones integrales o íntegrodiferenciales). Además, se entiende que representan derivadas fraccionarias del mismo tipo (Riemann-Liouville, Caputo, etc.) y que tienen un límite inferior de integración a común. Denición 43. Se denomina solución de la EDF anterior a cualquier función compleja de variable real y(x) que verique la igualdad (2.3.1). Denición 44. Una ecuación diferencial fraccionaria lineal de orden α∈C(<(α)>0) es una relación del tipo (Dα ay) (x) = a0(x)y(x) + m−1 X k=1 ak(x) (Dαk ay) (x) + b(x) (2.3.2) donde y(x) es una función compleja desconocida de dominio real, ak(x) (k= 1,2, ..., m −1) y b(x) son funciones complejas conocidas y Dαk a(k= 1,2, ..., m −1) son operadores diferenciales fraccionarios vericando 0<<(α1)<<(α2)<··· <<(αm−1)<<(α) y m≥2 . Cuando se cumple que las funciones ak(x) (k= 1,2, ..., m −1) y b(x) son constantes se dice que la EDF lineal es de coecientes constantes , y si se da que b(x)=0 se denomina EDF lineal homogénea . Problemas análogos a los de Cauchy y Dirichlet para ecuaciones diferenciales surgen también en el Cálculo Fraccionario. Aquí hemos de diferenciar explícitamente 25 y0(0) = b1 (2.3.16) Si aplicamos el operador integral Iα 0 sobre ambos miembros de la ecuación (2.3.14) obtenemos Iα 0CDα 0y(x) = (Iα 0f) (x) lo que, en virtud de (2.2.16), resulta en que la solución general de la EDF es y(x) = c0+c1x+ (Iα 0f) (x) (2.3.17) siendo c1, c2 constantes complejas dependientes de las condiciones iniciales. Si introducimos la condición (2.3.15) en la solución (2.3.17) y, dado que (Iα 0f) (0) = 0 , llegamos a que c0=b0 . Vamos a hacer uso de la segunda condición (2.3.16). Derivando (2.3.17) se obtiene y0(x) = c1+D1 0Iα 0f(x) que en base a (2.2.7) vale y0(x) = c1+Iα−1 0f(x) por lo que como Iα−1 0f(0) = 0 , se llega a c1=b1 . Por tanto, la solución particular del problema (2.3.14)-(2.3.15)-(2.3.16) será: y(x) = b0+b1x+ (Iα 0f) (x) =b0+b1x+1 Γ (α)Zx 0 (x−t)α−1f(t)dt (2.3.18) que como era de esperar haciendo α= 2 coincide con la conocida solución y(x) = b0+b1x+Zx 0 (x−t)f(t)dt de la EDO y00 (x) = f(x) con las mismas condiciones iniciales (2.3.15)-(2.3.16). 32 Método 2 (transformada de Laplace) Vamos a tratar ahora una ecuación algo más complicada, que será una generalización de la conocida ecuación diferencial de oscilación y00 (x) = λy (x) + f(x) . Sustituyendo la segunda derivada por una derivada de Caputo de orden α∈(1,2] obtenemos una EDF que denominaremos ecuación de oscilación fraccionaria . Para resolver el problema usaremos el método de la transformada de Laplace. Este sistema tiene la ventaja de que haciendo uso de la transformada podemos introducir las condiciones iniciales desde el primer momento, llegando directamente a la solución particular. Dada una función continua f: [0,∞)→C , λ∈C y α∈(1,2] vamos a resolver el problema tipo Cauchy con derivadas fraccionarias de Caputo dado por CDα 0y(x) = λy (x) + f(x) (x > 0) (2.3.19) Como 1< α ≤2 se tiene que n= [α] + 1 = 2 , por lo que necesitamos dos condiciones iniciales en el límite inferior de integración 0: y(0) = b0 (2.3.20) y0(0) = b1. (2.3.21) Aplicando la transformada de Laplace en ambos miembros de la ecuación 2.3.19, para lo cual utilizamos el resultado (2.3.8) para el cálculo de la transformada de una derivada de Caputo, se llega a sαY(s)−sα−1y(0) −sα−2y0(0) = λY (s) + F(s) donde Y(s) = L{y(x)} y F(s) = L{f(x)} . Introduciendo las condiciones iniciales (2.3.20)-(2.3.21) y despejando Y(s) obtenemos Y(s) = b0 sα−1 sα−λ+b1 sα−2 sα−λ+F(s) sα−λ. Se trata ahora de aplicar la transformada inversa de Laplace a los dos miembros de la ecuación, para recuperar y(x) . Utilizando la linealidad de L−1 , y(x) = b0L−1sα−1 sα−λ+b1L−1sα−2 sα−λ+L−1F(s) sα−λ. (2.3.22) 33 Calculemos por separado las tres transformadas inversas. El resultado (2.3.10) proporciona el valor de las dos primeras de forma directa: L−1sα−1 sα−λ=Eα,1(λxα) (2.3.23) L−1sα−2 sα−λ=xEα,2(λxα) (2.3.24) Para el cálculo de L−11 sα−λF(s) utilizamos la propiedad (2.3.12), L−11 sα−λF(s)=L−11 sα−λ∗L−1(F(s)) (2.3.25) con lo que reducimos el problema a la obtención de L−1(1 /sα−λ) . De nuevo volvemos a recurrir a (2.3.10), estableciendo los valores k= 0, β =α y a =λ , con lo que resulta L−11 sα−λ=xα−1Eα,α (λxα) , que junto con (2.3.25) proporciona L−11 sα−λF(s)=xα−1Eα,α (λxα)∗f(x) =Zx 0 (x−t)α−1Eα,α [λ(x−t)α]f(t)dt (2.3.26) En conclusión, (2.3.22), (2.3.23), (2.3.24) y (2.3.26) facilitan la solución particular al problema de Cauchy de partida: y(x) = b0Eα,1(λxα)+b1xEα,2(λxα)+Zx 0 (x−t)α−1Eα,α [λ(x−t)α]f(t)dt (2.3.27) Observación 65 . Las funciones Eα,1(λxα) , xEα,2(λxα) forman, por tanto, un sistema fundamental de soluciones de la correspondiente EDF homogénea CDα 0y(x)−λy (x) = 0 con α∈(1,2] . Observación 66 . El problema (2.3.14)-(2.3.15)-(2.3.16) estudiado con el método anterior es un caso particular, más sencillo, de este. Basta tomar λ= 0 en la ecuación 34 (2.3.19). Veamos que la solución obtenida aquí es coherente con la anterior. Haciendo λ= 0 en la solución (2.3.27) tenemos y(x) = b0Eα,1(0) + b1xEα,2(0) + Zx 0 (x−t)α−1Eα,α (0) f(t)dt y teniendo en cuenta que las funciones de Mittag-Leer verican Eα,1(0) = Eα,2(0) = 1 , Eα,α (0) = 1 /Γ(α) , llegamos a la solución (2.3.18) esperada: y(x) = b0+b1x+1 Γ (α)Zx 0 (x−t)α−1f(t)dt Observación 67 . Goreno y Mainardi [8] han obtenido una solución alternativa equivalente para este problema (con λ=−1 ). En su trabajo, la solución se expresa en términos de derivadas e integrales de funciones de Mittag-Leer de un solo parámetro Eα(z) : y(x) = b0Eα(−xα) + b1I1 0Eα(−xα) + Zx 0 d dtEα(−tα)f(x−t)dt 35 36 Capítulo 3 Aplicaciones dinámicas 3.1. Péndulo simple El péndulo simple es un sistema físico constituido por una partícula de masa m que, suspendida de un punto jo O por medio de una varilla de longitud L , puede oscilar en un plano vertical jo por efecto de la fuerza de la gravedad. La posición de la partícula en el instante t se especica mediante el ángulo θ que la varilla forma con la vertical en ese momento. El estudio del péndulo constituye uno de los problemas clásicos de la dinámica elemental, y todas sus componentes están perfectamente determinadas desde el punto de vista matemático, y son bien conocidas. Lo que se pretende con este apartado es realizar un análisis del movimiento pendular enfocado desde la óptica más extensa del Cálculo Fraccionario, y para ello se va a generalizar la idea clásica de aceleración (derivada segunda de la posición) a una derivada fraccionaria de Caputo de un orden comprendido entre 1.5 y 2. Obviamente el péndulo simple es un sistema idealizado, al que para su análisis vamos a exigir una serie de hipótesis simplicativas: La varilla que sujeta a la partícula carece de masa, es inextensible y siempre permanece rígida. El movimiento de la partícula se traza en dos dimensiones; es decir, la partícula no traza una elipse sino un arco. El sistema no pierde energía por efecto de la resistencia del aire ni por fricción alguna. 37 Figura 3.1.1: Fuerzas que intervienen en el péndulo 3.1.1. Ecuación diferencial fraccionaria Para determinar la función del movimiento del péndulo fraccionario empezamos buscando la ecuación diferencial fraccionaria que modeliza dicho movimiento. La ecuación buscada tiene por incógnita la función θ(t) , que determina en radianes el ángulo de la varilla en el instante t . La partícula se mueve sobre un arco de circunferencia bajo el efecto de dos fuerzas: su propio peso ( mg , donde g= 9,80665 m/ s2 es la fuerza gravitatoria ejercida por la Tierra cerca de su supercie) y la fuerza de tensión T ejercida por la varilla. Si descomponemos el peso en sus componentes tangencial y normal, observamos que la componente normal de esta fuerza se ve contrarrestada por la fuerza de tensión de la varilla (ver Figura 3.1.1). Por lo tanto, la única fuerza actuante en lo que concierne al movimiento del sistema será la componente tangencial del peso, que tendrá signo negativo por ir siempre en dirección opuesta al movimiento: F(t) = −mg sin (θ(t)) (3.1.1) En el modelo clásico, se introduce esta fuerza en la segunda Ley de Newton, F=ma, (3.1.2) siendo a la aceleración tangencial del movimiento, y se llega así a la ecuación tra38 dicional del péndulo. Aquí vamos a sustituir (3.1.2), en la que la aceleración a(t) es la derivada segunda de la posición r(t) , por una fórmula alternativa más general haciendo uso de derivadas fraccionarias: F(t) = mCDα 0r(t), (3.1.3) con 1,5< α ≤2 . Nos interesa una ecuación en la variable ángulo de la varilla θ(t) , y no en la variable posición de la partícula r(t) . Por ser un movimiento a lo largo de un arco de circunferencia de radio L tenemos que r(t) = Lθ (t) , con lo que CDα 0r(t) = LCDα 0θ(t) (3.1.4) donde CDα 0θ(t) será la aceleración angular. Finalmente, de (3.1.1), (3.1.3) y (3.1.4) obtenemos −mg sin (θ(t)) = mL CDα 0θ(t) que nos lleva a la ecuación diferencial del péndulo en su generalización fraccionaria: CDα 0θ(t) + g Lsin (θ(t)) = 0 (3.1.5) siendo 1,5< α ≤2 , que se reduce a la ecuación clásica d2θ dt2(t) + g Lsin (θ(t)) = 0 haciendo α= 2 . 3.1.2. Solución La EDF (3.1.5) no tiene fácil solución, al no tratarse de una ecuación lineal debido al término no lineal sin (θ(t)) . Ahora, recordemos que sin (θ)≈θ cuando θ es pequeño; en particular, θ y sin (θ) coinciden en las dos primeras cifras decimales cuando θ < π/12 (ver Cuadro 3.1). Por lo tanto, si nos restringimos a oscilaciones de pequeña amplitud (como las que describe, por ejemplo, el péndulo de un reloj), parece razonable simplicar nuestro modelo matemático sustituyendo sin (θ) por θ en la ecuación (3.1.5) 1 . La ecuación 1 Esto es lo que se hace en el estudio de la dinámica del péndulo clásico , y es razonable para valores pequeños del ángulo debido a que la solución constante cero es estable. En el caso fraccionario, la estabilidad de la solución cero, aunque parece obvia, es todavía un problema abierto (ver [11]). 39 resultante se reduce a CDα 0θ(t) + g Lθ(t) = 0 (3.1.6) Cuadro 3.1: Comparación de θ y sin (θ) θ ( º ) θ (rad) sin (θ) 0 0,00000 0,00000 2 0,03491 0,03490 5 0,08727 0,08716 10 0,17453 0,17365 15 0,26180 0,25882 20 0,34907 0,34202 25 0,43633 0,42262 30 0,52360 0,50000 Esta EDF encaja dentro del tipo de las estudiadas en el Método 2 en la página 33. Si establecemos las condiciones iniciales del péndulo, es decir, ángulo de partida θ0 y velocidad angular inicial v0 , θ(0) = θ0 (3.1.7) θ0(0) = v0 (3.1.8) entonces utilizando la fórmula (2.3.27), y dado que en este caso f= 0 , la solución particular, que determina la posición del péndulo en función del tiempo, es θ(t) = θ0Eα,1−g Ltα+v0tEα,2−g Ltα. (3.1.9) Por otra parte, en este sistema es habitual considerar nula la velocidad inicial de la partícula. Esto es, en origen el péndulo se deja oscilar sin ejercer sobre él ningún impulso. En esta situación la segunda condición inicial es θ0(0) = 0 (3.1.10) y la solución queda simplicada a θ(t) = θ0Eα,1−g Ltα=θ0Eα−g Ltα. (3.1.11) 40 En el caso particular α= 2 , teniendo en cuenta que E2,1−g Lt2= ∞ X k=0 −g Lt2k Γ (2k+ 1) = ∞ X k=0 (−1)k (2k)! rg Lt2k = cos rg Lt se obtiene la solución tradicional del péndulo linealizado d2θ dt2(t) + g Lθ(t)=0 con condiciones iniciales (3.1.7)-(3.1.10): θ(t) = θ0cos rg Lt 3.2. Proyectil Otro problema dinámico clásico es el del proyectil . Se trata de un sistema en el que una partícula es proyectada formando un determinado ángulo con la supercie de la Tierra. El proyectil, bajo el efecto de la fuerza de la gravedad, describirá una curva de tipo parabólico. Esta trayectoria está perfectamente determinada matemáticamente incluso para el caso, un poco más complicado, en el que consideremos otras fuerzas intervinientes en el sistema además de la gravedad, como son las de rozamiento. Áreas como la balística han estudiado en profundidad este tema. Aquí, siguiendo trabajos como los de Ebaid [15] u Otero [19], vamos a estudiar el movimiento de un proyectil desde el enfoque fraccionario. Para ello, tal y como hicimos con el péndulo, utilizaremos como hipótesis inicial una fórmula generalizada de la segunda Ley de Newton. Vamos a partir de unas hipótesis que simplican el problema, de forma que la única fuerza actuante sea la gravedad: El medio no ofrece oposición al avance del proyectil ni por resistencia del aire ni por ninguna otra fricción. No hay curvatura de la supercie terrestre, que es plana y sin rugosidades. El movimiento se traza dentro de un plano vertical jo OXY . La fuerza de la gravedad es uniforme, y no disminuye con la altura del proyectil. No se tiene en cuenta la fuerza de Colioris, debida al movimiento de rotación de la Tierra. 41 La constante de proporcionalidad c > 0 se denomina constante de amortiguamiento y, en el caso del aire, dependerá de las condiciones atmosféricas. Podría pensarse que no estamos teniendo en cuenta el peso mg del cuerpo como fuerza interviniente en el sistema. Por el contrario, el peso se ve contrarrestado permanentemente por una fuerza de igual módulo y signo contrario: la tensión ejercida por el muelle, que, jado en un punto, sujeta en todo momento el cuerpo. Por lo tanto ambas fuerzas se compensan mutuamente, y la fuerza neta resultante es nula. Tenemos así que la masa soporta la fuerza F=FH+FR . Para enfocar el problema del resorte desde el punto de vista fraccionario sustituimos de nuevo la fórmula de la segunda ley de Newton F=ma por la generalización fraccionaria F(t) = mCDα 0x(t), (3.3.3) con 1,5< α ≤2 . De esta manera obtenemos la nueva EDF que gobierna la dinámica del muelle: mCDα 0x(t) = FH(t) + FR(t)⇒ mCDα 0x(t) + cx0(t) + kx (t)=0 (3.3.4) a la que podemos añadir unas condiciones iniciales que determinen la posición inicial de la masa, x(0) = b0, (3.3.5) y su velocidad inicial, o el impulso que se le transmite en el inicio del experimento, x0(0) = b1. (3.3.6) 3.3.2. Soluciones El problema de Cauchy dado por (3.3.4)-(3.3.5)-(3.3.6) no tiene fácil solución, y en la bibliografía existente no se resuelve completamente para estos parámetros particulares. Kilbas, Srivastava y Trujillo [2] han encontrado dos soluciones particulares de la ecuación (3.3.4): x1(t) = ∞ X n=0 (−k/m)n n!tαn ∞ X j=0 Γ (n+j+ 1) Γ (αn +1+αj −j) (−c /m)jt(α−1)j j!+ +c m ∞ X n=0 (−k/m)n n!tαn+α−1 ∞ X j=0 Γ (n+j+ 1) Γ (αn +α+αj −j) (−c /m)jt(α−1)j j! (3.3.7) 48 x2(t) = ∞ X n=0 (−k/m)n n!tαn+1 ∞ X j=0 Γ (n+j+ 1) Γ (αn +2+αj −j) (−c /m)jt(α−1)j j!+ +c m ∞ X n=0 (−k/m)n n!tαn+α ∞ X j=0 Γ (n+j+ 1) Γ (αn +1+α+αj −j) (−c /m)jt(α−1)j j! (3.3.8) No hay garantía de que estas dos soluciones formen un sistema fundamental de soluciones de la ecuación. El encontrarlas se convierte en una cuestión bastante complicada. Por ejemplo, si α es un racional p /q con p < q , necesitaríamos p soluciones para congurar un sistema fundamental. Para un α real no racional, el proceso se complica más. 49 50 Bibliografía [1] SAMKO, S. G., KILBAS, A. A., y MARICHEV, O. I., Fractional Integrals and Derivatives: Theory and Applications , Gordon and Breach Science Publishers, Suiza, 1993. [2] KILBAS, A. A., SRIVASTAVA, H. M. y TRUJILLO, J. J., Theory and Applications of Fractional Dierential Equations , Elsevier, Amsterdam, 2006. [3] PODLUBNY, I., Fractional Dierential Equations , Academic Press, San Diego, 1999. [4] MILLER, K. S. y ROSS, B., An Introduction to the Fractional Calculus and Fractional Dierential Equations , Wiley and Sons, New York, 1993. [5] OLDHAM, K. B. y SPANIER, J., The Fractional Calculus , Academic Press, New York, 1974. [6] PODLUBNY, I., Geometric and physical interpretation of fractional integration and fractional dierentiation, Fract. Cal. Appl. Anal ., 5(4), (2002) 367386. [7] HEYMANS, N. y PODLUBNY, I., Physical interpretation of initial conditions for fractional dierential equations with Riemann-Liouville fractional derivatives, Rheologica Acta , vol. 45(5), (2006) 765-771 [8] GORENFLO, R. y MAINARDI, F., Fractional Calculus: Integral and dierential equations of fractional order, CISM courses and lectures, vol. 378, (1997) 223276. [9] VÁZQUEZ MARTÍNEZ, L., Una panorámica del cálculo fraccionario y sus aplicaciones, Rev. R. Acad. Cienc. Exact. Fís. Nat. , vol. 98(1), (2004) 17-25 [10] TENREIRO MACHADO, J., KIRYAKOVA, V. y MAINARDI, F., Recent history of fractional calculus, Commun Nonlinear Sci Numer Simulat , vol. 16(3), (2011) 11401153 51 [11] CHEN, F., NIETO, J.J. y ZHOU, Y., Global attractivity for nonlinear fractional dierential equations, Nonlinear Analysis: Real World Applications, 13, (2012) 287-298 [12] KILBAS, A., Some aspects of dierential equations of fractional order, Rev. R. Acad. Cienc. Exact. Fís. Nat. , vol. 98(1), (2004) 27-38 [13] AGRAWAL, O. P., A new lagrangian and a new Lagrange Equation of Motion for fractionally damped systems, Journal of Applied Mechanics , vol. 68, (2001) 339-341 [14] HAUBOLD, H. J., MATHAI, A. M. y SAXENA, R. K., Mittag-Leer functions and their applications, Journal of Applied Mathematics , vol. 2011, (2009) [15] EBAID, A., Analysis of projectile motion in view of fractional calculus, Applied Mathematical Modelling , 35, (2011) 1231-1239 [16] SHIMA, H., How far can Tarzan jump?, European Journal of Physics , vol. 33(6) (2012) 1687-1696 [17] BALEANU, D., ASAD, J. H. y PETRAS, I., Fractional-order two-electric pendulum, Romanian Reports in Physics , vol. 64(4), (2012) 907-914 [18] PIERANTOZZI, T., Estudio de generalizaciones fraccionarias de las ecuaciones estándar de difusión y de onda s (Tesis doctoral). Tutor: Vázquez, L., Universidad Complutense de Madrid, Facultad de Matemáticas (2006) [19] OTERO, O., Ecuaciones diferenciales de orden fraccionario (Trabajo de Fin de Máster). Tutor: Juan José Nieto Roig, Universidad de Santiago de Compostela, Facultad de Matemáticas (2012) [20] VELASCO CEBRIÁN, M. P., Modelos diferenciales y funciones especiales en el ámbito del cálculo fraccionario (Trabajo de Fin de Máster). Tutores: Trujillo, J. J. (Universidad de La Laguna) y Vázquez, L. (Universidad Complutense de Madrid) (2008) 52