Discretización Grünwald-Letnikov y Riemann-Liouville

Un operador no local obliga a revisar la historia en cada tiempo nuevo. La discretización debe aproximar la integral o la diferencia que define el operador y, al mismo tiempo, dejar visible cuánto cuesta conservar el pasado. Se desarrollan la suma de Grünwald-Letnikov y una cuadratura de la integral de Riemann-Liouville, junto con sus pesos, errores y costos.

Coeficientes de Grünwald-Letnikov

En una malla uniforme $t_n=a+nh$, la forma discreta básica es

$${}^{GL}D_{a+}^{q}x(t_n) \approx h^{-q}\sum_{k=0}^{n}w_k^{(q)}x_{n-k},$$

con

$$w_k^{(q)}=(-1)^k\binom{q}{k}.$$

Calcular gamma para cada $k$ es innecesario. Los pesos obedecen

$$w_0^{(q)}=1, \qquad w_k^{(q)}= \left(1-\frac{q+1}{k}\right)w_{k-1}^{(q)}.$$

La recurrencia se deriva del cociente entre dos binomios consecutivos. Para $0<q<1$, $w_0=1$ y los pesos restantes son negativos. La suma acumulada se aproxima a cero al crecer el historial, en concordancia con la derivada de una constante en el límite de historia apropiado; cerca de un terminal finito aparecen efectos de borde que deben conservarse.

Das desarrolla la aproximación GL y una recurrencia para sus coeficientes Das, 2020 Shantanu Das (2020) Kindergarten of Fractional Calculus Cambridge Scholars Publishing Ubicación consultada: cap. 3, §3.23, pp. 146-147 . Oliveira subraya que la suma puede truncarse para un cálculo numérico, lo que convierte el corte en parte explícita de la aproximación Oliveira, 2019 Edmundo Capelas de Oliveira (2019) Solved Exercises in Fractional Calculus Springer Ubicación consultada: cap. 5, §5.1, pp. 171-174; cap. 6, problema 25, p. 261 .

Semiderivada de una rampa

Tomemos $x(t)=t$, $q=1/2$, $a=0$, $t=1$ y cuatro subintervalos, de modo que $h=1/4$. Los primeros pesos son

$$w_0=1,\quad w_1=-\frac12,\quad w_2=-\frac18,\quad w_3=-\frac1{16},\quad w_4=-\frac5{128}.$$

Los valores históricos, desde el presente hacia atrás, son

$$x_4=1,\quad x_3=\frac34,\quad x_2=\frac12,\quad x_1=\frac14,\quad x_0=0.$$

La suma ponderada vale

$$1-\frac12\frac34-\frac18\frac12 -\frac1{16}\frac14=0.546875.$$

Como $h^{-1/2}=2$,

$${}^{GL}D^{1/2}x(1)\approx1.09375.$$

La referencia analítica RL —que coincide aquí con Caputo— es

$${}^{RL}D_{0+}^{1/2}t\big|_{t=1} =\frac{\Gamma(2)}{\Gamma(3/2)} =\frac{2}{\sqrt\pi} \approx1.128379.$$

El error absoluto de esta malla gruesa es aproximadamente $0.03463$. El número aislado no informa convergencia; se debe repetir con $h/2$, $h/4$ y una referencia evaluada con mayor precisión.

Cuadratura de la integral de Riemann-Liouville

Para $\alpha>0$,

$${}^{RL}I_{0+}^{\alpha}f(t_n) =\frac1{\Gamma(\alpha)} \sum_{j=0}^{n-1} \int_{t_j}^{t_{j+1}} (t_n-\tau)^{\alpha-1}f(\tau)\,d\tau.$$

Si $f$ se aproxima por el valor izquierdo $f_j$ en cada celda,

$${}^{RL}I_{0+}^{\alpha}f(t_n) \approx\frac{h^\alpha}{\Gamma(\alpha+1)} \sum_{j=0}^{n-1}b_{n-j-1}^{(\alpha)}f_j,$$

donde

$$b_m^{(\alpha)}=(m+1)^\alpha-m^\alpha.$$

