Método espectral
Los métodos espectrales son una clase de técnicas utilizadas en matemáticas aplicadas y computación científica para resolver numéricamente ciertas ecuaciones diferenciales. La idea es escribir la solución de la ecuación diferencial como una suma de ciertas "funciones base" (por ejemplo, como una serie de Fourier que es una suma de sinusoides) y luego elegir los coeficientes en la suma para satisfacer la ecuación diferencial lo mejor posible.
Los métodos espectrales y los métodos de elementos finitos están estrechamente relacionados y se basan en las mismas ideas; la principal diferencia entre ellos es que los métodos espectrales usan funciones de base que generalmente son distintas de cero en todo el dominio, mientras que los métodos de elementos finitos usan funciones de base que son distintas de cero solo en pequeños subdominios (soporte compacto). En consecuencia, los métodos espectrales conectan variables globalmente mientras que los elementos finitos lo hacen localmente. En parte por esta razón, los métodos espectrales tienen excelentes propiedades de error, con la llamada "convergencia exponencial" siendo lo más rápido posible, cuando la solución es fluida. Sin embargo, no se conocen resultados tridimensionales de captura de choque espectral de dominio único (las ondas de choque no son uniformes). En la comunidad de elementos finitos, un método en el que el grado de los elementos es muy alto o aumenta a medida que aumenta el parámetro de cuadrícula h a veces se denomina método de elementos espectrales.
Los métodos espectrales se pueden utilizar para resolver ecuaciones diferenciales (EDP, ODE, valores propios, etc.) y problemas de optimización. Cuando se aplican métodos espectrales a PDE dependientes del tiempo, la solución generalmente se escribe como una suma de funciones básicas con coeficientes dependientes del tiempo; sustituyendo esto en la EDO se obtiene un sistema de EDO en los coeficientes que se puede resolver usando cualquier método numérico para las EDO. Los problemas de valores propios para ODE se convierten de manera similar en problemas de valores propios de matrices.
Los métodos espectrales fueron desarrollados en una larga serie de artículos por Steven Orszag a partir de 1969, incluidos, entre otros, métodos de series de Fourier para problemas de geometría periódica, métodos espectrales polinómicos para problemas de geometría finita e ilimitada, métodos pseudoespectrales para problemas altamente no lineales y métodos de iteración espectral para una solución rápida de problemas de estado estacionario. La implementación del método espectral normalmente se logra ya sea con la colocación o un enfoque de Galerkin o Tau. Para problemas muy pequeños, el método espectral es único en el sentido de que las soluciones se pueden escribir simbólicamente, proporcionando una alternativa práctica a las soluciones en serie para ecuaciones diferenciales.
Los métodos espectrales pueden ser computacionalmente menos costosos y más fáciles de implementar que los métodos de elementos finitos; brillan mejor cuando se busca una alta precisión en dominios simples con soluciones fluidas. Sin embargo, debido a su naturaleza global, las matrices asociadas con el cálculo de pasos son densas y la eficiencia computacional se verá afectada rápidamente cuando haya muchos grados de libertad (con algunas excepciones, por ejemplo, si las aplicaciones de matrices se pueden escribir como transformadas de Fourier). Para problemas más grandes y soluciones no uniformes, los elementos finitos generalmente funcionarán mejor debido a las matrices dispersas y al mejor modelado de discontinuidades y curvas cerradas.
Ejemplos de métodos espectrales
Un ejemplo concreto y lineal
Aquí presumimos una comprensión del cálculo multivariado básico y la serie Fourier. Si g()x,Sí.){displaystyle g(x,y)} es una función conocida de valor complejo de dos variables reales, y g es periódica en x y y (es decir, g()x,Sí.)=g()x+2π π ,Sí.)=g()x,Sí.+2π π ){displaystyle g(x,y)=g(x+2piy)=g(x,y+2pi)}) entonces estamos interesados en encontrar una función f()x,Sí.Así que...
- ()∂ ∂ 2∂ ∂ x2+∂ ∂ 2∂ ∂ Sí.2)f()x,Sí.)=g()x,Sí.)para todosx,Sí.{displaystyle left {frac {partial ^{2}{partial x^{2}}}+{frac {partial ^{2}}{partial y^{2}}}right)f(x,y)=g(x,y)quad {text{for all }x,y}
donde la expresión de la izquierda denota las segundas derivadas parciales de f en x y y, respectivamente. Esta es la ecuación de Poisson, y puede interpretarse físicamente como algún tipo de problema de conducción de calor, o un problema de teoría potencial, entre otras posibilidades.
Si escribimos f y g en serie de Fourier:
- f=.. aj,kei()jx+kSí.){displaystyle f=:sum a_{j,k}e^{i(jx+ky)}
- g=.. bj,kei()jx+kSí.){displaystyle g=:sum b_{j,k}e^{i(jx+ky)}
y sustituimos en la ecuación diferencial, obtenemos esta ecuación:
- .. − − aj,k()j2+k2)ei()jx+kSí.)=.. bj,kei()jx+kSí.){displaystyle sum -a_{j,k}(j^{2}+k^{2})e^{i(jx+ky)}=sum b_{j,k}e^{i(jx+ky)}}
Hemos intercambiado diferenciación parcial con una suma infinita, lo cual es legítimo si asumimos por ejemplo que f tiene una segunda derivada continua. Por el teorema de unicidad para las expansiones de Fourier, debemos igualar los coeficientes de Fourier término por término, dando
- aj,k=− − bj,kj2+k2{displaystyle a_{j,k}=-{frac {b_{j,k} {j^{2}+k^{2}}} {}}} {cH}}}}}} {c}}}}}}}}} {cH}}}}}}}}}}} {cH}}}}}}}}}}}}} {cH}}}}}}}}}}}}}}}}
()*)
que es una fórmula explícita para los coeficientes de Fourier aj,k.
Con condiciones de frontera periódicas, la ecuación de Poisson posee una solución solo si b0,0 = 0. Por lo tanto, podemos elegir libremente a0,0 que será igual a la media de la resolución. Esto corresponde a elegir la constante de integración.
Para convertir esto en un algoritmo, sólo finitamente muchas frecuencias se resuelven. Esto introduce un error que puede demostrar ser proporcional a hn{displaystyle h^{n}, donde h:=1/n{displaystyle h:=1/n} y n{displaystyle n} es la frecuencia más alta tratada.
Algoritmo
- Computar la transformación Fourier (bj,k) de g.
- Computar la transformación Fourier (aj,k) de f a través de la fórmula (*).
- Computación f tomando una inversa transformación Fourier de (aj,k).
Dado que solo estamos interesados en una ventana finita de frecuencias (de tamaño n, digamos), esto se puede hacer usando un algoritmo de transformada rápida de Fourier. Por lo tanto, globalmente el algoritmo se ejecuta en tiempo O(n log n).
Ejemplo no lineal
Queremos resolver el problema forzado, transitorio y no lineal de Burgers' ecuación utilizando un enfoque espectral.
Dado u()x,0){displaystyle u(x,0)} sobre el dominio periódico x▪ ▪ [0,2π π ){displaystyle xin left[0,2piright]}, encontrar u▪ ▪ U{displaystyle uin {fn\\fnMitcal {U}} tales que
- 0}" xmlns="http://www.w3.org/1998/Math/MathML">∂ ∂ tu+u∂ ∂ xu=*** *** ∂ ∂ xxu+fО О x▪ ▪ [0,2π π ),О О t■0{displaystyle partial _{t}u+upartial _{x}u=rho partial _{xx}u+fquad forall xin left[0,2piright]
0" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/e0f0089fa0f063e144e70d33d2eaa65694929286" style="vertical-align: -0.838ex; width:46.235ex; height:2.843ex;"/>
donde ρ es el coeficiente de viscosidad. En forma conservadora débil esto se convierte en
- 0}" xmlns="http://www.w3.org/1998/Math/MathML">.∂ ∂ tu,v.=.∂ ∂ x()− − 12u2+*** *** ∂ ∂ xu),v.+.f,v.О О v▪ ▪ V,О О t■0{displaystyle leftlangle partial _{t}u,vrightrangle =leftlangle partial _{x}left(-{frac {1}{2}u^{2}+rho partial _{x}uright),vrightrangle +leftlangle f,vrightrangle quad forall vin {mathcal {V}},forall t confianza0}
0}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/4377e19bfda3863129cfea9ec78e67d7cc7ce061" style="vertical-align: -2.505ex; width:61.696ex; height:6.176ex;"/>
donde sigue la notación del producto interno. Integrando por partes y usando subvenciones de periodicidad
- 0.}" xmlns="http://www.w3.org/1998/Math/MathML">.. ∂ ∂ tu,v.. =.12u2− − *** *** ∂ ∂ xu,∂ ∂ xv.+.f,v.О О v▪ ▪ V,О О t■0.{displaystyle langle partial _{t}u,vrangle =leftlangle {frac {1}{2}u^{2}-rho partial _{x}u,partial ¿Qué? +leftlangle f,vrightrangle quad forall vin {mathcal {V}},forall t confianza0.}
0.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/60121e5fa4005afc74d83ef20e02a5f10cad5177" style="vertical-align: -2.505ex; width:56.339ex; height:6.176ex;"/>
Para aplicar el método de Fourier-Galerkin, elija ambos
- UN:={}u:u()x,t)=.. k=− − N/2N/2− − 1u^ ^ k()t)eikx}{displaystyle {mathcal {fnK}}=left{u:u(x,t)=sum ¿Por qué?
y
- VN:=lapso {}eikx:k▪ ▪ − − N/2,...... ,N/2− − 1}{displaystyle {mathcal {fn} {fn}fn}fnK}fnK}fn}fn} -N/2, 'dotsN/2-1 'right'
Donde u^ ^ k()t):=12π π .. u()x,t),eikx.. {displaystyle {hat {u}_{k}(t):={frac {1}{2pi}langle u(x,t),e^{ikx}rangle }. Esto reduce el problema para encontrar u▪ ▪ UN{displaystyle uin {fn\\fnMitcal {U}} {N}} tales que
- 0.}" xmlns="http://www.w3.org/1998/Math/MathML">.. ∂ ∂ tu,eikx.. =.12u2− − *** *** ∂ ∂ xu,∂ ∂ xeikx.+.f,eikx.О О k▪ ▪ {}− − N/2,...... ,N/2− − 1},О О t■0.{displaystyle langle partial ########################## {1}{2}}u^{2}-rho partial _{x}u,partial ¿Por qué?
0.}" aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/f4d0bd2dbe49c6577631ee066b01fcb3c990ce7b" style="vertical-align: -2.505ex; width:85.414ex; height:6.176ex;"/>
Usando la relación ortogonalidad .. eilx,eikx.. =2π π δ δ lk{displaystyle langle e^{ilx},e^{ikx}rangle =2pi delta ¿Qué? Donde δ δ lk{displaystyle delta _{lk} es el Kronecker delta, simplificamos los tres términos anteriores para cada uno k{displaystyle k} para ver
- .∂ ∂ tu,eikx.=.∂ ∂ t.. lu^ ^ leilx,eikx.=... l∂ ∂ tu^ ^ leilx,eikx.=2π π ∂ ∂ tu^ ^ k,.f,eikx.=... lf^ ^ leilx,eikx.=2π π f^ ^ k,y.12u2− − *** *** ∂ ∂ xu,∂ ∂ xeikx.=.12().. pu^ ^ peipx)().. qu^ ^ qeiqx)− − *** *** ∂ ∂ x.. lu^ ^ leilx,∂ ∂ xeikx.=.12.. p.. qu^ ^ pu^ ^ qei()p+q)x,ikeikx.− − .*** *** i.. llu^ ^ leilx,ikeikx.=− − ik2... p.. qu^ ^ pu^ ^ qei()p+q)x,eikx.− − *** *** k... llu^ ^ leilx,eikx.=− − iπ π k.. p+q=ku^ ^ pu^ ^ q− − 2π π *** *** k2u^ ^ k.{displaystyle {begin{aligned}leftlangle partial _{t}u,e^{ikx}rightrangle ¿Qué? =leftlangle sum _{l}partial ¿Qué? =2pi partial ¿Por qué? ¿Qué? ¿Por qué? {fnMicroc}left(sum _{p} {hat {u}_{p}e^{ipx}right)left(sum _{q}{hat {u}_{q}e^{iqx}right)-rho partial _{x}sum _{l} {hat {u}_{l}e^{ilx},partial ################################################################################################################################################################################################################################################################ ¿Qué? {fnK} {fnMicrosoft Sans Serif} -leftlangle rho isum _{l}l{hat {u}_{l}e^{ilx},ike^{ikx}rightrangle\\\fnt=-{frac {}{2}leftlangle {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} - ¿Qué? \fnc=-ipi ksum - ¿Qué? {fnh} {fnh}} {fnh}} {fnh}} {fnh}} {fn}}} {fn}} {fn}}} {fn}}}} {fnfn}} {fnf}}}} {fnfnf}}} {f}}}} {f}}}}} {f}}}}}} {f}}}}f}}}}}}}}}}}}}}}} {f} {f} {f} {f} {f}}} {f}}}}}}}} {f} {f} {f} {f}}}} {f}}}}} {f}f}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} {U}_{q}-2pi rho {}k}{2}{hat {fnMicrosoft Sans Serif}
Montar los tres términos para cada k{displaystyle k} para obtener
- 0.}" xmlns="http://www.w3.org/1998/Math/MathML">2π π ∂ ∂ tu^ ^ k=− − iπ π k.. p+q=ku^ ^ pu^ ^ q− − 2π π *** *** k2u^ ^ k+2π π f^ ^ kk▪ ▪ {}− − N/2,...... ,N/2− − 1},О О t■0.{displaystyle 2pi partial ¿Qué? {fnK}=-ipi} ksum - ¿Qué? {fnh} {fnK} {fnK}}} {fnK}}} {fn}} {fn}}} {fn}} {fn}}}}} {fn}} {fn}} {f}}}} {f}}}}}} {fnf}} {f}}}}} {f}}}} {f}}}}}}}} {f}}}}}} {f}}}}} {f}}} {f} {f}} {f}} {f}}}}}}}}}}}}}}}} {f}} {f} {f} {f}}}} {f}}}} {f} {f}}}f}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} {f}}} {U}_{q}-2pi rho {}k}{2}{hat {fnK}+2pi} {f} {f} {f}f}f}f}f} {f}f} {f}f}f}f}f} {f}}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f}f} kin left{-N/2,dotsN/2-1right},forall t título0.}
0. " aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/e4d577e7df77c99503d824db8de1b1800238dfc2" style="vertical-align: -3.338ex; width:81.973ex; height:5.843ex;"/>
Dividir a través de 2π π {displaystyle 2pi}, finalmente llegamos
- 0.}" xmlns="http://www.w3.org/1998/Math/MathML">∂ ∂ tu^ ^ k=− − ik2.. p+q=ku^ ^ pu^ ^ q− − *** *** k2u^ ^ k+f^ ^ kk▪ ▪ {}− − N/2,...... ,N/2− − 1},О О t■0.{displaystyle partial _{t}{hat {fnK}=- {fnMicroc} {ik}{2}sum - ¿Qué? {fnh} {fnK} {fnK}}} {fnK}}} {fn}} {fn}}} {fn}} {fn}}}}} {fn}} {fn}} {f}}}} {f}}}}}} {fnf}} {f}}}}} {f}}}} {f}}}}}}}} {f}}}}}} {f}}}}} {f}}} {f} {f}} {f}} {f}}}}}}}}}}}}}}}} {f}} {f} {f} {f}}}} {f}}}} {f} {f}}}f}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}} {f}}} {U}_{q}-rho {}k}{2}{hat {fnK} {fnK} {fnK}fnK} {fnK}}} {fnK}} {f}}fn}fnfnfnfn}fnfnfnf}}fnf}}\fnfnfnfnfnf}}}}}}}}}\\fn\fnfn\\\fnfnfnfn}}}}}}}}\\\\\\fnf}}}}}}}}}}}}\\\fn}}}}}}}}}}}}}\\\\\\\fn}}}}}}}}}}}}}}}}}}}}}}\\\\\\\\\ {f}_{k}quad kin left{-N/2,dotsN/2-1right},forall t título0.}
0. " aria-hidden="true" class="mwe-math-fallback-image-inline" src="https://wikimedia.org/api/rest_v1/media/math/render/svg/d86efb28aeeeaaf23f630933e3dee3ff7792beb6" style="vertical-align: -3.338ex; width:73.994ex; height:6.843ex;"/>
Con Fourier transformado las condiciones iniciales u^ ^ k()0){fnMicrosoft Sans Serif} y forzamiento f^ ^ k()t){displaystyle {hat {f}_{k}(t)}, este sistema acoplado de ecuaciones diferenciales ordinarias se puede integrar en el tiempo (usando, por ejemplo, una técnica Runge Kutta) para encontrar una solución. El término no lineal es una convolución, y hay varias técnicas basadas en la transformación para evaluarlo eficientemente. Vea las referencias de Boyd y Canuto et al. para más detalles.
Una relación con el método del elemento espectral
Uno puede demostrar que si g{displaystyle g} es infinitamente diferenciable, entonces el algoritmo numérico usando Fast Fourier Transforms convergerá más rápido que cualquier polinomio en el tamaño de la cuadrícula h. Es decir, para cualquier n título, hay un <math alttext="{displaystyle C_{n}Cn.JUEGO JUEGO {displaystyle C_{n}<img alt="{displaystyle C_{n} tal que el error es menos que Cnhn{displaystyle C_{n}h^{n} para todos los valores suficientemente pequeños h{displaystyle h}. Decimos que el método espectral es de orden n{displaystyle n}Por cada n título.
Debido a que un método de elementos espectrales es un método de elementos finitos de muy alto orden, existe una similitud en las propiedades de convergencia. Sin embargo, mientras que el método espectral se basa en la descomposición propia del problema de valor límite particular, el método de elementos finitos no utiliza esa información y funciona para problemas de valor límite elípticos arbitrarios.