Método de diferencias finitas

En análisis numérico , los métodos de diferencias finitas ( MDF ) son una clase de técnicas numéricas para resolver ecuaciones diferenciales mediante la aproximación de derivadas con diferencias finitas . Tanto el dominio espacial como el temporal (si corresponde) se discretizan , es decir, se dividen en un número finito de intervalos, y los valores de la solución en los extremos de los intervalos se aproximan resolviendo ecuaciones algebraicas que contienen diferencias finitas y valores de puntos cercanos.

Los métodos de diferencias finitas convierten ecuaciones diferenciales ordinarias (EDO) o ecuaciones diferenciales parciales (EDP), que pueden ser no lineales , en un sistema de ecuaciones lineales que se puede resolver mediante técnicas de álgebra matricial . Las computadoras modernas pueden realizar estos cálculos de álgebra lineal de manera eficiente, y esto, junto con su relativa facilidad de implementación, ha llevado al uso generalizado de los métodos de diferencias finitas en el análisis numérico moderno. ^{[ 1 ]} Hoy en día, los métodos de diferencias finitas son uno de los enfoques más comunes para la solución numérica de EDP, junto con los métodos de elementos finitos . ^{[ 1 ]}

Para una función n veces diferenciable, según el teorema de Taylor , el desarrollo en serie de Taylor viene dado por: $${\displaystyle f(x_{0}+h)=f(x_{0})+{\frac {f'(x_{0})}{1!}}h+{\frac {f^{(2)}(x_{0})}{2!}}h^{2}+\cdots +{\frac {f^{(n)}(x_{0})}{n!}}h^{n}+R_{n}(x),}$$

Donde n ! denota el factorial de n , y R _n ( x ) es un término restante, que denota la diferencia entre el polinomio de Taylor de grado n y la función original.

A continuación se describe el proceso para obtener una aproximación de la primera derivada de la función f truncando primero el polinomio de Taylor más el resto: $${\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x).}$$ Dividiendo por h obtenemos: $${\displaystyle {f(x_{0}+h) \over h}={f(x_{0}) \over h}+f'(x_{0})+{R_{1}(x) \over h}}$$ Resolver para${\displaystyle f'(x_{0})}$: $${\displaystyle f'(x_{0})={f(x_{0}+h)-f(x_{0}) \over h}-{R_{1}(x) \over h}.}$$

Suponiendo que${\displaystyle R_{1}(x)}$Si es suficientemente pequeño, la aproximación de la primera derivada de f es: $${\displaystyle f'(x_{0})\approx {f(x_{0}+h)-f(x_{0}) \over h}.}$$

Esto es similar a la definición de derivada, que es: $${\displaystyle f'(x_{0})=\lim _{h\to 0}{\frac {f(x_{0}+h)-f(x_{0})}{h}}.}$$ excepto por el límite hacia cero (el método recibe su nombre de este).

El error en la solución de un método se define como la diferencia entre la aproximación y la solución analítica exacta. Las dos fuentes de error en los métodos de diferencias finitas son el error de redondeo , que es la pérdida de precisión debida al redondeo informático de cantidades decimales, y el error de truncamiento o error de discretización , que es la diferencia entre la solución exacta de la ecuación diferencial original y la cantidad exacta suponiendo una aritmética perfecta (sin redondeo).

Para aproximar la solución de un problema mediante un método de diferencias finitas, primero se debe discretizar su dominio. Esto se suele hacer dividiendo el dominio en una cuadrícula uniforme (véase la imagen). De esta forma, los métodos de diferencias finitas generan conjuntos de aproximaciones numéricas discretas de la derivada, a menudo mediante un método de discretización temporal.

Una expresión de interés general es el error de truncamiento local de un método. Típicamente expresado usando la notación Big-O , el error de truncamiento local se refiere al error de una sola aplicación de un método. Es decir, es la cantidad${\displaystyle f'(x_{i})-f'_{i}}$si${\displaystyle f'(x_{i})}$se refiere al valor exacto y${\displaystyle f'_{i}}$a la aproximación numérica. El término restante del polinomio de Taylor se puede utilizar para analizar el error de truncamiento local . Utilizando la forma lagrangiana del resto del polinomio de Taylor para${\displaystyle f(x_{0}+h)}$, que es $${\displaystyle R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}(h)^{n+1}\,,\quad x_{0}<\xi <x_{0}+h,}$$ Se puede descubrir el término dominante del error de truncamiento local. Por ejemplo, utilizando nuevamente la fórmula de diferencias finitas hacia adelante para la primera derivada, sabiendo que${\displaystyle f(x_{i})=f(x_{0}+ih)}$, $${\displaystyle f(x_{0}+ih)=f(x_{0})+f'(x_{0})ih+{\frac {f''(\xi )}{2!}}(ih)^{2},}$$ y con alguna manipulación algebraica, esto conduce a $${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,}$$ y observando además que la cantidad de la izquierda es la aproximación obtenida mediante el método de diferencias finitas y que la cantidad de la derecha es la cantidad exacta de interés más un resto, que claramente corresponde al error de truncamiento local. Una expresión final de este ejemplo y su orden es: $${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+O(h).}$$

En este caso, el error de truncamiento local es proporcional al tamaño de los pasos. La calidad y la duración de la solución FDM simulada dependen de la selección de la ecuación de discretización y del tamaño de los pasos (de tiempo y espacio). La calidad de los datos y la duración de la simulación aumentan significativamente con un tamaño de paso menor. ^{[ 2 ]} Por lo tanto, es necesario un equilibrio razonable entre la calidad de los datos y la duración de la simulación para su uso práctico. Los pasos de tiempo grandes son útiles para aumentar la velocidad de simulación en la práctica. Sin embargo, los pasos de tiempo demasiado grandes pueden crear inestabilidades y afectar la calidad de los datos. ^{[ 3 ] [ 4 ]}

Los criterios de von Neumann y Courant-Friedrichs-Lewy se evalúan frecuentemente para determinar la estabilidad del modelo numérico. ^{[ 3 ] [ 4 ] [ 5 ] [ 6 ]}

Por ejemplo, consideremos la ecuación diferencial ordinaria. $${\displaystyle u'(x)=3u(x)+2.}$$ El método de Euler para resolver esta ecuación utiliza el cociente de diferencias finitas. $${\displaystyle {\frac {u(x+h)-u(x)}{h}}\approx u'(x)}$$ aproximar la ecuación diferencial sustituyéndola primero por${\displaystyle u'(x)}$luego aplicando un poco de álgebra (multiplicando ambos lados por h y luego sumando)${\displaystyle u(x)}$a ambos lados) para obtener $${\displaystyle u(x+h)\approx u(x)+h(3u(x)+2).}$$ La última ecuación es una ecuación de diferencias finitas, y al resolverla se obtiene una solución aproximada a la ecuación diferencial.

Consideremos la ecuación de calor normalizada en una dimensión, con condiciones de contorno de Dirichlet homogéneas.

$${\displaystyle {\begin{cases}U_{t}=U_{xx}\\U(0,t)=U(1,t)=0&{\text{(boundary condition)}}\\U(x,0)=U_{0}(x)&{\text{(initial condition)}}\end{cases}}}$$

Una forma de resolver numéricamente esta ecuación es aproximar todas las derivadas mediante diferencias finitas. Primero, se divide el dominio en el espacio utilizando una malla.${\displaystyle x_{0},\dots ,x_{J}}$y con el tiempo usando una malla${\displaystyle t_{0},\dots ,t_{N}}$Supongamos una partición uniforme tanto en el espacio como en el tiempo, de modo que la diferencia entre dos puntos espaciales consecutivos será h y entre dos puntos temporales consecutivos será k . Los puntos

$${\displaystyle u(x_{j},t_{n})=u_{j}^{n}}$$

representará la aproximación numérica de${\displaystyle u(x_{j},t_{n}).}$

Método explícito

Utilizando una diferencia hacia adelante en el tiempo${\displaystyle t_{n}}$y una diferencia central de segundo orden para la derivada espacial en la posición${\displaystyle x_{j}}$( FTCS ) proporciona la ecuación de recurrencia:

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}.}$$

Este es un método explícito para resolver la ecuación del calor unidimensional .

Uno puede obtener${\displaystyle u_{j}^{n+1}}$a partir de los otros valores de esta manera:

$${\displaystyle u_{j}^{n+1}=(1-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}$$

dónde${\displaystyle r=k/h^{2}.}$

Así pues, con esta relación de recurrencia, y conociendo los valores en el instante n , se pueden obtener los valores correspondientes en el instante n + 1.${\displaystyle u_{0}^{n}}$y${\displaystyle u_{J}^{n}}$deben ser reemplazados por las condiciones de contorno, en este ejemplo ambos son 0.

Se sabe que este método explícito es numéricamente estable y convergente siempre que...${\displaystyle r\leq 1/2}$. ^{[ 7 ]} Los errores numéricos son proporcionales al paso de tiempo y al cuadrado del paso espacial: $${\displaystyle \Delta u=O(k)+O(h^{2})}$$

Método implícito

Utilizando la diferencia hacia atrás en el tiempo${\displaystyle t_{n+1}}$y una diferencia central de segundo orden para la derivada espacial en la posición${\displaystyle x_{j}}$(El método de tiempo inverso y espacio centrado "BTCS") proporciona la ecuación de recurrencia:

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}.}$$

Este es un método implícito para resolver la ecuación del calor unidimensional .

Uno puede obtener${\displaystyle u_{j}^{n+1}}$a partir de la resolución de un sistema de ecuaciones lineales:

$${\displaystyle (1+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=u_{j}^{n}}$$

El esquema siempre es numéricamente estable y convergente, pero suele ser más intensivo computacionalmente que el método explícito, ya que requiere resolver un sistema de ecuaciones numéricas en cada paso de tiempo. Los errores son lineales con respecto al paso de tiempo y cuadráticos con respecto al paso espacial. $${\displaystyle \Delta u=O(k)+O(h^{2}).}$$

Método Crank - Nicolson

Finalmente, utilizando la diferencia central en el tiempo${\displaystyle t_{n+1/2}}$y una diferencia central de segundo orden para la derivada espacial en la posición${\displaystyle x_{j}}$("CTCS") proporciona la ecuación de recurrencia:

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k}}={\frac {1}{2}}\left({\frac {u_{j+1}^{n+1}-2u_{j}^{n+1}+u_{j-1}^{n+1}}{h^{2}}}+{\frac {u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{h^{2}}}\right).}$$

Esta fórmula se conoce como el método de Crank-Nicolson .

Uno puede obtener${\displaystyle u_{j}^{n+1}}$a partir de la resolución de un sistema de ecuaciones lineales:

$${\displaystyle (2+2r)u_{j}^{n+1}-ru_{j-1}^{n+1}-ru_{j+1}^{n+1}=(2-2r)u_{j}^{n}+ru_{j-1}^{n}+ru_{j+1}^{n}}$$

El esquema siempre es numéricamente estable y convergente, pero suele ser más intensivo computacionalmente, ya que requiere resolver un sistema de ecuaciones numéricas en cada paso de tiempo. Los errores son cuadráticos tanto en el paso de tiempo como en el paso espacial. $${\displaystyle \Delta u=O(k^{2})+O(h^{2}).}$$

Comparación

En resumen, el esquema de Crank - Nicolson suele ser el más preciso para pasos de tiempo pequeños. Para pasos de tiempo mayores, el esquema implícito funciona mejor, ya que requiere menos recursos computacionales. El esquema explícito es el menos preciso y puede ser inestable, pero también es el más fácil de implementar y el que requiere menos cálculos.

Aquí hay un ejemplo. Las figuras a continuación presentan las soluciones dadas por los métodos anteriores para aproximar la ecuación del calor.

$${\displaystyle U_{t}=\alpha U_{xx},\quad \alpha ={\frac {1}{\pi ^{2}}},}$$

con la condición de contorno

$${\displaystyle U(0,t)=U(1,t)=0.}$$

La solución exacta es

$${\displaystyle U(x,t)={\frac {1}{\pi ^{2}}}e^{-t}\sin(\pi x).}$$

Comparación de métodos de diferencias finitas

Método explícito ( no estable)

Método implícito (estable)

Método de Crank-Nicolson (estable)

El operador de Laplace (continuo) en${\displaystyle n}$-dimensiones viene dada por${\displaystyle \Delta u(x)=\sum _{i=1}^{n}\partial _{i}^{2}u(x)}$El operador de Laplace discreto${\displaystyle \Delta _{h}u}$depende de la dimensión${\displaystyle n}$.

En 1D el operador de Laplace se aproxima como $${\displaystyle \Delta u(x)=u''(x)\approx {\frac {u(x-h)-2u(x)+u(x+h)}{h^{2}}}=:\Delta _{h}u(x)\,.}$$ Esta aproximación se suele expresar mediante la siguiente plantilla.$${\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}1&-2&1\end{bmatrix}}}$$ y que representa una matriz tridiagonal simétrica. Para una cuadrícula equidistante se obtiene una matriz de Toeplitz .

