You are currently browsing the tag archive for the ‘diferencias finitas’ tag.

El caso mas sencillo es cuando tenemos cinco puntos:

(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3,y_3), (x_4,y_4),

de manera que:

x_1-x_0 = 2(x_2-x_1) = 2(x_3-x_2) = x_4-x_3.

Tal y como escribimos en el anterior post, el polinomio general para cinco puntos quedaria:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x) + y_3 l_3(x) + y_4 l_4(x),

donde:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4},

l_1(x) = \frac{x-x_0}{x_1-x_0} \frac{x-x_2}{x_1-x_2} \frac{x-x_3}{x_1-x_3} \frac{x-x_4}{x_1-x_4},

l_2(x) = \frac{x-x_0}{x_2-x_0} \frac{x-x_1}{x_2-x_1} \frac{x-x_3}{x_2-x_3} \frac{x-x_4}{x_2-x_4},

l_3(x) = \frac{x-x_0}{x_3-x_0} \frac{x-x_1}{x_3-x_1} \frac{x-x_2}{x_3-x_2} \frac{x-x_4}{x_3-x_4},

l_4(x) = \frac{x-x_0}{x_4-x_0} \frac{x-x_1}{x_4-x_1} \frac{x-x_2}{x_4-x_2} \frac{x-x_3}{x_4-x_3}.

Dado que tratamos de aproximar la segunda derivada, tenemos, en el caso de equidistancia quedaría:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12 h^2}.

Vamos a ver que queda ahora en nuestro caso. Reescribimos los l_i(x) en función de x_0 de manera que, por ejemplo, tenemos:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4} =

=\frac{x-x_0-2h}{x_0-x_0-2h} \frac{x-x_0-3h}{x_0-x_0-3h} \frac{x-x_0-4h}{x_0-x_0-4h} \frac{x-x_0-6h}{x_0-x_0-6h} =

=\frac{(x-x_0-2h)(x-x_0-3h)(x-x_0-4h)(x-x_0-6h)}{144h^4},

que para x=x_2=x_0+3h (2h+h), la segunda derivada queda:

\frac{d^2}{dx^2} l_0(x) |_{x=x_2} = -\frac{1}{72h^2}.

Haciendo lo mismo para todas las derivadas segundas de todos los l_i(x), obtenemos la siguiente formula en diferencias finitas:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-3} + 81 u_{u-1} - 160 u_i + 81 u_{i+1} - u_{i+3}}{72h^2},

donde los índices hacen referencia a una malla inicial equiespaciada.

Para el mismo denominador, antes teniamos:

\frac{d^2}{dx^2}u_i = \frac{-6u_{i-2}+96u_{i-1}-180u_i+96u_{i+1}-6u_{i+2}}{72 h^2}.

Podemos hacer lo mismo para cada uno de los puntos x=x_0, x=x_1=x_0+2h, x=x_3=x_0+4h y x=x_4=x_0+6h:

\frac{d^2}{dx^2}u_{i-3} = \frac{40 u_{i-3} -243 u_{i-1} + 352 u_i -162 u_{i+1} + 13 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i-2} = \frac{7 u_{i-3} - 32 u_i + 27 u_{i+1} - 2 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+1} =\frac{-2 u_{i-3} + 27 u_{i-1} - 32 u_i -+7 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+3} =\frac{13 u_{i-3} -162 u_{i-1} + 352 u_i -243 u_{i+1} + 40 u_{i+3}}{36 h^2}

 En el caso de x=x_0+4h (3h+h), tenemos:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-4} + 256_{u-1} - 510 u_i + 256 u_{i+1} - u_{i+4}}{140h^2}.

Anuncios

Dados n+1 puntos:

(x_0,y_0), \ldots, (x_n,y_n),

definimos el polinomio interpolador de Lagrange:

L(x) := \Sigma_{i=0}^n y_i l_i(x)

donde:

l_i(x) := \Pi_{0 \leq j \leq n,j \neq i} \frac{x-x_j}{x_i-x_j}.

Con esto, tenemos:

\frac{d^n}{dx^n}L(x) = \Sigma_{i=0}^n y_i \frac{d^n}{dx^n}l_i(x)

Supongamos que tenemos n=2:

(x_0,y_0), (x_1,y_1).

En este caso, tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x)

con:

l_0(x) = \frac{x-x_1}{x_0-x_1}

l_1(x) = \frac{x-x_0}{x_1-x_0}

Como tenemos dos puntos, únicamente podemos aproximar la primera derivada:

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x)

con:

\frac{d}{dx} l_0(x) = \frac{1}{x_0-x_1}

\frac{d}{dx} l_1(x) = \frac{1}{x_1-x_0},

de manera que:

\frac{d}{dx}L(x) = \frac{y_1 - y_0}{x_1 - x_0}.

Con esto, tenemos las aproximaciones de primer orden:

u_x(x_0) = u_x(x_1) = \frac{y_1 - y_0}{x_1-x_0},

que, en índices, queda:

\frac{d}{dx}u_i \approx \frac{u_{i+1}-u_i}{h_x},

\frac{d}{dx}u_{i+1} \approx \frac{u_{i+1}-u_i}{h_x},

donde h_x = x_1 - x_0 (y que es lógico, ya que la derivada será una constante).

Supongamos ahora que tenemos tres puntos:

(x_0,y_0),(x_1,y_1),(x_2,y_2).

