Spanish (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Cited by Google 
Similars in
SciELO 
Similars in Google 
Permalink
1. Introducción
El cálculo fraccionario surge para introducir derivadas de orden no entero y permite unificar los conceptos de integral y derivada. Ya desde los orígenes del cálculo clásico se especuló acerca del significado de una derivada de orden 1/2, de tal manera que D1/2D1/2 f = Df, como se evidencia en una carta de Gottfried Wilhelm Leibniz a Guillaume de L'Höpital en 1695 [1]. Sin embargo, su desarrollo se dió en los siglos XIX y XX, principalmente. Fue Niels Henrik Abel quien da un primer desarrollo de las nociones fundamentales de la integral y la derivada de orden no entero y su unificación [2, 3], trabajo análogo que desarrolló paralelamente Joseph Liouville [4, 5]. A lo largo de la historia se han dado múltiples definiciones para la integral y derivada de orden no entero, lo cual ha obstaculizado, en cierta medida la asimilación de estas técnicas en las diferentes áreas de las ciencias. Sin embargo, en la literatura moderna [6, 7] la mayoría de las definiciones suelen fundamentarse en la integral de Riemann-Liouville, y aunque hay diversas definiciones para la derivada [8, 9, 10] como: la derivada de Riemann-Liouville, de Caputo, de Grüwald-Letnikov, y de Riesz; la mayoría de las aplicaciones que usan ecuaciones diferenciales en derivadas fraccionarias lo hacen aplicando la derivada de Caputo pues las condiciones iniciales replican las que se tienen para las ecuaciones diferenciales en derivadas ordinarias. Es importante destacar, que además del interés teórico natural que recibe el cálculo fraccionario, se ha logrado aplicarlo en áreas tan diversas como: el estudio de la viscoelasticidad [11], la difusión anómala en sistemas biológicos [12], propagación ondulatoria [13], imágenes médicas por ultrasonido y elastografía [14], la espectroscopía de impedancia eléctrica [15], la teoría de control [16], el tratamiento de imágenes [17], entre muchos otros [18, 19]. Sin embargo, los cálculos para resolver modelos sencillos no son simples y con este fin se disponen de algunas ayudas computacionales. En el software MatLab se dispone de paquetes como [20] FOTF, NINTEGER, CRONE, FOMCON, DFOD, y los mencionados en [21], entre otros. En Fortran se tiene el código [22] denominado AstroFracTool. En cuanto a Python se dispone del paquete [23] FOMCOpy (extensión de FOMCO en Matlab) aplicado en el Internet de las Cosas, fracdiff para derivadas [24] y su extensión a integrales [25], como también differint [26] para el cálculo de derivadas. Sin embargo, no conocemos de un paquete más completo en Python que, además de calcular derivadas e integrales de orden no entero, pueda desarrollar la solución numérica de una ecuación en derivadas fraccionarias y un sistema de tales ecuaciones. Precisamente con este objetivo hemos desarrollado Numfracpy. A continuación, se describen los algoritmos empleados, se presentan los métodos implementados en la librería, seguido por su aplicación con ejemplos en varios contextos.
2. Metodología
A continuación indicamos los algoritmos implementados en la resolución numérica de integrales, derivadas y ecuaciones diferenciales.
2.1. Integral de Riemann-Liouville
La integral fraccionaria de Riemann-Liouville de orden α, con α > 0, de una función f (t), está definida por [6, 7, 28]
Esta definición es la noción más generalizada de integral en el cálculo fraccionario y se constituye en la piedra angular del mismo. Permite unificar los conceptos de derivada e integral del cálculo clásico y además, extenderlos a órdenes no enteros. En particular, nótese que rl D-1 a,t f (t ) = ∫ t a f (s) ds. Usualmente se toma el valor a = 0, lo cual seguiremos en este trabajo. También puede notarse que para a = 0 no está definida esta integral por la presencia del factor r(a) en el denominador. Obsérvese además que la misma notación indica que una integral puede verse como una derivada de orden negativo. Numéricamente la integral de Riemann-Liouville puede evaluarse usando algoritmos clásicos para evaluar una integral definida [29, 30, 31]. Por ejemplo, se puede hacer uso de una interpolación polinomial y usar fórmulas como la rectangular, trapezoidal, de Simpson y de Newton-Cotes [29]. También puede hacerse una interpolación por splines cúbicos, una interpolación gaussiana, y aún, usar el método lineal de multipaso [32]. Nosotros calculamos la integral de Riemann-Liouville haciendo uso de la función quad del paquete scipy [33]. Esta función basa su funcionamiento en rutinas de la librería QUADPACK de FORTRAN [34], la cual contiene numerosas funcionalidades para el cálculo de una integral definida y ha sido usada ampliamente durante décadas.
2.2. Derivada de Caputo
La derivada fraccionaria de Caputo de orden α de una función f(t), se fundamenta en la integral de Riemann-Liouville y está dada por [6, 7]
donde m es un entero positivo tal que m - 1 < α < m. Por tanto, para calcular la derivada de Caputo de una función f se deriva primero m veces y luego se aplica la integral de Riemann-Liouville de orden m -α. Asumimos, como es costumbre, que a = 0. En nuestra librería implementamos dos algoritmos denominados L1 y L2 [35] para calcular numéricamente la derivada de Caputo para los casos m = 1 y m = 2 respectivamente. Estos algoritmos tienen muy buena estabilidad [36, 37] y se aproxima de una manera simple la derivada de la función f.
2.2.1. Método L1, caso 0 < α <1
Para 0 < α < 1 la derivada fraccionaria de Caputo está dada por
En este método se realiza una discretización uniforme del intervalo temporal (0, t = t n ) en la forma [0, Δt, 2Δt,...,nΔt], de manera que tk = kΔt. La derivada f (t ) se aproxima por una diferencia finita hacia adelante y se obtiene
donde
2.2.2. Método L2, caso 1 < α < 2
Para 1 < α < 2 la derivada fraccionaria de Caputo está dada por
Haciendo una discretización similar a la del método L1 se cumple
En el último paso se hizo un cambio de variable y se tuvo en cuenta que el proceso se hace sobre todo el recorrido del sumatorio indicado. Aproximando la segunda derivada f´´ (t n -s) por
después de algo de álgebra elemental se obtiene
donde
2.3. Derivada de Riemann-Liouville
La derivada fraccionaria de Riemann-Liouville de orden α para una función f (t), está dada por [6, 7]
donde m = ⌈α⌉, y ⌈ ⌉ representa la función techo. En particular, para 0 < α < 1 se tiene m = 1 y
Para una función f suficientemente suave, es decir f e ∈C m [0,t], se satisface la siguiente relación [35]
Por tanto, los métodos L1 y L2 para los casos 0 < α< 1 y 1 < α< 2 respectivamente, para la derivada de Caputo se pueden utilizar directamente para aproximar la derivada de Riemann-Liouville. Esta es la manera en que calculamos la derivada de Riemann-Liouville en la librería Numfr cpy.
2.4. Derivada de Grünwald-Letnikov
La derivada de Grünwald-Letnikov [6, 7] de orden α > 0 de la función f (t), está definida por
Esta derivada puede verse como una aproximación por sumas de Riemann de a integra de Riemann-Liouvi e de orden -α [38]. Nótese que una integra de orden negativo se concibe como una derivada en e cálculo fraccionario. Se asume que e intervalo en la variable independiente se discretiza con un paso h y se toma α = 0 como de costumbre. En nuestra librería aproximamos la derivada de Grünwald-Letnikov directamente por su definición para 0 < α < 1. Esta definición puede escribirse en la forma
con ωj (α) = (−1) j(α j). Para el caso 1 < α < 2, el esquema anterior puede llevar a inestabilidades numéricas. Por lo cual, hacemos uso de la siguiente fórmula modificada de Grünwald-Letnikov, con p = 1 de acuerdo con [35]
2.5. Solución numérica de una ecuación diferencial
Consideramos ahora el siguiente problema de valor inicial
donde C D α 0,t es la derivada de Caputo y u 0 j es la derivada de u de orden j en t = 0. Este problema es equivalente a la siguiente ecuación integral de Volterra [39]
Para resolver este problema, se implementa el método fraccional de Adams basado en las siguientes relaciones [35].
donde
2.6. Sistema de ecuaciones diferenciales
Algunos métodos, como los de Runge-Kutta pueden extrapolarse en forma natural para considerar un sistema de ecuaciones diferenciales en el cálculo clásico. Para el caso de un sistema de ecuaciones diferenciales en derivadas fraccionarias podemos hacer algo similar con el método fraccional de Adams. En este punto, conviene hacer distinción entre una ecuación multi-término y un sistema multi-orden [40] como se explica a continuación.
2.6.1. Sistema multi-orden
Se define un sistema multi-orden a un sistema de ecuaciones diferenciales de la forma
Asumiendo 0 < α i < 1, las condiciones iniciales serán de la forma y i (0) = y ¡ ,0 con i = 1, 2,..., n. En caso que para alguna i, α i > 1, se deben proporcionar además los valores para y’ i (0), y’ i (0),..., y i m-1 (0) donde m = ⌈αi⌉. En este caso el método de Adams se extiende en forma natural y las relaciones para resolver el sistema de ecs. (22) está dado por
con b i,n y a i,n están dados por las ecs. (20) y (21) respectivamente.
2.6.2. Ecuación multi-término
Una ecuación multi-término en derivadas fraccionarias se define formalmente por la expresión
con condiciones iniciales
donde m = ⌈α n⌉. Se asume, sin pérdida de generalidad, que 0 < α1 < α2 … < α n y los enteros contenidos en el intervalo (0, α n ) se encuentra en la secuencia finita {αk}k n=1. Por tanto, se cumple 0 < αj+1 −αj ≤ 1para j = 1, 2,..., n - 1. Para resolver numéricamente una ecuación multi-orden tenemos en cuenta el siguiente teorema (ver [40]) que nos permite usar el sistema multi-orden asociado a tal ecuación.
Teorema 1 La ecuación multi-término (24) con condiciones iniciales (25) es equivalente al siguiente sistema multi-orden
con las siguientes condiciones iniciales
donde β 1 = α 1 y β j = α j −α j −1 para j = 2,3, . . . ,n. Notese que 0 < β j ≤ 1 para todo j.
3. Resultados y discusión
Numfracpy contiene los métodos que se enuncian en la tabla 2.
Es importante mencionar que los métodos asociados a las funciones de Mittag-Leffler fueron tomados de la fuente [41]. A su vez, este código se basa en el importante artículo de Garrappa [27]. Las funciones de Mittag-Leffler son esenciales para el cálculo fraccionario, y su papel es análogo al de la función exponencial en el cálculo clásico. La función de Mittag-Leffler de dos parámetros se define por la expresión [27]
La función de Mittag-Leffler de un parámetro está dada por E α (z) = E α 1 (z). Adicionalmente, la función de Mittag-Leffler de tres parámetros está dada por [27]
3.1. Comparación de resultados entre derivadas
Ilustraremos los resultados obtenidos para diversas derivadas en nuestra librería comparándolos con los siguientes resultados teóricos [26, 35]
Donde erf(x) es lafuncióndeerror.Recordemos quelarelación entre la derivada de Riemann-Liouville y Caputo está dada por (ec. (13) con m = 1)
En la tabla 1 se ilustran los resultados obtenidos para varias derivadas, en la última columna se muestra el error relativo. Se consideró unadiscretizacióndelintervalo[0,1] con120puntos y =1 2.
Estos resultados deben compararse con los obtenidos usando la librería differint [26] en Python.
3.2. Solución numérica de una ecuación diferencial
La siguiente es una ecuación básica que ha sido empleada para modelar distintos comportamientos [42]
Para 0 < α < 1, u(t) presenta un comportamiento monótono aproximándose a una asíntota horizontal como se ilustra en la figura 1. Este comportamiento puede modelar, por ejemplo, la evolución del esfuerzo en la viscoplasticidad [43].
Este resultado, en numfracpy, se obtiene usando el comando FODE( f, Initial, Interv, dx, alpha) con f (t, u) = 1 - u, Initial = [0], Interv = [0,5], dx = 0.01 y alpha recorre los valores 0.25, 0.50, 0.75 y 1.00 para obtener cada curva.
Análogamente se obtiene la solución a la ec. (36) para α= 1.00, 1.25, 1.50 y 1.75, como se ilustra en la figura 2. En este caso la situación física se equipara a la de un oscilador amortiguado.
Cabe resaltar que los resultados ilustrados en las gráficas 1 y 2 están en claro acuerdo con lo reportado en [42].
3.3. Solución de un sistema de ecuaciones diferenciales
Consideraremos ahora la ecuación de Bagley-Torvik, la cual ha sido usada ampliamente en la literatura [44, 45, 46] para ilustrar la aplicación de técnicas computacionales en una ecuación multi-término. Mediante la ecuación de Bagley-Torvik, que se ilustra a continuación, se modela el movimiento de una placa rígida conectada a un resorte con un punto fijo e inmersa en un fluido Newtoniano.
Donde y( x) representa el desplazamiento de la placa y las constantes A,B,C dependen de la viscosidad, la densidad del fluido, la rigidez del resorte y el área de la placa. Además, f (x) representa la fuerza externa. En esta ecuación, D 2 representa naturalmente la segunda derivada. Se cumple entonces, por la sección 2.6.1, que α1 = 1, α2 = 3/2 y α3 = 2, y por Teorema 1 β 1 = 1, β 2 = 1/2 y β 3 = 1/2. Por tanto, el sistema multi-orden asociado es pues estas son análogas a las condiciones iniciales clásicas e iguales a cero para las derivadas no enteras.
Las condiciones iniciales, conforme al Teorema 1, están dadas por
3.3.1. Ecuación Bagley-Torvik, caso A
Consideremos la ecuación de Bagley-Torvik con los parámetros A = B = C = 1, la función f (x) = x + 1, condiciones iniciales y(0) = 1, y'(0) = 1 y x ∈ [0,5]. La solución analítica está dada precisamente por [45]
Para obtener el resultado con numfracpy utilizamos el comando SystemFODEs(f,Initial,Interv,dx, alpha) para f = [f 1 , f 2 ,f 3], donde f 1 (x,y 1 ,y 2 ,y 3 ) = y2, f 2 (x,y 1 ,y 2 ,y 3 ) = y3 y f 3 (x,y1,y2,y3) = f (x) -y 3 -y1, Initial = [1,1,0], Interv = [0,5], dx = 0.001 y alpha = [1,1/2,1/2].
Nuestro resultado, ver figura 3, presenta una aceptable aproximación a la solución analítica con diferencias del orden de centésimas, este resultado puede contrastarse con el obtenido en [46], el cual desarrolla técnicas especiales para problemas con valor en la frontera. Cabe anotar que nuestro algoritmo utiliza un valor del paso dx bajo y puede mejorarse optimizando los algoritmos codificados. Sin embargo, son realmente escasas las librerías de software que pueden resolver un sistema de ecuaciones diferenciales en derivadas fraccionarias.
Figura 3: Solución de la ecuación Bagley-Torvik con A = B = C = 1, y(0) = 1, y'(0) = 1, x ∈ [0,5] y dx = 0.001. (a) Comparación entre las soluciones analítica y numérica. (b) Valor absoluto de la diferencia entre las soluciones numérica y analítica.
3.3.2. Ecuación Bagley-Torvik, caso B
Ahora consideraremos la ecuación de Bagley-Torvik con los parámetros A = B = C = 1, la función f (x) = x
3
+ 5x +--=
condiciones iniciales y(0) = 0, y'(0) = -1 y x ∈ [0,1]. En este caso la solución analítica está dada por [46]
El comando usado es análogo al que se empleó para el caso A y con el mismo valor para dx. La figura 4 (a) ilustra la comparación entre la solución analítica y la numérica, la cual muestra un muy buen ajuste entre ambos resultados. En la figura 4 (b) se ilustra la diferencia, en valor absoluto, entre los resultados analítico y numérico; la cual resulta ser del orden de milésimas.
4. Conclusiones
La librería de Python, Numfracpy, implementa la integral de Riemann-Liouville y tres diferentes derivadas fraccionarias: Derivada de Caputo, de Riemann-Liouville y de Grünwald-Letnikov; las cuales basan su definición en la primera integral. De acuerdo al desempeño mostrado para las funciones de la tabla 1 podemos observar diferencias relativas del orden de 10-4 con respecto al valor teórico. En comparación con la librería differint en [26] se obtienen mejores aproximaciones. Pocas librerías desarrolladas en Python existen actualmente [47], y es tal vez, differint la única otra librería que desarrolla diversas definiciones de la derivada fraccionaria.
Además, numfracpy logra desarrollar soluciones numéricas para ecuaciones diferenciales de la forma cD α u(t) = f (t,u(t)). Este proceso fue ilustrado con la solución numérica de la ecuación (36) que puede visualizarse en las figuras 1 y 2. Estos resultados están en claro acuerdo con lo reportado en la literatura, ver por ejemplo [42]. Otra característica notable de numfracpy es que puede resolver numéricamente un sistema de ecuaciones diferenciales en derivadas fraccionarias, lo cual no es nada común en los paquetes de software que implementan el cálculo fraccionario. En particular, analizamos la ecuación de Bagley-Torvik que permite modelar el movimiento de una placa conectada a un resorte bajo la acción de una fuerza externa e inmersa en un fluido newtoniano. Nuestro estudio consideró dos tipos de fuerza distintas y los resultados fueron contrastados con los analíticos en las figuras 3 y 4 obteniendo resultados aceptables y con diferencias respecto a la solución numérica en el orden las centésimas y milésimas. En este caso extrapolamos directamente el método fraccional de Adams que fue utilizado para el caso de una ecuación diferencial fraccionaria. Aunque en la literatura se han desarrollado algoritmos con mejor desempeño en casos particulares, ver por ejemplo [46], es notable tener una librería de acceso libre que pueda resolver un sistema de ecuaciones en derivadas fraccionarias. Es de anotar que nuestra librería solo implementa ecuaciones diferenciales usando la derivada de Caputo, como es costumbre en la literatura, pues las condiciones iniciales reflejan las de las derivadas clásicas con una relación física directa.
En resumen, Numfracpy es una librería en Python que desarrolla aproximaciones numéricas para calcular integrales, derivadas y ecuaciones diferenciales propias del cálculo fraccionario, y aunque ciertamente, los algoritmos implementados pueden optimizarse, es tal vez la única librería de Python multipropósito de su clase. Además, la adición de nuestra librería de cálculo fraccionario a las herramientas ya existentes en Python abre el camino a su uso como apoyo en cursos de matemáticas, ecuaciones diferenciales y ciencia en general, siendo una oportunidad para motivar el estudio del cálculo fraccionario y sus diferentes aplicaciones a nivel de estudiantes de pregrado, posgrado y su implementación en temas de investigación.
