Método de diferencias finitas

En el análisis numérico , los métodos de diferencias finitas ( FDM ) 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 dominio temporal (si corresponde) se discretizan , o se dividen en un número finito de intervalos, y los valores de la solución en los puntos finales de los intervalos se aproximan mediante la resolución de 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 pueden resolver mediante técnicas de álgebra matricial . Las computadoras modernas pueden realizar estos cálculos de álgebra lineal de manera eficiente, lo que, junto con su relativa facilidad de implementación, ha llevado al uso generalizado de FDM en el análisis numérico moderno. ^[1] Hoy en día, los FDM 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, por el teorema de Taylor la expansión en serie de Taylor se da como$${\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 muestra el proceso para derivar una aproximación para la primera derivada de la función f truncando primero el polinomio de Taylor más el resto: Dividiendo por h se obtiene: Resolviendo para :$${\displaystyle f(x_{0}+h)=f(x_{0})+f'(x_{0})h+R_{1}(x).}$$$${\displaystyle {f(x_{0}+h) \over h}={f(x_{0}) \over h}+f'(x_{0})+{R_{1}(x) \over h}}$$${\displaystyle f'(x_{0})}$$${\displaystyle f'(x_{0})={f(x_{0}+h)-f(x_{0}) \over h}-{R_{1}(x) \over h}.}$$

Suponiendo que es suficientemente pequeño, la aproximación de la primera derivada de f es:${\displaystyle R_{1}(x)}$$${\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: excepto por el límite hacia cero (el método recibe su nombre de esto).$${\displaystyle f'(x_{0})=\lim _{h\to 0}{\frac {f(x_{0}+h)-f(x_{0})}{h}}.}$$

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 , la pérdida de precisión debido al redondeo informático de cantidades decimales, y el error de truncamiento o error de discretización , 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 utilizar un método de diferencias finitas para aproximar la solución de un problema, primero se debe discretizar el dominio del problema. Esto se hace generalmente dividiendo el dominio en una cuadrícula uniforme (ver imagen). Esto significa que los métodos de diferencias finitas producen conjuntos de aproximaciones numéricas discretas a la derivada, a menudo de forma "escalonada en el tiempo".

Una expresión de interés general es el error de truncamiento local de un método. Expresado típicamente 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 si se refiere al valor exacto y a la aproximación numérica. El término restante del polinomio de Taylor se puede utilizar para analizar el error de truncamiento local . Usando la forma de Lagrange del resto del polinomio de Taylor para , se puede descubrir cuál es el término dominante del error de truncamiento local. Por ejemplo, usando nuevamente la fórmula de diferencia hacia adelante para la primera derivada, sabiendo que , y con alguna manipulación algebraica, esto lleva a y además notando que la cantidad de la izquierda es la aproximación del método de diferencia finita y que la cantidad de la derecha es la cantidad exacta de interés más un resto, claramente ese resto es el error de truncamiento local. Una expresión final de este ejemplo y su orden es:${\displaystyle f'(x_{i})-f'_{i}}$${\displaystyle f'(x_{i})}$${\displaystyle f'_{i}}$${\displaystyle f(x_{0}+h)}$$${\displaystyle R_{n}(x_{0}+h)={\frac {f^{(n+1)}(\xi )}{(n+1)!}}(h)^{n+1}\,,\quad x_{0}<\xi <x_{0}+h,}$$${\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},}$$$${\displaystyle {\frac {f(x_{0}+ih)-f(x_{0})}{ih}}=f'(x_{0})+{\frac {f''(\xi )}{2!}}ih,}$$$${\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 (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 el 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 a menudo para determinar la estabilidad del modelo numérico. ^{[3] [4] [5] [6]}

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

Considere 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 por diferencias finitas. Primero, particione el dominio en el espacio usando una malla y en el tiempo usando una malla . Suponga una partición uniforme tanto en el espacio como en el tiempo, por lo que la diferencia entre dos puntos espaciales consecutivos será h y entre dos puntos temporales consecutivos será k . Los puntos${\displaystyle x_{0},\dots ,x_{J}}$${\displaystyle t_{0},\dots ,t_{N}}$

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

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

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

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k\Delta t}}={\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 de calor unidimensional .

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

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

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

Entonces, con esta relación de recurrencia, y conociendo los valores en el tiempo n , se pueden obtener los valores correspondientes en el tiempo n +1. y deben reemplazarse por las condiciones de borde, en este ejemplo ambas son 0.${\displaystyle u_{0}^{n}}$${\displaystyle u_{J}^{n}}$

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

Usando 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 (el método de tiempo hacia atrás, espacio centrado "BTCS") se obtiene la ecuación de recurrencia:${\displaystyle x_{j}}$

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k\Delta t}}={\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 de calor unidimensional .

Resolviendo un sistema de ecuaciones lineales se puede obtener :${\displaystyle u_{j}^{n+1}}$

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

El esquema es siempre numéricamente estable y convergente, pero normalmente es más intensivo numéricamente 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 en el paso de tiempo y cuadráticos en el paso de espacio:$${\displaystyle \Delta u=O(k)+O(h^{2}).}$$

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

$${\displaystyle {\frac {u_{j}^{n+1}-u_{j}^{n}}{k\Delta t}}={\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 Crank-Nicolson .

Resolviendo un sistema de ecuaciones lineales se puede obtener :${\displaystyle u_{j}^{n+1}}$

$${\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 es siempre numéricamente estable y convergente, pero normalmente es más intensivo numéricamente, ya que requiere resolver un sistema de ecuaciones numéricas en cada paso temporal. Los errores son cuadráticos tanto en el paso temporal como en el paso espacial:$${\displaystyle \Delta u=O(k^{2})+O(h^{2}).}$$

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 exige 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 recursos numéricos.

A continuación se muestra un ejemplo. Las figuras siguientes presentan las soluciones obtenidas con 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).}$$

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

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

El caso 2D muestra todas las características del caso n-dimensional más general. Cada segunda derivada parcial debe aproximarse de manera similar al caso 1D , que generalmente se da mediante la siguiente plantilla$${\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}}}$$$${\displaystyle \Delta _{h}={\frac {1}{h^{2}}}{\begin{bmatrix}&1\\1&-4&1\\&1\end{bmatrix}}\,.}$$