Los pesos surgen al integrar exactamente el kernel y aproximar solo la función. Una interpolación lineal conduce a otra regla y otro orden de error. Das presenta aproximaciones por valores medios y valores ponderados para la integral RL Das, 2020 Shantanu Das (2020) Kindergarten of Fractional Calculus Cambridge Scholars Publishing Ubicación consultada: cap. 2, §2.23, pp. 96-99 .

Para comprobar la fórmula se toma $f(t)=1$. La suma telescópica da

$$\sum_{m=0}^{n-1}b_m^{(\alpha)}=n^\alpha,$$

y, por tanto,

$${}^{RL}I_{0+}^{\alpha}1(t_n) \approx\frac{h^\alpha n^\alpha}{\Gamma(\alpha+1)} =\frac{t_n^\alpha}{\Gamma(\alpha+1)},$$

que coincide exactamente con la solución analítica para esta función. Esa coincidencia es una prueba de fabricación útil, no demuestra el comportamiento para toda función.

Derivada de Riemann–Liouville mediante cuadratura

Para $0<q<1$,

$${}^{RL}D_{0+}^{q}x(t) =\frac{d}{dt}{}^{RL}I_{0+}^{1-q}x(t).$$

Una estrategia directa calcula primero la integral de orden $1-q$ con pesos de cuadratura y aproxima después la derivada exterior. Si

$$J_n\approx{}^{RL}I_{0+}^{1-q}x(t_n),$$

la diferencia regresiva produce

$${}^{RL}D_{0+}^{q}x(t_n) \approx\frac{J_n-J_{n-1}}{h}.$$

Este procedimiento hace visible la definición, pero combina error de cuadratura y error de diferenciación. Además, restar dos valores cercanos puede amplificar redondeo. Una fórmula de pesos obtenida al simplificar ambas operaciones puede ser más eficiente, pero debe comprobarse contra la construcción en dos etapas en una malla pequeña.

La derivada de Caputo también puede relacionarse con RL mediante la corrección inicial. Para $0<q<1$,

$${}^{C}D_{0+}^{q}x(t) = {}^{RL}D_{0+}^{q}[x(t)-x(0)].$$

Esta identidad proporciona una prueba de implementación: se calcula Caputo mediante incrementos y RL sobre la función corregida. Si no coinciden, se revisan el tratamiento de $x(0)$ y el primer intervalo.

Das compara fórmulas de cálculo GL, RL y Caputo y advierte que las condiciones en el punto inicial importan Das, 2020 Shantanu Das (2020) Kindergarten of Fractional Calculus Cambridge Scholars Publishing Ubicación consultada: cap. 4, §§4.15-4.18, pp. 183-197 .

Pérdida de orden por singularidades iniciales

Una solución de Mittag-Leffler puede ser continua y, al mismo tiempo, tener derivadas ordinarias poco regulares cerca de cero. Para una relajación,

$$x(t)=E_q(-\lambda t^q) =1-\frac{\lambda t^q}{\Gamma(q+1)}+\cdots.$$

Si $0<q<1$, la derivada ordinaria del término $t^q$ se comporta como $t^{q-1}$ y crece al acercarse al origen. Un análisis de convergencia que suponga derivadas ordinarias acotadas puede predecir un orden que la solución real no alcanza en una malla uniforme.

Una malla graduada conserva el intervalo $[0,T]$ y el número $N$ de subintervalos, pero coloca los nodos según

$$t_n=T\left(\frac{n}{N}\right)^\gamma, \qquad n=0,1,\ldots,N, \qquad \gamma\geq 1.$$

Con $\gamma=1$ se recupera la malla uniforme. Cuando $\gamma>1$, el paso

$$h_n=t_n-t_{n-1} =\frac{T}{N^\gamma}\left[n^\gamma-(n-1)^\gamma\right]$$

es pequeño cerca de $t=0$ y crece hacia $T$. Así se asignan más evaluaciones a la zona donde $t^{q-1}$ cambia con rapidez, sin aumentar $N$. Para $N=4$, la diferencia geométrica ya es visible:

$n$	uniforme: $t_n/T$, $\gamma=1$	graduada: $t_n/T$, $\gamma=2$
0	$0$	$0$
1	$0.25$	$0.0625$
2	$0.50$	$0.25$
3	$0.75$	$0.5625$
4	$1$	$1$

