Un método iterativo monótono para la

solución numérica de ecuaciones diferenciales

fraccionarias
Edgar Alejandro Guerrero Arroyo
24 de septiembre de 2012

CIMAT

1
”Desear algo es sólo anhelar. Querer conseguirlo es desearlo y
estar dispuesto a hacer lo necesario. Si realmente quieres
conseguir algo, encontrarás la manera. Si no, encontrarás el
pretexto.”

2
Índice
1. Introducción 5

2. Nociones de Cálculo Fraccionario 9

2.1. Aplicaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1. Ecuación fraccional de Schrödinger . . . . . . . . . . . 9
2.1.2. Ecuación integral de Abel . . . . . . . . . . . . . . . . 10
2.1.3. Propagación de ondas ultrasónicas en huesos esponjosos 11
2.1.4. Control lateral y longitudinal de vehı́culos autónomos . 11
2.1.5. Mecánica de fluı́dos . . . . . . . . . . . . . . . . . . . . 11
2.2. Integral de Riemann-Liouville . . . . . . . . . . . . . . . . . . 11
2.2.1. Propiedades básicas . . . . . . . . . . . . . . . . . . . . 12
2.3. Derivada Fraccional de Riemann-Liouville . . . . . . . . . . . 12
2.3.1. No localidad . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.2. Derivada de una función constante . . . . . . . . . . . 16
2.3.3. Funciones con derivada fraccional nula . . . . . . . . . 16
2.4. Otras derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3. Aproximación numérica de la función de Mittag-Leffler 19

3.1. La función de Mittag-Leffler . . . . . . . . . . . . . . . . . . . 19
3.2. División del plano complejo . . . . . . . . . . . . . . . . . . . 20
3.3. Zona G0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4. Zona G1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5. Zona G2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6. Cuadratura de Gauss-Laguerre . . . . . . . . . . . . . . . . . . 28
3.7. Observaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4. Solución numérica de una EDF 30

4.1. Aspectos teóricos de EDFs . . . . . . . . . . . . . . . . . . . . 30
4.2. Una discretización y aproximación numérica del método ite-
rativo monótono continuo . . . . . . . . . . . . . . . . . . . . 33
4.3. Aproximación numérica de la integral singular . . . . . . . . . 36

5. Resultados numéricos 38
5.1. Ejemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2. Ejemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.3. Ejemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.4. Ejemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.5. Ejemplo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.6. Observaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3
6. Estimación del orden de la derivada y sistemas Lotka-Volterra 48
6.1. Estimación del orden de derivación en EDFs . . . . . . . . . . 48
6.1.1. Formulación del problema inverso . . . . . . . . . . . . 48
6.1.2. Algoritmo Levenberg-Marquardt . . . . . . . . . . . . . 48
6.1.3. Algoritmo BFGS . . . . . . . . . . . . . . . . . . . . . 49
6.1.4. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.2. Sistema fraccional Lotka-Volterra . . . . . . . . . . . . . . . . 52
6.2.1. Sub y súper soluciones del sistema . . . . . . . . . . . 53
6.2.2. Iteración del sistema . . . . . . . . . . . . . . . . . . . 54
6.2.3. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . 54

7. Conclusiones 57

4
1. Introducción
Las ecuaciones diferenciales fraccionales, EDFs, han ganado una conside-
rable importancia debido a su variedad de aplicaciones en varios campos de
ciencias aplicadas e ingenierı́a, [11]. Se ha encontrado que el comportamiento
de muchos sistemas fı́sicos puede ser propiamente descrito mediante el uso
de la teorı́a de sistemas de orden fraccional.
Las derivadas fraccionales proporcionan un excelente instrumento para la
descripción de propiedades de materiales y procesos, tales como la memoria
o la herencia. Las ventajas reales de los sistemas de orden fraccional son que
tenemos más grados de libertad en el modelo y que una ”memoria”de estados
anteriores es incluı́da en dicho modelo.
En décadas recientes el campo del cálculo fraccional ha atraı́do el interés
de investigadores en varias áreas incluyendo matemáticas, fı́sica, quı́mica,
ingenierı́a e incluso finanzas y ciencias sociales. Pero, ¿porqué es tan impor-
tante el cálculo fraccional? Como se comenta en [12], hasta hace tiempos
recientes, el cálculo fraccional era considerado una mera teorı́a esotérica sin
aplicaciones, pero en la última década(s) ha habido una explosión de ac-
tividades de investigación en la aplicación del cálculo fraccional a muchos
campos cientı́ficos diversos que van desde la fı́sica del fenómeno de difusión y
advección, el control de sistemas hasta finanzas y economı́a. De hecho, en el
presente, aplicaciones y/o actividades relacionadas con el cálculo fraccional
han aparecido en por lo menos los siguientes campos:
Control fraccional de sistemas de ingenierı́a.
Avances en el cálculo de variaciones y el control óptimo de sistemas
dinámicos fraccionales.
Desarrollo de técnicas y herramientas analı́ticas y numéricas.
Exploraciones fundamentales de las relaciones constitutivas mecánicas,
eléctricas y térmicas, entre otras, de varios materiales de ingenierı́a tales
como polı́meros viscoelásticos, geles, hule espuma y tejidos animales,
ası́ como sus aplicaciones en ciencia e ingenierı́a.
Entendimiento fundamental de los fenómenos de onda y difusión, sus
medidas y verificaciones, incluyendo aplicaciones a la fı́sica de plasma
( tal como en la difusión en un Tokamak ).
Aplicaciones en bioingenierı́a y biomedicina.
Modelación térmica de sistemas de ingenierı́a tales como frenos y partes
de diversas maquinarias.

5
 Procesamiento de imágenes y señales.

Como se menciona en [11], son muchos los esquemas numéricos para resol-
ver EDFs que han sido introducidos en diversas publicaciones. Recientemen-
te, un gran esfuerzo se ha destinado a tratar de encontrar métodos robustos,
tanto analı́ticos como numéricos, para resolver EDFs de interés fı́sico.
En la década pasada, el desarrollo de métodos numéricos para la solu-
ción de ecuaciones que contienen derivada e integrales de orden fraccional
se han convertido en un tema “caliente”. Varios algoritmos simples han sido
publicados, produciendo soluciones aproximadas para EDFs. De hecho, en la
experiencia de [13], expertos en aplicaciones generalmente prefieren aplicar
alguno de los esquemas numéricos ahı́ considerados en vez de cualquier otro
esquema disponible de orden superior. De este modo se consideran usualmen-
te los siguiente métodos:

El método de la cuadratura implı́cita, introducido por Diethelm.

El método de predictor-corrector, ampliamente discutido por Diethelm,

Ford y Freed.

El método aproximado de Mittag-Leffler, considerado por Diethelm y

Luchko.

Un método de colocación, descrito por Blank.

El método de diferencias finitas, discutido por Gorenflo.

Es conocido que la solución numérica de EDFs es un problema muy costo-

so computacionalmente debido a la naturaleza de los operadores diferenciales
fraccionales. En [7], el autor Diethelm demuestra que la paralelización pue-
de ser utilizada para superar esas dificultades. Para este fin, se propone la
implementación de una versión fraccional de segundo orden del método de
Adams-Bashforth-Moulton en una computadora en paralelo.
Tal y como menciona Domenico Delbosco en [6], son muchas las publica-
ciones y libros sobre cálculo fraccional que han aparecido recientemente. La
mayorı́a de ellos son dedicados a la solución de EDFs lineales en términos
de funciones especiales. No existen contribuciones, al menos bien conocidas,
concernientes a EDFs no lineales de la forma

Dα u = f (t, u),
donde 0 < α < 1 y Dα es la derivada fraccional estándar de Riemann-
Liouville, considerada en R+ o en el intervalo (0, T ] con T > 0.

6
 En principio, uno puede reducir dicha ecuación a una ecuación integral
con una singularidad débil y aplicar a ésta técnicas básicas de análisis no
lineal (teoremas de puntos fijos, teorı́a de Leray-Schauder). En la práctica,
para obtener resultados explı́citos, uno debe tomar en cuenta la peculiaridad
del kernel.
En este trabajo consideraremos la existencia y unicidad de la solución
de un problema de condiciones iniciales para una EDF, usando el método
de sub y súper soluciones y su iterativo monótono asociado. Dicha técnica
monótona iterativa, combinada con el método de sub y súper soluciones, es
una herramienta poderosa para probar la existencia de soluciones de EDFs no
lineales. La técnica mencionada se reduce a calcular una secuencia iterativa
donde se deben resolver EDFs lineales no homogéneas en cada iteración. La
solución depende de la evaluación de la función especial de Mittag-Leffler.
En este trabajo se implementa una adaptación del algoritmo descrito en
[2]. Dicho algoritmo consiste en dividir el plano complejo en zonas disjuntas
donde se aproximará el valor de la función de una manera apropiada en cada
una de estas para evitar inestabilidades numéricas en su aproximación.
Mostraremos soluciones de ejemplos en la literatura de EDFs no lineales.
De este modo observaremos como nuestro método se autovalida al iniciar
nuestra secuencia iterativa con una sub y una súper solución.
En el capı́tulo 2 revisaremos algunos conceptos de cálculo fraccional que
nos serán de utilidad más adelante. Mostraremos algunas aplicaciones reales
de cálculo fraccional y revisaremos el concepto de derivada fraccional ası́ como
algunas definiciones de ésta. Pondremos especial énfasis en la definición de
Riemann-Liouville y mostraremos algunas de sus propiedades.
En el capı́tulo 3 expondremos la aproximación de la función de Mittag-
Leffler que implementamos. Describiremos la función y veremos algunas de
sus propiedades. De igual manera explicaremos cómo dividimos el dominio
de la función en zonas y explicaremos como aproximamos la función en cada
una de dichas zonas. De manera complementaria, mostraremos algoritmos de
la implementación que realizamos en cada una.
En el capı́tulo 4 analizaremos la solución numérica de una EDF. Vere-
mos los aspectos teóricos detrás del método propuesto ası́ como la solución
de EDFs lineales y la construcción de la secuencia iterativa que nos permi-
tirá resolver EDFs no lineales. Mostraremos los detalles de implementación
utilizados.
En el capı́tulo 5 ilustraremos el uso de nuestro método con diversos ejem-
plos de la literatura. Explicaremos la selección de sub y súper soluciones
apropiadas para cada caso. De igual manera mostraremos como la sub y la
súper solución convergen a la solución real de cada problema.
En el capı́tulo 6 expondremos la importancia de la estimación del orden

7
de la derivada en una EDF. De igual manera se analizará el uso de nuestro
método iterativo monótono en la solución de un sistema simple de EDFs del
tipo Lotka-Volterra. Se mostrarán resultados numéricos en ambos casos.
Finalmente remarcaremos conclusiones generales de nuestro trabajo en el
capı́tulo 7.

8
2. Nociones de Cálculo Fraccionario
El concepto de cálculo fraccionario ha sido popularizado últimamente por
su practicidad al modelar sistemas çon memoria”. A pesar de esto, no es un
tema completamente nuevo. Se tienen evidencias de que en 1695 L’Hopital
dn y 1
se cuestionó sobre el significado de dxn para n = 2 . En ese entonces, Leibniz