En este caso tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x),

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x) + y_2 \frac{d}{dx} l_2(x),

\frac{d^2}{dx^2} L(x) = y_0 \frac{d^2}{dx^2} l_0(x) + y_1 \frac{d^2}{dx^2} l_1(x) + y_2 \frac{d^2}{dx^2} l_2(x).

Ahora tenemos:

l_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{x^2 + (-x_1-x_2)x + x_1x_2}{(x_0-x_1)(x_0-x_2)},

l_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{x^2 + (-x_0-x_2)x + x_0x_2}{(x_1-x_0)(x_1-x_2)},

l_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{x^2 + (-x_0-x_1)x + x_0x_1}{(x_2-x_0)(x_2-x_1)}.

\frac{d}{dx} l_0(x) = \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)},

\frac{d}{dx} l_1(x) = \frac{2x+(-x_0-x_2)}{(x_1-x_0)(x_1-x_2)},

\frac{d}{dx} l_2(x) = \frac{2x+(-x_0-x_1)}{(x_2-x_0)(x_2-x_1)}.

\frac{d^2}{dx^2} l_0(x) = \frac{2}{(x_0-x_1)(x_0-x_2)},

\frac{d^2}{dx^2} l_1(x) = \frac{2}{(x_1-x_0)(x_1-x_2)},

\frac{d^2}{dx^2} l_2(x) = \frac{2}{(x_2-x_0)(x_2-x_1)}.

De esta manera, por ejemplo, tenemos que:

\frac{d^2}{dx^2}L(x) = y_0 \frac{2}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2}{(x_2-x_0)(x_2-x_1)},

que, en el caso de tener los puntos equiespaciados (h_x:=x_{i+1}-x_{i} con 0 \leq i \leq 1), queda la aproximación de segundo orden:

\frac{d^2}{dx^2}u(x) = \frac{y_0 - 2 y_1 + y_2}{h_x^2},

que además, como era de esperar, es independiente de x:

\frac{d^2}{dx^2}u_{i-1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_{i+1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}.

¿Qué pasa con la primera derivada? En este caso, el resultado si que depende de x (por lo que tendremos un resultado diferente en función de si la x vale x_0, x_1 o x_2) y lo que obtenemos es:

\frac{d}{dx} L(x) = y_0 \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

por lo que:

\frac{d}{dx} L(x_0) = y_0 \frac{2x_0 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_0 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_0 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_1) = y_0 \frac{2x_1 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_1 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_1 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_2) = y_0 \frac{2x_2 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_2 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_2 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

que, con equiespaciado, quedan:

L'(x_0) = \frac{3 y_0 -4 y_1 + y_2}{2 h_x},

L'(x_1) = \frac{y_0 -2 y_1 + y_2}{2 h_x} y

L'(x_2) = \frac{y_0 -4 y_1 + 3 y_2}{2 h_x}.

Laplaciano en cartesianas:

\Delta u = \Sigma_i \frac{\partial^2}{\partial x_i^2}u

1d

\frac{u_{i-1}-2u_i+u_{i+1}}{h^2} = f_i

\frac{1}{h^2}u_{i-1} + \frac{1}{h^2}u_{i+1} +\frac{-2}{h^2}u_i= f_i

2d

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

i fijo:

\frac{1}{h_y^2}u_{i,j-1} + \frac{1}{h_y^2}u_{i,j+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_x^2}u_{i-1,j} + \frac{-1}{h_x^2}u_{i+1,j})

j fijo:

\frac{1}{h_x^2}u_{i-1,j} + \frac{1}{h_x^2}u_{i+1,j} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_y^2}u_{i,j-1} + \frac{-1}{h_y^2}u_{i,j+1})

3d

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

i,j fijos:

\frac{1}{h_z^2}u_{i,j,k-1} + \frac{1}{h_z^2}u_{i,j,k+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k})

i,k fijos:

\frac{1}{h_y^2}u_{i,j-1,k} + \frac{1}{h_y^2}u_{i,j+1,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

j,k fijos:

\frac{1}{h_x^2}u_{i-1,j,k} + \frac{1}{h_x^2}u_{i+1,j,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

Aunque siempre podemos hacer cambios de coordenadas, vamos a ver como quedan los esquemas de diferencias finitas en sistemas no rectangulares: coordenadas cilíndricas, (\rho,\phi, z), y coordenadas esféricas, (r,\theta,\phi). Nos centraremos en la ecuación de Poisson aunque la técnica se puede extender de manera inmediata a cualquier tipo de PDE.

En coordenadas cilíndricas podemos escribir:

\nabla \cdot \nabla u = \frac{\partial^2}{\partial \rho^2}u + \frac{1}{\rho}\frac{\partial}{\partial \rho}u + \frac{1}{\rho^2}\frac{\partial^2}{\partial \phi^2}u + \frac{\partial^2}{\partial z^2} u = f,

que podemos discretizar como:

\frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{(\Delta \rho)^2} +

+ \frac{1}{\rho_{i,j,k}}\frac{u_{i+1,j,k}-u_{i-1,j,k}}{2\Delta \rho} +

+ \frac{1}{\rho_{i,j,k}^2} \frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{(\Delta \phi)^2} +

+ \frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{(\Delta z)^2} = f_{i,j,k}

En coordenadas esféricas tenemos:

\nabla \cdot \nabla u = \frac{\partial^2}{\partial r^2}u + \frac{2}{r} \frac{\partial}{\partial r}u + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}u + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial \theta}u + \frac{1}{r^2 \sin^2\theta} \frac{\partial^2}{\partial \phi^2}u = f