La comparación numérica mantiene fijos $q$, $\lambda$, $T$, $N$ y el método, y cambia únicamente $\gamma$. Para la relajación $x(t)=E_q(-\lambda t^q)$, los valores nodales de cada ejecución se interpolan linealmente en los mismos puntos $s_\ell=\ell T/K$. Allí se calculan

$$e_{\mathrm{global}}= \max_{0\leq \ell\leq K} \left|\widetilde{x}(s_\ell)-E_q(-\lambda s_\ell^q)\right|$$

y el mismo máximo restringido a $0\leq s_\ell\leq T/10$. Si la malla graduada reduce solo el error inicial, la tabla permite reconocer dónde se obtuvo la mejora; si aumenta el error lejos del origen, ese costo también queda visible. El valor de $\gamma$ y el número $K$ de puntos de comparación se fijan antes del experimento y se informan junto con $N$.

La singularidad también afecta la gráfica: omitir los primeros puntos puede ocultar el error dominante. Es mejor mostrar una ampliación del origen y explicar el criterio usado para incluir o excluir $t=0$.

Verificación mediante soluciones fabricadas

Probar únicamente el operador no detecta errores en el ensamblaje de una ecuación. Se elige una función exacta, por ejemplo

$$x_*(t)=t^\beta,$$

con $\beta$ compatible con Caputo. Para la ecuación

$${}^{C}D^qx(t)+\lambda x(t)=f(t),$$

se fabrica

$$f(t)= \frac{\Gamma(\beta+1)} {\Gamma(\beta-q+1)}t^{\beta-q} +\lambda t^\beta.$$

La solución exacta es $x_*$. La prueba recorre pesos, tratamiento del dato inicial, término lineal y fuerza. Se repite para varias $\beta$, incluyendo una función suave y otra con regularidad limitada. El forzamiento debe evaluarse con precisión suficiente para no convertirse en la principal fuente de error.

Una segunda prueba usa la relajación de Mittag-Leffler, donde $f=0$. Las dos se complementan: la solución fabricada aísla el orden algebraico de una potencia; Mittag-Leffler prueba una serie de potencias y una cola de memoria.

Coeficientes recursivos y cancelación

La recurrencia de los pesos evita evaluar cocientes de gamma grandes, pero una sucesión larga puede acumular redondeo. Se revisan tres invariantes prácticos:

1. $w_0=1$ y $w_1=-q$;
2. los primeros pesos coinciden con el binomio calculado a alta precisión;
3. las sumas parciales muestran el comportamiento esperado al crecer $n$.

Para órdenes cercanos a uno, muchos términos contribuyen mediante cancelación. Sumar en precisión insuficiente o en un orden desfavorable puede perder cifras. Una prueba con precisión extendida en una malla moderada permite saber si el error observado procede de la discretización o de la aritmética.

No se corrige una cancelación cambiando arbitrariamente el signo de pesos pequeños. Si se descartan por debajo de una tolerancia, el criterio debe medirse como una aproximación de memoria y compararse con la suma completa.

Paso variable y aproximación de memoria corta

Reducir $h$ mejora la resolución temporal; reducir el historial cambia el intervalo consultado. Un algoritmo puede usar paso pequeño y memoria insuficiente, o paso grueso e historia completa. Por eso el experimento separa dos ejes:

• refinamiento de $h$ con horizonte físico fijo;
• ampliación del horizonte de memoria con $h$ fijo.

El paso variable añade otra dificultad: los pesos uniformes ya no dependen solo de la edad expresada en número de celdas. Para verlo sin recurrir a una regla uniforme, se vuelve a la integral RL y se aproxima $f$ por su valor izquierdo $f_j=f(t_j)$ en cada celda. La contribución de $[t_j,t_{j+1}]$ al valor en $t_n$ es

$$\frac{f_j}{\Gamma(\alpha)} \int_{t_j}^{t_{j+1}}(t_n-\tau)^{\alpha-1}\,d\tau =a_{n,j}^{(\alpha)}f_j,$$

con el peso no uniforme

$$a_{n,j}^{(\alpha)}= \frac{(t_n-t_j)^\alpha-(t_n-t_{j+1})^\alpha} {\Gamma(\alpha+1)}.$$