expresó una respuesta para y = x, véase (1). A continuación daremos una

breve introdución sobre el concepto de derivada fracional y algunas de sus
aplicaciones.

2.1. Aplicaciones
Además de las ya mencionadas, mostramos aquı́ algunas aplicaciones más
del cálculo fraccionario y de las EDFs. Dichas aplicaciones fueron tomadas
de [10]. Más aplicaciones pueden ser consultadas en dicho artı́culo y en [4].

2.1.1. Ecuación fraccional de Schrödinger

La mecánica cuántica estandar puede ser estudiada a través de 3 enfoques
distintos:

Matrices mecánicas. Formulado por Werner Heisenberg, Max Born

y Pascual Jordan en 1925. Fue la primer formulación de la mecáni-
ca cuántica. Es una extensión del modelo de Bohr que describe como
ocurre un salto cuántico.

Ecuación de Schrödinger. Describe como el estado cuántico de un

sistema fı́sico cambia con respecto al tiempo. Fue formulado en 1925
por Erwin Schrödinger.La forma de la ecuación de Schrödinger depende
de la situación fı́sica en cuestión. En particular la ecuación dependiente
del tiempo no relativista de una sola partı́cula tiene la forma
d ~2
i~ Ψ(r, t) = − ∆Ψ(r, t) + V (r, t)Ψ(r, t),
dt 2m
donde m es la masa de la partı́cula, V es su energı́a potencial, ∆ es el
operador Laplaciano, Ψ es la función de onda de posición-espacio y ~
es la constante de Planck. Esta ecuación modela la energı́a total del sis-
tema como la suma de la energı́a cinética y la potencial de la partı́cula.
Es el enfoque análogo a los sistemas mecánicos de las ecuaciones de mo-
vimiento de la 2da ley de Newton en mecánica clásica. Las soluciones
de la ecuación de Schrödinger no sólo describen sistemas moleculares,
atómicos y subatómicos, sino también sistemas macroscópicos.

9
 Camino integral de Feynman. En mecánica cuántica, la formula-
ción del camino integral generaliza el principio de acción de la mecánica
clásica. Dicho enfoque remplaza la noción clásica de una única trayec-
toria de una partı́cula por un sistema con una infinidad de posibles
trayectorias para calcular la amplitud cuántica. La idea básica fue for-
mulada por Norbert Wiener pero el método completo fue desarrollado
por Richard Feynman en 1948. El camino integral de Feynman es el
camino integral sobre los caminos mecánico-cuánticos Brownianos.

El concepto de mecánica cuántica fraccional resulta de una generalización

de la mecánica cuántica estandar. Dicho concepto fue introducido por Nick
Laskin en 1999 como resultado de expandir el camino integral de Feynman,
del camino mecánico-cuántico Browniano a un camino mecánico-cuántico
de Lévy. Un camino integral sobre un camino mecánico-cuántico de Lévy
resulta en una generalización de la mecánica cuántica. La ecuación fraccional
de Schrödinger es una ecuación fundamental en mecánica cuántica fraccional.
Esta ecuación tiene la forma
d α
i~ Ψ(r, t) = Dα (−~2 ∆) 2 Ψ(r, t) + V (r, t)Ψ(r, t).
dt

2.1.2. Ecuación integral de Abel

Una de las primeras aplicaciones conocidas de la derivada fraccional fue
hecha por Abel en 1823, la cual consiste en encontrar la solución de una ecua-
ción integral deducida a partir de resolver el problema de la curva tautócrona.
Este problema intenta determinar la forma de una curva plana donde no hay
fricción tal que una partı́cula situada en dicha curva cae por efecto de la
gravedad en un tiempo fijo T que es independiente de su posición inicial. La
ecuación integral de Abel es
p Z η
2gT = (η − y)−0,5 f 0 (y)dy,
0

donde g es la aceleración gravitacional, (ξ, η) es la posición inicial y s = f (y)

es la ecuación de la curva buscada. Dicha ecuación es equivalente a la ecuación
integral fraccionaria
p
2gT = Γ(0,5)D−0,5 f 0 (η).

10
2.1.3. Propagación de ondas ultrasónicas en huesos esponjosos
El cálculo fraccional es usado para modelar la interacción de viscosidad
entre fluı́dos y una estructura sólida. Resultados experimentales son compa-
rados con predicciones teóricas para ondas lentas y rápidas transmitidas a
través de muestras de huesos esponjosos humanos.

2.1.4. Control lateral y longitudinal de vehı́culos autónomos

Esta aplicación utiliza controladores de orden fraccionarios (FOC), apli-
cados al problema de path-tracking en automóviles eléctricos autónomos. Un
modelo de dinámica lateral de vehı́culos industriales ha sido tomado en cuen-
ta para implementar controladores convencionales y FOC . Varios esquemas
de control con dichos controladores han sido simulados y comparados.

2.1.5. Mecánica de fluı́dos

El cálculo fraccional también ha sido utilizado en la solución de problemas
de dinámica de fluı́dos de difusión viscosa dependientes del tiempo. Como ya
hemos señalado con anterioridad, la derivada fraccional modela de manera
natural este tipo de fenómenos fı́sicos. Adicionalmente, también se ha utiliza-
do dicha derivada para modelar el problema de advección-difusión en fluı́dos.
Esto puede describir el movimiento de fluı́dos de diversas densidades a través
de otro fluı́do o medio poroso. Esto sirve para simular la difusión de conta-
minantes en el agua o para el movimiento de petróleo a través de las capas
del subsuelo, las cuales suelen ser no homogéneas.

2.2. Integral de Riemann-Liouville

Denotemos como C 0 (R+ ) el espacio de todas las funciones reales continuas
definidas en R+ = {x ∈ R, x > 0} y sea L1loc (R+ ) el espacio de todas las
funciones reales definidas en R+ tales que son Lebesgue integrables en cada
subintervalo acotado de R+ . Consideremos también el espacio C 0 (R+ 0 ) de
todas las funciones reales continuas en R0 = {x ∈ R, x ≥ 0}. Consideremos
las siguientes definiciones.

Definición. La primitiva fraccional de orden α > 0 de una función q : R+ →

R está dada por
Z t
α 1
I q(t) = (t − s)α−1 q(s)ds,
Γ(α) a
cuando el lado derecho está definido en cada punto en R+ .

11
 Por ejemplo, Iα q existe para todo α > 0, si q ∈ C 0 (R+ ) ∩ L1loc (R+ );
Notar también que cuando q ∈ C 0 (R+ α 0 +
0 ) entonces I q ∈ C (R0 ) y más aún
Iα q(0) = 0.

Observación. Si Iα2 q ∈ Dom(Iα1 ), entonces la igualdad

Iα1 Iα2 q = Iα1 +α2 q,

se cumple para todo α1 , α2 > 0.

Algunas propiedades básicas de esta integral se mencionan a continuación.

2.2.1. Propiedades básicas

Fijemos un intervalo (a, b). Para cada función integrable q en (a, b) el
operador Iα también es integrable en (a, b), convirtiéndolo en un operador
lineal en L1 (a, b). Iα resulta ser continuo respecto a la estructura de espacios
de Banach en L1 . Sea 1 ≤ p < α. En general, si q ∈ Lp (a, b), entonces
Iα q ∈ Lp (a, b) y se satisface la desigualdad

|b − a|α/p
||Iα q||p ≤ ||q||p .
αΓ(α)
Otra propiedad interesante de Iα es que su trasformada de Laplace queda
expresada de una manera particularmente simple

(LIα q)(t) = t−α Q(t),

donde Q denota la transformada de Laplace de q. De este modo Iα es un
multiplicador de Fourier.

2.3. Derivada Fraccional de Riemann-Liouville

Supongamos

f (t) = tk . (1)
Si tomamos la primera derivada de nuestra función obtenemos
d k
t = ktk−1 .
dt
Si repetimos este proceso tenemos en general que
dα k k!
α
t = tk−α .
dt (k − α)!

12
 Generalizando el factorial como la función Γ podemos definir la α-ésima
derivada de f para valores de α no necesariamente enteros de la siguiente
manera

dα k Γ(k + 1) k−α
α
t = t .
dt Γ(k − α + 1)
1
Si tomamos k = 1 y α = 2
obtenemos que
1
d2 Γ(1 + 1) 1− 1 1 1
1 t= 1 t 2 = 3 t2 ,
dt 2 Γ(1 − 2 + 1) Γ( 2 )
1
dándonos una expresión para la derivada de orden 2
de f (t) = t. Si nueva-
mente derivamos como antes obtenemos que
1 1
d2 1 1 1 d2 1 1 Γ(1 + 12 ) 1 − 1
3 t = 1 t = 3 t 2 2 = 1,
2 2
Γ( 23 ) dt 2 Γ( 2 ) Γ( 21 − 12 + 1)
1
dt 2 Γ( 2 )
conduciendo al resultado clásico conocido de
1 1
d2 d2 d
1 t = 1.
1 t=
dt dt 2 dt 2

En las figuras 1 a 3 se muestran derivadas de distintos órdenes fracciona-

rios para algunas funciones conocidas.

Figura 1: Derivada fraccional de diversos órdenes para la función f (t) = t

dα t t1−α
donde dt α = Γ(2−α) . Imagen tomada de [10].

Figura 2: Derivada fraccional de diversos órdenes para la función f (t) = et

donde ddtαe = ∞
α t P tk−α
k=0 Γ(k−α+1) . Imagen tomada de [10].

Figura 3: Derivada fraccional de diversos órdenes para la función f (t) =

α sin(t)
sin(t) donde d dt α = sin(t + απ
2
). Imagen tomada de [10].

En la literatura existen diversas definiciones de lo que representa una

derivada fraccional. Quizás la definición de derivada fraccional más amplia-
mente conocida es la de Riemann-Liouville, RL, expresada como se muestra
en la siguiente definición.

13
Definición. La derivada fraccional de Riemann-Liouville de orden 0 < α <
1 de una función continua u : R+ → R está dada por

d t
Z
1
α
D u(t) = (t − s)−α u(s)ds, (2)
Γ(1 − α) dt 0
asumiendo que el lado derecho está definida en cada punto sobre R+ .

Tenemos que

Dα Iα q = q, ∀q ∈ C 0 (R+ ) ∩ L1loc (R).

De la definición y el ejemplo arriba mencionados obtenemos el siguiente
lema.

Lema. Sea 0 < α < 1. Si q ∈ C 0 (R+ ) ∩ L1loc (R+ ), entonces la EDF

Dα q = 0,
tiene q(t) = ctα−1 , c ∈ R, como únicas soluciones.

De este lema deducimos la siguiente ley de composición.

Proposición. Supóngase que q ∈ C 0 (R+ ) ∩ L1loc (R+ ) con derivada fraccional

de orden 0 < α < 1 que pertenece a C 0 (R+ ) ∩ L1loc (R+ ). Entonces

Iα Dα q(t) = q(t) + ctα−1 ,

para algún c ∈ R. Cuando la función q está en C 0 (R+
0 ), entonces c = 0.

La definición de derivada fraccional de RL generaliza los ejemplos arriba

