Laplace

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;