que podemos discretizar como:

\frac{u_{i-1,j,k-2u_{i,j,k}+u_{i+1,j,k}}}{(\Delta r)^2} +

+ \frac{2}{r_{i,j,k}} \frac{u_{i+1,j,k}+u_{i-1,j,k}}{2\Delta r} +

+ \frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{(r_{i,j,k} \Delta \theta)^2} +

+ \frac{1}{r_{i,j,k}^2 \sin \phi_{i,j,k}} \frac{u_{i,j+1,k}-u_{i,j-1,k}}{2 \Delta \phi} +

+ \frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{(r_{i,j,k} \sin \phi_{i,j,k} \Delta \phi)^2} = f_{i,j,k}

Vamos a suponer n=3 para reducir el tamaño de las matrices.

Empezamos suponiendo que conocemos:

\frac{\partial}{\partial x}|_{0,0,}u, \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}u

\frac{\partial}{\partial y}|_{0,0}u, \frac{\partial}{\partial y}|_{1,0}u

\frac{\partial}{\partial y}|_{0,2}u, \frac{\partial}{\partial y}|_{1,2}u

u|_{2,0}, u|_{2,1}, u|_{2,2}

Discretizamos:

\frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{h^2} + \frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{h^2} = f_{0,0}

\frac{u_{-1,1}-2u_{0,1}+u_{1,1}}{h^2} + \frac{u_{0,0}-2u_{0,1}+u_{0,2}}{h^2} = f_{0,1}

\frac{u_{-1,2}-2u_{0,2}+u_{1,2}}{h^2} + \frac{u_{0,1}-2u_{0,2}+u_{0,3}}{h^2} = f_{0,2}

\frac{u_{0,0}-2u_{1,0}+u_{2,0}}{h^2} + \frac{u_{1,-1}-2u_{1,0}+u_{1,1}}{h^2} = f_{1,0}

\frac{u_{0,1}-2u_{1,1}+u_{2,1}}{h^2} + \frac{u_{1,0}-2u_{1,1}+u_{1,2}}{h^2} = f_{1,1}

\frac{u_{0,2}-2u_{1,2}+u_{2,2}}{h^2} + \frac{u_{1,1}-2u_{1,2}+u_{1,3}}{h^2} = f_{1,2}

En las fronteras, sabemos que:

\frac{u_{1,0}-u_{-1,0}}{2h} = \frac{\partial}{\partial x}|_{0,0}u \Leftrightarrow u_{-1,0}=u_{1,0}-2h \frac{\partial}{\partial x}|_{0,0}u

\frac{u_{1,1}-u_{-1,1}}{2h} = \frac{\partial}{\partial x}|_{0,1}u \Leftrightarrow u_{-1,1}=u_{1,1}-2h \frac{\partial}{\partial x}|_{0,1}u

\frac{u_{1,2}-u_{-1,2}}{2h} = \frac{\partial}{\partial x}|_{0,2}u \Leftrightarrow u_{-1,2}=u_{1,2}-2h \frac{\partial}{\partial x}|_{0,2}u

\frac{u_{0,1}-u_{0,-1}}{2h} = \frac{\partial}{\partial y}|_{0,0}u \Leftrightarrow u_{0,-1}=u_{0,1}-2h \frac{\partial}{\partial y}|_{0,0}u

\frac{u_{1,1}-u_{1,-1}}{2h} = \frac{\partial}{\partial y}|_{1,0}u \Leftrightarrow u_{1,-1}=u_{1,1}-2h \frac{\partial}{\partial y}|_{1,0}u

\frac{u_{0,3}-u_{0,1}}{2h} = \frac{\partial}{\partial y}|_{0,2}u \Leftrightarrow u_{0,3}=u_{0,1}+2h \frac{\partial}{\partial y}|_{0,2}u

\frac{u_{1,3}-u_{1,1}}{2h} = \frac{\partial}{\partial y}|_{1,2}u \Leftrightarrow u_{1,3}=u_{1,1}+2h \frac{\partial}{\partial y}|_{1,2}u

La matriz queda:

\left(  \begin{array}{ccc|ccc}  -4 & 2 & 0 & 2 & 0 & 0 \\  1 & -4 & 1 & 0 & 2 & 0 \\  0 & 2 & -4 & 0 & 0 & 2 \\ \hline  1 & 0 & 0 & -4 & 2 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 \\  0 & 0 & 1 & 0 & 2 & -4  \end{array}  \right)

Simetrizable como:

\left(  \begin{array}{ccc|ccc}  -1 & \frac{1}{2} & 0 & \frac{1}{2} & 0 & 0 \\  \frac{1}{2} & -2 & \frac{1}{2} & 0 & 1 & 0 \\  0 & \frac{1}{2} & -1 & 0 & 0 & \frac{1}{2} \\ \hline  \frac{1}{2} & 0 & 0 & -2 & 1 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 \\  0 & 0 & \frac{1}{2} & 0 & 1 & -2  \end{array}  \right)

Tenemos 6 ecuaciones con 6 incognitas y la matriz tiene rango 6, por lo que la solución es única.

