Ecuación de Laplace en Octave
Computación Científica
Facultad de Ciencias Matemáticas
Universidad Nacional Mayor de San Marcos
22 de octubre de 2010
Resumen
En este proyecto resolvemos numéricamente la ecuación de Laplace usando el método de diferencias finitas. Hacemos uso de
Octave para presentar los resultados.
Introducción
Las ecuaciones en derivadas parciales(PDE's) mas estudiadas se clasifican en hiperbólicas, parabólicas y elípticas. Dado que estas representan modelos muy aproximados de fénomenos en la naturaleza, es necesario hallar su solución, al menos de manera aproximada.
Presentaremos un esquema númerico de solución para una ecuacíon de tipo elíptico, la ecuación de Laplace.
Lo siquiente es el código desarrollado en Octave y un resultado obtenido. Éstare presentando mas detalles.
function V=Laplace(xi,yi,n,m,fr,f)
% laplaciano(u)=f(x,y)
% u(frontera)=h(x,y)
%n: numero de intervalos en x
%n+1 puntos x
%m: numero de intervalos en y
%m+1 puntos y
%------------------------
%tamaños de paso;
h=(xi(2)-xi(1))/n;
g=(yi(2)-yi(1))/m;
%-------------------------
%mallado
x=linspace(xi(1),xi(2),n+1);
y=linspace(yi(1),yi(2),m+1);
V=zeros(n+1,m+1);
%---------------------------
%funciones
F=inline(f,'x','y');
H=inline(fr,'x','y');
%-----------------------------
%construyendo la matriz A
%coeficientes
a=2*(1/(h^2)+1/(g^2));
b=-1/(h^2);
c=-1/(g^2);
%matrices de tamaño n-1
M_1=trigeneral(a,b,n-1);
M_2=trigeneral(c,0,n-1);
%m-1 bloques
A=trigeneral(M_1,M_2,m-1);
%--------------------------------
%en la frontera
for i=1:n+1
V(i,1)=H(x(i),yi(1));
V(i,m+1)=H(x(i),yi(2));
end
for i=1:m+1
V(1,i)=H(xi(1),y(i));
V(n+1,i)=H(xi(2),y(i));
end
%-------------------------
%los V(i,j) incognitas so los que estan en el interior
% 2<=i<=n ; 2<=j<=m
%Ahora formaremos elvector B del sistem : Av=B
B1=zeros((n-1)*(m-1),1);
B2=zeros((n-1)*(m-1),1);
B3=zeros((n-1)*(m-1),1);
B4=zeros((n-1)*(m-1),1);
%---------------------------
%la parte con f
for j=2:m
for i=2:n
B1( (j-2)*(n-1) +(i-1))=F(x(i),y(j));
end
end
B2(1:n-1:(n-1)*(m-1)-(n-2))=(-b)*V(1,2:m);
B3(n-1:n-1:(n-1)*(m-1))=(-b)*V(n+1,2:m);
B4(1:n-1)=(-c)*V(2:n,1);
B4((n-1)*(m-2)+1:(n-1)*(m-1))=(-c)*V(2:n,m+1);
%------------------
B=B1+B2+B3+B4;
w=A\B;
for j=2:m
for i=2:n
V(i,j)=w((j-2)*(n-1) +(i-1));
end
end
[X,Y]=meshgrid(x,y);
surf(X,Y,V')
function M=trigeneral(A,B,num)
m=size(A,1);TM=m*num;M=zeros(TM);
M(1:m,1:m)=A;
M(1:m,m+1:2*m)=B;
for i=m+1:m:TM-2*m+1,k=i+m-1;M(i:k,i:i+m-1)=A;M(i:k,i-m:i-1)=B;...
M(i:k,i+m:i+2*m-1)=B;end
k=TM+1-m;
M(k:TM,TM+1-m:TM)=A;
M(k:TM,TM-2*m+1:TM-m)=B;