Por tanto,

$${}^{RL}I_{0+}^{\alpha}f(t_n) \approx\sum_{j=0}^{n-1}a_{n,j}^{(\alpha)}f_j.$$

En una malla uniforme, $a_{n,j}^{(\alpha)}$ depende de $n-j$ y recupera los pesos $b_m^{(\alpha)}$ anteriores. En una malla no uniforme depende de los dos índices: cambiar un nodo modifica las distancias al punto de evaluación y obliga a recalcular las celdas afectadas.

Un cálculo corto comprueba la fórmula. Sean $\alpha=1/2$, $f(t)=1$ y los nodos $t_0=0$, $t_1=1/4$, $t_2=1$. En $t_2=1$,

$$a_{2,0}^{(1/2)}= \frac{1-\sqrt{3/4}}{\Gamma(3/2)}, \qquad a_{2,1}^{(1/2)}= \frac{\sqrt{3/4}}{\Gamma(3/2)}.$$

Los dos pesos suman $1/\Gamma(3/2)=2/\sqrt{\pi}$, exactamente el valor de ${}^{RL}I_{0+}^{1/2}1(1)$. La prueba también revela un error frecuente: sustituir ambos anchos por un único $h$ destruiría la telescopía.

Estudio de refinamiento

Una prueba de refinamiento mantiene fijos operador, terminal, orden, intervalo, datos iniciales y parámetros dimensionales. Se ejecutan mallas $h$, $h/2$ y $h/4$; en los nodos comunes se calcula

$$e_h=\max_n|x_h(t_n)-x_{ref}(t_n)|.$$

Si existe una referencia analítica, como una potencia o Mittag-Leffler, debe usarse. Si no existe, puede emplearse una solución de malla mucho más fina, indicando que es una referencia numérica. El orden observado entre dos refinamientos es

$$p_{obs}=\log_2\left(\frac{e_h}{e_{h/2}}\right).$$

El valor solo tiene sentido cuando los errores están en el régimen asintótico: mallas demasiado gruesas, redondeo, evaluación imprecisa de gamma o una historia truncada pueden deformarlo.

Costo de conservar toda la historia

En el paso $n$, una suma directa revisa aproximadamente $n$ valores. Para $N$ pasos, el número total de productos es

$$\sum_{n=1}^{N}n=\frac{N(N+1)}2,$$

de modo que el costo temporal crece como $O(N^2)$ y el almacenamiento de la trayectoria como $O(N)$. Esta cuenta se debe medir, no solo citar: registra tiempo, memoria y plataforma para varios $N$.

La suma GL con historia completa en el paso $n$ es

$$D_{\mathrm{full}}^q x_n =h^{-q}\sum_{k=0}^{n}w_k^{(q)}x_{n-k}.$$

Una memoria corta de $M$ pasos conserva solo los términos más recientes:

$$D_M^q x_n =h^{-q}\sum_{k=0}^{\min(n,M)}w_k^{(q)}x_{n-k}.$$

La ventana física es $L=Mh$. Si $M$ permanece fijo mientras $h$ disminuye, $L$ se acorta; para comparar refinamientos con el mismo horizonte físico se toma $M\approx L/h$. La diferencia frente al cálculo completo se mide, para cada paso, mediante

$$\varepsilon_{\mathrm{mem}}(n;M) =\left|D_{\mathrm{full}}^q x_n-D_M^q x_n\right|,$$

y puede resumirse con $\max_n\varepsilon_{\mathrm{mem}}(n;M)$. Si se necesita un error relativo, se divide ese máximo por $\max_n|D_{\mathrm{full}}^q x_n|$, siempre que el denominador no sea cero.

El efecto del corte se ve con cuatro datos. Para $q=1/2$, $h=1$ y $(x_0,x_1,x_2,x_3)=(1,2,3,4)$, el valor completo en $n=3$ es

$$D_{\mathrm{full}}^{1/2}x_3 =4-\frac12(3)-\frac18(2)-\frac1{16}(1) =2.1875.$$

ventana	términos retenidos	$D_M^{1/2}x_3$	$\varepsilon_{\mathrm{mem}}(3;M)$
$M=1$	2	$2.5000$	$0.3125$
$M=2$	3	$2.2500$	$0.0625$
$M=3$	4	$2.1875$	$0$