En el segundo caso, suponemos que todas las fronteras son Neumann:

\frac{\partial}{\partial x}|_{0,0}u, \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}u

\frac{\partial}{\partial y}|_{0,0}u, \frac{\partial}{\partial y}|_{1,0}u, \frac{\partial}{\partial y}|_{2,0}u

\frac{\partial}{\partial y}|_{0,2}u, \frac{\partial}{\partial y}|_{1,2}u, \frac{\partial}{\partial y}|_{2,2}u

\frac{\partial}{\partial x}|_{2,0}u, \frac{\partial}{\partial x}|_{2,1}u, \frac{\partial}{\partial x}|_{2,2}u

Si discretizamos:

\frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{h^2} + \frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{h^2} = f_{0,0}

\frac{u_{-1,1}-2u_{0,1}+u_{1,1}}{h^2} + \frac{u_{0,0}-2u_{0,1}+u_{0,2}}{h^2} = f_{0,1}

\frac{u_{-1,2}-2u_{0,2}+u_{1,2}}{h^2} + \frac{u_{0,1}-2u_{0,2}+u_{0,3}}{h^2} = f_{0,2}

\frac{u_{0,0}-2u_{1,0}+u_{2,0}}{h^2} + \frac{u_{1,-1}-2u_{1,0}+u_{1,1}}{h^2} = f_{1,0}

\frac{u_{0,1}-2u_{1,1}+u_{2,1}}{h^2} + \frac{u_{1,0}-2u_{1,1}+u_{1,2}}{h^2} = f_{1,1}

\frac{u_{0,2}-2u_{1,2}+u_{2,2}}{h^2} + \frac{u_{1,1}-2u_{1,2}+u_{1,3}}{h^2} = f_{1,2}

\frac{u_{1,0}-2u_{2,0}+u_{3,0}}{h^2} + \frac{u_{2,-1}-2u_{2,0}+u_{2,1}}{h^2} = f_{2,0}

\frac{u_{1,1}-2u_{2,1}+u_{3,1}}{h^2} + \frac{u_{2,0}-2u_{2,1}+u_{2,2}}{h^2} = f_{2,1}

\frac{u_{1,2}-2u_{2,2}+u_{3,2}}{h^2} + \frac{u_{2,1}-2u_{2,2}+u_{2,3}}{h^2} = f_{2,2}

En las fronteras, sabemos que:

\frac{u_{1,0}-u_{-1,0}}{2h} = \frac{\partial}{\partial x}|_{0,0}u \Leftrightarrow u_{-1,0}=u_{1,0}-2h \frac{\partial}{\partial x}|_{0,0}u

\frac{u_{1,1}-u_{-1,1}}{2h} = \frac{\partial}{\partial x}|_{0,1}u \Leftrightarrow u_{-1,1}=u_{1,1}-2h \frac{\partial}{\partial x}|_{0,1}u

\frac{u_{1,2}-u_{-1,2}}{2h} = \frac{\partial}{\partial x}|_{0,2}u \Leftrightarrow u_{-1,2}=u_{1,2}-2h \frac{\partial}{\partial x}|_{0,2}u

\frac{u_{0,1}-u_{0,-1}}{2h} = \frac{\partial}{\partial y}|_{0,0}u \Leftrightarrow u_{0,-1}=u_{0,1}-2h \frac{\partial}{\partial y}|_{0,0}u

\frac{u_{1,1}-u_{1,-1}}{2h} = \frac{\partial}{\partial y}|_{1,0}u \Leftrightarrow u_{1,-1}=u_{1,1}-2h \frac{\partial}{\partial y}|_{1,0}u

\frac{u_{2,1}-u_{2,-1}}{2h} = \frac{\partial}{\partial y}|_{2,0}u \Leftrightarrow u_{2,-1}=u_{2,1}-2h \frac{\partial}{\partial y}|_{2,0}u

\frac{u_{0,3}-u_{0,1}}{2h} = \frac{\partial}{\partial y}|_{0,2}u \Leftrightarrow u_{0,3}=u_{0,1}+2h \frac{\partial}{\partial y}|_{0,2}u

\frac{u_{1,3}-u_{1,1}}{2h} = \frac{\partial}{\partial y}|_{1,2}u \Leftrightarrow u_{1,3}=u_{1,1}+2h \frac{\partial}{\partial y}|_{1,2}u

\frac{u_{2,3}-u_{2,1}}{2h} = \frac{\partial}{\partial y}|_{2,2}u \Leftrightarrow u_{2,3}=u_{2,1}+2h \frac{\partial}{\partial y}|_{2,2}u

\frac{u_{3,0}-u_{1,0}}{2h} = \frac{\partial}{\partial x}|_{2,0}u \Leftrightarrow u_{3,0}=u_{1,0}+2h \frac{\partial}{\partial x}|_{2,0}u

\frac{u_{3,1}-u_{1,1}}{2h} = \frac{\partial}{\partial x}|_{2,1}u \Leftrightarrow u_{3,1}=u_{1,1}+2h \frac{\partial}{\partial x}|_{2,1}u

\frac{u_{3,2}-u_{1,2}}{2h} = \frac{\partial}{\partial x}|_{2,2}u \Leftrightarrow u_{3,2}=u_{1,2}+2h \frac{\partial}{\partial x}|_{2,2}u

La matriz, por tanto, queda:

\text{A6}=\left(  \begin{array}{ccc|ccc|ccc}  -4 & 2 & 0 & 2 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 \\  0 & 2 & -4 & 0 & 0 & 2 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & -4 & 2 & 0 & 1 & 0 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\  0 & 0 & 1 & 0 & 2 & -4 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 2 & 0 & 0 & -4 & 2 & 0 \\  0 & 0 & 0 & 0 & 2 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 2 & 0 & 2 & -4  \end{array}  \right)

Simetrizable como:

\text{A6s}=\left(  \begin{array}{ccc|ccc|ccc}  -1 & 1/2 & 0 & 1/2 & 0 & 0 & 0 & 0 & 0 \\  1/2 & -2 & 1/2 & 0 & 1 & 0 & 0 & 0 & 0 \\  0 & 1/2 & -1 & 0 & 0 & 1/2 & 0 & 0 & 0 \\ \hline  1/2 & 0 & 0 & -2 & 1 & 0 & 1/2 & 0 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\  0 & 0 & 1/2 & 0 & 1 & -2 & 0 & 0 & 1/2 \\ \hline  0 & 0 & 0 & 1/2 & 0 & 0 & -1 & 1/2 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 1/2 & -2 & 1/2 \\  0 & 0 & 0 & 0 & 0 & 1/2 & 0 & 1/2 & -1  \end{array}  \right)

En este caso, tenemos 9 ecuaciones con 9 incognitas pero la matriz tiene rango 8, por lo que tenemos infinitas soluciones. Hay que conservar.

Suponemos \Delta u = f en 2D, es decir,

\frac{\partial^2}{\partial x^2}u(x,y) + \frac{\partial^2}{\partial y^2}u(x,y) = f(x,y).

Miraremos como queda la matriz del sistema al discretizar, como simetrizarla y su rango en tres casos: condición Neuman respecto x en una frontera, condición Neumann respecto y en una frontera y condición Neumann respecto x e y en dos fronteras.

Discretizamos con n=5. Si todas las condiciones fueran Dirichlet, la matriz quedaría:

A_1 = \left(  \begin{array}{ccc|ccc|ccc}  -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\  0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\  0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4  \end{array}  \right) .

En este caso, A_1 \in \mathcal{M}(9 \times 9) y simétrica, lo que permite tratar de manera conjunta los problemas de existencia y unicidad de solución. Si calculamos su rango obtenemos 9 por lo que existe solución y es única. Desde el punto de vista algebraico, es el punto (u_{1,1},u_{1,2},u_{1,3},u_{2,1},u_{2,2},u_{2,3},u_{3,1},u_{3,2},u_{3,3}) intersección de 9 hiperplanos

-4x_{1,1} + x_{1,2} + x_{2,1} = f_{1,1},

x_{1,1}-4x_{1,2}+x_{1,3} + x_{2,2} = f_{1,2},

\ldots

en el espacio \mathbb{R}^9.

Si condiremos conocidos \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}u, \frac{\partial}{\partial x}|_{0,3}u en lugar de u_{0,1}, u_{0,2}, u_{0,3} (u_{0,0} y u_{0,4} son conocidos por las otras fronteras que son Dirichelt), tenemos:

A_2 = \left(  \begin{array}{ccc|ccc|ccc|ccc}  -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & -4 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \\ \hline  0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4  \end{array}  \right)

de manera que A_2 \in \mathcal{M}(12 \times 12) y no es simétrica. Sin embargo es facilmente simetrizable dividiendo las tres primera filas (hacemos lo mismo en el termino independiente) por 2:

A_2 = \left(  \begin{array}{ccc|ccc|ccc|ccc}  -2 & \frac{1}{2} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  \frac{1}{2} & -2 & \frac{1}{2} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & \frac{1}{2} & -2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \\ \hline  0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4  \end{array}  \right)

Tenemos 12 incognitas (u_{i,j} con i=0..3 y j=1..3) y el rango de A_2 es 12, por lo que la solución, nuevamente, es única.

Para el caso en el que conocemos \frac{\partial}{\partial y}|_{1,0}u, \frac{\partial}{\partial y}|_{2,0}u, \frac{\partial}{\partial y}|_{3,0}u en lugar de u_{1,0}, u_{2,0}, u_{3,0}, si el orden que tomamos es el contrario al tomado anteriormente llegaremos a la misma estructura de antes. Sin embargo, como en el siguiente caso nos veremos obligados a seleccionar uno de los dos, vamos a ver como queda este caso utilizando el mismo orden que antes:

A_3 = \left(  \begin{array}{cccc|cccc|cccc}  -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 \\  0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\  0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\  0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4  \end{array}  \right)

que podemos simetrizar facilmente y queda:

A_3 = \left(  \begin{array}{cccc|cccc|cccc}  -2 & 1 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \hline  \frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 \\  0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\  0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\  0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4  \end{array}  \right)

Tenemos 12 ecuaciones con 12 incognitas (u_{i,j} con i=1..3 y j=0..3) y el rango de A_3 es 12, por lo que la solución es única.

Finalmente, suponemos conocidos \frac{\partial}{\partial x}|_{0,0}u, \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}u, \frac{\partial}{\partial x}|_{0,3}u, \frac{\partial}{\partial y}|_{0,0}u, \frac{\partial}{\partial y}|_{1,0}u, \frac{\partial}{\partial y}|_{2,0}u, \frac{\partial}{\partial y}|_{3,0}u que incorpora 7 ecuaciones mas a las 9 que ya teniamos por lo que nos queda una matrix A_4 \in \mathcal{M}(16 \times 16):

