Método implícito de dirección alternada

En álgebra lineal numérica , el método implícito de dirección alternada (ADI) es un método iterativo utilizado para resolver ecuaciones matriciales de Sylvester . Es un método popular para resolver las grandes ecuaciones matriciales que surgen en la teoría y el control de sistemas , ^[1] y se puede formular para construir soluciones en una forma factorizada y de memoria eficiente. ^{[2] [3]} También se utiliza para resolver numéricamente ecuaciones diferenciales parciales parabólicas y elípticas , y es un método clásico utilizado para modelar la conducción de calor y resolver la ecuación de difusión en dos o más dimensiones. ^[4] Es un ejemplo de un método de división de operadores . ^[5]

ADI para ecuaciones matriciales

El método

El método ADI es un proceso de iteración de dos pasos que actualiza alternativamente los espacios de columnas y filas de una solución aproximada a . Una iteración ADI consta de los siguientes pasos: ^[6] ${\displaystyle AX-XB=C}$

1. Resuelva para , donde ${\displaystyle X^{(j+1/2)}}$ ${\displaystyle \left(A-\beta _{j+1}I\right)X^{(j+1/2)}=X^{(j)}\left(B-\beta _{j+1}I\right)+C.}$

2. Resuelva para , donde . ${\displaystyle X^{(j+1)}}$ ${\displaystyle X^{(j+1)}\left(B-\alpha _{j+1}I\right)=\left(A-\alpha _{j+1}I\right)X^{(j+1/2)}-C}$

Los números se denominan parámetros de desplazamiento y la convergencia depende en gran medida de la elección de estos parámetros. ^{[7] [8]} Para realizar iteraciones de ADI, se requiere una estimación inicial, así como parámetros de desplazamiento, . ${\displaystyle (\alpha _{j+1},\beta _{j+1})}$ ${\displaystyle K}$ ${\displaystyle X^{(0)}}$ ${\displaystyle K}$ ${\displaystyle \{(\alpha _{j},\beta _{j})\}_{j=1}^{K}}$

Cuándo utilizar ADI

Si y , entonces se pueden resolver directamente utilizando el método de Bartels-Stewart . ^[9] Por lo tanto, solo es beneficioso utilizar ADI cuando la multiplicación de matriz-vector y las resoluciones lineales que involucran y se pueden aplicar de manera económica. ${\displaystyle A\in \mathbb {C} ^{m\times m}}$ ${\displaystyle B\in \mathbb {C} ^{n\times n}}$ ${\displaystyle AX-XB=C}$ ${\displaystyle {\mathcal {O}}(m^{3}+n^{3})}$ ${\displaystyle A}$ ${\displaystyle B}$

La ecuación tiene una solución única si y solo si , donde es el espectro de . ^[1] Sin embargo, el método ADI funciona especialmente bien cuando y están bien separados, y y son matrices normales . Estos supuestos se cumplen, por ejemplo, mediante la ecuación de Lyapunov cuando es definida positiva . Bajo estos supuestos, se conocen parámetros de cambio casi óptimos para varias opciones de y . ^{[7] [8]} Además, se pueden calcular límites de error a priori, eliminando así la necesidad de monitorear el error residual en la implementación. ${\displaystyle AX-XB=C}$ ${\displaystyle \sigma (A)\cap \sigma (B)=\emptyset }$ ${\displaystyle \sigma (M)}$ ${\displaystyle M}$ ${\displaystyle \sigma (A)}$ ${\displaystyle \sigma (B)}$ ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle AX+XA^{*}=C}$ ${\displaystyle A}$ ${\displaystyle A}$ ${\displaystyle B}$

El método ADI todavía se puede aplicar cuando no se cumplen los supuestos anteriores. El uso de parámetros de desplazamiento subóptimos puede afectar negativamente a la convergencia, ^[1] y la convergencia también se ve afectada por la no normalidad de o (a veces de manera ventajosa). ^{[10] Se observa que los métodos} del subespacio de Krylov , como el método del subespacio racional de Krylov, ^[11] convergen típicamente más rápidamente que el ADI en este contexto, ^{[1] [3]} y esto ha llevado al desarrollo de métodos híbridos de proyección-ADI. ^[3] ${\displaystyle A}$ ${\displaystyle B}$