El caso 2D muestra todas las características del caso n-dimensional más general. Cada segunda derivada parcial necesita ser aproximada de forma similar al caso 1D. $${\displaystyle {\begin{aligned}\Delta u(x,y)&=u_{xx}(x,y)+u_{yy}(x,y)\\&\approx {\frac {u(x-h,y)-2u(x,y)+u(x+h,y)}{h^{2}}}+{\frac {u(x,y-h)-2u(x,y)+u(x,y+h)}{h^{2}}}\\&={\frac {u(x-h,y)+u(x+h,y)-4u(x,y)+u(x,y-h)+u(x,y+h)}{h^{2}}}\\&=:\Delta _{h}u(x,y)\,,\end{aligned}}}$$ que normalmente se proporciona mediante la siguiente plantilla.$${\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}&1\\1&-4&1\\&1\end{bmatrix}}\,.}$$

El método SBP-SAT ( suma por partes - término de aproximación simultánea ) es una técnica estable y precisa para discretizar e imponer condiciones de contorno de una ecuación diferencial parcial lineal bien planteada utilizando diferencias finitas de alto orden. ^{[ 8 ] [ 9 ]}

El método se basa en diferencias finitas donde los operadores de diferenciación exhiben propiedades de suma por partes . Típicamente, estos operadores consisten en matrices de diferenciación con plantillas de diferencias centrales en el interior y plantillas de frontera unilaterales cuidadosamente elegidas, diseñadas para imitar la integración por partes en el entorno discreto. Mediante la técnica SAT, las condiciones de frontera de la EDP se imponen débilmente, donde los valores de frontera se "ajustan" hacia las condiciones deseadas en lugar de cumplirse exactamente. Si los parámetros de ajuste (inherentes a la técnica SAT) se eligen adecuadamente, el sistema resultante de EDO exhibirá un comportamiento energético similar al de la EDP continua, es decir, el sistema no presenta crecimiento de energía no físico. Esto garantiza la estabilidad si se utiliza un esquema de integración con una región de estabilidad que incluya partes del eje imaginario, como el método de Runge-Kutta de cuarto orden . Esto convierte a la técnica SAT en un método atractivo para imponer condiciones de frontera para métodos de diferencias finitas de orden superior, a diferencia, por ejemplo, del método de inyección, que normalmente no será estable si se utilizan operadores de diferenciación de alto orden. La formulación aditiva de los términos SAT puede simplificar las implementaciones eficientes de este método en GPU u otras arquitecturas modernas de alto rendimiento al tratar las condiciones de contorno como términos de corrección separados, a diferencia de los métodos de diferencias finitas que imponen condiciones de contorno mediante la inyección directa de datos de contorno. ^{[ 10 ]}

