diferencias finitas

You are currently browsing articles tagged diferencias finitas.

El caso mas sencillo es cuando tenemos cinco puntos:

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

de manera que:

$latex 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:

$latex 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:

$latex 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}$,

$latex 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}$,

$latex 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}$,

$latex 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}$,

$latex 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:

$latex 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 $latex l_i(x)$ en función de $latex x_0$ de manera que, por ejemplo, tenemos:

$latex 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} =$

$latex =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} = $

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

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

$latex 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 $latex l_i(x)$, obtenemos la siguiente formula en diferencias finitas:

$latex 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:

$latex 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 $latex x=x_0$, $latex x=x_1=x_0+2h$, $latex x=x_3=x_0+4h$ y $latex x=x_4=x_0+6h$:

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

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

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

$latex 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 $latex x=x_0+4h (3h+h)$, tenemos:

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

Tags: , ,

Dados $latex n+1$ puntos:

$latex (x_0,y_0), ldots, (x_n,y_n)$,

definimos el polinomio interpolador de Lagrange:

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

donde:

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

Con esto, tenemos:

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

Supongamos que tenemos $latex n=2$:

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

En este caso, tenemos:

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

con:

$latex l_0(x) = frac{x-x_1}{x_0-x_1}$

$latex l_1(x) = frac{x-x_0}{x_1-x_0}$

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

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

con:

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

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

de manera que:

$latex frac{d}{dx}L(x) = frac{y_1 – y_0}{x_1 – x_0}$.

Con esto, tenemos las aproximaciones de primer orden:

$latex u_x(x_0) = u_x(x_1) = frac{y_1 – y_0}{x_1-x_0}$,

que, en índices, queda:

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

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

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

Supongamos ahora que tenemos tres puntos:

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

En este caso tenemos:

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

$latex 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)$,

$latex 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:

$latex 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)}$,

$latex 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)}$,

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

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

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

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

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

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

$latex 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:

$latex 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 ($latex h_x:=x_{i+1}-x_{i}$ con $latex 0 leq i leq 1$), queda la aproximación de segundo orden:

$latex 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 $latex x$:

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

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

$latex 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 $latex x$ (por lo que tendremos un resultado diferente en función de si la $latex x$ vale $latex x_0$, $latex x_1$ o $latex x_2$) y lo que obtenemos es:

$latex 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:

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

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

$latex 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:

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

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

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

Tags: , ,

Laplaciano en cartesianas:

$latex Delta u = Sigma_i frac{partial^2}{partial x_i^2}u$

1d

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

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

2d

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

$latex i$ fijo:

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

$latex j$ fijo:

$latex 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

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

$latex i,j$ fijos:

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

$latex (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})$

$latex i,k$ fijos:

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

$latex (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})$

$latex j,k$ fijos:

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

$latex (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})$

Tags: , , , , ,

Aunque siempre podemos hacer cambios de coordenadas, vamos a ver como quedan los esquemas de diferencias finitas en sistemas no rectangulares: coordenadas cilíndricas, $latex (rho,phi, z)$, y coordenadas esféricas, $latex (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:

$latex 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:

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

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

$latex + 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} + $

$latex + 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:

$latex 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^2sintheta} frac{partial}{partial theta}u + frac{1}{r^2 sin^2theta} frac{partial^2}{partial phi^2}u = f$

que podemos discretizar como:

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

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

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

$latex + 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} + $

$latex + 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}$

Tags: , , , , ,

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

Empezamos suponiendo que conocemos:

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

$latex frac{partial}{partial y}|_{0,0}u, frac{partial}{partial y}|_{1,0}u$

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

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

Discretizamos:

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

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

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

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

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

$latex 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:

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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:

$latex 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:

$latex 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 $latex 6$ ecuaciones con $latex 6$ incognitas y la matriz tiene rango $latex 6$, por lo que la solución es única.

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

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

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

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

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

Si discretizamos:

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

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

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

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

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

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

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

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

$latex 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:

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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$

$latex 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:

$latex 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:

$latex 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 $latex 9$ ecuaciones con $latex 9$ incognitas pero la matriz tiene rango $latex 8$, por lo que tenemos infinitas soluciones. Hay que conservar.

Tags: , , , , ,

Suponemos $latex Delta u = f$ en $latex 2D$, es decir,

$latex 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 $latex x$ en una frontera, condición Neumann respecto $latex y$ en una frontera y condición Neumann respecto $latex x$ e $latex y$ en dos fronteras.

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