\left(  \begin{array}{cccc|cccc|cccc|cccc}  -4 & 2 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  1 & -4 & 1 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & -4 & 1 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & -4 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline  1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \hline  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 & 1 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & -4 & 2 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4  \end{array}  \right),

simetrizable dividiendo la fila correspondiente a u_{0,0} por 4, y las correspondientes a u_{0,1}, u_{0,2}, u_{0,3}, u_{1,0},u_{2,0}, u_{3,0}  por 2, quedando:

\left(  \begin{array}{cccc|cccc|cccc|cccc}  -1 & \frac{1}{2} & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  \frac{1}{2} & -2 & \frac{1}{2} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & \frac{1}{2} & -2 & \frac{1}{2} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & \frac{1}{2} & -2 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline  \frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \hline  0 & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\ \hline  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & 0 & 0 & -2 & 1 & 0 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4  \end{array}  \right),

con lo que el sistema vuelve a ser compatible y determinado.

Suponemos n=5. En el caso de tener todas las fronteras con condiciones Dirichlet:

\frac{u_{0,1} -2u_{1,1} + u_{2,1}}{h^2} + \frac{u_{1,0} -2u_{1,1} + u_{1,2}}{h^2} = f_{1,1} para i,j=1,1,

\frac{u_{0,2} -2u_{1,2} + u_{2,2}}{h^2} + \frac{u_{1,1} -2u_{1,2} + u_{1,3}}{h^2} = f_{1,2} para i,j=1,2,

\frac{u_{0,3} -2u_{1,3} + u_{2,3}}{h^2} + \frac{u_{1,2} -2u_{1,3} + u_{1,4}}{h^2} = f_{1,3} para i,j=1,3,

\frac{u_{1,1} -2u_{2,1} + u_{3,1}}{h^2} + \frac{u_{2,0} -2u_{2,1} + u_{2,2}}{h^2} = f_{2,1} para i,j=2,1,

\frac{u_{1,2} -2u_{2,2} + u_{3,2}}{h^2} + \frac{u_{2,1} -2u_{2,2} + u_{2,3}}{h^2} = f_{2,2} para i,j=2,2,

\frac{u_{1,3} -2u_{2,3} + u_{3,3}}{h^2} + \frac{u_{2,2} -2u_{2,3} + u_{2,4}}{h^2} = f_{2,3} para i,j=2,3,

\frac{u_{2,1} -2u_{3,1} + u_{4,1}}{h^2} + \frac{u_{3,0} -2u_{3,1} + u_{3,2}}{h^2} = f_{3,1} para i,j=3,1,

\frac{u_{2,2} -2u_{3,2} + u_{4,2}}{h^2} + \frac{u_{3,1} -2u_{3,2} + u_{3,3}}{h^2} = f_{3,2} para i,j=3,2,

\frac{u_{2,3} -2u_{3,3} + u_{4,3}}{h^2} + \frac{u_{3,2} -2u_{3,3} + u_{3,4}}{h^2} = f_{3,3} para i,j=3,3,

de donde:

\begin{bmatrix} f_{1,1} -\frac{u_{1,0} + u_{0,1}}{h^2} & f_{1,2} - \frac{u_{0,2}}{h^2} & f_{1,3} - \frac{u_{0,3}+u_{1,4}}{h^2} \\ f_{2,1} -\frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} - \frac{u_{2,4}}{h^2} \\ f_{3,1} - \frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} - \frac{u_{4,2}}{h^2} & f_{3,3} - \frac{u_{4,3}+u_{3,4}}{h^2} \end{bmatrix}

En forma de matriz por bloques (para pensar en la simetrización):

\frac{1}{h^2} \begin{bmatrix} -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & -4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & -4 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & -4 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & -4 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 1 & -4 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & -4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & -4 \end{bmatrix} u_{i,j} = \begin{bmatrix} f_{1,1} -\frac{u_{1,0} + u_{0,1}}{h^2} \\ f_{1,2} - \frac{u_{0,2}}{h^2} \\ f_{1,3} - \frac{u_{0,3}+u_{1,4}}{h^2} \\ f_{2,1} -\frac{u_{2,0}}{h^2} \\ f_{2,2} \\ f_{2,3} - \frac{u_{2,4}}{h^2} \\ f_{3,1} - \frac{u_{3,0}+u_{4,1}}{h^2} \\ f_{3,2} - \frac{u_{4,2}}{h^2} \\ f_{3,3} - \frac{u_{4,3}+u_{3,4}}{h^2} \end{bmatrix}

¿Qué pasa ahora si en lugar de conocer u_{0,1}, u_{0,2}, u_{0,3} conocemos \frac{\partial}{\partial x}|_{0,1}u, \frac{\partial}{\partial x}|_{0,2}u, \frac{\partial}{\partial x}u|_{0,3}? Necesitamos tres ecuaciones mas:

\frac{u_{-1,1} -2u_{0,1} + u_{1,1}}{h^2} + \frac{u_{0,0} -2u_{0,1} + u_{0,2}}{h^2} = f_{0,1} para i,j=0,1

