6: Aproximación de diferencia finita

Presentamos aquí la diferenciación numérica, también llamada aproximación de diferencia finita. Esta técnica es comúnmente utilizada para discretizar y resolver ecuaciones diferenciales parciales.

Fórmulas de diferencia finita

Considere la aproximación de la serie Taylor para\(y(x+h)\) y\(y(x-h)\), dada por

\[\begin{aligned} &y(x+h)=y(x)+h y^{\prime}(x)+\frac{1}{2} h^{2} y^{\prime \prime}(x)+\frac{1}{6} h^{3} y^{\prime \prime \prime}(x)+\frac{1}{24} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots, \\ &y(x-h)=y(x)-h y^{\prime}(x)+\frac{1}{2} h^{2} y^{\prime \prime}(x)-\frac{1}{6} h^{3} y^{\prime \prime \prime}(x)+\frac{1}{24} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots . \end{aligned} \nonumber \]

Las definiciones estándar de los derivados dan las aproximaciones de primer orden

\[\begin{aligned} &y^{\prime}(x)=\frac{y(x+h)-y(x)}{h}+\mathrm{O}(h), \\ &y^{\prime}(x)=\frac{y(x)-y(x-h)}{h}+\mathrm{O}(h) . \end{aligned} \nonumber \]

La aproximación de segundo orden más utilizada se llama aproximación de diferencia central y viene dada por

