1% Nombres: -Rafael An�bal Calder�n Melara CM110992% -Estely Dolores M�ndez Moreno MM0920434% Fecha de creaci�n: 15 de Noviembre de 20165% Esta funci�n aproxima la soluci�n de la ecuaci�n de Poisson6% f_xx(x,y) + f_yy(x,y) = f(x,y); con x en [a,b], y en [c,d]7% sujeto a las condiciones iniciales8% u(x,y) = g(x,y) si x=a o x=b & y pertenece a [c,d]9% u(x,y) = g(x,y) si y=c o y=d & x pertenece a [a,b]10% El algoritmo utiliza el m�todo de diferencias finitas11% y a su vez el m�todo de Gauss Seidel.12% Notar que se utiliza la norma infinita en nuestro c�digo,13% pero se podr�a utilizar cualquier otra norma.14% INPUT: -- f: Funciones15% -- a,b,c,d: Valores frontera.16% -- n>=3, m>=3: Enteros que sirven para determinar las17% dimensiones de la cuadr�cula.18% -- tol: Tolerancia requerida.19% -- N: M�ximo n�mero de iteraciones.2021% OUTPUT: -- Tabla con las aproximaciones para u(x,y) y los puntos22% x & y.23% -- l: n�mero de iteraciones requeridas para obtener24% la soluci�n.2526function [] = PEFD(f,a,b,c,d,m,n,tol,N)2728% h y k representan el tama�o de paso para la construcci�n de los29% puntos x_i e y_j respectivamente.30h = (b-a)/n;31k = (d-c)/m;3233% Se crean x & y que son vectores que contienen las coordenadas34% para los distintos puntos de la cuadr�cula que utiliza35% el m�todo.36x = zeros(1,n-1);37for i = 1:n-138x(i) = a+i*h;39end40y = zeros(1,m-1);41for j = 1:m-142y(j) = c+j*k;43end4445% w es la matriz de ceros que luego se ir� actualizando con46% las aproximaciones w_ji de la soluci�n u(x_i,y_j).47w = zeros(n-1,m-1);4849% Se calculan lambda y mu para hacer m�s eficiente el c�digo en t�rmios de50% las operaciones que se realizan.51lambda = (h^2)/(k^2);52mu = 2*(1+lambda);5354% El bucle while se utiliza para realizar iteraciones del m�todo55% de Gauss-Seidel Optimizado56l = 1;57while l <= N5859% Aqu� se actualiza el elemento 1 de la �ltima columna de w.60z = (-(h^2)*f(x(1),y(m-1))+fung(a,y(m-1))+lambda*fung(x(1),d)+lambda*w(1,m-2)+w(2,m-1))/mu;61NORM = norm(z-w(1,m-1),inf);62w(1,m-1)=z;6364% El siguiente for actualiza los valores de w, desde el segundo65% al pen�ltimo elemento de w, de la �ltima columna.66for i = 2:n-267z = (-(h^2)*f(x(i),y(m-1))+lambda*fung(x(i),d)+w(i-1,m-1)+w(i+1,m-1)+lambda*w(i,m-2))/mu;68if norm(w(i,m-1)-z,inf) > NORM69NORM = norm(w(i,m-1)-z,inf);70end71w(i,m-1)=z;72end7374% Aqu� se actualiza el �ltimo valor de la �ltima columna de w.75z = (-(h^2)*f(x(n-1),y(m-1))+fung(b,y(m-1))+lambda*fung(x(n-1),d)+w(n-2,m-1)+lambda*w(n-1,m-2))/mu;76if norm(w(n-1,m-1)-z,inf) > NORM77NORM = norm(w(n-1,m-1)-z,inf);78end79w(n-1,m-1) = z;8081% El siguiente for actualiza los valores de w, desde el segundo82% al pen�ltimo, de la todas las fila.83for j = 2:m-28485% Aqu� actualizan los valores de w, desde el segundo al pen�ltimo,86% de la primera fila.87z = (-(h^2)*f(x(1),y(j))+fung(a,y(j))+lambda*w(1,j+1)+lambda*w(1,j-1)+w(2,j))/mu;88if norm(w(1,j)-z,inf) > NORM89NORM = norm(w(1,j)-z,inf);90end91w(1,j) = z;9293% Este for actualiza los valores de w, desde el segundo al pen�ltimo,94% de las filas 2 a la n-2.95for i = 2:n-296z = (-(h^2)*f(x(i),y(j))+w(i-1,j)+lambda*w(i,j+1)+w(i+1,j)+lambda*w(i,j-1))/mu;97if norm(w(i,j)-z,inf) > NORM98NORM = norm(w(i,j)-z,inf);99end100w(i,j) =z;101end102103% Aqu� se actualizan los valores de w, desde el segundo al pen�ltimo,104% de la �ltima fila105z = (-(h^2)*f(x(n-1),y(j))+fung(b,y(j))+w(n-2,j)+lambda*w(n-1,j+1)+lambda*w(n-1,j-1))/mu;106if norm(w(n-1,j)-z,inf) > NORM107NORM = norm(w(n-1,j)-z,inf);108end109w(n-1,j) =z;110end111112% Aqu� se actualiza el primer valor de la primera fila y primera columna de w. fila 1 y columna 1 de w.113z = (-(h^2)*f(x(1),y(1))+fung(a,y(1))+lambda*fung(x(1),c)+lambda*w(1,2)+w(2,1))/mu;114if norm(w(1,1)-z,inf) > NORM115NORM = norm(w(1,1)-z,inf);116end117w(1,1) =z;118119% El siguiente for actualiza los valores de w, desde el segundo al pen�ltimo120% de la primera columna.121for i = 2:n-2122z = (-(h^2)*f(x(i),y(1))+lambda*fung(x(i),c)+w(i-1,1)+lambda*w(i,2)+w(i+1,1))/mu;123if norm(w(i,1)-z,inf) > NORM124NORM = norm(w(i,1)-z,inf);125end126w(i,1) =z;127end128129% Aqu� se actualiza el �ltimo valor de la ultima columna de la matriz w.130z = (-(h^2)*f(x(n-1),y(1))+fung(b,y(1))+lambda*fung(x(n-1),c)+w(n-2,1)+lambda*w(n-1,2))/mu;131if norm(w(n-1,1)-z,inf) > NORM132NORM = norm(w(n-1,1)-z,inf);133end134w(n-1,1) =z;135136% Este if verifica que se ha alcanzado la tolerancia deseada.137if NORM <= tol138break139end140141% Aqu� se actualiza el contador del bucle while.142l = l+1;143end144145% Este �ltimo if env�a un mensaje de error si se ha revasado el n�mero146% de iteraciones deseado.147if l > N148error('El numero maximo de iteraciones fue excedido')149end150I=zeros(m-1,1);151J=zeros(m-1,1);152%imprime los t�tulos de cada columna de la tabla resultante153fprintf(' i j xi xj w(i,j) \n')154for j=1:m-1 %se llena los valores para los �ndices j en la tabla resultante155J(j)=j;156end157for i=1:n-1 %se llenan los vectores con las coordenadas de x e y158x1=zeros(m-1,1);159for j=1:m-1160I(j)=i; %se llena los valores para los indices i en la tabla resultante161x1(j)=x(i);162end163A= [ I J x1 transpose(y) transpose(w(i,:)) ]; %guarda la matriz con164%las aproximaciones para cada i165%imprime los elemento de la tabla166fprintf([repmat('%10.5f\t', 1, size(A, 2)) '\n'], A');167end168% muestra el n�mero de iteraciones169fprintf('\n Y el numero de iteraciones necesarias fue: %d ', l)170end % termina el codigo171172% Ejecuci�n del ejemplo 2 , de la secci�n 12.1 del libro de Burden173% f = @(x,y)x*exp(y);174% a=0; b=2; c=0; d=1; m=5; n=6; tol=1e-10; N=100;175% PEFD(f,a,b,c,d,m,n,tol,N)176177178179