\frac{u_{-1,2} -2u_{0,2} + u_{1,2}}{h^2} + \frac{u_{0,1} -2u_{0,2} + u_{0,3}}{h^2} = f_{0,2} para i,j=0,2

\frac{u_{-1,3} -2u_{0,3} + u_{1,3}}{h^2} + \frac{u_{0,2} -2u_{0,3} + u_{0,4}}{h^2} = f_{0,3} para i,j=0,3

y

\frac{u_{1,1}-u_{-1,1}}{2h} = \frac{\partial}{\partial x}|_{0,1}u \Leftrightarrow u_{-1,1} = u_{1,1} - 2h \, \frac{\partial}{\partial x}|_{0,1}u

\frac{u_{1,2}-u_{-1,2}}{2h} = \frac{\partial}{\partial x}|_{0,2}u \Leftrightarrow u_{-1,2} = u_{1,2} - 2h \, \frac{\partial}{\partial x}|_{0,2}u

\frac{u_{1,3}-u_{-1,3}}{2h} = \frac{\partial}{\partial x}|_{0,3}u \Leftrightarrow u_{-1,3} = u_{1,3} - 2h \, \frac{\partial}{\partial x}|_{0,3}u

por lo que:

\begin{bmatrix} f_{0,1} +\frac{2h \, \frac{\partial}{\partial x}|_{0,1}u - u_{0,0}}{h^2} & f_{0,2} + \frac{2h \, \frac{\partial}{\partial x}|_{0,2}u}{h^2} & f_{0,3} + \frac{ 2h \, \frac{\partial}{\partial x}|_{0,3}u - u_{0,4}}{h^2} \\ f_{1,1} -\frac{u_{1,0}}{h^2} & f_{1,2} & f_{1,3} - \frac{u_{1,4}}{h^2} \\ f_{2,1} -\frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} - \frac{u_{2,4}}{h^2} \\ f_{3,1} - \frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} - \frac{u_{4,2}}{h^2} & f_{3,3} - \frac{u_{4,3}+u_{3,4}}{h^2} \end{bmatrix}

La matriz queda:

\frac{1}{h^2} \begin{bmatrix} -4 & 1 & 0 & 2 & 0 & 0 & \ldots \\ 1 & -4 & 1 & 0 & 2 & 0 & \ldots \\ 0 & 1 & -4 & 0 & 0 & 2 & \ldots \\ 1 & 0 & 0 & -4 & 1 & 0 & \ldots \\ 0 & 1 & 0 & 1 & -4 & 1 & \ldots \\ 0 & 0 & 1 & 0 & 1 & -4 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}

Que podemos simetrizar:

\frac{1}{h^2} \begin{bmatrix} -2 & \frac{1}{2} & 0 & 1 & 0 & 0 & \ldots \\ \frac{1}{2} & -2 & \frac{1}{2} & 0 & 1 & 0 & \ldots \\ 0 & \frac{1}{2} & -2 & 0 & 0 & 1 & \ldots \\ 1 & 0 & 0 & -4 & 1 & 0 & \ldots \\ 0 & 1 & 0 & 1 & -4 & 1 & \ldots \\ 0 & 0 & 1 & 0 & 1 & -4 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}

con:

\begin{bmatrix} \frac{1}{2}(f_{0,1} +\frac{2h \, \frac{\partial}{\partial x}|_{0,1}u - u_{0,0}}{h^2}) & \frac{1}{2}(f_{0,2} + \frac{2h \, \frac{\partial}{\partial x}|_{0,2}u}{h^2}) & \frac{1}{2}(f_{0,3} + \frac{ 2h \, \frac{\partial}{\partial x}|_{0,3}u - u_{0,4}}{h^2}) \\ f_{1,1} -\frac{u_{1,0}}{h^2} & f_{1,2} & f_{1,3} - \frac{u_{1,4}}{h^2} \\ f_{2,1} -\frac{u_{2,0}}{h^2} & f_{2,2} & f_{2,3} - \frac{u_{2,4}}{h^2} \\ f_{3,1} - \frac{u_{3,0}+u_{4,1}}{h^2} & f_{3,2} - \frac{u_{4,2}}{h^2} & f_{3,3} - \frac{u_{4,3}+u_{3,4}}{h^2} \end{bmatrix}

Si las condiciones las tenemos sobre la derivada en el extremo opuesto llegaremos a la misma estructura pero en la parte inferior de la frontera y de la matriz.

Si las condiciones las tenemos sobre derivadas en la otra dirección, podemos llegar también a estas estructuras tomando el orden de variables donde tiene prioridad la variable contraria a la tomada en los casos anteriores.

En el post anterior hablamos sobre condiciones de frontera y su transferencia entre mallas pero no comentamos en el caso de que las condición haga referencia al valor de la derivada y no al de la función: condición de Neumann.

En 1D supongamos que ahora tenemos \frac{\partial^2}{\partial x^2}u = f en [a,b] con u(a)=u_a pero \frac{\partial}{\partial x} = du_b. Suponiendo de nuevo n=8, las ecuaciones nos quedan:

\frac{u_0 -2u_1 + u_2}{h^2} = f_1 para i=1,

\frac{u_1 -2u_2 + u_3}{h^2} = f_2 para i=2,

\frac{u_2 -2u_3 + u_4}{h^2} = f_3 para i=3,

\frac{u_3 -2u_4 + u_5}{h^2} = f_4 para i=4,