mostrados de manera natural.
Como mencionamos anteriormente, hay ciertas propiedades conocidas de
las derivadas de orden entero que ya no se cumplen cuando el orden de deriva-
ción es fraccionario. Algunas de estas propiedades se indican a continuación
para el caso 0 < α < 1. En general se considera sólo el caso en el interva-
lo (0, 1) pues para otros intervalos, los resultados se puede representar en
términos del intervalo antes mencionado.

2.3.1. No localidad
En general, para conocer la el valor de dicha derivada en un momento
t se necesita conocer las evaluaciones de Iα q en todos los tiempos anterio-
res. Esto ocurre incluso para el esquema propuesto de diferencias finitas (7).
En general, el número de términos utilizados aumenta conforme t lo hace.

14
Esto representa un problema al calcular la derivada para tiempos grandes
mediante aproximaciones numéricas. Una manera de disminuir el tiempo de
cómputo utilizado es utilizar estrategias de cómputo paralelo. En particu-
lar, podemos consultar [7] donde se explica una estrategia en paralelo para
reducir el tiempo de cómputo de nuestro método iterativo.

2.3.2. Derivada de una función constante

Recordemos que en cálculo clásico tenemos que la derivada de una función
constante es 0. Veamos qué pasa si derivamos una función constante con (2).
Vemos que

d t
Z
k
α
D k= (t − s)−α ds
Γ(1 − α) dt 0
k d  (t − s)1−α t
= − )
Γ(1 − α) dt 1−α 0
k d  t1−α 
=
Γ(1 − α) dt 1 − α
k
= t−α . (3)
Γ(1 − α)
Como vemos, la derivada fraccional de una función constante k es 0 y en
general depende de t. Notamos también que no está definida en t = 0.

2.3.3. Funciones con derivada fraccional nula

Como vimos anteriormente, la derivada de una función constante no es
necesariamente 0. Propongamos funciones del tipo u(t) = tα−k . Vemos que

d t
Z
1
α α−k
D t = (t − s)−α tα−k ds.
Γ(1 − α) dt 0
Recordando el siguiente resultado de integración

Z t
Γ(M )Γ(N )
(t − s)M −1 (s − h)N −1 ds = (t − h)M +N −1 , M, N > 0,
h Γ(M + N )
tenemos que

1 d Γ(1 − α)Γ(1 + α − k) 1−k

Dα tα−k = t
Γ(1 − α) dt Γ(2 − k)

15
 Γ(1 + α − k)
= (1 − k)t−k . (4)
Γ(2 − k)
En particular observamos que lo obtenido bajo k = α coincide con (3)
con constante 1 y que Dα tα = Γ(1 + α) es constante. Adicionalmente, vemos
Dα tα−1 = 0. De lo que hemos discutido anteriormente no es difı́cil ver que,
en general

cΓ(1 + α − k)
Dα ctα−k = (1 − k)t−k .
Γ(2 − k)
El conocer derivadas fraccionales de funciones particulares será útil al
momento de construir el punto de partida de nuestro algoritmo iterativo.

2.4. Otras derivadas

A manera de complemento podemos citar otras definiciones de derivada
fraccional de la literatura que también son utilizadas. Una en particular es
la definición de M. Caputo dada por
Z t
α 1 u(τ )dτ
D u(t) = , n − 1 < α < n. (5)
Γ(α − n) a (t − τ )α+1−n
Usualmente es la segunda definición más popular de derivada fraccional.
En la literatura también aparece la definición de derivada fraccional de
Canavati de orden v dada por

1 d t d(n)
Z
α
Dv u(t) = (t − s)α−1 (n) u(s)ds, t ∈ (0, T ], (6)
Γ(α) dt 0 dt
donde v > 0, n es la parte entera de v y α = v−n. Esta derivada procede de las
propiedades de la integral de Riemann-Liouville. De hecho, es precisamente
la inversa de dicho operador integral.
También hay un esquema de diferencias finitas conocido como la definición
de Grünwald-Letnikove dada por
[ t−α
h
]  
α −α
X α
j
D u(t) = lı́m h (−1) u(t − jh). (7)
h→0
j=0
j

El método iterativo monótono que desarrollamos está basado en un teore-

ma que parte de la definición de RL. Por tanto, utilizaremos dicha definición
para tener resultados teóricos en los cuales apoyarnos. Aun ası́, el método
númerico tiene potencial de aplicación a las demás definiciones de derivadas

16
fraccionales. Más definiciones de la derivada fraccional pueden ser consulta-
das en [10].

17
3. Aproximación numérica de la función de
Mittag-Leffler
La función (generalizada) de Mittag-Leffler Eα,β (z) es una función entera
con dos parámetros α, β ∈ C, [2]. La función Eα,1 lleva el nombre de Mittag-
Leffler, quien la introdujera en 1903 en una publicación en integrales de
Laplace-Abel. Poco después de su introducción dicha función fue generalizada
por Wiman.
La función generalizada de ML con un parámetro α no negativo es monóto-
na si y sólo si 0 < α ≤ 1, Re(β) ≥ α. Esta función juega un papel central
en cálculo fraccional y en sus aplicaciones pues está ı́ntimamente relacionada
con la función propia del operador de derivación fraccional.
El avance en este campo requiere del cálculo de los valores numéricos de la
función generalizada de ML para argumentos complejos arbitrarios, ası́ como
el estudio de sus propiedades.
La función de ML muestra comportamientos muy diferentes en el plano
complejo con respecto a los parámetros α y β, por lo que se deben utilizar
diferentes aproximaciones para diversas regiones en el plano. En [2] se presen-
ta un método numérico estable basado en relaciones de recursión, ası́ntotas
mejoradas exponencialmente, y representaciones integrales para calcular la
función generalizada de ML para parámetros reales α > 0 y β, y argumentos
complejos arbitrarios z.
En sección repasamos la función de ML y mostramos su implementación
para un dominio real. Dividiremos el eje real en zonas donde aproximaremos
numéricamente y de manera apropiada la ML en cada zona basados en [2].

3.1. La función de Mittag-Leffler

La función de ML es una función compleja con 2 parámetros, α, β ∈ C.
Cuando Re(α) > 0 podemos definir dicha función mediante la siguiente serie
de potencias
∞
X zk
Eα,β (z) = . (8)
k=0
Γ(αk + β)
En este caso la serie converge ∀z ∈ C por lo cual Eα,β es una función
entera.
En particular tenemos que si α = β = 1, entonces
∞ ∞
X zk X zk
E1,1 (z) = = = ez , (9)
k=0
Γ(k + 1) k=0 k!

18
que es la función exponencial compleja. Por otro lado, si consideramos la
EDF más simple

Dα u = u,
entonces sabemos de [3] que para β = 0 la solución queda expresada como

u(t) = tα−1 Eα,α (t).

Observamos una similitud al caso entero donde dtd et = et . En este sentido,
la ML es una generalización de la exponencial. Resulta de especial interés
para aplicaciones en fı́sica el caso donde β = 1 donde la solución queda como
u(t) = Eα,1 (t).

3.2. División del plano complejo

Para nuestros fines suponemos z ∈ R y α = β ∈ (0, 1). Para aproximar
la función de ML dividimos el eje real en secciones donde la función muestra
diferentes comportamientos de estabilidad. Luego entonces, se deben de uti-
lizar diferentes esquemas de aproximación para las diferentes particiones del
mismo.
Denotemos la bola abierta de radio r centrado en el origen y su cerradura
correspondiente como

D(r) = {z ∈ R : |z| < r}

y

D(r) = {z ∈ R : |z| ≤ r},

respectivamente. De este modo dividiremos el eje real en 3 zonas disjuntas en
las cuales la ML será aproximada de manera distinta. Esto para procurar una
estabilidad numérica que nos permita obtener una aproximación aceptable
de la función en cada una de las zonas.
La primer zona la denotaremos por G0 = D(r0 ) para un r0 < 1. La zona
−
2 la denotaremos por G1 como la unión de 2 subzonas, G+ 1 y G1 , definidas
respectivamente por

G+
1 = {z ∈ R \ D(r1 ) : z > 0},

G−
1 = {z ∈ R \ D(r1 ) : z < 0},

19
para un r1 > r0 . Finalmente, la tercer zona, denotada por G2 , corresponde a
la zona entre G0 y G1 y la definiremos como la unión de dos subzonas G+ 2 y
−
G2 definidas por

G+
2 = {z ∈ D(r1 ) \ G0 : z > 0},

G−
2 = {z ∈ D(r1 ) \ G0 : z < 0}.

A continuación daremos un “paseo” a través de cada una de las zonas en

cuestión para dar detalles sobre la aproximación utilizada según sea el caso
ası́ como algoritmos para realizar dicha aproximación.

3.3. Zona G0
En esta zona la aproximación de ML por medio de series de potencias
M −1
X zk
Eα,α (z) = , (10)
k=0
Γ(α(k + 1))
resulta adecuada. Para el algoritmo optamos por tomar r0 = 0,95 basándonos
en [2]. El error cometido al aproximar ML por (10) lo podemos controlar
mediante el siguiente teorema.

Teorema. Sea  > 0 y z ∈ C. Si |z| < 1 y

    
2−β ln((1 − |z|))
M ≥ max + 1, +1 ,
α ln(|z|)

entonces el error cometido queda acotado como

∞
X zk
RG0 (z) = ≤ ,
k=M
Γ(αk + β)

donde [a] denota el mı́nimo entero mayor que a. Tomamos un valor de  =

1e−8 al momento de acotar el error. El algoritmo 1 refleja como realizar dicha
aproximación en G0 .
En las figuras 4 y 5 se muestran las aproximaciones de la ML en G0
usando el algoritmo 1 para r0 = 0,95 y varios valores de α.

20
Algorithm 1 Aproximando la ML en G0
Require: z ∈ G0 ; α, β ∈ (0, 1),  > 0.
Ensure: aprox
 = Eα,β (z) con precisión
 .

2−β ln((1−|z|))
M ← max α
+ 1, ln(|z|)
+1
aprox = 0
for k = 0 → M do
zk
aprox ← aprox + Γ(αk+β)
end for
return aprox

(c)
(a)
(b)
α=
0,1
0,3
0,5

Figura 4: Aproximaciones de la ML en G0 .

(b)
(a)
α=
0,7
0,9

Figura 5: Aproximaciones de la ML en G0 .

3.4. Zona G1
Definimos la notación
 
1 1−β 1 x
A(z, α, β, x) = z α exp z α cos .
α α
Luego entonces, aproximaremos la ML en esta zona mediante ası́ntotas
exponencialmente mejoradas donde, para z ∈ G1 , tendremos que
 PM −1 −k
z
 − k=1 Γ(β−αk)
 z ∈ G−
1
Eα,β (z) = (11)
 A(z, α, β, 0) − PM −1 z−k
 +
z ∈ G1 .
k=1 Γ(β−αk)

El error al aproximar ası́ la ML está controlado por el siguiente teorema

1
Teorema. Sea α ∈ (0, 1),  > 0 y z ∈ C. Para M ≈ α1 |z| α y |z| ≥ (−2log C )α

21
el error cometido por (11) satisface
M −1
z −k
 X 
RG1 (z) = Eα,β (z) − − ≤ ,
k=1
Γ(β − αk)

donde C es una constante que solo depende de α y β.

Dicha constante es aproximada como

1
C≈ . (12)
π sin(πα)

