Diff finitas calculo 3

´ ´ ´ PRACTICAS DE CALCULO NUMERICO III ´ PRACTICA 6: la ecuaci´n de Laplace en un dominio rectangular mediante o diferencias finitas
Vamos a comparar los m´todos de diferencias finitas y elementos finitos (con la fore mulaci´n del m´todo de Galerkin) para resolver la ecuaci´n de Laplace con condiciones o e o de contorno de tipo Dirichlet. Este es el problema correspondiente a la obtenci´n de o ladistribuci´n de temperaturas en una placa rectangular con unas distribuci´n fija de o o temperatura en los bordes de la placa. Comenzaremos con la aproximaci´n de diferencias finitas. o Nuestro problema de prueba ser´: a

y 1
φ(x)

u=0

ux x + uy y =0

ψ(y)

0
φ(x) ψ(y)

u=0

1

x

= 100sinh(pi/10)sin(pi x/10)/sinh(pi) = 100sin(pi/10)sinh(pi y/10)/sinh(pi)

cuya soluci´n anal´o ıtica es: uan (x, y) = 100 sinh(πy/10) sin(πx/10) sinh(π)

1

Diferencias finitas:

La aproximaci´n de diferencias finitas (con un esquema de 5 puntos para el laplaciano) de o la ecuaci´n de Poisson o ∇2 u ≡ uxx + uyy = f (x, y) responde, como sabemos, al siguiente m´todo: e u(x + h, y) + u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y) ≈ h2 f (x, y) En esta expresi´n asumimos, porsimplicidad, un mismo espaciado h en las dos direco ciones. Tomaremos el origen de coordenadas en el v´rtice inferior izquierdo de la placa. En e (1)

esta pr´ctica resolveremos la ecuaci´n de Laplace tanto mediante esta regla de 5 puntos a o como mediante la de 9 puntos. Nuestra rutina de diferencias finitas (dfelip.m) tendr´ la siguiente sintaxis: a
function [x,y,u,lap]=dfelip(h,a,b,fx0,fxa,fy0,fyb)Donde el significado de los inputs es el siguiente:
% % % % % % % h: a: b: fx0: fxa: fy0: fyb: espaciado del reticulo en cada una de las dos direcciones. longitud de la placa en la direccion x longitud de la placa en la direccion y funcion correspondiente a la condicion de Dirichlet para x=0 Idem para x=a Idem para y=0 Idem para y=b

y el significado de los outputs es:
% % % % x: y: u: lap:vector vector matriz matriz de coordenadas en la direccion x de cordenadas en la direccion y del campo de temperaturas en los puntos del reticulo laplaciana obtenida

Damos el programa dfelip expl´ ıcitamente:
function [x,y,u,L]=dfelip(h,a,b,fx0,fxa,fy0,fyb,rule) [x,y,u,IP,nv,in,jn,aa]=gridrD(h,a,b,fx0,fxa,fy0,fyb); [L,B]=feval(rule,u,IP,nv,in,jn,aa); v=L\B; for a=1:nv u(in(a),jn(a))=v(a); end;Vemos que el programa llama a dos funciones. La primera de ellas, gridrD.m, crea los vectores x e y, inicializa la funci´n soluci´n considerando las condiciones de contorno o o y hace los mismo con otras variables necesarias para la construcci´n del sistema lineal. o La segunda se utiliza para la construcci´n de la matriz laplaciana L y la de t´rminos o e independientes B; cuandorule=’lapla5’, la l´ ınea [L,B]=feval(rule,u,IP,nv,in,jn,aa); es, por supesto, equivalente a [L,B]=lapla5(u,IP,nv,in,jn,aa); Describiremos m´s adelante en qu´ consiste la rutina lapla5.m. Antes, es necesario ena e tender la utilidad de gridrD.m. La rutina gridrD.m crea la matriz u(i, j), i = 1...nx, j = 1...ny que va a representar a nuestro dominio rectangular y donde se almacenar´n los valores de la soluci´nnum´rica. a o e Inicializaremos la matriz u(i, j) rellenando todas sus posiciones con un valor num´rico e (IP) que no se vaya a alcanzar en la frontera del dominio; por ejemplo IP=10000. A continuaci´n especificaremos las condiciones de contorno (Dirichlet) del problema fijando o los valores u(1, j), u(nx, j), u(i, 1), u(i, ny) (nx y ny el n´mero en las direcciones x e y). u Con una la elecci´n de IPconveniente, los valores asignados distinguir´n a los puntos (i, j) o a interiores de los que no lo son, pues s´lo los interiores verificar´n que u(i, j) =IP. o a

Como sabemos, la resoluci´n mediante diferencias finitas desemboca en la resoluci´n de o o un sistema lineal en el que las inc´gnitas ser´n los valores num´ricos de la soluci´n, u(i, j), o a e o en los puntos interiores. Parte de...