$latex 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, $latex 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 $latex 9$ por lo que existe solución y es única. Desde el punto de vista algebraico, es el punto $latex (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 $latex 9$ hiperplanos

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

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

$latex ldots$

en el espacio $latex mathbb{R}^9$.

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

$latex 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 $latex 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 $latex 2$:

$latex 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 $latex 12$ incognitas ($latex u_{i,j}$ con $latex i=0..3$ y $latex j=1..3$) y el rango de $latex A_2$ es $latex 12$, por lo que la solución, nuevamente, es única.

Para el caso en el que conocemos $latex frac{partial}{partial y}|_{1,0}u, frac{partial}{partial y}|_{2,0}u, frac{partial}{partial y}|_{3,0}u$ en lugar de $latex 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:

$latex 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:

$latex 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 $latex 12$ ecuaciones con $latex 12$ incognitas ($latex u_{i,j}$ con $latex i=1..3$ y $latex j=0..3$) y el rango de $latex A_3$ es $latex 12$, por lo que la solución es única.

Finalmente, suponemos conocidos $latex 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 $latex 7$ ecuaciones mas a las $latex 9$ que ya teniamos por lo que nos queda una matrix $latex A_4 in mathcal{M}(16 times 16)$:

$latex 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 $latex u_{0,0}$ por $latex 4$, y las correspondientes a $latex u_{0,1}, u_{0,2}, u_{0,3}, u_{1,0},u_{2,0}, u_{3,0}$  por $latex 2$, quedando:

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

Tags: , , , , , ,

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

$latex 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 $latex i,j=1,1$,

$latex 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 $latex i,j=1,2$,

$latex 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 $latex i,j=1,3$,

$latex 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 $latex i,j=2,1$,

$latex 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 $latex i,j=2,2$,

$latex 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 $latex i,j=2,3$,

$latex 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 $latex i,j=3,1$,

$latex 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 $latex i,j=3,2$,

$latex 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 $latex i,j=3,3$,

de donde:

$latex 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):

$latex 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 $latex u_{0,1}, u_{0,2}, u_{0,3}$ conocemos $latex 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:

$latex 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 $latex i,j=0,1$

$latex 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 $latex i,j=0,2$

$latex 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 $latex i,j=0,3$

y

$latex 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$

$latex 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$

$latex 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:

$latex 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:

$latex 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:

$latex 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:

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

Tags: , , , , ,

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 $latex 1D$ supongamos que ahora tenemos $latex frac{partial^2}{partial x^2}u = f$ en $latex [a,b]$ con $latex u(a)=u_a$ pero $latex frac{partial}{partial x} = du_b$. Suponiendo de nuevo $latex n=8$, las ecuaciones nos quedan:

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

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

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

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

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

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

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

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

$latex frac{u_1 – u_{-1}}{2h} = frac{d}{dx}u_{0} = du_0$,

que, despejando, nos da:

$latex 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:

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

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

$latex 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:

$latex 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:

$latex 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 $latex -2$ y $latex -1$ en la diagonal en función de si es Dirichlet o Neumann.

Tags: , , , , ,

En $latex 1D$ tenemos:

$latex u_{xx} = f(x)$,

con $latex 0 < x < 1$ y $latex u(0)= u(1)=0$. Particionamos el dominio en $latex n$ subintervalos introduciendo los puntos $latex x_i = ih$ con $latex 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 $latex x_i$):

$latex 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 $latex xi_i in (x_{i-1}, x_{i+1})$ por lo que:

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

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

$latex frac{1}{h^2} begin{bmatrix} -2 & 1 & & & & & \ 1 & -2 & 1 & & & & \ & 1 & -2 & 1 & & & \ & & cdot & cdot & cdot & & \ & & & cdot & cdot & cdot & \ & & & & 1 & -2 & 1 \ & & & & & 1 & -2end{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:

$latex 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 $latex xi_i in (x_{i-2},x_{i+2})$.

En $latex 2D$ tenemos:

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

con $latex 0 < x < 1$, $latex 0 < y < 1$ y $latex u=0$ en el cuadrado unidad. Particionamos el dominio en $latex n_x times n_y$ subintervalos introduciendo los puntos $latex (x_i,y_j) = (ih_x, jh_y)$ con $latex h_x=frac{1}{n_x}$ y $latex 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:

$latex 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 $latex v_{0,j} = v_{n_x,j} = v_{i,0} = u_{i,n_y} = 0$. Tenemos un error local de truncamiento del orden de $latex O(h_x^2+h_y^2)$.

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

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

En $latex 3D$ tenemos:

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

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

$latex 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:

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

si $latex h_x = h_y = h_z = h$.

Tags: , ,

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

La función kernel, $latex 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: $latex W(r-r’,h) geq 0$ dentro del dominio.
  2. Soporte compacto: $latex W(r-r’,h) = 0$ fuera del dominio.
  3. Normalizada: $latex int W(r-r’,h) dr’ = 1$.
  4. Comportamiento de función delta: $latex 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.

$latex nabla f(vec{r}) = sum_b frac{m_b}{rho_b} f_b W(vec{r}-vec{r’},h)$

Tags: , , , ,

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.

Tags: , , , , , , , ,

FireStats icon Powered by FireStats