Para nuestro algoritmo tomamos, basados en el teorema anterior

   α
1 1 
M= |z| α + 1, r1 = − 2 log , (13)
α C

donde [x] denota el mı́nimo entero mayor que x. El algoritmo 2 refleja como
realizar dicha aproximación en G1 .

Algorithm 2 Aproximando la ML en G1
Require: z ∈ G1 ; α, β ∈ (0, 1),  > 0.
Ensure: aprox = Eα,β (z) con precisión .
1
M ← α1 |z| α + 1

1
C ← π sin(πα)
 α
r1 ← − 2 log C
aprox = 0
for k = 1 → M do
z −k
aprox ← aprox − Γ(β−αk)
end for
if z ∈ G+
1 then
approx ← approx + A(z, α, β, 0)
end if
return aprox

En las figuras 6 y 7 se muestran las aproximaciones de la ML en G− 1

mientras que en 8 y 9 se muestran aproximaciones en G+1 usando (2) para r0 =
0,95 y varios valores de α. Acotamos los valores extremos de z dependiendo
del valor de α pues de este último depende el tamaño de r1 . Dicho lı́mite se
menciona en las figuras correspondientes.

22
 (a)
(b)
(c)
−3 ≤≤
−7
−10
z≤
−r1 , α =
0,1
0,3
0,5

Figura 6: Aproximaciones de la ML en G−
1.

(a)
(b)
−15 ≤
−25
z≤
−r1 , α =
0,7
0,9

Figura 7: Aproximaciones de la ML en G−
1.

(c)
(b)
(a)
r1 ≤
z≤
10,
3, αα==
7,
0,5
0,3
0,1

Figura 8: Aproximaciones de la ML en G+
1.

(a)
(b)
r1 ≤
z≤
15, α =
25,
0,7
0,9

Figura 9: Aproximaciones de la ML en G+
1.

3.5. Zona G2
Introducimos las siguientes notaciones
1 y 1−β
ω(x, y, α, β) = x α sin + y(1 + ),
α α

1 r sin[ω(r, φ, α, β) − φ] − z sin[ω(r, φ, α, β)]

B(r, α, β, z, φ) = A(r, α, β, φ) .
π r2 − 2rz cos φ + z 2
Luego entonces, aproximaremos ML en G2 mediante las siguientes repre-

23
sentaciones integrales

 R∞
 0 B(r, α, β, z, 2πα
3
)dr z ∈ G−
2
Eα,β (z) = R∞ (14)
A(z, α, β, 0) + B(r, α, β, z, πα)dr z ∈ G+
2.

0

El error cometido al aproximar la ML con (14) queda acotado por

Z ∞
RG2 (Rmax , α, β, z, φ) ≤ B(r; α, β, z, φ)dr < ,
Rmax

que claramente depende la la elección de Rmax . De [2] sabemos que

β
α

 max 2α , 2|z|, − 2 ln π212  z ∈ G−

2
Rmax ≥
α (15)
max 1, 2|z|, − ln π z ∈ G+

2,

6

para 0 ≤ β ≤ 1. La evaluación numérica de la parte integral sugerida en [2]

es Gauss-Kronrod o Gauss-Lobatto. Sin embargo, para nuestros propósitos
requerimos evaluar ML muchas veces y rápidamente por lo cual optamos por
utilizar una cuadratura de Gaussiana con nodos de Gauss-Laguerre 1 que
resulta más apropiada en nuestro caso. En la siguiente sección veremos esta
cuadratura un poco más a detalle.
El algoritmo 3 refleja como realizar dicha aproximación en G2 . En nuestro
caso aproximamos la integral utilizando 15 nodos de Laguerre y definimos
ψ(r) = f (r)er para tener f (r) = ψ(r)e−r .

Algorithm 3 Aproximando la ML en G2
Require: z ∈ G2 ; α, β ∈ (0, 1),  > 0.
Ensure: aprox = Eα,β (z) con precisión .
if z ∈ G−
2 then
R∞
aprox ← 0 B(r, α, β, z, 2πα 3
)dr {Integrando con Gauss-Laguerre}
else R∞
aprox ← A(z, α, β, 0) + 0 B(r, α, β, z, πα)dr
end if
return aprox

En las figuras 10 y 11 se muestran las aproximaciones de la ML en G− 2

mientras que en 12 y 13 se muestran aproximaciones en G+2 usando (2) para
r0 = 0.95, r1 como en (13) y varios valores de α.
1
Ver [15]

24
 (a)
(b)
(c)
α=
0,1
0,3
0,5

Figura 10: Aproximaciones de la ML en G−

2.

(a)
(b)
α=
0,7
0,9

Figura 11: Aproximaciones de la ML en G−

2.

(c)
(a)
(b)
α=
0,1
0,3
0,5

Figura 12: Aproximaciones de la ML en G+

2.

(b)
(a)
α=
0,9
0,7

Figura 13: Aproximaciones de la ML en G+

2.

3.6. Cuadratura de Gauss-Laguerre

Los polinomios ortogonales juegan un papel crucial en la elaboración de
fórmulas de cuadratura de máximo grado de exactitud. Sean x0 , ..., xn puntos
distintos Ren el intervalo [−1, 1]. Para la aproximación de la integral pesada
1
Iw (q) = −1 q(t)w(t)dt, siendo q ∈ C 0 ([−1, 1]), consideramos una regla de
cuadratura del tipo
n
X
In,m (q) = wk q(xk ), (16)
k=0

donde los pesos wk son coeficientes por determinar de un modo apropiado y

xk son los nodos donde se evaluará la función que no son más que los ceros de
los polinomios ortogonales en el intervalo [−1, 1] de la cuadratura utilizada.
De [15] obtenemos el siguiente resultado.

Lema. El máximo grado de exactitud de la fórmula de cuadratura (16) es

2n+1.

25
 Nótese que podemos mapear los nodos y pesos en [−1, 1] a un intervalo
arbitrario [a, b] mediante
a+b b−a
t = Φ(ξ) = ξ+ ,
2 2
obteniendo de este modo que
Z b
a+b 1
Z
q(t)dt = (q ◦ Φ)(ξ)dξ.
a 2 −1

Luego entonces, podemos usar la cuadratura con pesos wk0 = a+b 2

wi y
0
nodos xi = Φ(xi ) en el intervalo [a, b]. Notar que esta fórmula mantiene en
el intervalo [a, b] el mismo grado de exactitud que el de la fórmula generada
en [−1, 1].
Sin embargo, en nuestro caso el intervalo no está acotado superiormente.
Para este caso utilizamos una fórmula de cuadratura Gaussiana cuyos nodos
son ceros de los polinomios ortogonales de Laguerre. Estos polinomios, orto-
gonales en el intervalo [0, +∞) con respecto a la función de peso w(t) = e−t ,
están definidos como
dn −t n
Ln (t) = et (e t ), n ≥ 0,
dtn
y satisfacen la relación de recursividad

Ln+1 (t) = (2n + 1 − t)Ln (t) − n2 Ln−1 (t), n ≥ 0,
L−1 = 0, L0 = 1.
t
R ∞ Para cualquier
R ∞ −t función q, definimos φ(t) = q(t)e . Entonces, I(q) =
0
q(t)dt = 0 e φ(t)dt, obteniendo de este modo que
n
X (n!)2 (2n)
I(q) = wk φ(xk ) + φ (ξ), 0 < ξ < +∞. (17)
k=1
(2n!)
Denotemos por P2n−1 el espacio de polinomios de grado 2n−1. De la ecua-
ción anterior uno concluye que la fórmula de cuadratura de Gauss-Laguerre
es exacta para funciones q del tipo φe−t , donde φ ∈ P2n−1 . En un sentido
generalizado, podemos decir que se tiene un grado óptimo de exactitud de
2n − 1.

3.7. Observaciones
Al implementar la ML en el eje real notamos que hay algunas discon-
tinuidades numéricas al pasar de una zona a otra. Dichas discontinuidades

26
varı́an en base a diversos factores como son el número de nodos de Laguerre
tomados y la manera de construir la ψ en G2 ası́ como el valor que toma α.
Sin embargo dichas discontinuidades no representan una fuente considerable
de error por lo cual aproximaremos la ML del modo descrito en esta sección.
También notamos que para valores de α cercanos a 0 se requieren más su-
mandos para aproximar las series de potencias en G0 y G1 . Análogamente,
el número requerido de dichos sumandos disminuye al tener α cercano a 1.
Con esta sección damos por sentado la aproximación numérica de la ML para
valores reales de z y α = β ∈ (0, 1).

27
4. Solución numérica de una EDF
En este capı́tulo introduciremos la teorı́a necesaria para desarrollar nues-
tro método iterativo monótono, también estableceremos dicho método en el
caso continuo, su discretización y detallaremos su implementación. La de-
mostración de los resultados puede ser consultada en [1].

4.1. Aspectos teóricos de EDFs

Nos enfocaremos en considerar la existencia y unicidad de la solución del
siguiente problema de valores iniciales para una EDF, usando el método de
sub y súper soluciones y su método monótono iterativo asociado
 α
D u(t) = f (t, u), t ∈ (0, T ].
(18)
t1−α u(t)|t=0 = u0 ,

donde 0 < T < +∞ y Dα representa la derivada fraccional de Riemann-

Liouville de orden 0 < α < 1.
EDFs aparecen frecuentemente en diferentes areas de investigación e in-
genierı́a, tales como quı́mica, fı́sica, control de sistemas dinámicos etc. Re-
cientemente, mucha gente ha puesto atención en el resultado de existencia
de la solución del problema de valores iniciales para EDFs. A continuación
enlistaremos algunos resultados importantes.

Observación. El problema de valores iniciales

 α
D u(t) = f (u), t ∈ (0, T ].
u(a) = u0 ,
donde 0 < a < T , tiene sentido.

Observación. Sea 0 < α < 1. Asumamos f (x) ∈ C(R+ ) ∩ L1loc (R+ ). Enton-
ces para todo (a, b) ∈ R+ × R, el problema de valores iniciales
 α
D u(x) = f (x)
u(a) = b,
tiene una única solución en C(R+ ) ∩ L1loc (R+ ) dada por

Z a  xα−1 Z x
 1 α−1 1
u(x) = b − (a − t) f (t)dt α−1 + (x − t)α−1 f (t)dt.
Γ(α) 0 a Γ(α) 0

28
Observación. Sea C1−α ([0, T ]) = {u ∈ C(0, T ]; t1−α u ∈ C([0, T ])}. Supon-
gamos 0 < α < 1 y y(x) ∈ C1−α ([0, b]).
(a) Si

lı́m [x1−α y(x)] = c, c ∈ R,

x→0+
entonces

I1−α y(0+) := lı́m+ I1−α y(x) = cΓ(α).

x→0

(b) Si

lı́m I1−α y(x) = b, b∈R

x→0+

y existe el lı́mite lı́mx→0+ [x1−α y(x)], entonces

b
lı́m+ [x1−α y(x)] = .
x→0 Γ(α)
Definición. Llamaremos a una función u(t) una solución clásica del proble-
ma (18), si:
(i) u(t) es continua en (0, T ], t1−α u(t) es continua en [0, T ], y si la integral
fraccional I1−α u(t) es continuamente diferenciable para t ∈ (0, T ];
(ii) u(t) satisface (18).
Para el problema (18), tenemos las siguientes definiciones de sub y súper
soluciones.
Definición. Una función u ∈ C1−α ([0, T ]) es llamada una súper solución
del problema (18), si satisface
 α
