Convolución adaptativa, rápida y con poca memoria para ecuaciones de evolución con poca memoria

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Convoluci´on adaptativa, r´apida y con poca memoria para ecuaciones de evoluci´on con memoria M. L´ opez-Fern´ andez 1, C. Lubich2, A. Sch¨ adle3 1Dpto. de Matem´aticas, Universidad Aut´onoma de Madrid, Campus de Cantoblanco, Ctra. de Colmenar Viejo Km. 15, E-28049 Madrid. E-mail: [email protected]. 2Mathematisches Institut, Universit¨at T¨ubingen, Auf der Morgenstelle 10, D–72076 T¨ubingen, Germany. E-mail: [email protected]. 3ZIB Berlin, Takustr. 7, D-14195 Berlin, Germany. E-mail: [email protected]. Palabras clave: convolution quadrature, adaptivity, Volterra integral equations, numerical inverse Laplace transform, anomalous diffusion. Resumen Presentamos un algoritmo de paso variable para aproximar convoluciones que aparecen como t´erminos con memoria en ecuaciones de evoluci´on. Para avanzar Npasos, el algoritmo requiere s´olo O(Nlog N) operaciones y O(log N) datos en la memoria activa, en lugar de las O(N2) operaciones y O(N) datos en la memoria necesarios para una implementaci´on m´as directa. Una caracter´ıstica b´asica del algoritmo es la reducci´on, v´ıa integrales de contorno, a ecuaciones diferenciales escalares, que se integran num´ericamente con paso variable. En lugar del n´ucleo de convoluci´on, el m´etodo utiliza los valores de su transformada de Laplace. Ilustramos el algoritmo aplic´andolo a un modelo para reacciones qu´ımicas con difusi´on inhibida. 1. Introducci´on Consideramos el c´alculo de convoluciones Zt 0 f(t−τ)g(τ)dτ , 0≤t≤T, (1) posiblemente con n´ucleo matricial fy funci´on vectorial g, y la resoluci´on de ecuaciones de evoluci´on con memoria conteniendo integrales de convoluci´on de este tipo donde gno es una funci´on conocida a priori, sino que g(τ) depende de la soluci´on en tiempo τde la ecuaci´on integral o integro-diferencial. 1 M. L´opez-Fern´andez, C. Lubich, Achim Sch¨adle En trabajos anteriores [7, 3, 8], se han desarrollado algoritmos que son r´apidos y con poca memoria: para aproximar (1) en una secuencia t=nh (n= 0,1, . . . , N) con tama˜no de paso constante hyNh =T, estos algoritmos necesitan O(Nlog N) operaciones, O(log N) evaluaciones de la tansformada de Laplace F=Lf, ninguna de f, y O(log N) entradas en la memoria activa. En el paso n-´esimo, gse eval´ua en tn=nh, pero el algoritmo “olvida”la historia g(tj) para j < n y s´olo guarda en memoria unas pocas combinaciones lineales de dichos valores de g. Esto contrasta con las O(N2) operaciones, O(N) evaluaciones del n´ucleo f, y O(N) datos en memoria necesarios para una implementaci´on m´as directa de una f´ormula de cuadratura para (5). Se˜nalamos adem´as que en muchas aplicaciones es la transformada de Laplace Fla que se conoce de antemano, en lugar del n´ucleo f. Una caracter´ıstica b´asica del algoritmo r´apido que presentamos aqu´ı es la reducci´on, v´ıa la representaci´on mediante integrales de contorno, a ecuaciones diferenciales de la forma y0=λy +gpara ciertos λ∈C, que son aproximadas num´ericamente. La integraci´on de estas ecuaciones diferenciales no tiene que hacerse necesariamente con paso temporal fijo h, como se hace en [3, 7, 8], sino que el tama˜no del paso se puede adaptar al comportamiento de g. Esta observaci´on es el punto de partida del algoritmo “adaptativo”, r´apido y con poca memoria que presentamos aqu´ı. La necesidad de utilizar paso de tiempo variable en las aplicaciones para resolver ecuaciones de evoluci´on integro-diferenciales aparece en diversas referencias de la literatura [1, 2]. Sin embargo, ninguno de los algoritmos propuestos en estas referencias aprovecha la estructura de convoluci´on para reducir las O(N2) operaciones yO(N) datos en la memoria necesarios para dar Npasos de tiempo. El algoritmo de convoluci´on propuesto aqu´ı funciona en caso de que la transformada de Laplace Fsea sectorial: Fes anal´ıtica en un sector |arg(s−σ)|< π −ϕcon ϕ < 1 2π, y en este sector, |F(s)| ≤ M|s|−νpara alg´un Mreal y ν > 0,(2) siendo un ingrediente esencial la inversi´on num´erica de la transformada F. Para ello, discretizamos la f´ormula de inversi´on f(t) = 1 2πi ZΓ etλ F(λ)dλ, t > 0,(3) tomando como Γ un contorno apropiado en el sector de holomorf´ıa de F. Con un s´olo contorno Γ no podemos obtener una aproximaci´on uniformemente buena para todo t∈ (0, T], as´ı que utilizamos diferentes hip´erbolas Γ`para intervalos de aproximaci´on uniforme que podemos tomar con crecimiento geom´etrico, del tipo [B`−1h∗, B`+1h∗] con B∈N, B≥1, digamos B= 5, y con un paso m´ınimo h∗. El n´umero necesario de contornos est´a por tanto acotado por logB(T/h∗), en lugar del logB(N) para el algoritmo de paso fijo. Para m´as ejemplos y detalles v´ease [4]. 2. Preparaci´on: La inversi´on num´erica de la transformada de Laplace En la f´ormula de inversi´on (3), elegimos Γ como la rama izquierda de una hip´erbola parametrizada por γ(x) = µ(1 −sin(a+ix)) + σ, donde µ > 0 es un par´ametro de escala, 2 Convoluci´on adaptativa, r´apida y con poca memoria para ecuaciones de evoluci´on σes el desplazamiento en (2) y 0 <a<π/2−ϕ. De este modo Γ es la rama izquierda de la hip´erbola con centro en (µ, 0), focos en (0,0), (2µ, 0) y as´ıntotas formando ´angulos ±(π/2 + a) con el eje real, de modo que Γ permanece dentro del sector (2) de holomorf´ıa de F. Tras parametrizar (3), aproximamos faplicando a la integral a lo largo del eje real que resulta la regla del trapecio truncada, esto es f(t) = 1 2πi ZΓ etλ F(λ)dλ ≈ K X k=−K wketλkF(λk),(4) donde los pesos wky los nodos de cuadratura λkest´an dados por wk=τ 2πi γ0(kτ) y λk=γ(kτ), con τ > 0 un par´ametro apropiado. El buen comportamiento de esta f´ormula de cuadratura para aproximar (3) se debe a las buenas propiedades de la regla del trapecio cuando el integrando se puede prolongar anal´ıticamente a una banda horizontal en torno al eje real [9, 10]. Enunciamos a continuaci´on la siguiente cota de error, que decae exponencialmente con respecto al n´umero de nodos de cuadratura. Para m´as detalles, nos referimos a [5] y [6]. Teorema 1 [6] Supongamos que la transformada de Laplace Fsatisface la condici´on de sectorialidad (2). Para T > 0,Λ≥1,0<a<π/2−ϕ, y K≥1fijos, existen cantidades positivas C1, C2, C, c que dependen de ayΛ(Ctambi´en depende de Ta menos que σ < 0 en (2)) tales que la selecci´on de par´ametros τ=C1/K yµ=C2K/(Λt0)da lugar a un error de cuadratura en (4) acotado por |EK(t)| ≤ C tν−1¡²+e−cK¢, uniformemente para t∈[t0,Λt0], con 0< t0≤T/Λarbitrario. En la cota, νes el exponente en (2) y²la precisi´on en las evaluaciones de la transformada de Laplace Fy en las operaciones elementales involucradas en (4). De este modo, K=O(log 1 ε) nodos de cuadratura son suficientes para obtener una precisi´on O(ε) en la aproximaci´on de la integral de contorno (3), uniformemente para t∈[t0,Λt0]. En la pr´actica, tomamos a≈1 2(π 2−ϕ) y calculamos los valores de C1yC2 siguiendo el proceso de optimizaci´on descrito en [6]. 3. El algoritmo de convoluci´on de paso variable, r´apido y con poca memoria 3.1. Reducci´on local a ecuaciones diferenciales El objetivo es aproximar u(t) = Zt 0 f(t−τ)g(τ)dτ (5) en una secuencia de tiempos 0 < t1< . . . < tN, espaciados arbitrariamente. Para un tn dado, insertamos la f´ormula de inversi´on de la transformada de Laplace en (5) y escribimos Ztn 0 f(tn−τ)g(τ)dτ =Ztn 0 1 2πi ZΓ e(tn−τ)λF(λ)dλ g(τ)dτ . (6) 3 M. L´opez-Fern´andez, C. Lubich, Achim Sch¨adle La inversion num´erica de la transformada de Laplace se lleva a cabo aplicando la cuadratura (4). Como hemos visto en la Secci´on 2, se puede utilizar el mismo contorno Γ para aproximar fen diferentes valores de t∈[t0,Λt0], para un ratio dado Λ ≥1. En (6) necesitamos aproximar f(tn−τ) para tn−τ∈[0, tn], de modo que no podemos utilizar un ´unico contorno Γ y necesitamos dividir la integral en (5) en varias partes. Para tiempos intermedios apropiados 0 < t−< t+< tn, con (tn−t−)/(tn−t+)≤Λ, seleccionamos un contorno adecuado Γ para cada intervalo [tn−t+, tn−t−] y aproximamos Zt+ t− f(tn−τ)g(τ)dτ =Zt+ t− 1 2πi ZΓ e(tn−τ)λF(λ)dλ g(τ)dτ ≈Zt+ t− K X k=−K wke(tn−τ)λkF(λk)g(τ)dτ = K X k=−K wkF(λk)e(tn−t+)λky(t+, t−, λk), donde y(t+, t−, λk) es la soluci´on en t+de la EDO lineal no homog´enea y0=λky+g, con y(t−) = 0, para −K≤k≤K. Aproximamos y(t+, t−, λk) interpolando glinealmente a trozos en los tj,j= 0, . . . , N−1, e integrando de forma exacta. Denotando por ¯y(t+, t−, λk) a la aproximaci´on resultante, calculamos: Zt+ t− f(tn−τ)g(τ)dτ ≈ K X k=−K wkF(λk)e(tn−t+)λk¯y(t+, t−, λk).(7) 3.2. Llenado del mosaico La clave del algoritmo es la forma que en que se seleccionan los tiempos para dividir la integral en (5). Para ello, seguimos el esquema representado por el mosaico en el tri´angulo {(t, τ):0≤τ≤t≤T}que mostramos en la Figura 1, donde las piezas crecen geom´etricamente a medida que aumenta la distancia a la diagonal. Por el momento fijamos un paso m´ınimo h∗. Si en el interior del tri´angulo la l´ınea vertical de (tn,0) a (tn, tn) corta a L piezas distintas del mosaico, obtenemos los valores t− `yt+ `, para 1 ≤`≤L, como el menor y el mayor de los tiempos tj, 1 ≤j≤n, contenidos en la pieza `-´esima respectivamente. Cuando la pieza `-´esima contiene alguno de los tiempos del mallado, el ´ındice `se recoge en un conjunto de ´ındices J. Los valores t± `=t± n,` dependen de naunque por sencillez omitimos nen la notaci´on. Las piezas del mismo tama˜no forman una clase que representa la distancia a la diagonal en el mosaico y corresponden por tanto a un intervalo de aproximaci´on distinto y un contorno diferente para realizar la inversi´on de la transformada, lo que se traduce en un conjunto distinto de 2K+ 1 ecuaciones diferenciales escalares. Los intervalos de aproximaci´on para los valores tn−t± `son del tipo I`= [B`−1h∗, B`+1h∗], 2≤`≤L, de modo que el ratio Λ es B2. Dado que los tiempos tjno est´an equiespaciados, en esta partici´on aparecen “huecos” entre los valores t+ `+myt− `, que corresponden en la Figura 1 a pares de l´ıneas horizontales con alguna frontera de las piezas del mosaico enmedio. De este modo, dividimos (5) en 2|J|+ 1 partes u(tn) = ˜u(0)(tn) + X `∈J u(`)(tn) + X `∈J ˜u(`)(tn) (8) donde u(`)(tn) = Zt+ ` t− ` f(tn−τ)g(τ)dτ (9) 4 Convoluci´on adaptativa, r´apida y con poca memoria para ecuaciones de evoluci´on 0 1 2 3 4 0 1 2 3 4 t15 t3 +=t11 t1 −=t12 t1 +=t13 t τ Figura 1: Mosaico en el tri´angulo τ≤tpara B= 3 con los tiempos tjindicados por puntos. Cada rect´angulo monocrom´atico completamente cerrado por l´ıneas negras corresponde a un conjunto de soluciones y(t+ `, t− `, λk),1≤k≤2K+ 1, de EDOs lineales. se aproxima siguiendo (7) y ˜u(0)(tn) = Ztn tn−1 f(tn−τ)g(τ)dτ y ˜u(`)(tn) = Zt− ` t+ `+m f(tn−τ)g(τ)dτ (10) corresponden al paso de tn−1atny a los “huecos” entre t+ `+myt− `, respectivamente (v´ease la Figura 1). Estas partes se calculan mediante “pasos directos” que explicamos en la siguiente secci´on. El procedimiento formal es el siguiente: dado un paso m´ınimo h∗y una base B, para cada tntomamos Lel menor entero tal que podemos representar dtn/h∗e= 2+PL `=1 b`B`−1con b`∈ {1,2, . . . , B}, donde dxedenota el menor de lo enteros mayores que x.t+ `es el mayor y t− `el menor de los valores en {tj:j= 0, . . . , n}, tales que h∗PL k=`+1 bkBk−1≤t− `≤t+ `≤h∗PL k=`bkBk−1. Los t± `son tales que los l´ımites de integraci´on en u(`)(tn) verifican tn−t± `∈I`. El llenado del mosaico representado en la Figura 1 se lleva a cabo avanzando todas las ecuaciones diferenciales para todos los contornos en cada paso de tiempo, de modo que el algoritmo va “olvidando” los valores pasados de g, a excepci´on los valores en t± `, necesarios para para los pasos directos descritos a contiuaci´on. En total se guardan tres copias de las soluciones a las EDOs. Pseudoc´odigos para la organizaci´on de esta descomposici´on y m´as detalles se pueden encontrar en [4]. 3.3. Pasos directos Los valores ˜u(`)(tn) que cubren los “huecos” [t+ `+m, t− `] est´an dados por ˜u(`)(tn) = Zt− ` t+ `+m f(tn−τ)g(τ)dτ =Ztj+1 tj f(tn−τ)g(τ)dτ, (11) 5 M. L´opez-Fern´andez, C. Lubich, Achim Sch¨adle para alg´un j∈ {0,1, . . . , N −1}. En [tj, tj+1], aproximamos g(t) por la funci´on lineal g(t)≈¯g(t) = gj+gj+1 −gj hj+1 (t−tj), hj+1 =tj+1 −tj, con gj=g(tj), j= 0,1, . . . , N. Prolongando ¯g(t) a todo [0, tn], calculamos Ztj+1 tj f(tn−τ) ¯g(τ)dτ =Ztn tj f(tn−τ) ¯g(τ)dτ −Ztn tj+1 f(tn−τ) ¯g(τ)dτ =L−1[F· L¯g(·+tj)](tn−tj)− L−1[F· L¯g(·+tj+1)](tn−tj+1) =L−1·F1gj+F2 gj+1 −gj hj+1 ¸(tn−tj)− L−1·F1gj+1 +F2 gj+1 −gj hj+1 ¸(tn−tj+1), donde F1(s) = F(s)/s yF2(s) = F(s)/s2. Aproximamos las transformadas inversas f1(t) = ¡L−1F1¢(t) y f2(t) = ¡L−1F2¢(t) en t=tn−tjyt=tn−tj+1 utilizando de nuevo (4) con los contornos correspondientes a los intervalos de aproximaci´on I`1yI`2con tn−tj+1 ∈I`2 ytn−tj∈I`1(notamos que f1(0) = f2(0) = 0). El paso directo se aproxima entonces por la combinaci´on lineal Ztj+1 tj f(tn−τ) ¯g(τ)dτ =f1(tn−tj)gj+f2(tn−tj)gj+1 −gj hj+1 −f1(tn−tj+1)gj+1 −f2(tn−tj+1)gj+1 −gj hj+1 . (12) Nota acerca de la complejidad. Dada una secuencia de tiempos cualquiera 0 < t1< . . . < tN=Tcon paso m´ınimo h∗= m´ınj(tj+1 −tj), el algoritmo aproxima (5) con error menor que εrealizando un total de O(NLK) = O(Nlog T h∗log 1 ε) operaciones y guardando en memoria O(LK) = O(log T h∗log 1 ε) vectores del tama˜no de g. Notemos que el algoritmo es independiente del m´etodo elegido para seleccionar el tama˜no de los pasos en la integraci´on num´erica de las ecuaciones integrales e integrodiferenciales. En el ejemplo siguiente la estrategia seguida se basa en controlar el error en la interpolaci´on lineal a trozos de gen (5), que por otra parte es la ´unica fuente de error en el esquema aparte de la discretizaci´on de las integrales de contorno, pero otras alternativas son posibles. Para m´as ejemplos v´ease [4]. 4. Din´amica de reacciones qu´ımicas con difusi´on inhibida Consideramos tres especies moleculares A,ByC, con concentraciones u1,u2yu3 respectivamente, que reaccionan en la forma A+Bk1 −→ C,Ck2 −→ A+B,Ck3 −→ A+P, donde Pes el producto resultante. La difusi´on de las especies es an´omala, dando lugar a una ecuaci´on de reacci´on difusi´on con un t´ermino de memoria. Siguiendo el modelo propuesto en [11], (cf. [2]), consideramos el siguiente sistema de ecuaciones integro-diferenciales: u1(t)−u1(0) =∂−α t£K∆u1(t)−k1u1(t)u2(t)+(k2+k3)u3(t)¤ u2(t)−u2(0) =∂−α t£K∆u2(t)−k1u1(t)u2(t) + k2u3(t)¤ u3(t)−u3(0) =∂−α t£K∆u3(t) + k1u1(t)u2(t)−(k2+k3)u3(t)¤, (13) 6 Convoluci´on adaptativa, r´apida y con poca memoria para ecuaciones de evoluci´on donde ∆ = ∂xx es el Laplaciano en 1D con condiciones de contorno peri´odicas en [−5,5] y∂−α tdenota la integral fraccionaria de orden 0 < α < 1, definida por el operador de Riemann Liouville ∂−α tg(t) = 1 Γ(α)Zt 0 (t−τ)α−1g(τ)dτ. (14) El n´ucleo de convoluci´on f(t) = tα−1/Γ(α) es d´ebilmente singular y tiene transformada de Laplace F(s) = s−α. Elegimos K= 0.5, k1= 1, k2= 2 y k3= 3, e integramos hasta T= 30. Fijamos α= 0.5 y consideramos funciones escal´on suavizadas como dato inicial. Denotando u= [u1, u2, u3]T,I3la matriz identidad 3 ×3, y siguiendo la notaci´on de la Secci´on 3.3 para los pasos directos, obtenemos la ecuaci´on discreta ³I3⊗I−f2(hn) hn (K I3⊗S+R)´un=u0+1 Γ(α)Ztn−1 0 g(¯ u(τ)) (t−τ)1−αdτ +³f1(hn)−f2(hn) hn´¡K I3⊗S+R¢un−1+k1f1(hn)e⊗un−1 1un−1 2, donde ¯ udenota el interpolante lineal a trozos de uen los tiempos t0, t1, . . . , tny, para M nodos en la discretizaci´on espacial y vun vector columna de longitud 3M, definimos g(v) = ¡K I3⊗S+R¢v+k1¡e⊗v1v2¢,e= [−1,−1,1]T, con Sla aproximaci´on de ∂xx por diferencias finitas de segundo orden y condiciones de contorno peri´odicas y Rla matriz 3M×3M R=  0 0 (k2+k3)IM 0 0 k2IM 0 0 −(k2+k3)IM  , con IMla matriz identidad M×M. Dada una tolerancia Tol, la estrategia para elegir el nuevo paso hn+1 es Ch2 n+1γ0 n= 0.8·Tol ,(15) con γ0 n=k¯g0(tn)kyC≈RT 0|f(t)|dt. Adicionalmente imponemos hn/4< hn+1 <4hn. El test para el nuevo paso propuesto hn+1 es Ch2 n+1γ0 n+1 ≤Tol , donde el nuevo valor g(tn+hn+1) se utiliza en el c´alculo de γ0 n+1. Si se cumple el test, aceptamos hn+1 y tomamos tn+1 =tn+hn+1. Si no, repetimos el procedimiento con hn+1 obtenido a partir de (15) con γ0 n+1 en lugar de γ0 n. Este procedimiento se repite si es necesario hasta que se verifique el test. La Figura 2 muestra los resultados num´ericos obtenidos al aplicar el algoritmo. La inversi´on de la transformada de Laplace se ha llevado a cabo con a= 1, d= 0.5, y K= 40 nodos de cuadratura, lo que da C1= 6.036 y C2= 0.0739. El algoritmo se ha aplicado con base B= 5, lo que proporciona intervalos de aproximaci´on del tipo [t0,25t0]. Agradecimientos El primer autor ha sido respaldado por la DGI-MCYT a trav´es del proyecto MTM 2004-07194, cofinanciado con fondos FEDER. El segundo autor por DFG, SFB 382 y el tercero por DFG Research Center Matheon “Mathematics for key technologies”, Berl´ın. 7 M. L´opez-Fern´andez, C. Lubich, Achim Sch¨adle −5 0 5 0 0.2 0.4 0.6 0.8 1 t 0 t 1 t 4 t 8 t 16 t 27 x A −5 0 5 0 0.2 0.4 0.6 0.8 1 t 0 t 1 t 4 t 8 t 16 t 27 x B −5 0 5 0 2 4 6 8 x 10−3 t 0 t 1 t 4 t 8 t 16 t 27 x C 0 10 20 30 40 10−5 10−4 10−3 10−2 10−1 100 tol = 1×10−6 tol = 1×10−2 tol = 1×10−3 tol = 1×10−4 tol = 1×10−5 t stepsize 10−5 10−4 10−3 10−2 10−5 10−4 10−3 tol absolute error 101102103104 10−5 10−4 10−3 number of steps absolute error Figura 2: Arriba: Soluciones para las tres especies A,B,C, en tiempos diferentes. Abajo: Evoluci´on del tama˜no de paso, error vs tolerancias y error vs n´umero de pasos. Referencias [1] K. Adolfsson, M. Enelund, and S. Larsson,Adaptive discretization of fractional order viscoelasticity using sparse time history, Computer Meth. Appl. Mech. Eng., 193 (2004), pp. 4567–4590. [2] X. Cao, K. Burrage, and F. Abdullah,A variable stepsize implementation for fractional differential equations, Report, 2006. [3] R. Hiptmair and A. Sch¨ adle,Non-reflecting boundary conditions for Maxwell’s equations, Computing, 71 (2003), pp. 265–292. [4] M. L´ opez-Fern´ andez, C. Lubich, and A. Sch¨ adle,Adaptive, fast and oblivious convolution in evolution equations with memory, Preprint 2006. [5] M. L´ opez-Fern´ andez and C. Palencia,On the numerical inversion of the Laplace transform of certain holomorphic mappings, Appl. Numer. Math., 51 (2004), pp. 289–303. [6] M. L´ opez-Fern´ andez, C. Palencia, and A. Sch¨ adle,A spectral order method for inverting sectorial Laplace transforms, SIAM J. Numer. Anal., 44 (2006), pp. 1332–1350. [7] C. Lubich and A. Sch¨ adle,Fast convolution for nonreflecting boundary conditions, SIAM J. Sci. Comput., 24 (2002), pp. 161–182. [8] A. Sch¨ adle, M. L´ opez-Fern´ andez, and C. Lubich,Fast and oblivious convolution quadrature, SIAM J. Sci. Comput. 28 (2006), pp. 421–438. [9] F. Stenger,Approximations via Whittaker’s Cardinal Function, J. Approx. Theory, 17 (1976), pp. 222–240. [10] F. Stenger,Numerical methods based on Whittaker Cardinal, or sinc Functions, SIAM Rev., 23 (1981), pp. 165–224. [11] S. B. Yuste, L. Acedo, K. Lindenberg,Reaction front in an A+B →C reaction-subdiffusion process, Phys. Rev. E 69, (2004), pp. 036126. 8