Selección de parámetros de desplazamiento y ecuación de error ADI

El problema de encontrar buenos parámetros de cambio no es trivial. Este problema se puede entender examinando la ecuación de error de ADI. Después de las iteraciones, el error viene dado por ${\displaystyle K}$

${\displaystyle X-X^{(K)}=\prod _{j=1}^{K}{\frac {(A-\alpha _{j}I)}{(A-\beta _{j}I)}}\left(X-X^{(0)}\right)\prod _{j=1}^{K}{\frac {(B-\beta _{j}I)}{(B-\alpha _{j}I)}}.}$

La elección da como resultado el siguiente límite en el error relativo: ${\displaystyle X^{(0)}=0}$

${\displaystyle {\frac {\left\|X-X^{(K)}\right\|_{2}}{\|X\|_{2}}}\leq \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2},\quad r_{K}(M)=\prod _{j=1}^{K}{\frac {(M-\alpha _{j}I)}{(M-\beta _{j}I)}}.}$

donde es la norma del operador . El conjunto ideal de parámetros de desplazamiento define una función racional que minimiza la cantidad . Si y son matrices normales y tienen descomposiciones propias y , entonces ${\displaystyle \|\cdot \|_{2}}$ ${\displaystyle \{(\alpha _{j},\beta _{j})\}_{j=1}^{K}}$ ${\displaystyle r_{K}}$ ${\displaystyle \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2}}$ ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle A=V_{A}\Lambda _{A}V_{A}^{*}}$ ${\displaystyle B=V_{B}\Lambda _{B}V_{B}^{*}}$

${\displaystyle \|r_{K}(A)\|_{2}\|r_{K}(B)^{-1}\|_{2}=\|r_{K}(\Lambda _{A})\|_{2}\|r_{K}(\Lambda _{B})^{-1}\|_{2}}$.

Parámetros de cambio casi óptimos

Los parámetros de desplazamiento casi óptimos se conocen en ciertos casos, como cuando y , donde y son intervalos disjuntos en la línea real. ^{[7] [8]} La ecuación de Lyapunov , por ejemplo, satisface estos supuestos cuando es definida positiva . En este caso, los parámetros de desplazamiento se pueden expresar en forma cerrada utilizando integrales elípticas , y se pueden calcular numéricamente fácilmente. ${\displaystyle \Lambda _{A}\subset [a,b]}$ ${\displaystyle \Lambda _{B}\subset [c,d]}$ ${\displaystyle [a,b]}$ ${\displaystyle [c,d]}$ ${\displaystyle AX+XA^{*}=C}$ ${\displaystyle A}$

De manera más general, si se conocen conjuntos cerrados y disjuntos y , donde y , el problema de selección del parámetro de desplazamiento óptimo se resuelve aproximadamente encontrando una función racional extremal que alcance el valor ${\displaystyle E}$ ${\displaystyle F}$ ${\displaystyle \Lambda _{A}\subset E}$ ${\displaystyle \Lambda _{B}\subset F}$

${\displaystyle Z_{K}(E,F):=\inf _{r}{\frac {\sup _{z\in E}|r(z)|}{\inf _{z\in F}|r(z)|}},}$

donde el ínfimo se toma sobre todas las funciones racionales de grado . ^[8] Este problema de aproximación está relacionado con varios resultados en la teoría del potencial , ^{[12] [13]} y fue resuelto por Zolotarev en 1877 para = [a, b] y ^[14] La solución también se conoce cuando y son discos disjuntos en el plano complejo. ^[15] ${\displaystyle (K,K)}$ ${\displaystyle E}$ ${\displaystyle F=-E.}$ ${\displaystyle E}$ ${\displaystyle F}$

Estrategias heurísticas de cambio de parámetros