D u(t) ≥ f (t, u), t ∈ (0, T ].
(19)
t1−α u(t)|t=0 ≥ u0 .
Definición. Una función u ∈ C1−α ([0, T ]) es llamada una sub solución
del problema (18), si satisface
 α
D u(t) ≤ f (t, u), t ∈ (0, T ].
(20)
t1−α u(t)|t=0 ≤ u0 .
De aquı́ en adelante, supondremos que

u(t) ≥ u(t), t ∈ (0, T ]; t1−α u(t)|t=0 ≥ t1−α u(t)|t=0 , (21)

y definimos el sector

29
 hu, ui = hu ∈ C1−α ([0, T ]); u(t) ≥ u(t) ≥ u(t), t ∈ (0, T ],

t1−α u(t)|t=0 ≥ t1−α u(t)|t=0 ≥ t1−α u(t)|t=0 i.

El siguiente es un resultado de existencia de la solución para el proble-
ma lineal de valores iniciales para EDFs y una propiedad de la definición
de Riemann-Liouville de cálculo fraccional, los cuales son importantes para
obtener resultados de existencia y unicidad de las soluciones del problema
(18).

Lema. El problema lineal de valores iniciales

 α
D u + γu = q(t), t ∈ (0, T ].
(22)
t1−α u(t)|t=0 = u0 ,
donde γ es una constante y q ∈ C([0, T ] × R) tiene por solución la siguiente
representación integral

Z t
α−1 α
u(t) = Γ(α)u0 t Eα,α (−γt ) + (t − s)α−1 Eα,α (−γ(t − s)α )q(s)ds, (23)
0

donde Eα,α denota la función de ML. En particular, cuando γ = 0, el

problema de valores iniciales (22) tiene por solución
Z t
α−1 1
u(t) = u0 t + (t − s)α−1 q(s)ds.
Γ(α) 0
El siguiente resultado desempeña un papel muy importante en nuestro
análisis restante.

Lema. Si w ∈ C1−α ([0, T ]) y satisface las relaciones

 α
D w + γ0 w ≥ 0, t ∈ (0, T ],
t1−α w(t)|t=0 ≥ 0,
donde γ0 > − Γ(1+α)
Tα
es una constante. Entonces w ≥ 0 para t ∈ (0, T ].

Suponemos que f satisface la siguiente condición.

f (t, u1 ) − f (t, u2 ) ≥ −γ(u1 − u2 ), u ≤ u1 ≤ u2 ≤ u, (24)

donde γ > − Γ(1+α)
Tα
es una constante y u, u ∈ C1−α ([0, T ]) son respectivamen-
te sub y súper soluciones de (18). Claramente, esta condición se satisface con

30
γ = 0, cuando f es monótona no decreciente en u. En términos de (24), la
función

F (t, u) = γu + f (t, u) (25)

es monótona no decreciente en u para u ∈ hu, ui.
Γ(1+α)
Suponemos también la existencia de una constante Tα
> γ ≥ −γ, tal
que

f (t, u1 ) − f (t, u2 ) ≤ γ(u1 − u2 ), u ≤ u1 ≤ u2 ≤ u, (26)

donde u, u ∈ C1−α ([0, T ]) son respectivamente sub y súper soluciones de (18).
El siguiente es un teorema de existencia y unicidad de la solución del
problema (18).

Teorema. Supongamos que u, u ∈ C1−α ([0, T ]) son respectivamente sub y

súper soluciones de (18), tales que (21) se cumple, f ∈ C([0, T ] × R), y
(24) se satisface. Entonces el problema (18) tiene una solución en el sector
hu, ui. Adicionalmente, si (26) se satisface, entonces el problema (18) tiene
una única solución en el sector hu, ui.

El siguiente corolario se desprende de los resultados anteriores.

Corolario 1. Sea Lipγ (hu, ui) el espacio de las funciones f que satisfacen la
condición |f (t, u1 ) − f (t, u2 )| ≤ γ|u1 − u1 |, u1 , u2 ∈ hu, ui. Si f ∈ Lipγ (hu, ui)
satisface
Γ(1 + α)
0<γ< , (27)
Tα
entonces (18) tiene una única solución en el sector hu, ui.

4.2. Una discretización y aproximación numérica del

método iterativo monótono continuo
Partiendo de la demostración del teorema de existencia y unicidad ante-
rior de Zhang Shuqin en [1] podemos ver que el problema (18) es equivalente
al siguiente problema
 α
D v + γv = γv + f (t, v), t ∈ (0, T ].
1−α (28)
t v(t)|t=0 = u0 ,
donde γ es una constante conocida.
Para poder aplicar nuestro método de sub y súper soluciones y su asociado
monótono iterativo para considerar la existencia de la solución del problema

31
(18), añadimos γ (γ es la constante de (24)) en ambos lados de la EDF
del problema (18), y escogemos un iterando adecuado v (0) ; construimos la
secuencia {v (k) }, k = 1, 2, ..., del siguiente proceso iterativo
 α (k)
D v + γv (k) = γv (k−1) + f (t, v (k−1) ), t ∈ (0, T ].
1−α (k) (29)
t v (t)|t=0 = u0 .
Como para cada k el lado derecho de (29) es conocido, el lema (23) implica
que la secuencia {v (k) } está bien definida. Resulta de particular interés la
secuencia obtenida de (29) con una sub o súper solución del problema (18)
como iteración inicial. Denotemos la secuencia con la iteración inicial v (0) = u
como {v (k) } y la secuencia con v (0) = u como {v (k) }. Las secuencias {v (k) } y
{v (k) } construidas con (29) mantienen las propiedades de monotonı́a

u ≤ v (k) ≤ v (k+1) ≤ v (k+1) ≤ v (k) ≤ u, t ∈ (0, T ], (30)

para todo k = 1, 2, .... De hecho, sea r = v (0) − v (1) . De (21), (19) y (29), y
v (0) = u, tenemos que

Dα r + γr = (Dα u + γu) − (γv (0) + f (t, v (0) ))

= Dα u − f (t, u) ≥ 0, t ∈ (0, T ],

t1−α r(t)|t=0 ≥ u0 − u0 = 0.
En vista del lema (4.1), r ≥ 0 para t ∈ (0, T ], lo que implica que
v ≤ v (0) = u, t ∈ (0, T ]. Un argumento similar usando la propiedad de
(1)

sub solución del problema (18) conduce a que v (1) ≥ v (0) = u, t ∈ (0, T ]. Sea
r(1) = v (1) − v (1) . Por (24), (25) y (21), tenemos que

Dα r(1) + γr(1) = γv (0) + f (t, v (0) ) − (γv (0) + f (t, v (0) ))

= γ(u − u) + f (t, u) − f (t, u)

≥ 0, t ∈ (0, T ],

t1−α r(1) (t)|t=0 = v (1) (0) − v (1) (0) = 0.

De nuevo, por el lema (4.1), r(1) ≥ 0 para t ∈ (0, T ], la conclusión de
arriba muestra que

32
 u = v (0) ≤ v (1) ≤ v (1) ≤ v (0) = u, t ∈ (0, T ].
Supongamos, por inducción

u ≤ v (k−1) ≤ v (k) ≤ v (k) ≤ v (k−1) ≤ u, t ∈ (0, T ]. (31)

Entonces por (29), (25), (26) y (31), la función r(k) = v (k) −v (k+1) satisface
las relaciones

Dα r(k) + γr(k) = γv (k−1) + f (t, v (k−1) ) − (γv (k) + f (t, v (k) ))

≥ 0, t ∈ (0, T ],

t1−α r(k) (t)|t=0 = 0.

Por el lema (4.1), r(k) ≥ 0, esto es v (k+1) ≤ v (k) para t ∈ (0, T ]. Un
razonamiento similar conlleva a v (k) ≤ v (k+1) y v (k+1) ≤ v (k+1) para t ∈ (0, T ].
Ası́, la propiedad de monotonı́a (30) continua por el principio de inducción.
De este modo podemos resolver una EDF no lineal mediante el método
monótono iterativo partiendo de una sub o súper solución del problema (18).
En cada iteración el lado derecho de (29) es conocido por lo cual podemos
utilizar el resultado (23) para resolver la EDF lineal asociada. En este con-
texto, nos resta mencionar como podemos construir sub y súper soluciones
adecuadas del problema (18) para comenzar nuestro proceso iterativo. Para
este fin el siguiente corolario será de gran utilidad.
Corolario 2. Supongamos que el corolario 1 se cumple y que

f (t, 0) ≥ 0, u0 ≥ 0. (32)
Supongamos que existe una constante ρ > 0 tal que

ρt−α
f (t, ρ) ≤ , t ∈ (0, T ], u0 ≤ 0. (33)
Γ(1 − α)
Luego entonces, el problema (18) tiene una única solución u ∈ h0, ρi.
Con esto en mente podemos construir sub y súper soluciones de (18) cons-
tantes de manera relativamente simple para tiempos t ∈ (0, T ] adecuados.
Computacionalmente, lo más trascendente de calcular en la solución del caso
lineal es la evaluación numérica de la función de ML y aproximación de la
integral singular en (23). Luego entonces, nos resta realizar un análisis de
esta última.

33
4.3. Aproximación numérica de la integral singular
La parte integral en la solución de (23) debe ser tratada con cuidado pues
contiene una singularidad en el lı́mite de integración superior. En este caso
un cambio de variable es usualmente la solución. Partimos de

ψ(s)ds = (t − s)α−1 Eα,α (−d(t − s)α )q(s)ds. (34)

Denotemos por ϕ(s) = Eα,α (−d(t − s)α )q(s). Realizando un cambio de
variable del tipo
1
s = t − xα , (35)
tenemos que

Z t Z t Z tα
α−1 1 1−α 1 1
ψ(s)ds = (t − s) ϕ(s)ds = x α (t − t + x α )α−1 ϕ(t − x α )dx
0 0 α 0

Z tα Z tα
1 1 1 1 1
= ϕ(t − x )dx =
α Eα,α (−d(t − t + x α )α )q(t − x α )dx
α 0 α 0

Z tα Z tα
1 1 1
= Eα,α (−dx)q(t − x )dx =
α ψ 0 (x)dx.
α 0 α 0
De este modo obtenemos una integral no singular. La idea del cambio de
variable se puede analizar más a detalle en [16]. Para aproximar esta última
integral optamos por una cuadratura con nodos de Gauss-Legendre. De este
modo, la aproximación de dicha integral queda expresada como
Z t
1X
ψ(s)ds = wi ψ 0 (xi ), (36)
0 α i

siendo wi y xi los pesos y nodos de la cuadratura. Para mayor detalle de las

cuadraturas se puede consultar [15]. De este modo, ya estamos en condiciones
de calcular (23).
Con esto tenemos el camino libre para calcular la solución de una EDF
no lineal de manera numérica. Partiendo de una sub o súper solución de (18)
conocida podemos calcular la secuencia (29) numéricamente. Para ilustrar
esto mostramos el algoritmo 4 para la solución de una EDF no lineal. En
este algoritmo suponemos f ∈ Lipγ (hu, ui) por lo que la solución es única en
el sector. Por el momento asumiremos que conocemos una sub o súper solu-
ción de (18). También conocemos el intervalo de tiempo (0, T ] y el grado de

