En 1D tenemos:

u_{xx} = f(x),

con 0 < x < 1 y u(0)= u(1)=0. Particionamos el dominio en n subintervalos introduciendo los puntos x_i = ih con h=\frac{1}{n}. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas (a partir de Taylor alrededor de x_i):

\frac{\partial^2}{\partial x^2} u(x_i) = \frac{u_{i-1} -2 u_i + u_{i+1}}{h^2} - \frac{h^2}{12} \frac{\partial^4}{\partial x^4}u(\xi_i)

con \xi_i \in (x_{i-1}, x_{i+1}) por lo que:

\frac{v_{i-1}-2v_i+v_{i+1}}{h^2} = f_i.

con un error de truncamiento local del orden de O(h^2) y donde utilizamos la v para hablar de la aproximación a la u. Si lo escribimos de manera matricial, tenemos:

\frac{1}{h^2} \begin{bmatrix} -2 & 1 & & & & & \\ 1 & -2 & 1 & & & & \\ & 1 & -2 & 1 & & & \\ & & \cdot & \cdot & \cdot & & \\ & & & \cdot & \cdot & \cdot & \\ & & & & 1 & -2 & 1 \\ & & & & & 1 & -2\end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \\ v_3 \\ \cdot \\ \cdot \\ v_{n-2} \\ v_{n-1} \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \\ f_3 \\ \cdot \\ \cdot \\ f_{n-2} \\ f_{n-1} \end{bmatrix}

Existen otros esquemas de mayor orden para la aproximación, por ejemplo el esquema central de cuarto orden:

\frac{\partial^2}{\partial x^2} u(x_i) = \frac{-u_{i-2}+16u_{i-1}-30u_{i}+16u_{i+1}-u_{i+2}}{12 h^2} + \frac{h^4}{90} \frac{\partial^6}{\partial x^6}u(\xi_i)

con \xi_i \in (x_{i-2},x_{i+2}).

En 2D tenemos:

u_{xx}+u_{yy} = f(x,y),

con 0 < x < 1, 0 < y < 1 y u=0 en el cuadrado unidad. Particionamos el dominio en n_x \times n_y subintervalos introduciendo los puntos (x_i,y_j) = (ih_x, jh_y) con h_x=\frac{1}{n_x} y h_y = \frac{1}{n_y}. De esta manera, para cada uno de los puntos interiores reemplazamos las derivadas por su aproximación de segundo orden en diferencias finitas:

\frac{v_{i-1\,j}-2v_{i\,j}+v_{i+1\,j}}{h_x^2}+\frac{v_{i\,j-1}-2v_{i\,j}+v_{i\,j+1}}{h_y^2} = f_{i\,j}

con v_{0\,j} = v_{n_x\,j} = v_{i\,0} = u_{i\,n_y} = 0. Tenemos un error local de truncamiento del orden de O(h_x^2+h_y^2).

Si hacemos h_x = h_y = h nos queda una plantilla, que representa un operador que interactua localmente en la malla, como esta:

A = \frac{1}{h^2} \begin{bmatrix} & 1 & \\ 1 & -4 & 1 \\ & 1 & \end{bmatrix}

En 3D tenemos:

u_{xx} + u_{yy} + u_{zz} = f(x,y,z),

con 0 < x, y, z < 1 y u=0 en la frontera del cubo unitario. Procediento de la misma manera que antes tenemos:

\frac{v_{i-1\,j\,k}-2v_{i\,j\,k}+v_{i+1\,j\,k}}{h_x^2}+\frac{v_{i\,j-1\,k}-2v_{i\,j\,k}+v_{i\,j+1\,k}}{h_y^2} + \frac{v_{i\,j\,k-1}-2v_{i\,j\,k}+v_{i\,j\,k+1}}{h_z^2}= f_{i\,j\,k}

En este caso, la manera mas sencilla de representar el operador local es mediante el tensor:

A = \frac{1}{h^2} a_{ijk} con a_{0jk}=\begin{bmatrix} & & \\ & 1 & \\ & & \end{bmatrix}, a_{1jk}=\begin{bmatrix} & 1 & \\ 1 & -6 & 1 \\ & 1 & \end{bmatrix} y a_{2jk} = \begin{bmatrix} & & \\ & 1 & \\ & & \end{bmatrix}

si h_x = h_y = h_z = h.

Anuncios