Cuando se sabe menos acerca de y , o cuando o son matrices no normales, puede que no sea posible encontrar parámetros de desplazamiento casi óptimos. En este contexto, se puede utilizar una variedad de estrategias para generar buenos parámetros de desplazamiento. Estas incluyen estrategias basadas en resultados asintóticos en la teoría del potencial, ^[16] utilizando los valores de Ritz de las matrices , , , y para formular un enfoque voraz, ^[17] y métodos cíclicos, donde se reutiliza la misma pequeña colección de parámetros de desplazamiento hasta que se cumple una tolerancia de convergencia. ^{[17] [10]} Cuando se utiliza el mismo parámetro de desplazamiento en cada iteración, ADI es equivalente a un algoritmo llamado método de Smith. ^[18] ${\displaystyle \sigma (A)}$ ${\displaystyle \sigma (B)}$ ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle A}$ ${\displaystyle A^{-1}}$ ${\displaystyle B}$ ${\displaystyle B^{-1}}$

ADI factorizada

En muchas aplicaciones, y son matrices dispersas muy grandes, y se pueden factorizar como , donde , con . ^[1] En tal entorno, puede que no sea factible almacenar la matriz potencialmente densa de forma explícita. Se puede utilizar una variante de ADI, llamada ADI factorizada, ^{[3] [2]} para calcular , donde . La eficacia de la ADI factorizada depende de si se aproxima bien mediante una matriz de bajo rango. Se sabe que esto es cierto bajo varios supuestos sobre y . ^{[10] [8]} ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle C}$ ${\displaystyle C=C_{1}C_{2}^{*}}$ ${\displaystyle C_{1}\in \mathbb {C} ^{m\times r},C_{2}\in \mathbb {C} ^{n\times r}}$ ${\displaystyle r=1,2}$ ${\displaystyle X}$ ${\displaystyle ZY^{*}}$ ${\displaystyle X\approx ZY^{*}}$ ${\displaystyle X}$ ${\displaystyle A}$ ${\displaystyle B}$

ADI para ecuaciones parabólicas

Históricamente, el método ADI se desarrolló para resolver la ecuación de difusión 2D en un dominio cuadrado utilizando diferencias finitas. ^[4] A diferencia del ADI para ecuaciones matriciales, el ADI para ecuaciones parabólicas no requiere la selección de parámetros de desplazamiento, ya que el desplazamiento que aparece en cada iteración está determinado por parámetros como el paso de tiempo, el coeficiente de difusión y el espaciado de la cuadrícula. La conexión con el ADI en ecuaciones matriciales se puede observar cuando se considera la acción de la iteración ADI en el sistema en estado estable.

Ejemplo: ecuación de difusión 2D

Figura de plantilla para el método implícito de dirección alternada en ecuaciones de diferencias finitas

El método tradicional para resolver numéricamente la ecuación de conducción de calor es el método Crank-Nicolson . Este método da como resultado un conjunto muy complicado de ecuaciones en múltiples dimensiones, que son costosas de resolver. La ventaja del método ADI es que las ecuaciones que deben resolverse en cada paso tienen una estructura más simple y se pueden resolver de manera eficiente con el algoritmo de matriz tridiagonal .

Consideremos la ecuación de difusión lineal en dos dimensiones,

${\displaystyle {\partial u \over \partial t}=\left({\partial ^{2}u \over \partial x^{2}}+{\partial ^{2}u \over \partial y^{2}}\right)=(u_{xx}+u_{yy})}$

El método implícito de Crank-Nicolson produce la siguiente ecuación de diferencias finitas:

${\displaystyle {u_{ij}^{n+1}-u_{ij}^{n} \over \Delta t}={1 \over 2(\Delta x)^{2}}\left(\delta _{x}^{2}+\delta _{y}^{2}\right)\left(u_{ij}^{n+1}+u_{ij}^{n}\right)}$

dónde:

${\displaystyle \Delta x=\Delta y}$

y es el segundo operador de diferencia central para la coordenada p -ésima ${\displaystyle \delta _{p}^{2}}$

${\displaystyle \delta _{p}^{2}u_{ij}=u_{ij+e_{p}}-2u_{ij}+u_{ij-e_{p}}}$

con o para o respectivamente (y una abreviatura de puntos reticulares ). ${\displaystyle e_{p}=(1,0)}$ ${\displaystyle (0,1)}$ ${\displaystyle p=x}$ ${\displaystyle y}$ ${\displaystyle ij}$ ${\displaystyle (i,j)}$