34
discretización P rec de las funciones. En este último punto, P rec denotará la
cantidad de subtiempos entre dos unidades de tiempo. Posteriormente, ilus-
traremos en la siguiente sección el proceso completo con varios ejemplos.

Algorithm 4 Aproximando la solución de una EDF no lineal

Require: Sub o súper solución v 0 de (18), α ∈ (0, 1), tiempo de calculo T ,
Tamaño de la discretización P rec, condición inicial u0 , constante Lipschitz
γ de f en el sector, número máximo de iteraciones permitidas maxiter.
Ensure: Discretización de u(t).
Obtener una discretización de tamaño P rec de v 0 .
for k = 1 : maxiter do
Obtener v k+1 de solucionar el problema lineal
 α (k)
D v + γv (k) = γv (k−1) + f (t, v (k−1) ), t ∈ (0, T ].
1−α (k)
t v (t)|t=0 = u0 .

end for
return u = v k .

Para obtener una aproximación a la solución de una EDF lineal podemos

utilizar el algoritmo 5.

Algorithm 5 Aproximando la solución de una EDF lineal

Require: α ∈ (0, 1), tiempo de cálculo T , Tamaño de la discretización P rec,
condición inicial u0 , constante Lipschitz γ de f en el sector, discretización
de v k−1 .
Ensure: Discretización de v k .
Obtener una discretización de tamaño P rec de q(t) = γv (k−1) +f (t, v (k−1) ).
for t ∈ (0, T ] do
Calcular la solución vk dada por (23) utilizando (36) y aproximando Eα,α
como se ilustró anteriormente.
end for
return v k .

Con esto en mente estamos en condiciones de pasar a ilustrar nuestro

método con ejemplos prácticos.

35
5. Resultados numéricos
A continuación ilustramos el uso de nuestro método con varios ejemplos
de EDF no lineales. En dichos ejemplos se muestra a detalle la selección de
sub y súper soluciones y demás consideraciones pertinentes.

5.1. Ejemplo 1
Sea f (t, u) dada por

f (t, u) = u2 , u0 = 0, t ∈ (0, 100], (37)
2
donde  > 0 es una constante. Podemos verificar las hipótesis necesarias para
utilizar el corolario 2. Inicialmente vemos que (32) se cumple pues f (t, 0) = 0.
Para no tener conflictos con (33) asumimos u0 = 0. Fijemos la constante
ρ = 1 para tener el sector h0, 1i. Es simple verificar que f ∈ Lip (h0, 1i) pues
∀u1 , u2 ∈ h0, 1i
  
| u21 − u22 | = |u1 − u2 ||u1 + u2 |
2 2 2
 
|u1 + u2 | ≤ 2 = .
2 2
Obtenemos que f ∈ Lip (h0, 1i). Para encontrar una constante  apropia-
da consideremos que se deben satisfacer (27) y (33). Para esta última debe
pasar que

 t−α
≤ , t ∈ (0, 100].
2 Γ(1 − α)
Si t es cercano a 0, entonces (33) se cumple sin problemas. Si t = T ,
entonces debe pasar que
 1
≤ ,
2 Γ(1 − α)T α
por lo que tomamos
1
= .
Γ(1 − α)T α
Esta elección de  satisface (27) pues

1 Γ(1 + α)
0<= α
≤ ⇐⇒
Γ(1 − α)T Tα

36
 1
≤ Γ(1 + α).
Γ(1 − α)
Esta última condición se cumple para varios valores de α ∈ (0, 1). En
particular para α = 0,5. En este ejemplo tomaremos este valor.
Habiendo cumplido las hipótesis del corolario 1 tenemos que u = 0 y
u = 1 son respectivamente sub y súper soluciones de (18). En las figuras 14
y 15 mostramos algunas iteraciones de (29) para las sub y súper soluciones
arriba citadas.

(a)
(b)
(c)
Ite
210

Figura 14: Secuencia iterativa (iteraciones 0-2) para la sub y súper soluciones
de (37).

(c)
(b)
(a)
Ite
543

Figura 15: Secuencia iterativa (iteraciones 3-5) para la sub y súper soluciones
de (37).

Notemos cómo la escala de las figuras cambia conforme itera nuestro

método; en la iteración 0 la escala es de 1e2 mientras que en la iteración 5
es de 1e−5 . En la tabla siguiente mostramos cómo decae la norma entre las
discretizaciones de u(t) y u(t) denotadas respectivamente como u[t] y u[t]
dada por ||u[t] − u[t]||.

Iteración ||u[t] − u[t]||

0 71009.253
1 152.2688
2 8.2745
3 0.41154
4 0.018951
5 0.00081568
Notamos una rápida convergencia al obtener un error no significativo en
tan sólo 5 iteraciones del algoritmo. Esto refleja la pronta convergencia entre
sub y súper soluciones a la solución real, que en este caso es u(t) = 0.

37
5.2. Ejemplo 2
Variando ligeramente la función del ejemplo anterior obtenemos

f (t, u) = u2 + t, u0 = 0, t ∈ (0, 4]. (38)
2
Al igual que antes tomamos  = γ. Usando el corolario 2 como en el
ejemplo anterior, variamos el valor de la súper solución a ρ = 25Γ(1 − α).
En este contexto, y de manera similar al ejemplo (37), tenemos que f ∈
Lipγ (h0, ρi). Tomamos  = γ = ρ22 . Variando la subsolución, vemos que
u = −1 es sub solución de (18) para el ejemplo (38), pues satisface (23).
Tomando α = 0,5 mostramos el comportamiento de la secuencia iterativa
(29) del problema (38) para unas cuantas iteraciones en las figuras 16 y 17.

(c)
(b)
(a)
Ite
210

Figura 16: Secuencia iterativa (iteraciones 0-2) para la sub y súper soluciones
de (38).

(c)
(a)
(b)
Ite
543

Figura 17: Secuencia iterativa (iteraciones 3-5) para la sub y súper soluciones
de (38).

Igual que antes, mostramos como decae la norma ||u[t] − u[t]|| en la si-
guiente tabla comparativa.

Iteración ||u[t] − u[t]||

0 2865.7412
1 6.9772
2 0.015409
3 3,0645e−5
4 5,562e−8
5 9,3318e−11

38
5.3. Ejemplo 3
Tratemos un ejemplo un poco más complicado tomado de [5]. Considere-
mos

40320 8−α Γ(5 − α2 ) 4− α α

4 3 3
f (t, u) = t −3 α t
2 + 2,25Γ(1 + α) + (1,5t 2 − t ) − u 2 ,
Γ(9 − α) Γ(5 + 2 )
(39)
1
para t ∈ (0, 20 ], u0 = 0. Buscando satisfacer las hipótesis del corolario 2
tomamos ρ = 1. Nos interesa que f ∈ Lipγ (h0, 1i). Para encontrar un
γ apropiado utilizamos el teorema del valor medio para obtener que γ =
d
sup{| du f (t, u)|h0,1i }. Luego entonces
n d o n 3 1 o 3
γ = sup f (t, u) = sup u2 = .
du h0,1i 2 h0,1i 2
No es difı́cil ver que f (t, 0) ≥, por lo que (32) se cumple. El tiempo máxi-
1
mo T = 20 fue escogido para satisfacer las hipótesis restantes más fácilmente.
Para el mismo fin tomaremos α = 0,5. Claramente vemos que

3 Γ(1 + α) Γ(1,5) √
0< < = = 20 Γ(1,5) ≈ 3,9633,
2 Tα √1
20

por lo cual se satisface (27). Finalmente vemos que

40320 8−α Γ(5 − α2 ) 4− α α

4 3
f (t, 1) = t −3 α t
2 + 2,25Γ(1 + α) − 1 + (1,5t 2 − t ) .
Γ(9 − α) Γ(5 + 2 )

No es complicado de observar numéricamente que para el valor escogido

de α se cumple (33) en el intervalo (0, T ]. Mostramos el comportamiento de
la secuencia iterativa (29) del problema (39) para unas cuantas iteraciones
en las figuras 18 y 19.

(c)
(a)
(b)
Ite
210

Figura 18: Secuencia iterativa (iteraciones 0-2) para la sub y súper soluciones
de (39).

En la tabla siguiente mostramos como decae la norma ||u[t] − u[t]||.

39
 (a)
(b)
(c)
Ite
543

Figura 19: Secuencia iterativa (iteraciones 3-5) para la sub y súper soluciones
de (39).

Iteración ||u[t] − u[t]||

0 7.0711
1 0.51222
2 0.042302
3 0.0032694
4 2,3627e−4
5 1,6553e−5

5.4. Ejemplo 4
Consideremos la función tomada de [6] dada por
M tα
f (t, u) = (u + M )2 , M > 0, u0 = 0, t ∈ (0, 1]. (40)
πα
Supongamos que M = 41 y α = 0,5. De nueva cuenta utilizaremos el
corolario 2 para obtener sub y súper soluciones de (18) para (40). Tomemos
nuevamente ρ = 1. Utilizando el teorema del valor medio obtenemos que
n d o n 2M tα o
γ = sup f (t, u) = sup (u + M )
du h0,1i πα h0,1i

2M (M + 1)T α 5T α
= = .
πα 8π α
Vemos que

5T α Γ(1 + α)
0< α
≈ 0,3526 < ≈ 0,8862,
8π Tα
por lo que se cumple (27). Por otro lado, observamos que

M 3 tα
f (t, 0) = > 0 ∀t ∈ (0, T ],
πα
por lo que (32) se satisface. Finalmente, no es difı́cil ver que numéricamente

40
 25tα tα
f (t, 1) = < ∀t ∈ (0, T ],
64π α Γ(1 − α)
por lo cual se satisface (33). En las figuras 20 y 21 mostramos el comporta-
miento de la secuencia iterativa (29) para el problema (40).

(c)
(a)
(b)
Ite
310

Figura 20: Secuencia iterativa (iteraciones 0-2) para la sub y súper soluciones
de (40).

(c)
(b)
(a)
Ite
975

Figura 21: Secuencia iterativa (iteraciones 3-5) para la sub y súper soluciones
de (40).

En la tabla siguiente mostramos como decae la norma ||u[t] − u[t]||.

Iteración ||u[t] − u[t]||

0 10
1 3.1218
2 0.18166
3 0.0071331
4 0.00022292
5 5,8087e−6

5.5. Ejemplo 5
Consideremos la función no lineal f dada por

f (t, u) = au − bu2 , u0 = 1, t ∈ (0, 200]. (41)

Consideremos los valores constantes de a = 0,1 y b = 1e−4 . Esta ecuación
modela el crecimiento de una especie en el caso de derivadas enteras y su
solución es representada por la función logı́stica. Dicha función resulta estar
asintóticamente acotada por 0 y ab . De este modo tomaremos un α = 0,95

