""
8.5 Ejercicio. Resolución de la ecuación de Laplace en un dominio bidimensional
Resolver la ecuación de Laplace en un dominio rectangular por diferencias finitas con condiciones de contorno Dirichlet. La ecuación de Laplace es la ecuación elíptica más sencilla: ∇2φ=0, que formulada en dos dimensiones se convierte en:
Si se toman diferencias finitas centradas de segundo orden para puntos dependiendo de i y j llegamos a la ecuación en diferencias siguiente:
φ(i+1,j)-2φ(i,j)+φ(i-1,) |
|
Δ x2 |
|
+ |
φ(i,j+1)-2φ(i,j)+φ(i,j-1) |
|
Δ y2 |
|
=0 |
Esta ecuación puede ser expresada en forma de sistema de ecuaciones lineales Aϕ=b, donde b es un término independiente que aparece al aplicar las condiciones de contorno y ϕ es la matriz de incógnitas φ expresada en forma de vector columna. La traslación del problema bidimensional a un vector de incógnitas es un paso que nos puede costar de entender pero si tenemos una matriz de incógnitas el tensor del sistema tendrá tres dimensiones con lo que ya no tendremos rutinas escritas para resolver el sistema.
Usaremos como parámetros dx=2, dy=3, n=50, m=50. No podemos utilizar muchos más elementos porque se nos quedaríamos sin memoria. Esta matriz de incógnitas se va a convertir en un vector nנm, lo que significa que la matriz del sistema va a tener (nנm)2 elementos. Para un dominio de 100 por 100 puntos llegamos a 108 puntos. Sería un buen caso de aplicación de matrices sparse
8.5.1 Guía para la resolución del ejercicio
- Escribir una función place que implemente la transformación de la matriz al vector. La entrada serán los índices i y j y el número de elementos por columna n. La salida de la función será la posición en el vector posterior: place=i+n(j-1).
- Crear la matriz del sistema con la ecuación en diferencias y la función creada en el apartado anterior
- Poner las condiciones de contorno al gusto en el vector del término independiente y resolver el sistema lineal
- Para resolver el sistema lineal del modo usual basta con hacer A.
- Para resolver el sistema con matrices sparse primero creamos la matriz sparse con:
spA=sparse(A)
y luego resolveremos el sistema del modo usual. Es una buena idea eliminar la matriz del sistema de la memoria.
8.5.2 Solución del ejercicio
Primero definimos las constantes del problema dx=2;dy=3;n=50;m=25;
A continuación creamos una función que reordena cualquier elemento de una matriz bidimensional en un vector en el que se concatenan las columnas. function place=place(i,j,n)
place=i+n*(j-1);
end
Ahora definimos la matriz del sistema teniendo en cuenta que la incógnita no va a ser una matriz de dos dimensiones sino un vector. Esto hace que la matriz del sistema tenga en realidad nmנnm elementos. A=zeros(n*m,n*m);
for i=1:n*m
A(i,i)=1;
end
for i=2:n-1
for j=2:m-1
A(place(i,j,n),place(i,j,n))=-2*(1/(dx*dx)+1/(dy*dy));
A(place(i,j,n),place(i-1,j,n))=1/(dx*dx);
A(place(i,j,n),place(i+1,j,n))=1/(dx*dx);
A(place(i,j,n),place(i,j-1,n))=1/(dy*dy);
A(place(i,j,n),place(i,j+1,n))=1/(dy*dy);
end
end
Una vez definida la matriz del sistema creamos el vector del término independiente que contiene las condiciones de contorno. i=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
end
i=n;
for j=1:m
b(place(i,j,n))=-sin(pi*(j-1)/(m-1));
end
j=1;
for i=1:n
b(place(i,j,n))=sin(pi*(i-1)/(n-1));
end
j=n;
for i=1:n
b(place(i,j,n))=0;
end
Y finalmente resolvemos el sistema lineal. T=A \b';
T=reshape(T,n,m);
contour(T);
El resultado del script son las figuras 8.3 y 8.4.
Figure 8.3: Superficie solución
En el capítulo dedicado a la representación gráfica de soluciones hemos hablado sobre lo poco útiles que suelen ser las superfícies coloreadas para dar un resultado. Aunque una de estas representaciones pueda parecer muy explicativa e intuitiva nunca debemos perder de vista que lo que nos importa es mostrar un resultado y una superfície no lo consigue. Por muy bien que coloquemos los colores o las barras no conocemos el valor de cada punto con precisión y si encima colocamos líneas de nivel dotamos al gráfico de un exceso de información.
Aunque sea mucho más espartano y pueda parecer inapropiado la solución al problema son las curvas de nivel, en inglés contours . Aunque en un principio cueste un poco más entender la solución estamos ofreciendo muchísima más información de un modo netamente más simple. Se podría decir que no hay motivos para no utilizar las curvas de nivel a parte de intentar demostrar nuestra pericia en el uso de Matlab.
Figure 8.4: Superficie solución
Los gráficos de curvas de nivel pueden ajustarse perfectamente a nuestras necesidades. Una práctica muy común y beneficiosa es la de embeber el valor de cada curva dentro de la misma para no obligar al lector a mirar la leyenda.
8.5.2.1 Mejorar la solución. Cálculo de tiempos.
Si analizamos el código de la solución es evidente que estamos utilizando pocas funciones de la biblioteca. Además estamos haciendo una de las pocas cosas casi prohibidas en Matlab que es crear una matriz por fuerza bruta.
Un primer paso para mejorar la ejecución de un código es comprobar qué partes consumen más recursos e intentar mejorarlas aisladamente. Si no somos capaces de conseguirlo nos replantearemos la solución de un modo global. Para calcular los tiempos llenaremos el código con las funciones tic y toc. dx=2;
dy=3;
n=50;
m=50;
function place=place(i,j,n)
place=i+n*(j-1);
end
A=zeros(n*m,n*m);
tic
for i=1:n*m
...
end
toc;disp('Creacion de la matriz del sistema'),disp(toc);tic;
i=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
...
b(place(i,j,n))=0;
end
toc;disp('Creacion del vector b'),disp(toc);tic;
T=A \b';
toc;disp('Resolucion del sistema'),disp(toc);
T=reshape(T,n,m);
Los tiempos en cada paso son: >> ejercicio4
Creacion de la matriz del sistema
1.0611
Creacion del vector b
0.017038
Resolucion del sistema
2.1457
Parece que debemos mejorar la construcción de la matriz del sistema y la resolución del mismo. ¿Cómo? El primer paso suele ser analizar la forma de la matriz con la función spy. >> spy(A)
Que tiene como resultado el siguiente patrón (figura 8.5):
Figure 8.5: Patrón de elementos no nulos de la matriz del sistema
Esto nos demuestra que es una matriz n-diagonal por bandas lo que entra en la definición de matriz sparse. Parece que la solución a nuestros ligeros problemas de velocidad pueden solucionarse con el uso de este tipo de matrices.
8.5.2.2 Resolución del problema mediante matrices sparse(+)
Uno de los errores cometidos en la primera solución propuesta es que no utiliza las rutinas de creación de matrices disponibles. La alternativa es utilizar bucles for para crearla por fuerza bruta, solución de ningún modo recomendable. La mejor opción suele ser utilizar la función diag pero como estamos tratando con matrices sparse utilizaremos spdiags.
8.5.2.3 Resolución del problema con un método iterativo.
Para resolver el problema por un método directo es estrictamente necesario romper la matriz para convertirla en un vector. Esto hace que no podamos utilizar los operadores diferenciales en forma de matriz tal como veremos en el ejercicio 8.6. Un truco para mantener la integridad de la matriz y poder utilizar estos operadores es resolver el problema con un método iterativo. Aunque estemos multiplicando una matriz llena por matrices casi vacías el resultado se acelera considerablemente.
La teoría de este método es bastante sencilla. Si aplicamos el operador laplaciano a la matriz bidimensional el problema numérico se reescribe del siguiente modo.
Dxφ+φDy⊤=0
Este planteamiento carece de solución analítica pero sirve para plantear un problema iterativo en el que sólo hay que evaluar la expresión. Empezaremos la subrutina del mismo modo que la anterior, definiendo la función que recoloca los elementos de la matriz bidimensional. dx=2;dy=3;global n=50;global m=50;
% resolucion del problema con un solver iterativo construccion de los
% operadores en diferencias
function place=place(i,j,n)
place=i+n*(j-1);
end
Ahora definiremos las matrices de los operadores diferenciales. Las definimos como variables globales al igual que las dimensiones de la matriz porque son necesarios dentro de la función que calcula el paso de iteración y ésta puede tener sólo un argumento de entrada. global DX=diag([1,-2/(dx*dx).*ones(1,n-2),1],0).+...
diag([0,1/(dx*dx).*ones(1,n-2)],1).+...
diag([1/(dx*dx).*ones(1,n-2),0],-1);
global DY=diag([1,-2/(dy*dy).*ones(1,n-2),1],0).+...
diag([0,1/(dy*dy).*ones(1,n-2)],1).+...
diag([1/(dy*dy).*ones(1,n-2),0],-1);
A continuación la función que calcula el término en cada iteración. Notemos el hecho de que el argumento que recibe la función debe ser un vector por requerimientos de la rutina de resolución. Para calcular sobre las matrices utilizamos la función reshape. function [AT]=calcstuff(T)
global n
global m
global DX
global DY
rT=reshape(T,n,m);
AT=reshape((DX*rT)+(rT*DY'),n*m,1);
end
Ahora es el turno del término independiente que se define gracias a la función place definida anteriormente. i=1;j=1;b=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
end
i=n;
for j=1:m
b(place(i,j,n))=-sin(pi*(j-1)/(m-1));
end
j=1;
for i=1:n
b(place(i,j,n))=sin(pi*(i-1)/(n-1));
end
j=n;
for i=1:n
b(place(i,j,n))=0;
end
Y finalmente resolvemos el sistema. tic;T=pcg('calcstuff',b',tol=1e-6,maxit=100);
disp('Tiempo de calculo'),disp(toc)
El método utilizado para resolver el problema es el del ``Gradiente Conjugado Precondicionado''. Para utilizarlo debemos tener en cuenta que el número de iteraciones por defecto no es sufciente; lo corregimos y lo situamos en 100. El tiempo de proceso ha pasado de más de tres segundos a 0.249292. Este programa no sería posible sin el uso inteligente de las variables globales tal como las hemos utilizado en otros ejemplos.
8.6 Ejercicio. Un problema de calor unidimensional
Este ejercicio requiere la función trisolve que está disponible en Octave. Matlab no tiene ninguna rutina para resolver directamente sistemas tridiagonales porque trata este tipo de matrices como sparse. Para resolver el mismo problema en Matlab utilizaremos la función spdiags para montar la matriz sparse y luego resolveremos el sistema de ecuaciones del modo usual.
Este ejercicio es una pequeña demostración de hasta dónde puede llegar el ahorro de código con Matlab. Se trata de resolver la ecuación del calor unidimensional con condiciones de contorno Dirichlet. La discretización espacial serán diferencias finitas de sevundo órden y la temporal será un esquema Crank Nicholson. Tomaremos este esquema porque permite usar saltos de tiempo más grandes, de hecho nos permitirá llegar al resultado sin tener que calcular ninguna iteración temporal. La solución propuesta tiene sólo unas pocas líneas.
Una solución tan breve requiere algo de preparación matemática de modo que hay que estar muy atento al planteamiento analítico del problema.
Discretización de la ecuación completa
La ecuación del calor unidimensional es la EDP parabólica más sencilla:
∂tφ=∂xxφ
Como cualquier EDP parabólica se puede formular de la siguientes manera:
Si utilizamos un esquema Crank Nicholson la discretización temporal será de la forma:
Discretizando también el lado derecho de la ecuación con diferencias finitas de segundo orden llegamos a la ecuación discretizada:
φin+1- |
|
漂R> 缂R> 缂R> 輯FONT> |
φi+1n+1-2φin+1+φi-1n+1 |
|
Δ x2 |
|
?> ?> ?> ?NT> |
=φin+ |
|
漂R> 缂R> 缂R> 輯FONT> |
|
?> ?> ?> ?NT> |
Formulación matricial del problema.
La ecuación anterior puede escribirse del modo siguiente:
Agrupando términos, diciendo que B?=I?Δ t/2A y utilizando un Δ x=1 (en el fondo es sólo adimensionalizar):
B-φn+1=B+φn
donde la matriz A tiene la forma:
A= |
漂R> 缂R> 缂R> 缂R> 缂R> 缂R> 缂R> 輯FONT> |
-2 |
1 |
0 |
0 |
⋯ |
0 |
1 |
-2 |
1 |
0 |
|
0 |
0 |
1 |
-2 |
1 |
|
0 |
⋮ |
|
|
⋰ |
|
⋮ |
0 |
|
0 |
1 |
-2 |
1 |
0 |
⋯ |
0 |
0 |
1 |
-2 |
|
?> ?> ?> ?> ?> ?> ?> ?NT> |
Condiciones de contorno e iniciales.
- Condiciones de contorno: φ=0 en x=0 y φ=1 en x=10
- Condiciones iniciales: φ(x)=0 en los puntos interiores.
- Dar una solución al problema para cualquier tiempo inferior a t=10.
8.6.1 Guía para la resolución del ejercicio
La dificultad de la resolución de este problema radica en la correcta situación de las condiciones de contorno en la matriz, el uso correcto de las submatrices y entender el funcionamiento de trisolve.
8.6.2 Solución del ejercicio (Octave)
Parte del ejercicio es entender qué estamos haciendo para llegar a la solución. Lo bueno de abreviar en la escritura es que además el código resultante suele ser mucho más rápido. phi=zeros(10,1); phi(10,1)=1;
dt=0.1;
for i=1:50 #Activate loop in case of precision leak
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=trisolve(-[0.5*dt*ones(8,1);0],[1;(1+dt)*ones(8,1);1],...
-[0;0.5*dt*ones(8,1)],phi);
end
De este modo llegamos al resultado del problema con una precisión aceptable. La figura del perfil de temperaturas es:
Figure 8.6: Figura solución
8.6.3 Solución del ejercicio (Matlab)
Aunque la resolución de sistemas tridiagonales es un clásico del cálculo numérico Matlab no cuenta con ninguna función para ello. El motivo es que decidieron incluir las matrices tridiagonales dentro de las matrices sparse. Las matrices tridiagonales o con pocas bandas se resuelven eficientemente mediante un método directo con lo que en este caso Matlab tiene un error de diseño.
No tendría ningún sentido programarse una función que resolviera sistemas tridiagonales en Matlab porque todas estas rutinas son en el fondo interfaces a bibliotecas en Fortran y en C por motivos de velocidad.
En este caso la dimensión de la matriz no justifica el uso del almacenamiento sparse de modo que utilizaremos la función diag para crear una matriz según sus diagonales: phi=zeros(10,1); phi(10,1)=1; dt=0.1;
B=diag([1,(1+dt)*ones(1,8),1],0)+...
diag(-[0.5*dt*ones(1,8),0],-1)+...
diag(-[0,0.5*dt*ones(1,8)],1);
for i=1:50
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=B \phi;
end
8.6.4 Comprobación de la evolución temporal
La siguiente modificación del programa permite comprobar la evolución de los pasos temporales sacando por pantalla varios estados temporales de la solución. Aunque el problema es estable para saltos temporales muy grandes preferiremos acortarlo por motivos de precisión. Representaremos la solución en los tiempos t=[0.1,1,2,5,10,20,50,100]. Muy atentos al modo de definir la condición lógica que indica cuándo pintar la curva: phi=zeros(10,1); phi(10,1)=1;dt=0.1;
# times to output
tout=[0.1,1,2,5,10,20,50,100];
hold on
for i=1:500
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=trisolve(-[0.5*dt*ones(8,1);0],[1;(1+dt)*ones(8,1);1],...
-[0;0.5*dt*ones(8,1)],phi);
if sum((dt*i)*ones(1,length(tout))==tout)==1
plot(phi)
end
end
hold off
title('Solucion de la ecuacion del calor unidimensional')
xlabel('Longitud (adimensional)')
ylabel('Temperatura (adimensional)')
legend({'0.1','1','2','5','10','20','50','100'},2)
Llegamos finalmente a la siguiente gráfica:
Figure 8.7: Evolución del perfil de temperaturas
- 1
- El uso de lenguajes ``pegamento'' es una práctica cada día más habitual. Escribir un código enteramente en C o Fortran es bastante costoso tanto por el tiempo de desarrollo como por el esfuerzo necesario. Muchas veces se empieza con un prototipo escrito con un lenguaje de RAD (Rapid Application Development) y posteriormente se reescriben sus partes.
El mayor obstáculo en los proyectos grandes es que deben participar varios desarrolladores. Cada uno puede tener sus preferencias y las guías de estilo no siempre son una solución efectiva porque sólo arreglan el formato del código escrito. Todas estas dificultades unidas al aumento de velocidad de lenguajes como Python o Matlab ha llevado a la creación del concepto de ``lenguaje pegamento''.
Se trata de dejar a medias el código entre el prototipo y el programa definitivo. Partiendo del prototipo se van analizando las partes que consumen más recursos y, sin cambiar las cabeceras ni el esquema de variables, se reescriben en algún lenguaje rápido. Se para cuando se consigue un equilibrio entre nivel de interactividad, velocidad y coste de mantenimiento de código.
Probablemente el mejor ``lenguaje pegamento'' sea Python gracias a Pyrex, SWIG y F2Py. El primero es capaz de convertir el código Python en C++ automáticamente (rendimiento máximo y coste cero) y el segudo y el tercero son generadores automáticos de interfaces, SWIG para C y C++ y F2Py para Fortran 77 y Fortran 95.
- 2
- Este ejercicio se ha resuelto en un AMD Athlon 3000+.
""
