Año 2 · Módulo 14

135.1. Evolución temporal de una densidad

Las ecuaciones de Langevin especifican cómo evoluciona una realización del estado cuando se conocen la condición inicial y los incrementos aleatorios. Muchos problemas experimentales, sin embargo, preguntan por la probabilidad de encontrar el sistema en una región. La ecuación de Fokker–Planck estudia esa distribución sin seguir una por una todas las trayectorias. Su derivación combina la regla de Itô con un balance de probabilidad y deja explícito el papel de las fronteras.

Pensemos en partículas suspendidas en un líquido. La posición de una partícula se modifica por una fuerza sistemática y por choques microscópicos que no seguimos individualmente. Al observarla repetidas veces encontramos una trayectoria irregular. Si repetimos el experimento con numerosas partículas independientes, cada una con las mismas condiciones estadísticas iniciales, podremos construir un histograma de posiciones para cada instante. El histograma puede ensancharse, desplazarse o concentrarse aunque ninguna partícula represente por sí sola la evolución completa.

La ecuación de Fokker–Planck, también conocida en determinados contextos como ecuación de Kolmogórov hacia delante, describe cómo evoluciona una densidad de probabilidad bajo un modelo de difusión. Debe su nombre a los trabajos de Adriaan Fokker y Max Planck a comienzos del siglo XX; la formulación general se enmarca en el desarrollo posterior de los procesos de Markov. Su empleo se extiende desde física estadística y química hasta finanzas cuantitativas, ecología y neurociencia. En todos esos ámbitos, la ecuación adquiere significado solo después de definir qué variable representa el estado y qué hipótesis justifican tratar sus cambios como una difusión continua.

Consideraremos primero el modelo escalar de Itô

dXt=a(Xt,t) dt+b(Xt,t) dWt.dX_t=a(X_t,t)\,dt+b(X_t,t)\,dW_t.

El coeficiente aa es la deriva, que gobierna, bajo las condiciones técnicas apropiadas, el cambio condicional medio por unidad de tiempo. El coeficiente bb determina la amplitud de los incrementos aleatorios. WtW_t es un movimiento browniano estándar, cuyos incrementos en un intervalo corto de longitud hh tienen varianza hh. Salvo indicación contraria, los términos de ruido de este capítulo se interpretan en el sentido de Itô.

Si p(x,t)p(x,t) es la densidad de XtX_t, su evolución obedecerá

∂tp(x,t)=−∂x[a(x,t)p(x,t)]+12∂x2[b(x,t)2p(x,t)].\boxed{\partial_t p(x,t)=-\partial_x[a(x,t)p(x,t)] +\frac12\partial_x^2[b(x,t)^2p(x,t)].}

Las trayectorias impulsadas directamente por ruido browniano son, en general, no diferenciables; lo que puede ser suficientemente regular para satisfacer una ecuación diferencial parcial es la densidad que resulta de promediar sobre realizaciones. Tampoco existe siempre una densidad clásica, porque pueden presentarse estados discretos, masas puntuales o degeneraciones que obliguen a trabajar con medidas o soluciones débiles. Empezaremos con coeficientes regulares y difusión no degenerada para comprender el mecanismo antes de examinar sus límites.

A la izquierda, trayectorias individuales irregulares con una media que se aproxima al equilibrio; a la derecha, un histograma con la curva de densidad gaussiana exacta en el instante uno.
Figura 135.1. Veintiocho trayectorias generadas con transiciones exactas de Ornstein–Uhlenbeck y un histograma independiente de 28 000 estados en t=1. La línea continua del primer panel representa la media condicional; la curva del segundo es la densidad gaussiana exacta. Semilla 13501. Ampliar figura ↗

135.2. Forma débil de Fokker–Planck y ley de conservación

La derivación más instructiva comienza sin manipular directamente una densidad. Elegimos una función suave f(x)f(x) que represente una magnitud observable del estado. La fórmula de Itô establece que