La tabla no promete que el error disminuya de forma monótona para cualquier señal, porque los pesos y los datos pueden cancelarse; por eso la comparación se repite en toda la trayectoria. Con historia completa, cada paso consulta un número creciente de estados y el costo acumulado es $O(N^2)$. Con la ventana, cada paso usa a lo sumo $M+1$ términos, de modo que el costo es $O(NM)$ y el historial de trabajo ocupa $O(M)$; guardar además toda la trayectoria de salida sigue requiriendo $O(N)$. Una elección de $M$ queda justificada por la tabla conjunta de $L$, error frente a la historia completa, tiempo y memoria.

Pruebas mínimas de discretización

Una implementación de GL o RL debe superar varios controles:

1. Funciones constantes y potencias. Comparar con fórmulas gamma y revisar el comportamiento en el terminal.
2. Límite entero. Recuperar diferencia o integral ordinaria al aproximarse al orden correspondiente.
3. Solución de Mittag-Leffler. Resolver una relajación lineal y comparar en un intervalo fijo.
4. Refinamiento. Separar error de paso, evaluación de funciones especiales e historia truncada.
5. Costo. Medir crecimiento temporal y memoria con el mismo código.
6. Reproducibilidad. Guardar convención de pesos, precisión, terminal y versión del algoritmo.

Herrmann resuelve un problema físico con la derivada de Riesz y hace explícitos el operador, las condiciones del problema y la construcción numérica Herrmann, 2014 Richard Herrmann (2014) Fractional Calculus: An Introduction for Physicists 2.ª ed. · World Scientific Ubicación consultada: cap. 12, §§12.1-12.2, pp. 161-175 . Esos elementos forman el contrato mínimo para que otra persona pueda reproducir una aplicación.

Convergencia, estabilidad y acotación responden preguntas diferentes. Una tabla de refinamiento estudia si el error disminuye al reducir el paso. Un análisis de estabilidad determina cómo evolucionan las perturbaciones del problema continuo o del esquema discreto. Una trayectoria acotada en un intervalo finito solo describe esa ejecución. El informe numérico debe nombrar la propiedad que realmente se midió y acompañarla con su prueba correspondiente.

Los pesos de Grünwald–Letnikov y las cuadraturas de Riemann–Liouville provienen de identidades analíticas distintas, aunque ambas produzcan sumas sobre el pasado Das, 2020 Shantanu Das (2020) Kindergarten of Fractional Calculus Cambridge Scholars Publishing Ubicación consultada: cap. 2, §2.23, pp. 96-99; cap. 3, §§3.23-3.24, pp. 146-149 . Una implementación debe separar la identidad continua, la discretización, el error medido y el eventual truncamiento de historia. Mezclar estos elementos puede hacer que un ahorro de memoria parezca, erróneamente, una propiedad del operador.

Ejercicios: pesos, error y tiempo de ejecución

1. Genera los primeros veinte pesos GL para tres órdenes y comprueba la recurrencia.
2. Repite la semiderivada de $t$ con $N=8$, $16$ y $32$.
3. Deriva la regla RL por tramos constantes y pruébala sobre $1$ y $t$.
4. Mide $p_{obs}$ con una referencia analítica evaluada a mayor precisión.
5. Resuelve la relajación de Mittag-Leffler con $\gamma=1$, $2$ y $3$. Mantén fijo $N$ y compara el error global, el error en $[0,T/10]$ y la distribución de nodos.
6. Reproduce los pesos no uniformes de la malla $(0,1/4,1)$; después cambia $f(t)=1$ por $f(t)=t$ y compara con la integral exacta.
7. Compara historia completa y una ventana física fija. Para al menos tres valores de $M$, reporta $L=Mh$, $\max_n\varepsilon_{\mathrm{mem}}(n;M)$ y el tiempo; no mantengas solo $M$ fijo al refinar.
8. Grafica tiempo contra $N$ para la historia completa y para dos ventanas. Contrasta el crecimiento observado con $N^2$ y $NM$.

Una trayectoria calculada aporta evidencia numérica cuando se contrasta con una identidad, una prueba de refinamiento y una medición de costo.