La consistencia de la aproximación antes mencionada se puede demostrar para funciones altamente regulares, como . La afirmación es${\displaystyle u\in C^{4}(\Omega )}$$${\displaystyle \Delta u-\Delta _{h}u={\mathcal {O}}(h^{2})\,.}$$

Para demostrar esto, es necesario sustituir las expansiones de la serie de Taylor hasta el orden 3 en el operador de Laplace discreto.

De manera similar a las funciones subarmónicas continuas, se pueden definir funciones subarmónicas para aproximaciones de diferencias finitas.${\displaystyle u_{h}}$$${\displaystyle -\Delta _{h}u_{h}\leq 0\,.}$$

Se puede definir una plantilla general de tipo positivo mediante$${\displaystyle {\begin{bmatrix}&\alpha _{N}\\\alpha _{W}&-\alpha _{C}&\alpha _{E}\\&\alpha _{S}\end{bmatrix}}\,,\quad \alpha _{i}>0\,,\quad \alpha _{C}=\sum _{i\in \{N,E,S,W\}}\alpha _{i}\,.}$$

Si es subarmónico (discreto), entonces se cumple la siguiente propiedad de valor medio donde la aproximación se evalúa en puntos de la cuadrícula y se supone que la plantilla es de tipo positivo.${\displaystyle u_{h}}$$${\displaystyle u_{h}(x_{C})\leq {\frac {\sum _{i\in \{N,E,S,W\}}\alpha _{i}u_{h}(x_{i})}{\sum _{i\in \{N,E,S,W\}}\alpha _{i}}}\,,}$$

Una propiedad de valor medio similar también se aplica al caso continuo.

Para una función subarmónica (discreta) se cumple lo siguiente , donde son discretizaciones del dominio continuo , respectivamente el límite .${\displaystyle u_{h}}$$${\displaystyle \max _{\Omega _{h}}u_{h}\leq \max _{\partial \Omega _{h}}u_{h}\,,}$$${\displaystyle \Omega _{h},\partial \Omega _{h}}$${\displaystyle \Omega }$${\displaystyle \partial \Omega }$

Un principio de máximo similar también se aplica al caso continuo.

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 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 diferencia central en el interior con plantillas de contorno unilaterales cuidadosamente elegidas diseñadas para imitar la integración por partes en el entorno discreto. Usando la técnica SAT, las condiciones de contorno de la EDP se imponen débilmente, donde los valores de contorno son "arrastrados" hacia las condiciones deseadas en lugar de cumplirse exactamente. Si los parámetros de ajuste (inherentes a la técnica SAT) se eligen correctamente, el sistema resultante de EDO exhibirá un comportamiento de energía similar al de la EDP continua, es decir, el sistema no tiene 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 incluye partes del eje imaginario, como el método de Runge-Kutta de cuarto orden . Esto hace que la técnica SAT sea un método atractivo para imponer condiciones de contorno para métodos de diferencias finitas de orden superior, en contraste, por ejemplo, con el método de inyección, que normalmente no será estable si se utilizan operadores de diferenciación de orden superior.

• 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) [1]. Contiene una breve introducción orientada a la ingeniería a 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.