41
para mostrar el comportamiento parecido al caso clásico. Al saber esto de
antemano, podemos proponer u(t) = 0 y u(t) = ab = 1000. Utilizamos γ =
0,05.
En las figuras 22 y 23 mostramos el comportamiento de la secuencia
iterativa(29) para el problema (40).

(a)
(b)
(c)
Ite
950

Figura 22: Secuencia iterativa (iteraciones 0-2) para la sub y súper soluciones
de (41).

(c)
(a)
(b)
Ite
15
20
30

Figura 23: Secuencia iterativa (iteraciones 3-5) para la sub y súper soluciones
de (41).

En la tabla siguiente mostramos como decae la norma ||u[t] − u[t]||.

Iteración ||u[t] − u[t]||

0 141421.3562
6 119054.6834
9 60967.4806
15 12119.3008
20 1804.0748
25 0.52813
Notamos como requerimos de más iteraciones en este caso al utilizar tiem-
pos grandes. Como fue pronosticado, la cota superior de la solución es ab .
Otros experimentos realizados muestran como la curva solución se modifica
en escala al disminuir el parámetro α. En general, observamos que mientras
α disminuye, la cota superior disminuye también mientras que el tiempo re-
querido para obtener la solución debe aumentar. Por citar un ejemplo, para
α = 0,5 la cota superior es de 820 y se requiere calcular la solución para
T = 2200.

42
5.6. Observaciones
El método iterativo propuesto aproxima una solución de manera rápida
al comenzar a iterar con una sub o una súper solución del problema. Esto
refleja la rapidez y robustez del método propuesto en la solución de EDFs no
lineales.
Un inconveniente al momento de implementar este algoritmo es el tiempo
elevado de cómputo utilizado. Sin embargo, este puede ser dramáticamente
reducido sin necesidad de paralelizar el algoritmo. Notemos primero que el
mayor tiempo de cómputo se dedica a calcular la aproximación de la ML.
Los puntos a evaluar en dicha función cambian al momento de cambiar el
intervalo de tiempo. Sin embargo, también notamos que en cada iteración
los puntos a evaluar son exactamente los mismos para cada iteración pues
la ML sólo depende del tiempo actual de cálculo el cual se repite en cada
iteración. De este modo podemos ahorrar bastante tiempo al guardar los
valores calculados de la ML durante la primera iteración para no volver a
calcularlos en las iteraciones restantes. Ası́, solamente la primera iteración
calculará los valores de la ML que usaremos a lo largo de todas las iteraciones.
Otro detalle que nos permite ahorrar tiempo es el de usar pocos nodos de
Legendre y de Laguerre en la aproximación de (36) y de la ML respectiva-
mente. En nuestro caso usamos sólo 15 sumandos lo cual nos da una precisión
bastante aceptable para prácticamente cualquier fin. También es bueno que
dichos valores de los sitios y pesos de la cuadratura sean precalculados una
sola vez y solamente ser leı́dos durante su uso en el programa.

43
6. Estimación del orden de la derivada y sis-
temas Lotka-Volterra
En este capı́tulo exploramos la sensibilidad de EDFs al orden fraccional
de diferenciación y la aplicabilidad del método iterativo monótono a sistemas
de segundo orden del tipo Lotka-Volterra.

6.1. Estimación del orden de derivación en EDFs

En contraste con modelos clásicos, se puede suponer que un modelo es
fraccional aunque se desconozca el orden de la derivada. Luego entonces, de
manera natural hay que estimarlo. En este contexto surgen dos problemas.
En primer lugar la sensibilidad de la solución de una EDF al orden de deri-
vación fraccional asociado y en segundo lugar la factibilidad de estimar dicho
orden. Un problema inverso natural en este contexto es el utilizar la función
f y estimaciones de û obtenidas con nuestro método para estimar el coefi-
ciente fraccionario de derivación utilizado. A continuación describiremos el
problema a detalle y mostraremos algunos resultados numéricos.

6.1.1. Formulación del problema inverso

Para este problema generaremos datos sintéticos û para una f y α dados.
Podemos notar que este problema inverso carece de solución única por lo que
tendremos que conformarnos con la de norma mı́nima. Denotemos por F(α)
el operador de nuestro método iterativo para resolver (18) con el α y f dados.
De este modo F(α) = û. Podemos aproximar la solución a nuestro problema
inverso al minimizar la siguiente función
1X
f (α) = (ûk − Fk (α))2 . (42)
2 k=1
Para minimizar f compararemos los métodos BFGS y Levenberg-Marquardt,
(LM). El primero fue elegido por ser un método iterativo muy robusto. La
estructura de la función a minimizar nos exige aprovechar la estructura de
mı́nimos cuadrados por lo que el LM fue elegido también.

6.1.2. Algoritmo Levenberg-Marquardt

El LM es un algoritmo de optimización que provee una solución numéri-
ca al problema de minimizar una función, generalmente no lineal, dentro de
un espacio de parámetros de la función. Dichos problemas de minimización
aparecen especialmente al ajustar mı́nimos cuadrados lineales y no lineales.

44
El LM interpola entre el algortimo de Gauss-Newton, GN, y un descenso de
gradiente. LM es más robusto que el GN en general por esta modificación
al momento de resolver las ecuaciones normales. Esencialmente, el algorit-
mo se comporta como un descenso de gradiente lejos de la solución y como
un GN cerca de ésta. El secreto del LM es alterar la diagonal de la matriz
cuadrada del producto de Jacobianos en las ecuaciones normales del GN.
Hay varias maneras de resolver las ecuaciones normales. Podemos usar Cho-
lesky si cond(J T ) no es muy grande pues en otro caso el sistema será mal
condicionado conduciendo a errores numéricos en la solución sin mencionar
que el Jacobiano puede tener muchas escalas. También podemos utilizar QR
que es más estable que Cholesky en este sentido. Idealmente usarı́amos SVD
pero la complejidad de programación de este ultimo nos orillá a optar por
QR. Adicionalmente es una buena idea comprobar que la solución de las
ecuaciones normales sı́ mejora el punto donde nos encontramos. Si no es ası́,
entonces podemos hacer un backtracking sobre el parámetro que multiplica
a la identidad hasta que obtengamos una mejora. El algoritmo del LM se
puede encontrar en [17].

6.1.3. Algoritmo BFGS

El algortimo BFGS es un método para resolver problemas de optimización
no lineal. Este método aproxima el método de Newton por lo cual puede no
converger a menos que la función sea cuadrática en su expansión de Taylor
alrededor del punto óptimo. Este método utiliza solamente el gradiente de la
función.
La idea del algoritmo se puede sintetizar en lo siguiente: partiendo de
un punto inicial x0 obtenemos la evaluación del gradiente en dicho punto
∇f (x0 ) := ∇f0 . Inicializamos una matriz H0 como una matriz simétrica y
definida positiva. Luego calculamos la dirección de Newton pk = −H0 ∇f0 .
Usando dicha dirección calculamos un valor α que satisfaga suficiente descen-
so. Luego actualizamos x1 = x0 + αp0 . Posteriormente comenzamos a iterar
como se muestra a continuación en el algoritmo 6.
Podemos utilizar un esquema de diferencias finitas para calcular el gra-
diente numérico de la función. Como condición de suficiente descenso pode-
mos usar la condición de Armijo utilizando backtracking para encontrar el
αk . Al inicializar H0 se pueden considerar las siguientes opciones:

H0 = ∇2 f0 .

H0 = J0T J0 + τ I (Levenberg-Marquardt).

H0 = I.

45
Algorithm 6 Algoritmo del BFGS
Require: x0 ∈ Rn , f ∈ C 2 , tol ∈ R;
k = 0;
while ∇fk > tol do
sk = xk − xk − 1;
yk = ∇fk − ∇fk−1 ;
sk = xk − xk − 1;
ρk = yt1sk ;
k
Hk = (I − ρk sk ykT )Hk−1 (I − ρk yk sTk ) + ρk sk sTk ;
pk = −Hk ∇fk ;
Calculamos αk tal que satisfaga suficiente descenso en pk ;
xk+1 = xk + αk pk ;
end while

Resaltamos en particular la tercer opción pues no tenemos que calcular

segundas derivadas (como en la primera) y tampoco tenemos que preocupar-
nos por un posible número de condición grande del Jacobiano J0 (como en
la segunda) además de lucir más simple e intuitivo.

6.1.4. Resultados
Para ejemplificar nuestro análisis utilizaremos la función dada en (39)
y un valor de α = 0,5. En las siguientes tablas mostramos los resultados
obtenidos.

BFGS
LM
Iteración α f (α)
Iteración α f (α)
0 0.8 1.60717
0 0.8 1.60717
1 0.578122 0.367011
1 0.738178 9.70465
2 0.469005 0.0855874
2 0.656725 7.29619
3 0.532405 0.0741657
3 0.503453 0.211426
4 0.510856 0.00899683
4 0.500205 0.0126457
5 0.494708 0.00226772
5 0.500014 0.000851371
6 0.500661 3.46344e−5
6 0.500001 5.76922e−5
7 0.500038 1.14161e−7
7 0.5 3.91115e−6
8 0.5 6.19775e−12
8 0.5 2.65157e−7
9 0.5 6.19775e−12

46
 LM vs BFGS
α = 0.1 BFGS LM α = 0.2 BFGS LM α = 0.3 BFGS LM
Iters - - Iters - 9 Iters - 8
Time - - T - 17.4 T - 15
α = 0.4 BFGS LM α = 0.45 BFGS LM α = 0.55 BFGS LM
Iters - 7 Iters 7 7 Iters 6 7
T - 14 T 28.2 13. T 23.6 13.8
α = 0.6 BFGS LM α = 0.7 BFGS LM α = 0.75 BFGS LM
Iters 8 8 Iters 9 8 Iters - 9
T 29.95 15.8 T 71.1 15.1 T - 17.3
α = 0.8 BFGS LM α = 0.9 BFGS LM α = 0.99 BFGS LM
Iters - 8 Iters - 8 Iters - -
T - 15.5 T - 16.4 T - -
Apreciamos como el LM resulta ser más robusto para los puntos iniciales
utilizados a lo largo de todo el dominio de α, además de ser más rápido
también que el BFGS. Esto último se debe a que el BFGS requiere de más
evaluaciones de función por las aproximaciones de las derivadas y por el
backtracking que contiene. Mientras α está más cercano a 0, el cálculo de la
Mittag-Leffler es más costoso.
Resulta interesante el conocer como es la dependencia de α y F(α). Para
este fin mostramos dicha gráfica para varios valores de α en (0,1) en la figura
24. En dicho intervalo, variamos el tamaño de la discretización utilizada.

(c)
(a)
(b)
10
20
30
pun-
tos

Figura 24: Gráfica de α vs F(α) para distintos tamaños de discretización.

Al graficar α vs F(α) notamos que el comportamiento engañosamente

estable de solo utilizar una discretización de 10 puntos en (0,1) se desenmas-
cara si aumentamos el número de puntos. Aunado a esto está el hecho de que
F(α) se dispara cuando α se aproxima a 0. Este comportamiento se ilustra
en la figura anteriomente mencionada.

6.2. Sistema fraccional Lotka-Volterra

Ejemplificaremos ahora el uso de nuestro método para aproximar la so-
lución de un sistema tipo Lotka-Volterra de EDFs tomado de [9]. Partimos
del sistema tradicional de presa-depredador de Lotka-Volterra dado por