\frac{u_4 -2u_5 + u_6}{h^2} = f_5 para i=5,

\frac{u_5 -2u_6 + u_7}{h^2} = f_6 para i=6,

\frac{u_6 -2u_7 + u_8}{h^2} = f_7 para i=7,

La única diferencia con respecto al caso anterior es que, en la primera ecuación, desconocemos el valor de u_0 pero  conocemos el de su primera derivada. Sabemos que:

\frac{u_1 - u_{-1}}{2h} = \frac{d}{dx}u_{0} = du_0,

que, despejando, nos da:

u_1 - u_{-1} = 2h \, du_0 \Leftrightarrow u_{-1} = u_1 - 2h \, du_0,

Como tenemos una incognita mas por determinar, añadimos una nueva ecuación:

\frac{u_{-1} -2u_0 + u_1}{h^2} = f_0 para i=0,

donde reescribimos el valor de u_{-1} según acabamos de determinar:

\frac{ u_1 - 2h \, du_0 -2u_0 + u_1}{h^2} = \frac{ -2u_0 + 2u_1 - 2h \, du_0}{h^2} = f_0.

Por lo tanto,  en forma matricial tenemos:

\frac{1}{h^2} \begin{bmatrix} -1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 &0 \\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0 &0 \\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{bmatrix} \begin{bmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ u_6 \\u_7 \end{bmatrix} = \begin{bmatrix} \frac{1}{2} (f_0 + \frac{2h \, du_0}{h^2}) \\ f_1 \\ f_2 \\ f_3 \\f_4 \\ f_5 \\ f_6 \\ f_7 - \frac{u_8}{h^2}\end{bmatrix}.

De la misma manera, en el caso por el otro extremo, llegariamos a:

\frac{1}{h^2} \begin{bmatrix} -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 &0 \\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0 &0 \\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ u_6 \\ u_7 \\u_8 \end{bmatrix} = \begin{bmatrix} \ f_1 - \frac{u_0}{h^2} \\ f_2 \\ f_3 \\ f_4 \\f_5 \\ f_6 \\ f_7 \\ \frac{1}{2}(f_8 + \frac{2h \, du_8}{h^2})\end{bmatrix}.

En resumen, básicamente hay que hacer dos trabajos: en primer lugar, construir el termino independiente de manera apropiada para incorporar la información de las fronteras; en segundo, llegados a los extremos, escoger entre -2 y -1 en la diagonal en función de si es Dirichlet o Neumann.

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.

En el artíclo [Gingold & Monaghan, 1977] se introducen

La función kernel, W(\vec{r}-\vec{r'},h), es una función que permite interpolar los valores de cualquier propiedad del fluido en función del valor de las partículas del entorno. Su papel es similar al de los diferentes esquemas de diferencias en el ámbito de las Diferencias Finitas o las funciones de forma en los Elementos Finitos.

Existen diferentes funciones kernel: Gaussiana, cuadrática, spline cúbico, quíntica, etc.

La función kernel debe cumplir:

  1. Positiva: W(r-r',h) \geq 0 dentro del dominio.
  2. Soporte compacto: W(r-r',h) = 0 fuera del dominio.
  3. Normalizada: \int W(r-r',h) dr' = 1.
  4. Comportamiento de función delta: \lim_{h \rightarrow 0} W(r-r',h) dr' = \delta(r-r').
  5. Monotona decreciente.

Para derivar, tomamos la derivada analítica de la suma aproximada. De esta manera, como la derivada de la función kernel es conocida, no necesitamos diferencias finitas y el conjunto de ecuaciones PDE pasa a ser ODE.

\nabla f(\vec{r}) = \sum_b \frac{m_b}{\rho_b} f_b W(\vec{r}-\vec{r'},h)

En el artículo  [Gingold & Monaghan, 1977] se presenta por primera vez el método Smoothed Particle Hydrodynamics. Los autores, originalmente, buscaban un método que permitiera tratar problemas en astrofísica asimétricos (sin simetría esférica, sin simetría axial, etc.) . En estos casos, los métodos de diferencias finitas no se adaptan bien, pues requieren elevar el número de puntos en la malla para seguir con la precisión deseada su evolución, lo cual complica enormemente la evolución de las integrales múltiples.

Lo que pensaron es utilizar la descripción Lagrangiana del flujo del fluido que centra su atención en los elementos del fluido. En la discretización, estos elementos se mueven siguiendo las leyes de Newton con fuerzas debidas a los gradientes de presión y a otras fuerzas de cuerpo como la gravedad, rotación o magnéticas. ¿Qué método utilizar para determinar las fuerzas que actuan en un momento determinado sobre un elemento de fluido?

Para empezar, para elementos de fluido de igual masa, el número de elementos por unidad de volumen debe ser proporcional a la densidad. Además, sin ningún tipo especial de simetría, la posicion de los elementos será aleatoria de acuerdo con la densidad. Para recuperar la densidad de la distribución conocida de los elementos es equivalente a recuperar una distribución de probabilidad a partir de una muestra. Existen dos métodos para conseguir esto que funcionan bien con problemas de fluidos: el smoothing kernel method y  la técnica del spline delta. Ambos métodos se pueden pensar como la aproximación de una integral por el procedimiento de Monte Carlo.

octubre 2017
L M X J V S D
« Ago    
 1
2345678
9101112131415
16171819202122
23242526272829
3031