noviembre 2013

You are currently browsing the monthly archive for noviembre 2013.

El blog, que empezó para servir de mero canal de comunicación asíncrono entre dos o tres personas, acaba de superar las mil visitas mensuales. Aunque no conozco a casi ninguno de mis lectores,

¡GRACIAS A TODOS! 🙂

blogStats   blogMap

Tags:

psiAlpha

phiAlpha1

psiAlpha2

psiAlpha3

Tags: , , , ,

Seguimos con el post anterior. Vamos a utilizar ahora siete puntos separados de la siguiente manera:

$latex x_0$

$latex x_1:=x_0+n^2 h$

$latex x_2:=x_0 + (n^2+n) h$

$latex x_3:=x_0 + (n^2+n+1) h$

$latex x_4:=x_0 + (n^2+n+2) h$

$latex x_5:=x_0 + (n^2+2n+2) h$

$latex x_6:=x_0 + (2n^2+2n+2) h$

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: , , , , ,

Ya hace casi dos meses que no escribo nada en el blog :-(. Es hora de añadir algo.

El Block Cyclic Reduction es un algoritmo rápido para resolver sistemas de ecuaciones de manera directa cuya matriz tiene estructura tridiagonal por bloques.

Supongamos que tenemos una discretización de un volumen mediante cinco puntos para las $latex x$, tres puntos para las $latex y$ y siete puntos para las $latex z$ (necesitamos que sean números impares).

Para la dimensión $latex z$ tenemos:

$latex A_3 = left(
begin{array}{ccccccc}
6 & -1 & 0 & 0 & 0 & 0 & 0 \
-1 & 6 & -1 & 0 & 0 & 0 & 0 \
0 & -1 & 6 & -1 & 0 & 0 & 0 \
0 & 0 & -1 & 6 & -1 & 0 & 0\
0 & 0 & 0 & -1 & 6 & -1 & 0 \
0 & 0 & 0 & 0 & -1 & 6 & -1\
0 & 0 & 0 & 0 & 0 & -1 & 6
end{array}
right)$

$latex -I_3 =left(
begin{array}{ccccccc}
-1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & -1 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & -1 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & -1 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & -1 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & -1 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & -1
end{array}
right)$

Para la dimensión $latex y$ no queda:

$latex A_2 = left(
begin{array}{ccc}
A_3 & -I_3 & 0 \
-I_3 & A_3 & -I_3 \
0 & -I_3 & A_3
end{array}
right); ,,, -I_2 = left(
begin{array}{ccc}
-I_3 & 0 & 0 \
0 & -I_3 & 0 \
0 & 0 & -I_3
end{array}
right)$

Y, finalmente, para la dimensión $latex x$ obtenemos:

$latex A_1 = left(
begin{array}{ccccc}
A_2 & -I_2 & 0 & 0 & 0 \
-I_2 & A_2 & -I_2 & 0 & 0 \
0 & -I_2 & A_2 & -I_2 & 0 \
0 & 0 & -I_2 & A_2 & -I_2 \
0 & 0 & 0 & -I_2 & A_2
end{array}
right)$.

Ahora todo es fácilmente generalizable a $latex n$ dimensiones $latex x_1, x_2, ldots, x_n$ con $latex t_1, t_2, ldots, t_n$ puntos en cada dimensión, ya que nos queda $latex A_i$ una matriz tridiagonal por bloques con $latex t_i times t_i$ bloques donde cada uno está formado por matrices $latex A_{i-1}$ en la diagonal, $latex -I_{i-1}$ arriba y debajo de la diagonal y $latex 0$ en el resto, todas de dimensión $latex t_{i-1} times t_{i-1}$ con $latex i=1,ldots, n$. En $latex A_n$ tenemos el stencil $latex (-1,2n,-1)$.

En el caso que estamos tratando, para cada fila par de $latex A_1$ podemos multiplicarla por $latex A_2$ y sumarla a las filas impares anterior y posterior, de manera que eliminamos las variables pares:

$latex -v^1_{j-2} + ((A_2)^2 – 2I_2)v^1_j – v^1_{j+2} = b^1_{j-1} + A_2 b^1_j + b^1_{j+1}$,

(el $latex ^1$ hace referencia a las $latex x$ y no corresponde a potencias… Para resolver las $latex k$ impares simplemente:

$latex A_2 v^1_k = b^1_{k} + v^1_{k-1} + v^1_{k+1}$. )

Y resulta que definiendo $latex B_2:=(A_2)^2 – 2 I_2$, nos queda una matriz $latex B_1$ de la misma estructura que $latex A_1$, por lo que podemos aplicar este metodo de manera recursiva hasta llegar a tener una sola ecuación (para esto, no basta con que tengamos un número impar de nodos, necesitamos algo como $latex N=2^{k+1}-1$).

Llegados a este punto, si estamos con $latex v^n$, resolvemos como queramos, pero si estamos en otra variable $latex v^i$ con $latex ineq n$, la matriz vuelve a tener la misma estructura, pero con una dimensión menos, y volvemos a aplicar el mismo algoritmo.

Tenemos dos lugares distintos donde aplicar recursividad: el primero por ir agrupando de tres en tres las ecuaciones de una dimensión de trabajo ($latex A_i rightarrow B_i rightarrow ldots rightarrow F_i$, $latex F$ por poner algo :-)); el otro cuando, llegados a una sola ecuación aplicando la recursividad anterior, nos quedan mas dimensiones por resolver ($latex B_{j+1}, B_{j+2}, ldots $).

Tags: ,

FireStats icon Powered by FireStats