Integrador exponencial

Los integradores exponenciales son una clase de métodos numéricos para la resolución de ecuaciones diferenciales ordinarias , específicamente problemas de valor inicial . Esta amplia clase de métodos del análisis numérico se basa en la integración exacta de la parte lineal del problema de valor inicial. Dado que la parte lineal se integra exactamente, esto puede ayudar a mitigar la rigidez de una ecuación diferencial. Los integradores exponenciales pueden construirse de forma explícita o implícita para ecuaciones diferenciales ordinarias numéricas , o bien servir como integrador temporal para ecuaciones diferenciales parciales numéricas .

Estos métodos, que datan al menos de la década de 1960, fueron reconocidos por Certaine ^{[ 1 ]} y Pope. ^{[ 2 ]} Últimamente, los integradores exponenciales se han convertido en un área activa de investigación, véase Hochbruck y Ostermann (2010). ^{[ 3 ]} Desarrollados originalmente para resolver ecuaciones diferenciales rígidas , los métodos se han utilizado para resolver ecuaciones diferenciales parciales , incluyendo problemas hiperbólicos y parabólicos ^{[ 4 ]} como la ecuación del calor .

Consideramos problemas de valor inicial de la forma,

dónde${\displaystyle L}$está compuesto por términos lineales y${\displaystyle N}$está compuesto por términos no lineales . Estos problemas pueden surgir de un problema de valor inicial más típico. $${\displaystyle u'(t)=f(u(t)),\qquad u(t_{0})=u_{0},}$$ después de linealizar localmente alrededor de un estado fijo o local${\displaystyle u^{*}}$:$${\displaystyle L={\frac {\partial f}{\partial u}}(u^{*});\qquad N=f(u)-Lu.}$$ Aquí,${\displaystyle {\frac {\partial f}{\partial u}}}$se refiere a la derivada parcial de${\displaystyle f}$con respecto a${\displaystyle u}$(el jacobiano de f ).

Integración exacta de este problema desde el tiempo 0 hasta un tiempo posterior.${\displaystyle t}$se puede realizar utilizando exponenciales matriciales para definir una ecuación integral para la solución exacta: ^{[ 3 ]}

Esto es similar a la integral exacta utilizada en el teorema de Picard-Lindelöf . En el caso de${\displaystyle N\equiv 0}$, esta formulación es la solución exacta de la ecuación diferencial lineal .

Los métodos numéricos requieren una discretización de la ecuación (2). Pueden basarse en discretizaciones de Runge-Kutta , ^{[ 5 ] [ 6 ] [ 7 ]} métodos lineales multipaso o una variedad de otras opciones.

Se ha demostrado que los métodos exponenciales de Rosenbrock son muy eficientes para resolver grandes sistemas de ecuaciones diferenciales ordinarias rígidas, generalmente resultantes de la discretización espacial de EDP dependientes del tiempo (parabólicas). Estos integradores se construyen a partir de una linealización continua de (1) a lo largo de la solución numérica.${\displaystyle u_{n}}$

dónde${\displaystyle L_{n}={\frac {\partial f}{\partial u}}(u_{n})}$,${\ Displaystyle N_ {n} (u) = f (u) -L_ {n} u}$Este procedimiento goza de la ventaja, en cada paso, de que${\displaystyle {\frac {\partial N_{n}}{\partial u}}(u_{n})=0.}$ Esto simplifica considerablemente la derivación de las condiciones de orden y mejora la estabilidad al integrar la no linealidad.${\displaystyle N(u(t))}$Nuevamente, al aplicar la fórmula de variación de constantes (2) se obtiene la solución exacta en el instante${\displaystyle t_{n+1}}$como

La idea ahora es aproximar la integral en (4) mediante alguna regla de cuadratura con nodos${\displaystyle c_{i}}$y pesos${\displaystyle b_{i}(h_{n}L_{n})}$(${\displaystyle 1\leq i\leq s}$). Esto produce la siguiente clase de${\displaystyle s}$-Métodos de Rosenbrock exponenciales explícitos en etapas, véanse Hochbruck y Ostermann (2006), Hochbruck, Ostermann y Schweitzer (2009): $${\displaystyle {\begin{aligned}U_{ni}&=e^{c_{i}h_{n}L_{n}}u_{n}+h_{n}\sum _{j=1}^{i-1}a_{ij}(h_{n}L_{n})N_{n}(U_{nj}),\\u_{n+1}&=e^{h_{n}L_{n}}u_{n}+h_{n}\sum _{i=1}^{s}b_{i}(h_{n}L_{n})N_{n}(U_{ni})\end{aligned}}}$$ con${\displaystyle u_{n}\approx u(t_{n})}$,${\displaystyle U_{ni}\approx u(t_{n}+c_{i}h_{n})}$,${\displaystyle h_{n}=t_{n+1}-t_{n}}$. Los coeficientes${\ Displaystyle a_ {ij} (z), b_ {i} (z)}$Por lo general, se eligen como combinaciones lineales de las funciones completas.${\displaystyle \varphi _{k}(c_{i}z),\varphi _{k}(z)}$, respectivamente, donde $${\displaystyle \varphi _{0}(z)=e^{z},\quad \varphi _{k}(z)=\int _{0}^{1}e^{(1-\theta )z}{\frac {\theta ^{k-1}}{(k-1)!}}d\theta ,\quad k\geq 1.}$$ Estas funciones satisfacen la relación de recursión. $${\displaystyle \varphi _{k+1}(z)={\frac {\varphi _{k}(z)-\varphi _{k}(0)}{z}},\ k\geq 0.}$$ Al introducir la diferencia${\ Displaystyle D_ {ni} = N_ {n} (U_ {ni}) -N_ {n} (u_ {n})}$, pueden ser reformulados de una manera más eficiente para su implementación (ver también ^{[ 3 ]} ) como $${\displaystyle {\begin{aligned}U_{ni}&=u_{n}+c_{i}h_{n}\varphi _{1}(c_{i}h_{n}L_{n})f(u_{n})+h_{n}\sum _{j=2}^{i-1}a_{ij}(h_{n}L_{n})D_{nj},\\u_{n+1}&=u_{n}+h_{n}\varphi _{1}(h_{n}L_{n})f(u_{n})+h_{n}\sum _{i=2}^{s}b_{i}(h_{n}L_{n})D_{ni}.\end{aligned}}}$$

Para implementar este esquema con tamaño de paso adaptativo, se pueden considerar, para el propósito de la estimación de errores locales, los siguientes métodos integrados. $${\displaystyle {\bar {u}}_{n+1}=u_{n}+h_{n}\varphi _{1}(h_{n}L_{n})f(u_{n})+h_{n}\sum _{i=2}^{s}{\bar {b}}_{i}(h_{n}L_{n})D_{ni},}$$ que utilizan las mismas etapas${\displaystyle U_{ni}}$pero con pesas${\displaystyle {\bar {b}}_{i}}$.

Para mayor comodidad, los coeficientes de los métodos exponenciales explícitos de Rosenbrock, junto con sus métodos integrados, pueden representarse utilizando el llamado cuadro de Butcher reducido de la siguiente manera:

Cox y Matthews ^{[ 5 ]} describen un método de diferenciación temporal exponencial (ETD) de cuarto orden que derivaron utilizando Maple .

Usamos su notación y asumimos que la función desconocida es${\displaystyle u}$y que tenemos una solución conocida${\displaystyle u_{n}}$en ese momento${\displaystyle t_{n}}$Además, haremos un uso explícito de un lado derecho que posiblemente dependa del tiempo:${\displaystyle {\mathcal {N}}={\mathcal {N}}(t,u)}$.

Primero se construyen tres valores de etapa: $${\displaystyle {\begin{aligned}a_{n}&=e^{Lh/2}u_{n}+L^{-1}\left(e^{Lh/2}-I\right){\mathcal {N}}(t_{n},u_{n})\\b_{n}&=e^{Lh/2}u_{n}+L^{-1}\left(e^{Lh/2}-I\right){\mathcal {N}}(t_{n}{+}{\tfrac {h}{2}},a_{n})\\c_{n}&=e^{Lh/2}a_{n}+L^{-1}\left(e^{Lh/2}-I\right)\left(2{\mathcal {N}}(t_{n}{+}{\tfrac {h}{2}},b_{n})-{\mathcal {N}}(t_{n},u_{n})\right)\end{aligned}}}$$ La actualización final la proporciona, $${\displaystyle {\begin{aligned}u_{n+1}=e^{Lh}u_{n}+h^{-2}L^{-3}{\Bigl \{}{\hphantom {{}+2}}&\left[-4-Lh+e^{Lh}\left(4-3Lh+(Lh)^{2}\right)\right]{\mathcal {N}}(t_{n},u_{n})\\{}+2&\left[2+Lh+e^{Lh}\left(-2+Lh\right)\right]\left[{\mathcal {N}}(t_{n}{+}{\tfrac {h}{2}},a_{n})+{\mathcal {N}}(t_{n}{+}{\tfrac {h}{2}},b_{n})\right]\\{}+{\hphantom {2}}&\left[-4-3Lh-(Lh)^{2}+e^{Lh}\left(4-Lh\right)\right]{\mathcal {N}}(t_{n}{+}h,c_{n}){\Bigr \}}.\end{aligned}}}$$

Si se implementa de forma ingenua, el algoritmo anterior sufre inestabilidades numéricas debido a errores de redondeo de punto flotante . ^{[ 10 ]} Para ver por qué, considere la primera función, $${\displaystyle \varphi _{1}(z)={\frac {e^{z}-1}{z}},}$$ que está presente en el método de Euler de primer orden, así como en las tres etapas de ETDRK4. Para valores pequeños de${\displaystyle z}$, esta función sufre errores de cancelación numérica. Sin embargo, estos problemas numéricos pueden evitarse evaluando la${\displaystyle \varphi _{1}}$función mediante un enfoque de integral de contorno ^{[ 10 ]} o mediante una aproximación de Padé . ^{[ 11 ]}

Los integradores exponenciales se utilizan para la simulación de escenarios rígidos en computación científica y visual , por ejemplo, en dinámica molecular , ^{[ 12 ]} para simulación de circuitos VLSI , ^{[ 13 ] [ 14 ]} y en gráficos por computadora . ^{[ 15 ]} También se aplican en el contexto de métodos híbridos de Monte Carlo . ^{[ 16 ]} En estas aplicaciones, los integradores exponenciales muestran la ventaja de una gran capacidad de paso de tiempo y alta precisión. Para acelerar la evaluación de funciones matriciales en escenarios tan complejos, los integradores exponenciales se combinan a menudo con métodos de proyección de subespacios de Krylov.

• métodos lineales generales

• integradores en GPGPUs
• código para un integrador exponencial sin malla