\[y^{\prime}(x)=\frac{y(x+h)-y(x-h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

La aproximación de diferencia finita a la segunda derivada se puede encontrar considerando

\[y(x+h)+y(x-h)=2 y(x)+h^{2} y^{\prime \prime}(x)+\frac{1}{12} h^{4} y^{\prime \prime \prime \prime}(x)+\ldots, \nonumber \]

de donde encontramos

\[y^{\prime \prime}(x)=\frac{y(x+h)-2 y(x)+y(x-h)}{h^{2}}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

A menudo se requiere un método de segundo orden\(x\) en los límites del dominio. Para un punto límite a la izquierda, un método de diferencia hacia adelante de segundo orden requiere la serie Taylor adicional

\[y(x+2 h)=y(x)+2 h y^{\prime}(x)+2 h^{2} y^{\prime \prime}(x)+\frac{4}{3} h^{3} y^{\prime \prime \prime}(x)+\ldots \nonumber \]

Combinamos la serie Taylor para\(y(x+h)\) y\(y(x+2 h)\) para eliminar el término proporcional a\(h^{2}\):

\[y(x+2 h)-4 y(x+h)=-3 y(x)-2 h y^{\prime}(x)+\mathrm{O}\left(\mathrm{h}^{3}\right) \text {. } \nonumber \]

Por lo tanto,

\[y^{\prime}(x)=\frac{-3 y(x)+4 y(x+h)-y(x+2 h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

Para un punto límite a la derecha, enviamos\(h \rightarrow-h\) a encontrar

\[y^{\prime}(x)=\frac{3 y(x)-4 y(x-h)+y(x-2 h)}{2 h}+\mathrm{O}\left(h^{2}\right) . \nonumber \]

Ejemplo: la ecuación de Laplace

Como ejemplo de la técnica de diferencia finita, consideremos cómo discretizar la ecuación bidimensional de Laplace

\[\left(\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}}\right) \Phi=0 \nonumber \]

en el dominio rectangular\([0,1] \times[0,1]\). Asumimos aquí que los valores de\(\Phi\) son conocidos en los límites del dominio. Formamos una cuadrícula bidimensional con\(N\) intervalos en cada dirección junto con puntos finales que se encuentran en los límites escribiendo

\[\begin{array}{ll} x_{i}=i h, & i=0,1, \ldots, N, \\ y_{j}=j h, & j=0,1, \ldots, N, \end{array} \nonumber \]

donde\(h=1 / N\). También denotaremos el valor de\(\Phi\left(x_{i}, y_{j}\right)\) by\(\Phi_{i, j}\). La aproximación de diferencia finita para las segundas derivadas en el punto interior\(\left(x_{i}, y_{j}\right)\) da como resultado una ecuación que escribimos en la forma

\[4 \Phi_{i, j}-\Phi_{i+1, j}-\Phi_{i-1, j}-\Phi_{i, j+1}-\Phi_{i, j-1}=0, \nonumber \]

válido para\(i=1,2, \ldots, N-1\) y\(j=1,2, \ldots, N-1\). Los valores límite\(\Phi_{0, j}, \Phi_{N, j}, \Phi_{i, 0}\), y\(\Phi_{i, N}\) se supone que se dan.

Se puede observar que (6.10) representa un sistema de ecuaciones\((N-1)^{2}\) lineales para\((N-1)^{2}\) incógnitas. Podemos escribir esto como una ecuación matricial si decidimos cómo ordenar las incógnitas\(\Phi_{i, j}\) en un vector. El orden estándar se llama orden natural, y procede desde la parte inferior de la cuadrícula a lo largo de las filas que se mueven hacia la parte superior. Esto ordena las incógnitas como

\[\Phi=\left[\Phi_{1,1}, \Phi_{2,1}, \ldots, \Phi_{(N-1), 1}, \ldots, \Phi_{1,(N-1)}, \Phi_{2(N-1)}, \ldots, \Phi_{(N-1),(N-1)}\right]^{T} . \nonumber \]

Para ilustrar la construcción de la ecuación matricial, consideramos el caso\(N=3\), con dos puntos interiores en cada dirección. Las cuatro ecuaciones lineales resultantes con los términos límite escritos en el lado derecho son

\[\begin{aligned} &4 \Phi_{1,1}-\Phi_{2,1}-\Phi_{1,2}=\Phi_{0,1}+\Phi_{1,0} \\ &4 \Phi_{2,1}-\Phi_{1,1}-\Phi_{2,2}=\Phi_{3,1}+\Phi_{2,0} \\ &4 \Phi_{1,2}-\Phi_{2,2}-\Phi_{1,1}=\Phi_{0,2}+\Phi_{1,3} \\ &4 \Phi_{2,2}-\Phi_{1,2}-\Phi_{2,1}=\Phi_{3,2}+\Phi_{2,3} \end{aligned} \nonumber \]

y la ecuación matricial correspondiente es

\[\left(\begin{array}{rrrr} 4 & -1 & -1 & 0 \\ -1 & 4 & 0 & -1 \\ -1 & 0 & 4 & -1 \\ 0 & -1 & -1 & 4 \end{array}\right)\left(\begin{array}{l} \Phi_{1,1} \\ \Phi_{2,1} \\ \Phi_{1,2} \\ \Phi_{2,2} \end{array}\right)=\left(\begin{array}{l} \Phi_{0,1}+\Phi_{1,0} \\ \Phi_{3,1}+\Phi_{2,0} \\ \Phi_{0,2}+\Phi_{1,3} \\ \Phi_{3,2}+\Phi_{2,3} \end{array}\right) . \nonumber \]

El patrón aquí puede no ser obvio, pero la matriz laplaciana se descompone en matrices de 2 -por- 2 bloques.

Con un poco de pensamiento, se puede obtener un resultado general para una\(N_{y}\) cuadrícula\(N_{x}\) -by- de puntos internos (con\(N_{x}+1\) e\(N_{y}+1\) intervalos en las direcciones\(x\) - e y, respectivamente). La matriz Laplaciana luego se descompone en matrices\(N_{x}\) -por-\(N_{x}\) bloque. La diagonal contiene\(N_{y}\) de estas matrices\(N_{x}\) -por-\(N_{x}\) bloque, cada una de las cuales son tridiagonales con un 4 en la diagonal y a\(-1\) en las diagonales fuera. Inmediatamente por encima y por debajo de estas matrices de bloques en la diagonal se encuentran matrices de\(N_{y}-1\) bloques también de tamaño\(N_{x}\) -by-\(N_{x}\), cada una de las cuales son diagonales\(-1\) a lo largo de la diagonal.

El código MATLAB para la matriz Laplaciana se puede encontrar en la web en la función sp_laplace.m Este código fue escrito por Bordner y Saied en 1995, y he escrito una versión más moderna y rápida de este código en sp_laplace_new.m.

Un método alternativo de solución, que posteriormente utilizaremos in\(\$ 17.2 .4\), incluye los valores límite en el vector de solución. Si una de las condiciones de contorno está en posición\(n\) en este vector, entonces la fila\(n\) de la matriz del lado izquierdo tendrá solo una en la diagonal con todos los demás elementos iguales a cero (es decir, una fila de la matriz de identidad), y el elemento correspondiente en la derecha- lado vector tendrá el valor en el límite. Aquí, con el número total de puntos de cuadrícula (incluidos los puntos de límite) en las\(x\)\(y\) direcciones y dadas por\(N_{x}\) y\(N_{y}\), la matriz del lado izquierdo se genera usando\(A=\) sp_laplace_new (N_X, N_Y ), y todas las filas correspondientes a los valores de límite se sustituyen por las filas correspondientes de la matriz de identidad. Esta formulación puede ser un poco más fácil de codificar, y también más fácil de incorporar otras condiciones de contorno más generales.