Método de diferencias finitas

Clase de técnicas numéricas

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]

Derive el cociente de diferencias del polinomio de Taylor

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}}.}$$

Precisión y orden

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).

El método de diferencias finitas se basa en la discretización de una función en una cuadrícula.

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 pasos más pequeños. ^[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]}

Ejemplo: ecuación diferencial ordinaria

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).}$$

Ejemplo: La ecuación del calor

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}).}$

Método explícito

La plantilla para el método explícito más común para la ecuación de calor.

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}}={\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/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})}$$

Método implícito

La plantilla del método implícito.

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}}={\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}).}$$

Método de Crank-Nicolson

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}}={\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 .

La plantilla de 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 de tiempo. Los errores son cuadráticos tanto en el paso de tiempo como en el paso de espacio: $${\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 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).}$$

Comparación de métodos de diferencias finitas

Ejemplo: El operador de Laplace

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}}\,.}$$

Consistencia

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.

Propiedades

Subarmónico

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\,.}$$

Valor medio

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.

Principio máximo

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

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 "tirados" 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.

Véase también

• Método de elementos finitos
• Diferencia finita
• Dominio del tiempo de diferencias finitas
• Método de diferencias infinitas
• Plantilla (análisis numérico)
• Coeficientes de diferencias finitas
• Plantilla de cinco puntos
• Teorema de Lax-Richtmyer
• Métodos de diferencias finitas para la determinación del precio de opciones
• Esquema de diferenciación en contra del viento para convección
• Esquema de diferenciación central
• Ecuación de Poisson discreta
• Operador discreto de Laplace

Referencias

1. Christian Grossmann; Hans-G. Roos; Martin Stynes ​​(2007). Tratamiento numérico de ecuaciones diferenciales parciales . Springer Science & Business Media. pág. 23. ISBN 978-3-540-71584-9.
2. Arieh Iserles (2008). Un primer curso sobre el análisis numérico de ecuaciones diferenciales . Cambridge University Press. pág. 23. ISBN 9780521734905.
3. Hoffman JD; Frankel S (2001). Métodos numéricos para ingenieros y científicos . CRC Press, Boca Raton.
4. Jaluria Y; Atluri S (1994). "Transferencia de calor computacional". Mecánica computacional . 14 (5): 385–386. Código Bibliográfico :1994CompM..14..385J. doi :10.1007/BF00377593. S2CID  119502676.
5. Majumdar P (2005). Métodos computacionales para transferencia de calor y masa (1.ª ed.). Taylor y Francis, Nueva York.
6. Smith GD (1985). Solución numérica de ecuaciones diferenciales parciales: métodos de diferencias finitas (3.ª ed.). Oxford University Press.
7. Crank, J. Las matemáticas de la difusión . 2.ª edición, Oxford, 1975, pág. 143.
8. Bo Strand (1994). "Suma por partes para aproximaciones de diferencias finitas para d/dx". Journal of Computational Physics . 110 (1): 47–67. Bibcode :1994JCoPh.110...47S. doi :10.1006/jcph.1994.1005.
9. Mark H. Carpenter; David I. Gottlieb; Saul S. Abarbanel (1994). "Condiciones de contorno estables en el tiempo para esquemas de diferencias finitas que resuelven sistemas hiperbólicos: metodología y aplicación a esquemas compactos de alto orden". Journal of Computational Physics . 111 (2): 220–236. Bibcode :1994JCoPh.111..220C. doi :10.1006/jcph.1994.1057. hdl : 2060/19930013937 .

Lectura adicional

• 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.