df(Xt)=[a(Xt,t)f′(Xt)+12b(Xt,t)2f′′(Xt)]dt+b(Xt,t)f′(Xt) dWt.df(X_t)=\left[a(X_t,t)f'(X_t) +\frac12b(X_t,t)^2f''(X_t)\right]dt +b(X_t,t)f'(X_t)\,dW_t.

Bajo condiciones de integrabilidad y adaptación adecuadas, la integral estocástica del último término tiene esperanza cero. Por tanto,

ddtE[f(Xt)]=E ⁣[a(Xt,t)f′(Xt)+12b(Xt,t)2f′′(Xt)].\frac{d}{dt}\mathbb E[f(X_t)] =\mathbb E\!\left[a(X_t,t)f'(X_t) +\frac12b(X_t,t)^2f''(X_t)\right].

Si el proceso admite densidad, podemos escribir E[f(Xt)]=∫f(x)p(x,t) dx\mathbb E[f(X_t)]=\int f(x)p(x,t)\,dx. La expresión anterior se convierte en

∫f(x)∂tp(x,t) dx=∫[a(x,t)f′(x)+12b(x,t)2f′′(x)]p(x,t) dx.\int f(x)\partial_t p(x,t)\,dx =\int\left[a(x,t)f'(x)+\frac12b(x,t)^2f''(x)\right]p(x,t)\,dx.

Aplicando integración por partes una vez al término de deriva y dos veces al término de difusión, obtenemos

∫f ∂tp dx=∫f[−∂x(ap)+12∂x2(b2p)]dx,\int f\,\partial_t p\,dx =\int f\left[-\partial_x(ap)+\frac12\partial_x^2(b^2p)\right]dx,

si los términos de frontera desaparecen, por ejemplo porque ff tiene soporte compacto dentro del dominio. Como la igualdad se cumple para todas esas funciones de prueba, identificamos la ecuación de Fokker–Planck en sentido débil. Bajo hipótesis adicionales de regularidad, la misma identidad se interpreta punto por punto y obtenemos la ecuación diferencial parcial clásica.

La construcción aclara la relación entre dos operadores. Sobre las funciones observables actúa el generador infinitesimal hacia atrás,

Ltf=af′+12b2f′′.\mathcal L_t f=af'+\tfrac12b^2f''.

Sobre densidades actúa su adjunto formal

Lt∗p=−∂x(ap)+12∂x2(b2p).\mathcal L_t^*p=-\partial_x(ap)+\tfrac12\partial_x^2(b^2p).

Las derivadas aparecen en lugares diferentes en uno y otro operador. El primero indica cómo cambia la esperanza de una función del estado; el adjunto, cómo se redistribuye la probabilidad entre estados. Confundirlos conduce a ecuaciones incorrectas cuando los coeficientes dependen de la posición.

Existe además una conexión directa con C126. Allí una cadena de Markov en tiempo continuo obedecía una ecuación maestra que transfería masa probabilística entre estados mediante tasas de salto. Aquí el espacio de estados puede ser continuo y los pequeños incrementos se aproximan mediante un término local de transporte más otro de segundo orden. Ambas descripciones son ecuaciones hacia delante para una ley de probabilidad, aunque su estructura y sus hipótesis microscópicas sean diferentes. No toda ecuación maestra se reduce a Fokker–Planck; saltos grandes, distribuciones de cola pesada o mecanismos con memoria pueden requerir operadores no locales o ecuaciones de otro tipo.

Fokker–Planck puede leerse también como una ecuación de balance. Definimos el coeficiente de difusión

D(x,t)=12b(x,t)2D(x,t)=\frac12b(x,t)^2

y la corriente de probabilidad

J(x,t)=a(x,t)p(x,t)−∂x[D(x,t)p(x,t)].J(x,t)=a(x,t)p(x,t)-\partial_x[D(x,t)p(x,t)].

La ecuación se escribe entonces

∂tp+∂xJ=0.\boxed{\partial_t p+\partial_xJ=0.}

Esta forma indica que, dentro del dominio, la probabilidad solo se transporta entre regiones. Una corriente positiva indica movimiento neto de probabilidad hacia la derecha, aunque haya realizaciones individuales que se mueven en ambos sentidos. Del mismo modo, una densidad local constante es compatible con movimiento, porque la probabilidad que entra en una región puede compensar exactamente la que sale.

Integrando la ecuación sobre un intervalo [ℓ,r][\ell,r] obtenemos

ddt∫ℓrp(x,t) dx=J(ℓ,t)−J(r,t).\frac{d}{dt}\int_{\ell}^{r}p(x,t)\,dx=J(\ell,t)-J(r,t).

La variación de probabilidad dentro del intervalo depende solo del flujo en sus extremos. Este resultado es útil para detectar errores en un modelo numérico. Si el sistema tiene fronteras impermeables y la probabilidad total aumenta o disminuye de forma apreciable, el algoritmo ha violado una propiedad estructural. En cambio, si existe una frontera absorbente que retira trayectorias, la masa de la densidad de los supervivientes puede disminuir legítimamente.

Conviene detenerse en una peculiaridad del ruido multiplicativo. Si DD depende de xx, entonces

J=ap−D ∂xp−(∂xD)p.J=ap-D\,\partial_xp-(\partial_xD)p.

El término adicional con ∂xD\partial_xD aparece aunque la deriva de Itô sea aa. Por eso identificar automáticamente apap con transporte determinista y −D∂xp-D\partial_xp con difusión puede ser incompleto. Esa separación sencilla solo es válida para un coeficiente de difusión constante; en el caso general hay que derivar el producto DpDp.

Tampoco sería correcto corregir esta fórmula con intuiciones procedentes de Stratonovich sin revisar qué ecuación se está describiendo. Si empezamos con dXt=α(Xt)dt+b(Xt)∘dWtdX_t=\alpha(X_t)dt+b(X_t)\circ dW_t, primero convertimos a Itô usando a=α+12bb′a=\alpha+\tfrac12bb' en una dimensión. Después aplicamos la ecuación hacia delante. El proceso físico y su densidad no cambian por esa conversión correcta, pero la letra que ocupa el lugar de la deriva sí cambia. Esta es una razón para mantener explícita la convención de ruido en artículos, modelos y código.

Finalmente, la corriente tiene unidades verificables. Si xx tiene unidades de longitud, pp tiene unidades de longitud inversa, aa de longitud por tiempo y DD de longitud cuadrada por tiempo. En consecuencia, JJ tiene unidades de tiempo inverso. Una interpretación física que dé dimensiones incompatibles señala un error más básico que la elección del método numérico.

Densidad con intervalo sombreado y corriente superpuesta; líneas verticales marcan dos fronteras donde la diferencia de flujos determina el cambio de probabilidad interior.
Figura 135.2. Densidad gaussiana y corriente de una difusión con deriva constante sobre la recta. El área sombreada muestra la masa en un intervalo; su variación temporal coincide exactamente con la corriente que entra menos la que sale por sus fronteras. Curvas y tasas calculadas analíticamente. Ampliar figura ↗

135.3. Difusión libre y transporte de probabilidad

El caso más sencillo elimina la deriva y conserva una amplitud de ruido constante

dXt=2D dWt,D>0.dX_t=\sqrt{2D}\,dW_t,\qquad D>0.

Entonces Fokker–Planck se reduce a la ecuación del calor

∂tp=D∂x2p.\partial_t p=D\partial_x^2p.

Si todas las partículas comienzan exactamente en x0x_0, la distribución inicial es una masa puntual, que representamos formalmente con una delta de Dirac. Para tiempos positivos, la solución fundamental es

p(x,t∣x0)=14πDtexp⁡ ⁣[−(x−x0)24Dt].p(x,t\mid x_0)=\frac{1}{\sqrt{4\pi Dt}} \exp\!\left[-\frac{(x-x_0)^2}{4Dt}\right].

Se trata de una densidad normal con media x0x_0 y varianza 2Dt2Dt. Su desviación estándar crece como 2Dt\sqrt{2Dt}, resultado coherente con el escalamiento de las caminatas aleatorias de C129. Ese crecimiento con la raíz cuadrada del tiempo no es universal en todo movimiento irregular; depende del régimen difusivo, del carácter browniano de los incrementos y de la ausencia de mecanismos capaces de producir escalamiento anómalo.

La solución permite separar dos preguntas. Si se observa una partícula, su desplazamiento puede ser positivo o negativo. Si se observan miles, el centro del conjunto permanece en x0x_0 cuando no existe deriva, mientras la nube se ensancha. No hay contradicción, porque el promedio no describe el movimiento de cada integrante. En un modelo continuo regular, la probabilidad de encontrar la partícula en un punto exacto es cero; para hacer predicciones debemos integrar la densidad sobre intervalos observables.

Supongamos que D=0,4D=0{,}4 unidades cuadradas por segundo y que partimos de x0=0x_0=0. Tras cinco segundos, la varianza es 44 unidades cuadradas y la desviación estándar vale 22 unidades. Aproximadamente el 68 % de la probabilidad se encuentra entre −2-2 y 22, puesto que la distribución es normal. Si duplicamos el tiempo, la varianza se duplica y la desviación estándar solo aumenta en un factor 2\sqrt2.

Para una condición inicial distribuida, la densidad posterior se obtiene por convolución con la solución fundamental

p(x,t)=∫Rexp⁡[−(x−y)2/(4Dt)]4πDtp0(y) dy.p(x,t)=\int_{\mathbb R}\frac{\exp[-(x-y)^2/(4Dt)]} {\sqrt{4\pi Dt}}p_0(y)\,dy.

Cada punto de la distribución inicial contribuye con una pequeña campana gaussiana a la densidad final. La difusión suaviza los contrastes espaciales, pero sin confinamiento la masa sigue extendiéndose y no existe una distribución estacionaria normalizable sobre R\mathbb R. La forma y la existencia de una distribución estacionaria dependerán del dominio, las fuerzas y las fronteras, cuestión central en C136.

Incorporamos ahora una velocidad de deriva constante vv,

dXt=v dt+2D dWt.dX_t=v\,dt+\sqrt{2D}\,dW_t.

La ecuación para la densidad queda

∂tp=−v∂xp+D∂x2p.\partial_t p=-v\partial_xp+D\partial_x^2p.

Para una posición inicial precisa, la solución vuelve a ser gaussiana, con un centro que se desplaza,

p(x,t∣x0)=14πDtexp⁡ ⁣[−(x−x0−vt)24Dt].p(x,t\mid x_0)=\frac{1}{\sqrt{4\pi Dt}} \exp\!\left[-\frac{(x-x_0-vt)^2}{4Dt}\right].

La media es x0+vtx_0+vt y la varianza, 2Dt2Dt; transporte y dispersión dependen de parámetros distintos. En datos experimentales, atribuir el desplazamiento medio a la difusión, o el ensanchamiento a la deriva, puede producir estimaciones incorrectas de fuerzas o velocidades de transporte.

La corriente es

J=vp−D∂xp.J=vp-D\partial_xp.

El término vpvp transporta el perfil en bloque. El término −D∂xp-D\partial_xp tiende a trasladar probabilidad desde regiones más densas hacia regiones menos densas. Ambas contribuciones se suman, aunque no siempre apuntan en la misma dirección; en la parte posterior de una nube en movimiento, la difusión puede actuar contra el desplazamiento promedio.

Una magnitud adimensional útil es el número de Péclet, que compara transporte dirigido y difusión en una escala espacial LL,

Pe=∣v∣LD.\mathrm{Pe}=\frac{|v|L}{D}.

Cuando Pe\mathrm{Pe} es grande, el desplazamiento por deriva domina frente a la dispersión difusiva sobre esa escala. Cuando es pequeño, los gradientes producidos por difusión tienen mayor peso relativo. La interpretación exige declarar LL, porque cambiar la longitud observada cambia el número de Péclet aunque el sistema microscópico sea el mismo. En sistemas multiescala, un mismo mecanismo puede parecer dominado por transporte a una escala y por difusión a otra.

En un intervalo con fronteras, la solución gaussiana sobre la recta deja de ser suficiente. Una frontera reflectante redistribuye trayectorias que intentarían salir, mientras una frontera absorbente retira del conjunto a quienes llegan a ella. La distribución debe satisfacer condiciones matemáticas compatibles con ese mecanismo. Un ejemplo clínico o logístico de tiempos de espera puede requerir además estados de llegada y salida discretos; observar fluctuaciones no justifica sustituirlos por una ecuación de difusión.

Una ecuación como la anterior es, por tanto, una hipótesis dinámica precisa. Separa desplazamiento y ensanchamiento, especifica la ley de incrementos, permite predecir distribuciones futuras y expone qué medidas serían pertinentes para estimar sus parámetros. Que una nube de datos sea aproximadamente gaussiana, en cambio, no prueba que el sistema siga una difusión browniana; inferir el mecanismo requiere datos temporales y comparación con alternativas.

Dos paneles de distribuciones gaussianas en tres instantes. Las curvas sin deriva se ensanchan sin cambiar de centro; con deriva, su centro avanza a la derecha.
Figura 135.3. Núcleos gaussianos analíticos a tres tiempos para difusión libre y difusión con deriva constante. En ambos casos, la varianza es 2Dt. La deriva cambia la media en vt, mientras que la difusión determina el ensanchamiento. Ampliar figura ↗

135.4. El proceso de Ornstein–Uhlenbeck como equilibrio dinámico

La difusión libre se dispersa indefinidamente. Para comprender una distribución que se aproxima a un perfil estable necesitamos introducir una fuerza restauradora. El ejemplo canónico es el proceso de Ornstein–Uhlenbeck, surgido del estudio del movimiento browniano y empleado después en modelos de fluctuaciones alrededor de un estado de referencia. Escribimos

dXt=−κ(Xt−μ) dt+2D dWt,κ,D>0.dX_t=-\kappa(X_t-\mu)\,dt+\sqrt{2D}\,dW_t, \qquad \kappa,D>0.

El parámetro μ\mu determina el centro hacia el cual apunta la deriva. El parámetro κ\kappa gobierna la rapidez de restauración y DD, la intensidad de difusión. La ecuación de Fokker–Planck correspondiente es

∂tp=κ∂x[(x−μ)p]+D∂x2p.\partial_t p=\kappa\partial_x[(x-\mu)p] +D\partial_x^2p.

La fuerza restauradora compite con el ensanchamiento producido por las fluctuaciones. Si comenzamos en un valor preciso x0x_0, el proceso sigue una distribución normal cuya media y varianza son

m(t)=μ+(x0−μ)e−κt,m(t)=\mu+(x_0-\mu)e^{-\kappa t},
s2(t)=Dκ(1−e−2κt).s^2(t)=\frac{D}{\kappa}(1-e^{-2\kappa t}).

En tiempos muy cortos, como e−2κt≈1−2κte^{-2\kappa t}\approx1-2\kappa t, la varianza crece aproximadamente como 2Dt2Dt, igual que en difusión libre. A tiempos largos, el ensanchamiento se detiene y la varianza se aproxima a D/κD/\kappa. Incluso en la situación estacionaria, una trayectoria individual sigue variando; lo que se estabiliza es la ley de probabilidad.

Podemos obtener las mismas ecuaciones desde el generador, sin resolver primero la ecuación diferencial parcial. Para f(x)=xf(x)=x, tenemos f′=1f'=1 y f′′=0f''=0, por lo que

ddtE[Xt]=−κ(E[Xt]−μ).\frac{d}{dt}\mathbb E[X_t]=-\kappa(\mathbb E[X_t]-\mu).

Para f(x)=x2f(x)=x^2, la fórmula de Itô da una ecuación de segundo momento. Combinando ambas obtenemos

ddtVar⁡(Xt)=−2κVar⁡(Xt)+2D.\frac{d}{dt}\operatorname{Var}(X_t) =-2\kappa\operatorname{Var}(X_t)+2D.

Si la varianza inicial es v0v_0, la solución general es v(t)=D/κ+(v0−D/κ)e−2κtv(t)=D/\kappa+(v_0-D/\kappa)e^{-2\kappa t}. Para una condición inicial puntual, v0=0v_0=0 y recuperamos la expresión anterior. Los momentos pueden calcularse así directamente desde el generador, aun cuando la densidad completa sea difícil de obtener.

El carácter gaussiano del ejemplo se debe a que la ecuación estocástica es lineal y el ruido es gaussiano. Si el estado inicial tiene una distribución arbitraria, su distribución a un tiempo finito es la convolución de una transformación lineal de la inicial con un ruido gaussiano independiente. No tiene por qué ser exactamente gaussiana de inmediato, aunque converge hacia la distribución estacionaria normal en las condiciones habituales del modelo. Que las ecuaciones de momentos se cierren no prueba, en general, la gaussianidad.

La escala característica de relajación de la media es 1/κ1/\kappa, mientras la diferencia de varianza frente al valor estacionario decae con la escala 1/(2κ)1/(2\kappa); la relajación no tiene una única escala temporal. Esta distinción reaparecerá al estudiar corrientes, sistemas multivariados y aproximaciones al equilibrio.

La media se aproxima al centro mu, la banda de dispersión se estabiliza y las distribuciones a tiempos crecientes se acercan a una gaussiana estacionaria.
Figura 135.4. Media condicional, banda de una desviación estándar y densidades gaussianas exactas a distintos tiempos para Ornstein–Uhlenbeck. La distribución converge hacia N(μ,D/κ); la banda expresa dispersión de estados, no un intervalo de confianza de un estimador. Ampliar figura ↗

135.5. Cuando hay varias variables y cuando el ruido depende del estado

La mayor parte de los sistemas complejos requiere varias variables de estado. Sea Xt∈RdX_t\in\mathbb R^d, con deriva vectorial a(x,t)a(x,t) y una matriz B(x,t)B(x,t) que multiplica un movimiento browniano de mm componentes independientes. La ecuación de Itô se escribe

dXt=a(Xt,t) dt+B(Xt,t) dWt.dX_t=a(X_t,t)\,dt+B(X_t,t)\,dW_t.

Definimos la matriz de covarianza infinitesimal A=BBTA=BB^{\mathsf T}. Sus componentes indican qué variables fluctúan conjuntamente por efecto de los mismos canales de ruido. La ecuación hacia delante es

∂tp=−∑i=1d∂xi(aip)+12∑i,j=1d∂xi∂xj(Aijp).\boxed{\partial_t p=-\sum_{i=1}^{d}\partial_{x_i}(a_ip) +\frac12\sum_{i,j=1}^{d}\partial_{x_i}\partial_{x_j}(A_{ij}p).}

Los términos mixtos son relevantes cuando hay correlación instantánea entre incrementos de distintas componentes. Por ejemplo, dos recursos que reciben perturbaciones de una fuente ambiental común pueden fluctuar de manera correlacionada. Sustituir AA por una matriz diagonal en ese caso modifica la ley del proceso. La simetría y semidefinición positiva de AA son condiciones estructurales, porque la matriz debe describir una covarianza posible.

En varias dimensiones podemos definir una corriente vectorial cuyos componentes satisfacen

Ji=aip−12∑j∂xj(Aijp),∂tp=−∇⋅J.J_i=a_ip-\frac12\sum_j\partial_{x_j}(A_{ij}p), \qquad \partial_t p=-\nabla\cdot J.

Lo que cambia la densidad es la divergencia de la corriente; un flujo puede circular localmente sin acumular ni vaciar probabilidad en ningún punto. La estacionariedad requiere ∇⋅J=0\nabla\cdot J=0, mientras que la ausencia de corriente exige la condición más fuerte J=0J=0. En una dimensión, dentro de un intervalo sin fuentes, estacionariedad significa una corriente espacialmente constante; en varias dimensiones pueden existir campos de corriente cerrados y estructuras de circulación. Ese será el eje de C136.

El coeficiente de ruido dependiente del estado introduce otra dificultad. Consideremos el proceso geométrico dXt=rXtdt+σXtdWtdX_t=rX_tdt+\sigma X_tdW_t, iniciado en X0>0X_0>0. Su ecuación de Fokker–Planck es

∂tp=−∂x(rxp)+σ22∂x2(x2p),x>0.\partial_t p=-\partial_x(rxp) +\frac{\sigma^2}{2}\partial_x^2(x^2p),\qquad x>0.

La difusión local es proporcional a x2x^2, de modo que las regiones con valores grandes tienen fluctuaciones absolutas más intensas. La solución se obtiene transformando Yt=log⁡XtY_t=\log X_t mediante Itô, lo que da dYt=(r−σ2/2)dt+σdWtdY_t=(r-\sigma^2/2)dt+\sigma dW_t. Por eso XtX_t sigue una distribución lognormal en la variable original. Una trayectoria que parte de un valor positivo permanece, en este modelo, positiva en todo tiempo finito.

Este ejemplo muestra que observar la forma de una densidad requiere elegir correctamente la variable. Una distribución aproximadamente normal después de aplicar logaritmos puede corresponder a una distribución claramente asimétrica en la escala original. El cambio de variable transforma tanto la deriva como la difusión y altera la densidad mediante el jacobiano correspondiente. El modelo matemático debe respetar esas transformaciones.

La existencia de coeficientes dependientes del estado obliga, por último, a verificar regularidad, crecimiento, conservación del dominio y tratamiento de posibles degeneraciones. La ecuación de Fokker–Planck puede seguir siendo válida en sentido débil, aunque una densidad suave y única no esté garantizada en todos los casos. Las afirmaciones sobre positividad o sobre la accesibilidad de una frontera requieren analizar el proceso concreto.

135.6. Las fronteras forman parte del mecanismo

Al derivar una ecuación hacia delante sobre toda la recta solemos ocultar los términos de frontera. Un modelo real rara vez dispone de un dominio ilimitado sin restricciones. Una concentración no puede tomar valores negativos; la proporción de una población está acotada entre cero y uno; una partícula puede encontrarse dentro de una caja. La forma diferencial de la ecuación no dice qué ocurre al llegar a un límite, así que debemos especificarlo y traducirlo a una condición sobre densidad o corriente.

Una frontera reflectante es un límite que no deja escapar probabilidad. Su condición natural es J⋅n=0J\cdot n=0, donde nn es la normal exterior. En una dimensión se reduce a J=0J=0 en cada extremo reflectante. Si DD es constante y no hay deriva normal, la condición se expresa como ∂xp=0\partial_xp=0; con deriva, la condición correcta es ap−D∂xp=0ap-D\partial_xp=0. Reemplazar cualquier frontera reflectante por derivada espacial nula puede violar la conservación.

Una frontera absorbente retira del conjunto activo las trayectorias que llegan a ella. En difusiones regulares no degeneradas con absorción por primera llegada, suele imponerse p=0p=0 en la frontera para la densidad de trayectorias todavía no absorbidas. La masa total interior disminuye porque sale una corriente hacia el límite. El flujo de salida permite calcular distribuciones de tiempos de llegada, siempre que se atienda a las convenciones de signo y a las características del dominio. En una ecuación cuyo estado de frontera absorbido se incorpora como masa puntual separada, la probabilidad total del proceso completo sí permanece normalizada; disminuye únicamente la masa que permanece en el interior.

Las fronteras periódicas identifican extremos opuestos de un intervalo, como una posición angular sobre una circunferencia. Entonces densidad y corriente han de concordar al cerrar el dominio; una corriente neta puede circular indefinidamente sin acumulación local. Las fronteras llamadas naturales o inaccesibles exigen otro tratamiento; no se comportan necesariamente como paredes físicas, porque algunas trayectorias no pueden alcanzarlas en tiempo finito según los coeficientes del modelo.

A partir del balance integrado sobre el dominio Ω\Omega tenemos

ddt∫Ωp dx=−∫∂ΩJ⋅n dS.\frac{d}{dt}\int_{\Omega}p\,dx =-\int_{\partial\Omega}J\cdot n\,dS.

Esta igualdad es una prueba de consistencia previa a cualquier solución. Si un algoritmo devuelve masa creciente en un sistema pasivo con única frontera absorbente, probablemente hay una incompatibilidad entre el flujo, el signo de la normal y la actualización numérica. Si una condición reflectante deja escapar probabilidad, la implementación no representa el mecanismo anunciado. Verificar estas propiedades es más informativo que obtener una curva aparentemente suave.

También conviene separar tres objetos: la distribución de XtX_t sobre todos los estados, la distribución restringida a quienes no fueron absorbidos y la distribución condicionada a la supervivencia. La segunda es una subprobabilidad y no integra necesariamente uno; la tercera se normaliza dividiendo por la probabilidad de supervivencia. Mezclar esos objetos lleva a conclusiones falsas acerca de estacionariedad, supervivencia o riesgo acumulado. Volveremos sobre ellos al estudiar escape y cuasiestacionariedad.

Tres perfiles de densidad: reflexión y periodicidad conservan masa uno, mientras la absorción presenta masa interior menor que uno y densidad nula en ambos extremos.
Figura 135.5. Soluciones analíticas por el método de imágenes para difusión en [0,1] con frontera reflectante, absorbente y periódica. Las tres parten de la misma posición y usan el mismo D y tiempo. Solo la condición absorbente pierde masa interior; la curva no se renormaliza. Ampliar figura ↗

135.7. Experimento reproducible con una densidad que se forma en el tiempo

Podemos verificar Fokker–Planck sin programar de inmediato un solucionador de ecuaciones diferenciales parciales. Usaremos la solución exacta del proceso de Ornstein–Uhlenbeck y generaremos muchas realizaciones independientes en varios instantes. Después compararemos cantidades derivadas matemáticamente con estimaciones Monte Carlo y calcularemos cuánto difiere un histograma de su densidad teórica integrada por intervalos.

Tomaremos κ=1,3\kappa=1{,}3, μ=0,7\mu=0{,}7, D=0,6D=0{,}6, X0=−1X_0=-1 y un conjunto de tiempos positivos. Para cada tiempo podemos simular directamente la distribución exacta, puesto que

Xt=m(t)+s(t)Z,Z∼N(0,1).X_t=m(t)+s(t)Z,\qquad Z\sim N(0,1).

Esto evita confundir un error de discretización temporal con uno de muestreo. La aproximación mediante pasos de Euler–Maruyama también podría usarse, pero incorporaría un sesgo dependiente de dtdt que convendría analizar por separado. La comprobación de este capítulo se limita a la propagación de la ley teórica y al comportamiento del histograma bajo muestras finitas.

El siguiente código requiere NumPy; emplea solamente su generador seudoaleatorio y operaciones vectorizadas. Para las probabilidades por intervalos usa la función de distribución normal a través de math.erf, sin depender de bibliotecas adicionales. La prueba numérica tiene tolerancias explícitas, elegidas para un tamaño de muestra suficientemente grande y una semilla reproducible.

import math
import numpy as np

rng = np.random.default_rng(1352026)
kappa, mu, D, x0 = 1.3, 0.7, 0.6, -1.0
n = 100_000

def cdf_normal(x, mean, sd):
    z = (x - mean) / (sd * math.sqrt(2.0))
    return 0.5 * (1.0 + math.erf(z))

for t in (0.25, 1.0, 2.5):
    mean = mu + (x0 - mu) * math.exp(-kappa * t)
    var = (D / kappa) * (1.0 - math.exp(-2.0 * kappa * t))
    sd = math.sqrt(var)
    sample = mean + sd * rng.standard_normal(n)

    # Todos los intervalos cubren la recta: se agregan colas infinitas.
    edges = np.r_[[-np.inf], np.linspace(mean-4*sd, mean+4*sd, 61),
                  [np.inf]]
    counts, _ = np.histogram(sample, bins=edges)
    empirical = counts / n
    expected = np.array([
        cdf_normal(float(edges[j+1]), mean, sd)
        - cdf_normal(float(edges[j]), mean, sd)
        for j in range(len(edges)-1)
    ])
    l1 = np.abs(empirical - expected).sum()
    assert abs(sample.mean() - mean) < 0.015
    assert abs(sample.var() - var) < 0.018
    assert abs(empirical.sum() - 1.0) < 1e-12
    assert l1 < 0.045
    print(f"t={t:4.2f} media={sample.mean():.4f} "
          f"varianza={sample.var():.4f} distancia_L1={l1:.4f}")

La distancia L1L^1 discreta calculada aquí es la suma de las diferencias absolutas entre las probabilidades empíricas y teóricas de cada intervalo, una medida distinta del error puntual de la densidad. Con intervalos demasiado estrechos, el histograma puede parecer más irregular aun cuando las masas por intervalo sean coherentes. Una comparación correcta debe indicar el tamaño de la muestra, los límites de los intervalos y qué propiedad se pretende contrastar.

La comprobación muestra que el código genera observaciones de la distribución analítica anunciada y que los estadísticos elegidos son consistentes con ella en la muestra realizada. La validez física del proceso de Ornstein–Uhlenbeck es otra cuestión; para contrastarlo con un sistema empírico habría que comparar autocorrelaciones temporales, escalas de relajación, distribución de incrementos y estabilidad de los parámetros entre contextos. Una coincidencia entre histogramas aislados no identifica una dinámica única.

135.8. Qué informa una solución para la densidad

La ecuación de Fokker–Planck convierte una descripción por trayectorias en una descripción por distribuciones. Hemos partido de una ecuación de Itô, aplicado la fórmula de cambio de variables a funciones observables y obtenido el operador adjunto que gobierna la densidad. La forma de conservación permitió interpretar la corriente como transporte neto de probabilidad, mientras los ejemplos mostraron dos comportamientos distintos: difusión sin confinamiento y difusión sometida a una deriva restauradora.

Esta traducción descansa en varios supuestos. Primero, la densidad de un instante no determina las trayectorias que la produjeron. Dos modelos pueden compartir la misma distribución marginal y diferir en correlación temporal, circulación de probabilidades o respuesta a intervenciones. La inferencia dinámica exige observaciones temporales y una estructura de estados que conserve la información relevante. La misma precaución vale para poblaciones heterogéneas, en las que el promedio de comportamientos individuales puede ocultar mecanismos distintos.

Segundo, no siempre está justificado sustituir un proceso de saltos por una difusión. La aproximación puede funcionar cuando abundan incrementos pequeños y el escalamiento permite retener esencialmente los primeros momentos locales pertinentes. Puede fallar con acontecimientos grandes y raros, sistemas de tamaño reducido o dinámicas no markovianas. Incluso cuando un límite difusivo es matemáticamente correcto, sus predicciones para colas extremas y tiempos de primera llegada pueden requerir una verificación independiente. Conviene juzgar la aproximación por las cantidades que necesitamos estimar, más allá de su ajuste a un histograma central.

Tercero, el movimiento browniano y la ecuación diferencial parcial están ligados mediante un convenio preciso de coeficientes. En nuestro modelo, D=b2/2D=b^2/2. Si el ruido se expresa como 2D dWt\sqrt{2D}\,dW_t, el término de difusión es D∂x2pD\partial_x^2p cuando DD es constante. Perder el factor dos produce errores de varianza, escalas de tiempo y distribuciones estacionarias. Con ruido multiplicativo, además, hay que derivar el producto DpDp completo. La convención Itô–Stratonovich debe quedar resuelta antes de deducir el operador hacia delante.

Cuarto, la elección de frontera modifica la predicción. El mismo operador local puede describir un proceso que refleja, uno que se absorbe o uno que recorre periódicamente un dominio. Ninguna gráfica obtenida en el interior del intervalo autoriza a omitir esas condiciones. Una probabilidad que abandona el dominio puede representar una pérdida real, una transición a otro estado o una restricción artificial del análisis. Debemos declarar cuál de estas situaciones está modelada.

La ecuación de Fokker–Planck convierte la información local de deriva y difusión en una ley de evolución de densidades, sujeta a condiciones iniciales y de frontera. Su forma conservativa permite comprobar balances de masa, mientras que la solución gaussiana de Ornstein–Uhlenbeck es un caso exacto para verificar integradores y estimadores de densidad.

Los acontecimientos que dependen del recorrido, como el máximo de una trayectoria o su primer tiempo de llegada, no se recuperan en general a partir de una colección de marginales instantáneas. Cuando existe absorción, la densidad interior pierde masa; al condicionarla en supervivencia se obtiene una distribución diferente. Precisar cuál de esos objetos se presenta evita atribuir a un gráfico probabilidades que no contiene.

Referencias

  • Gardiner, C. W. (2009). Stochastic Methods: A Handbook for the Natural and Social Sciences (4.ª ed.). Springer. Generadores, ecuaciones hacia delante y aproximaciones de difusión.
  • Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications (6.ª ed.). Springer. Fórmula de Itô y procesos de difusión.
  • Pavliotis, G. A. (2014). Stochastic Processes and Applications: Diffusion Processes, the Fokker–Planck and Langevin Equations. Springer. DOI: 10.1007/978-1-4939-1323-7.
  • Risken, H. (1989). The Fokker–Planck Equation: Methods of Solution and Applications (2.ª ed.). Springer. Soluciones, corrientes y métodos de análisis.
  • van Kampen, N. G. (2007). Stochastic Processes in Physics and Chemistry (3.ª ed.). Elsevier. Ecuaciones maestras y aproximaciones continuas.