• KW Morton y DF Mayers, Solución numérica de ecuaciones diferenciales parciales: una introducción . Cambridge University Press, 2005.
• Autar Kaw y E. Eric Kalu, Métodos numéricos con aplicaciones , (2008). Contiene una breve introducción orientada a la ingeniería al FDM (para EDO) en el Capítulo 08.07 .
• John Strikwerda (2004). Esquemas de diferencias finitas y ecuaciones diferenciales parciales (2.ª  ed.). SIAM. ISBN 978-0-89871-639-9.
• Smith, GD (1985), Solución numérica de ecuaciones diferenciales parciales: métodos de diferencias finitas, 3.ª ed. , Oxford University Press
• Peter Olver (2013). Introducción a las ecuaciones diferenciales parciales . Springer. Capítulo 5: Diferencias finitas. ISBN 978-3-319-02099-0..
• Randall J. LeVeque , Métodos de diferencias finitas para ecuaciones diferenciales ordinarias y parciales , SIAM, 2007.
• Sergey Lemeshevsky, Piotr Matus, Dmitriy Poliakov (Eds): "Esquemas exactos de diferencias finitas", De Gruyter (2016). DOI: https://doi.org/10.1515/9783110491326 .
• Mikhail Shashkov: Métodos de diferencias finitas conservativos en cuadrículas generales , CRC Press, ISBN 0-8493-7375-1 (1996).