47
 x0 = x(A − Bx − Cy) (43)
y 0 = y(−D + Ex), (44)
para ciertas condiciones iniciales. En este contexto, x e y son funciones de
t ∈ (0, T ] y A, B, C, D, E escalares. En una forma más general podemos
expresar el sistema anterior en su forma fraccional como
 α
D 1 x = x(A − Bx − Cy)
(45)
t1−α x(t)|t=0 = x0 ,
 α
D 2 y = y(−D + Ex)
(46)
t1−α y(t)|t=0 = y0 ,
para t ∈ (0, T ].
Supongamos α1 = α2 = 0,9 y los valores escalares A = 0,12, B = 0,
C = E = 6e−6 , D = 0,18 en un itervalo de tiempo (0, 60]. Dichos valores son
tomados para mantener una relación de 1000 a 1 entre el término lineal y el
cuadrático. Esto representa un modelo realista como se ilustró en el ejemplo
5 con (41). Tomemos condiciones iniciales x0 = 1 y y0 = 2. De este modo
nuestro sistema queda expresado como
 α
D x + γ1 x = −Exy
(47)
t1−α x(t)|t=0 = 1,
 α
D y + γ2 y = Exy
(48)
t1−α y(t)|t=0 = 2,
donde γ1 = −0,12, γ2 = 0,18 y E = 6e−6 ayudan a denotar los coeficientes
lineales y cuadráticos respectivamente. Nos basaremos en la teorı́a básica co-
nocida para tratar sistemas bajo derivadas enteras. Dicha información puede
ser consultada en [14]. Para esto construiremos sub y súper soluciones del
sistema. Posteriormente, partiremos de dichos sistemas iniciales e iteraremos
con un predictor-corrector.

6.2.1. Sub y súper soluciones del sistema

Si de la ecuación (47) eliminamos el término cuadrático nos quedamos
con la ecuación
 α
D x + γ1 x = 0
t1−α x(t)|t=0 = 1,
la cual tiene una solución analı́tica dada por (23). Dicha solución queda
expresada como

48
 x(t) = Γ(α)tα−1 Eα,α (−γ1 tα ). (49)
No es difı́cil ver que x es súper solución de (47) pues
 α
D x + γ1 x = 0 ≥ −Exy
t1−α x(t)|t=0 = 1.
De este modo, podemos usar nuestro método para calcular la solución de
la ecuación restante. En este caso no hace falta sumar una constante del lado
derecho como se hizo anteriormente en (29). En vez de eso, podemos iterar
directamente como
 α k+1
D y + γ2 y k+1 = Exy k
1−α k+1
t y (t)|t=0 = 2.
Luego, x y y son súper soluciones de nuestro sistema de ecuaciones. Del
mismo modo que antes, podemos considerar la ecuación (48) sin la parte
cuadrática. De este modo tenemos
 α
D y + γ2 y = 0
t1−α y(t)|t=0 = 2,
la cual tiene una solución analı́tica expresada como

y(t) = Γ(α)tα−1 Eα,α (−γ2 tα ). (50)

Podemos notar que

Dα y + γ2 y = 0 ≤ Exy


t1−α y(t)|t=0 = 2,
por lo cual y es sub solución de (48). Utilizando nuestra propuesta iterativa
podemos obtener la solución de
 α k+1
D x + γ1 xk+1 = −Exk y
1−α k+1
t x (t)|t=0 = 1.
Análogamente, x y y son sub soluciones del sistema de ecuaciones. Este
será nuestro punto de partida para poder comenzar nuestro proceso iterativo
de una solución del sistema formado por (47) y (48).

49
6.2.2. Iteración del sistema
Reflejando nuevamente la teorı́a establecida para el caso con derivadas
de orden entero, hacemos uso de un predictor corrector como iterador. De
este modo, inicializamos x0 , y 0 con las sub o súper soluciones antes descritas
y calculamos un predictor-corrector donde en cada iteración se resuelven las
siguientes ecuaciones
 α k
D X + γ1 X k = −Exk y k
(51)
t1−α X k (t)|t=0 = 1,
 α k
D Y + γ2 Y k = EX k y k
(52)
t1−α Y k (t)|t=0 = 2,
 α k+1
D x + γ1 xk+1 = −EX k Y k
1−α k+1 (53)
t x (t)|t=0 = 1,
 α k+1
D y + γ2 y k+1 = EX k Y k
(54)
t1−α y k+1 (t)|t=0 = 2.
En cada iteración se resuelven 4 EDFs. Las primeras 2 sirven de predictor
y las ultimas 2 de corrector. Con esto tenemos un algoritmo iterativo para
resolver sistemas de este tipo.

6.2.3. Resultados
En las figuras 25 y 26 se muestran las iteraciones del sistema acorde al
predictor corrector antes descrito e iniciando en las sub y súper soluciones
del sistema respectivamente.

(c)
(b)
(a)
Ite
210

Figura 25: Iterados de la solución del sistema de ecuaciones (47) comenzando

en la sub solución del mismo.

50
 (a)
(b)
(c)
Ite
210

Figura 26: Iterados de la solución del sistema de ecuaciones (47) comenzando

en la súper solución del mismo.

Apreciamos como la convergencia del sistema es inmediata y ambos pun-

tos iniciales conducen a la misma solución. Adicionalmente, mostramos en la
figura 27 la gráfica de puntos de la forma (x(t), y(t)).

51
Figura 27: Gráfica de la solución (x(t)), y(t)).

52
7. Conclusiones
Hemos demostrado como el método de solución propuesto es útil al re-
solver EDFs no lineales en un enfoque general y de manera robusta. Este
método es eficiente y aproxima la solución esperada en un número pequeño
de iteraciones (alrededor de 5). Los algoritmos mostrados para el cálculo
numérico utilizado fueron implementados sin problema de manera directa.
De este modo también se proporciona una guı́a para implementar el algorit-
mo propuesto sin contratiempos teniendo especial cuidado en especificar los
detalles importantes al momento de realizar aproximaciones numéricas.
También se ilustró, mediante diversos ejemplou os, como construir sub y
súper soluciones de nuestro problema ası́ como ciertos parámetros necesarios
para comenzar a iterar con nuestro método. Se hizo particular énfasis en la
utilidad del corolario 2, pues permite la construcción de sub y súper soluciones
lineales de una manera relativamente simple.
Mediante el problema inverso ilustramos la relación de nuestro método
con el ı́ndice fraccionario. Dicha relación arrojó resultados que reflejan un
comportamiento amigable para valores cercanos a 1 y problemáticos en el
extremo opuesto (valores cercanos a 0). Esto resalta la conveniencia de nues-
tro método para problemas prácticos pues en la mayorı́a de las aplicaciones
reales se trabaja con ı́ndices fraccionarios en el intervalo [0.5,1).
Vimos como podemos utilizar nuestro método iterativo para resolver sis-
temas de EDFs, en especial del tipo Lotka-Volterra. Además de mostrar la
convergencia de nuestra propuesta de predictor-corrector, hemos proporcio-
nado un algoritmo para que su implementación no tenga mayores complica-
ciones.
Con todo lo anterior ponemos en evidencia la robustez de nuestro método
al resolver diversas EDFs no lineales, su rápida convergencia al utilizar solo
pocas iteraciones para obtener las soluciones para diversos tamaños de tiempo
y su rapidez al momento de ejecución de nuestro problema. Las primeras dos
caracterı́sticas resultan muy novedosas al no haber en la literatura métodos de
solución para EDFs no lineales en una perspectiva general. El tercer punto
resalta una caracterı́stica muy deseable en este tipo de algoritmos pues es
bien conocido que la mayor complicación al tratar con EDFs es su efecto de
memoria lo cual se refleja en un alto tiempo de cómputo.
De este modo, esperamos que las contribuciones hechas en esta tesis sirvan
para el desarrollo de mejores estrategias de solución de EDFs ası́ como de
poder contribuir a la solución real de una amplia gama de problemas reales
modelados con EDFs.

53
Referencias
[1] Zhang Shuqin, Monotone iterative method for initial value problem in-
volving Riemann-Liouville fractional derivates, Nonlinear Analysis 71
(2009), 2087-2093.
[2] Hansjörg Seybold and Rudolf Hilfer, Numerical algorithm for calculating
the generalized Mittag-Leffler function, SIAM J. Number. Anal Vol. 47,
No. 1, 69-88,2008.
[3] H. J. Seybold, R. Hilfer, Numerical results for the generalized Mittag-
Leffler function, Fractional calculus & applied analysis, Vol 8, No 2,
2005.
[4] Mehdi Rahimy, Applications of Fractional Differential Equations, Ap-
plied Mathematical Sciences, Vol. 4, no. 50, 2010, 2453 - 2461.
[5] F. Farokhi, M. Haeri, M. S. Tavazoei, Comparing Numerical Methods
for Solving Nonlinear Fractional Order Differential Equations, Springer
Science+Business Media B.V. 2010.
[6] Domenico Delbosco, Luigi Rondino, Existence and Uniqueness for a
Nonlinear Fractional Differential Equation, Journal of Mathematical
Analysis and Applications 204, 1996, 609-625.
[7] Kai Diethelm, An Efficient Parallel Algorithm for the Numerical Solu-
tion of Fractional Differential Equations, Fractional Calculus & Applied
Analysis, Vol 14, Num 3, 2011.
[8] Igor Podlubny, Fractional Differential Equations, Academic Press, 1999.
[9] Ivo Petrás, An Effective Numerical Method and Its Utilization to Solu-
tion of Fractional Models Used in Bioengineering Applications, Advances
in Differential Equations, 2011.
[10] Mehdi Dalir, Majid Bashour, Applications of Fractional Calculus, Ap-
plied Mathematical Science, Vol. 4, no. 21,1021-1032, 2010.
[11] Zaid M. Odibat, Analytic study on linear systems of fractional differen-
tial equations, Computers and Mathematics with Applications 59, 2010,
1171-1183.
[12] J. Tenreiro Machado, Virginia Kiryakova, Francesco Mainardi, Recent
history of fractional calculus, Commun Nonlinear Sci Numer Simulat,
2010.

54
[13] Neville J. Ford, Joseph A. Connolly, Comparison of numerical methods
for differential equations, Communications on Pure an Applied Analysis,
Volume 5, Number 2, June 2006.

[14] G. S. Ladde, V. Lakshmikantham, A. S. Vatsala, Monotone Iterative

Techniques for Nonlinear Differential Equations, Pitman Advances Pu-
blishing Program, 1985.

[15] Alfio Quarteroni, Riccardo Sacco, Fausto Saleri, Numerical Mathema-

tics, Texts in Applied Mathematics 37, Springer, pag 428.

[16] William H. Press, Brian P. Flannery, Saul A. Teukolsky, William T.

Vetterling, Numerical Recipes in C: The Art of Scientific Computing,
Second Edition, Cambridge University Press, 2002.

[17] A Brief Description of the Levenberg-Marquardt Algorithm Implemen-

ted by levmar, Manolis I. A. Lourakis Institute of Computer Scien-
ce Heraklion, Crete, GREECE February 11, 2005, www.ics.forth.gr/
~lourakis/levmar/levmar.pdf.

55