Después de realizar un análisis de estabilidad , se puede demostrar que este método será estable para cualquier . ${\displaystyle \Delta t}$

Una desventaja del método de Crank-Nicolson es que la matriz de la ecuación anterior tiene un ancho de banda que generalmente es bastante grande. Esto hace que la solución directa del sistema de ecuaciones lineales sea bastante costosa (aunque existen soluciones aproximadas eficientes, por ejemplo, el uso del método del gradiente conjugado preacondicionado con factorización de Cholesky incompleta ).

La idea detrás del método ADI es dividir las ecuaciones de diferencias finitas en dos, una con la derivada de x tomada implícitamente y la siguiente con la derivada de y tomada implícitamente.

${\displaystyle {u_{ij}^{n+1/2}-u_{ij}^{n} \over \Delta t/2}={\left(\delta _{x}^{2}u_{ij}^{n+1/2}+\delta _{y}^{2}u_{ij}^{n}\right) \over \Delta x^{2}}}$

${\displaystyle {u_{ij}^{n+1}-u_{ij}^{n+1/2} \over \Delta t/2}={\left(\delta _{x}^{2}u_{ij}^{n+1/2}+\delta _{y}^{2}u_{ij}^{n+1}\right) \over \Delta y^{2}}}$

El sistema de ecuaciones involucrado es simétrico y tridiagonal (con ancho de banda de 3), y normalmente se resuelve utilizando el algoritmo de matriz tridiagonal .

Se puede demostrar que este método es incondicionalmente estable y de segundo orden en el tiempo y el espacio. ^[19] Existen métodos ADI más refinados, como los métodos de Douglas, ^[20] o el método del factor f ^[21], que se pueden utilizar para tres o más dimensiones.

Generalizaciones

El uso del método ADI como esquema de división de operadores se puede generalizar, es decir, podemos considerar ecuaciones de evolución generales

${\displaystyle {\dot {u}}=F_{1}u+F_{2}u,}$

donde y son operadores (posiblemente no lineales) definidos en un espacio de Banach . ^{[22] [23]} En el ejemplo de difusión anterior tenemos y . ${\displaystyle F_{1}}$ ${\displaystyle F_{2}}$ ${\displaystyle F_{1}={\partial ^{2} \over \partial x^{2}}}$ ${\displaystyle F_{2}={\partial ^{2} \over \partial y^{2}}}$

ADI fundamental (FADI)

Simplificación de ADI a FADI

Es posible simplificar el método ADI convencional en el método ADI fundamental, que solo tiene operadores similares en los lados izquierdos mientras que no tiene operadores en los lados derechos. Esto puede considerarse como el esquema fundamental (básico) del método ADI, ^{[24] [25]} sin más operadores (a reducir) en los lados derechos, a diferencia de la mayoría de los métodos implícitos tradicionales que generalmente consisten en operadores en ambos lados de las ecuaciones. El método FADI conduce a ecuaciones de actualización más simples, más concisas y eficientes sin degradar la precisión del método ADI convencional.

Relaciones con otros métodos implícitos

Muchos métodos implícitos clásicos de Peaceman-Rachford, Douglas-Gunn, D'Yakonov, Beam-Warming, Crank-Nicolson, etc., pueden simplificarse a esquemas implícitos fundamentales con lados derechos libres de operadores. ^[25] En sus formas fundamentales, el método FADI de precisión temporal de segundo orden puede relacionarse estrechamente con el método fundamental localmente unidimensional (FLOD), que puede actualizarse a una precisión temporal de segundo orden, como para las ecuaciones de Maxwell tridimensionales ^{[26] [27]} en electromagnetismo computacional . Para ecuaciones de conducción y difusión de calor bidimensionales y tridimensionales, tanto los métodos FADI como FLOD pueden implementarse de manera más simple, más eficiente y estable en comparación con sus métodos convencionales. ^{[28] [29]}

Referencias

1. Simoncini, V. (2016). "Métodos computacionales para ecuaciones matriciales lineales". SIAM Review . 58 (3): 377–441. doi :10.1137/130912839. hdl : 11585/586011 . ISSN  0036-1445. S2CID  17271167.
2. Li, Jing-Rebecca ; White, Jacob (2002). "Solución de bajo rango de ecuaciones de Lyapunov". Revista SIAM sobre análisis de matrices y aplicaciones . 24 (1): 260–280. doi :10.1137/s0895479801384937. ISSN  0895-4798.
3. Benner, Peter; Li, Ren-Cang; Truhar, Ninoslav (2009). "Sobre el método ADI para ecuaciones de Sylvester". Revista de Matemática Computacional y Aplicada . 233 (4): 1035–1045. Bibcode :2009JCoAM.233.1035B. doi : 10.1016/j.cam.2009.08.108 . ISSN  0377-0427.
4. Peaceman, DW; Rachford Jr., HH (1955), "La solución numérica de ecuaciones diferenciales parabólicas y elípticas", Revista de la Sociedad de Matemáticas Industriales y Aplicadas , 3 (1): 28–41, doi :10.1137/0103003, hdl : 10338.dmlcz/135399 , MR  0071874.
5. * Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Sección 20.3.3. Métodos de división de operadores en general". Recetas numéricas: el arte de la computación científica (3.ª ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8Archivado desde el original el 11-08-2011 . Consultado el 18-08-2011 .
6. Wachspress, Eugene L. (2008). "El camino hacia un solucionador de ecuaciones de Lyapunov". Computers & Mathematics with Applications . 55 (8): 1653–1659. doi : 10.1016/j.camwa.2007.04.048 . ISSN  0898-1221.
7. Lu, An; Wachspress, EL (1991). "Solución de ecuaciones de Lyapunov mediante iteración implícita de dirección alternada". Computers & Mathematics with Applications . 21 (9): 43–58. doi :10.1016/0898-1221(91)90124-m. ISSN  0898-1221.
8. Beckermann, Bernhard; Townsend, Alex (2017). "Sobre los valores singulares de matrices con estructura de desplazamiento". Revista SIAM sobre análisis de matrices y aplicaciones . 38 (4): 1227–1248. arXiv : 1609.09494 . doi :10.1137/16m1096426. ISSN  0895-4798. S2CID  3828461.
9. Golub, G.; Van Loan, C (1989). Cálculos matriciales (Cuarta edición). Baltimore: Johns Hopkins University. ISBN 1421407949.OCLC 824733531  .
10. Sabino, J (2007). Solución de ecuaciones de Lyapunov a gran escala mediante el método de Smith modificado por bloques . Tesis doctoral, Universidad Rice . hdl :1911/20641.
11. Druskin, V.; Simoncini, V. (2011). "Subespacios de Krylov racionales adaptativos para sistemas dinámicos de gran escala". Systems & Control Letters . 60 (8): 546–560. doi :10.1016/j.sysconle.2011.04.013. ISSN  0167-6911.
12. Saff, EB; Totik, V. (11 de noviembre de 2013). Potenciales logarítmicos con campos externos . Berlín. ISBN 9783662033296.OCLC 883382758  .{{cite book}}: CS1 maint: location missing publisher (link)
13. Gonchar, AA (1969). "Problemas de Zolotarev relacionados con funciones racionales". Matemáticas de la URSS-Sbornik . 7 (4): 623–635. Código Bibliográfico :1969SbMat...7..623G. doi :10.1070/SM1969v007n04ABEH001107.
14. Zolotarev, DI (1877). "Aplicación de funciones elípticas a cuestiones de funciones que se desvían mínima y máximamente de cero". Zap. Imp. Akad. Nauk. San Petersburgo . 30 : 1–59.
15. Starke, Gerhard (julio de 1992). "Casi circularidad para el problema racional de Zolotarev en el plano complejo". Journal of Approximation Theory . 70 (1): 115–130. doi : 10.1016/0021-9045(92)90059-w . ISSN  0021-9045.
16. Starke, Gerhard (junio de 1993). "Puntos de Fejér-Walsh para funciones racionales y su uso en el método iterativo ADI". Journal of Computational and Applied Mathematics . 46 (1–2): 129–141. doi : 10.1016/0377-0427(93)90291-i . ISSN  0377-0427.
17. Penzl, Thilo (enero de 1999). "Un método Smith cíclico de bajo rango para ecuaciones de Lyapunov dispersas y de gran tamaño". Revista SIAM de informática científica . 21 (4): 1401–1418. Bibcode :1999SJSC...21.1401P. doi :10.1137/s1064827598347666. ISSN  1064-8275.
18. Smith, RA (enero de 1968). "Matriz de ecuaciones XA + BX = C". Revista SIAM de Matemáticas Aplicadas . 16 (1): 198–201. doi :10.1137/0116017. ISSN  0036-1399.
19. Douglas, J. Jr. (1955), "Sobre la integración numérica de u _xx + u _yy = u _t mediante métodos implícitos", Journal of the Society for Industrial and Applied Mathematics , 3 : 42–65, MR  0071875.
20. Douglas, Jim Jr. (1962), "Métodos de dirección alterna para tres variables espaciales", Numerische Mathematik , 4 (1): 41–63, doi :10.1007/BF01386295, ISSN  0029-599X, S2CID  121455963.
21. Chang, MJ; Chow, LC; Chang, WS (1991), "Método implícito mejorado de dirección alterna para resolver problemas transitorios de difusión de calor tridimensional", Transferencia numérica de calor, Parte B: Fundamentos , 19 (1): 69–84, Bibcode :1991NHTB...19...69C, doi :10.1080/10407799108944957, ISSN  1040-7790.
22. Hundsdorfer, Willem; Verwer, Jan (2003). Solución numérica de ecuaciones de advección-difusión-reacción dependientes del tiempo . Berlín, Heidelberg: Springer Berlin Heidelberg. ISBN 978-3-662-09017-6.
23. Lions, PL; Mercier, B. (diciembre de 1979). "Algoritmos de división para la suma de dos operadores no lineales". Revista SIAM sobre análisis numérico . 16 (6): 964–979. Bibcode :1979SJNA...16..964L. doi :10.1137/0716071.
24. Tan, EL (2007). "Algoritmo eficiente para el método ADI-FDTD tridimensional incondicionalmente estable" (PDF) . IEEE Microwave and Wireless Components Letters . 17 (1): 7–9. doi :10.1109/LMWC.2006.887239. hdl :10356/138245. S2CID  29025478.
25. Tan, EL (2008). "Esquemas fundamentales para métodos de dominio temporal de diferencias finitas implícitos, incondicionalmente estables y eficientes" (PDF) . IEEE Transactions on Antennas and Propagation . 56 (1): 170–177. arXiv : 2011.14043 . Bibcode :2008ITAP...56..170T. doi :10.1109/TAP.2007.913089. hdl :10356/138249. S2CID  37135325.
26. Tan, EL (2007). "Método LOD-FDTD incondicionalmente estable para ecuaciones de Maxwell en 3D" (PDF) . IEEE Microwave and Wireless Components Letters . 17 (2): 85–87. doi :10.1109/LMWC.2006.890166. hdl :10356/138296. S2CID  22940993.
27. Gan, TH; Tan, EL (2013). "Método LOD-FDTD fundamental incondicionalmente estable con precisión temporal de segundo orden y divergencia compatible" (PDF) . IEEE Transactions on Antennas and Propagation . 61 (5): 2630–2638. Bibcode :2013ITAP...61.2630G. doi :10.1109/TAP.2013.2242036. S2CID  7578037.
28. Tay, WC; Tan, EL; Heh, DY (2014). "Método fundamental localmente unidimensional para simulación térmica 3-D". IEICE Transactions on Electronics . E-97-C (7): 636–644. Bibcode :2014IEITE..97..636T. doi :10.1587/transele.E97.C.636. hdl : 10220/20410 .
29. Heh, DY; Tan, EL; Tay, WC (2016). "Método implícito de dirección alternante rápida para simulación térmica transitoria eficiente de circuitos integrados". Revista internacional de modelado numérico: redes electrónicas, dispositivos y campos . 29 (1): 93–108. doi :10.1002/jnm.2049. hdl : 10356/137201 . S2CID  61039449.