1

2
% Nombres: -Rafael An�bal Calder�n Melara  CM11099
3
%          -Estely Dolores M�ndez Moreno   MM09204
4

5
% Fecha de creaci�n: 15 de Noviembre de 2016
6
% Esta funci�n aproxima la soluci�n de la ecuaci�n de Poisson
7
%      f_xx(x,y) + f_yy(x,y) = f(x,y); con x en [a,b], y en [c,d]
8
% sujeto a las condiciones iniciales
9
%      u(x,y) = g(x,y) si x=a o x=b & y pertenece a [c,d]
10
%      u(x,y) = g(x,y) si y=c o y=d & x pertenece a [a,b]
11
% El algoritmo utiliza el m�todo de diferencias finitas 
12
% y a su vez el m�todo de Gauss Seidel.
13
% Notar que se utiliza la norma infinita en nuestro c�digo, 
14
% pero se podr�a utilizar cualquier otra norma.
15
% INPUT:  -- f: Funciones
16
%        -- a,b,c,d: Valores frontera.
17
%        -- n>=3, m>=3: Enteros que sirven para determinar las 
18
%                   dimensiones de la cuadr�cula.
19
%        -- tol: Tolerancia  requerida.
20
%        -- N: M�ximo n�mero de iteraciones.
21

22
% OUTPUT: -- Tabla con las aproximaciones para u(x,y) y los puntos 
23
%            x & y.
24
%         -- l: n�mero de iteraciones requeridas para obtener 
25
%              la soluci�n.
26

27
function [] = PEFD(f,a,b,c,d,m,n,tol,N)
28

29
% h y k representan el tama�o de paso para la construcci�n de los
30
% puntos x_i e y_j respectivamente.
31
h = (b-a)/n;
32
k = (d-c)/m;
33

34
% Se crean x & y que son vectores que contienen las coordenadas  
35
% para los distintos puntos de la cuadr�cula que utiliza 
36
% el m�todo.
37
x = zeros(1,n-1);
38
for i = 1:n-1
39
    x(i) = a+i*h;
40
end
41
y = zeros(1,m-1);
42
for j = 1:m-1
43
    y(j) = c+j*k;
44
end
45

46
% w es la matriz de ceros que luego se ir� actualizando con 
47
% las aproximaciones w_ji de la soluci�n u(x_i,y_j).
48
w = zeros(n-1,m-1);
49

50
% Se calculan lambda y mu para hacer m�s eficiente el c�digo en t�rmios de
51
% las operaciones que se realizan. 
52
lambda = (h^2)/(k^2);
53
mu = 2*(1+lambda);
54

55
% El bucle while se utiliza para realizar iteraciones del m�todo
56
% de Gauss-Seidel Optimizado
57
l = 1;
58
while l <= N 
59
    
60
    % Aqu� se actualiza el elemento 1 de la �ltima columna de w.
61
    z = (-(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;
62
    NORM = norm(z-w(1,m-1),inf);
63
    w(1,m-1)=z;
64
    
65
    % El siguiente for actualiza los valores de w, desde el segundo
66
    % al pen�ltimo elemento de w, de la �ltima columna.
67
    for i = 2:n-2
68
        z = (-(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;
69
        if norm(w(i,m-1)-z,inf) > NORM
70
            NORM = norm(w(i,m-1)-z,inf);
71
        end
72
        w(i,m-1)=z;
73
    end
74
    
75
    % Aqu� se actualiza el �ltimo valor de la �ltima columna de w.
76
    z = (-(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;
77
    if norm(w(n-1,m-1)-z,inf) > NORM
78
        NORM = norm(w(n-1,m-1)-z,inf);
79
    end
80
    w(n-1,m-1) = z;
81
    
82
    % El siguiente for actualiza los valores de w, desde el segundo 
83
    % al pen�ltimo, de la todas las fila.
84
    for j = 2:m-2
85
        
86
        % Aqu� actualizan los valores de w, desde el segundo al pen�ltimo, 
87
        % de la primera fila.
88
        z = (-(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;
89
        if norm(w(1,j)-z,inf) > NORM
90
            NORM = norm(w(1,j)-z,inf);
91
        end
92
        w(1,j) = z;
93
        
94
        % Este for actualiza los valores de w, desde el segundo al pen�ltimo, 
95
        % de las filas 2 a la n-2.
96
        for i = 2:n-2
97
            z = (-(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;
98
            if norm(w(i,j)-z,inf) > NORM
99
                NORM = norm(w(i,j)-z,inf);
100
            end
101
            w(i,j) =z;
102
        end
103
        
104
        % Aqu� se actualizan los valores de w, desde el segundo al pen�ltimo,
105
        % de la �ltima fila
106
        z = (-(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;
107
        if norm(w(n-1,j)-z,inf) > NORM
108
            NORM = norm(w(n-1,j)-z,inf);
109
        end
110
        w(n-1,j) =z;
111
    end
112
    
113
    % Aqu� se actualiza el primer valor de la primera fila y primera columna de w. fila 1 y columna 1 de w.
114
    z = (-(h^2)*f(x(1),y(1))+fung(a,y(1))+lambda*fung(x(1),c)+lambda*w(1,2)+w(2,1))/mu;
115
    if norm(w(1,1)-z,inf) > NORM
116
        NORM = norm(w(1,1)-z,inf);
117
    end
118
    w(1,1) =z;
119
    
120
    % El siguiente for actualiza los valores de w, desde el segundo al pen�ltimo
121
    % de la primera columna.
122
    for i = 2:n-2
123
        z = (-(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;
124
        if norm(w(i,1)-z,inf) > NORM
125
            NORM = norm(w(i,1)-z,inf);
126
        end    
127
        w(i,1) =z;
128
    end
129
    
130
    % Aqu� se actualiza el �ltimo valor de la ultima columna de la matriz w.
131
    z = (-(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;
132
    if norm(w(n-1,1)-z,inf) > NORM
133
        NORM = norm(w(n-1,1)-z,inf);
134
    end
135
    w(n-1,1) =z;
136
    
137
    % Este if verifica que se ha alcanzado la tolerancia deseada.
138
    if NORM <= tol
139
        break
140
    end
141
    
142
    % Aqu� se actualiza el contador del bucle while.
143
    l = l+1;
144
end
145

146
% Este �ltimo if env�a un mensaje de error si se ha revasado el n�mero
147
% de iteraciones deseado.
148
if l > N
149
error('El numero maximo de iteraciones fue excedido') 
150
end
151
I=zeros(m-1,1);
152
J=zeros(m-1,1);
153
%imprime los t�tulos de cada columna de la tabla resultante
154
fprintf('      i              j              xi              xj              w(i,j) \n')
155
for j=1:m-1 %se llena los valores para los �ndices j en la tabla resultante
156
    J(j)=j;
157
end
158
for i=1:n-1 %se llenan los vectores con las coordenadas de x e y
159
    x1=zeros(m-1,1);
160
    for j=1:m-1
161
        I(j)=i; %se llena los valores para los indices i en la tabla resultante
162
        x1(j)=x(i);
163
    end
164
    A= [ I  J  x1 transpose(y) transpose(w(i,:)) ]; %guarda la matriz con 
165
    %las aproximaciones para cada i
166
    %imprime los elemento de la tabla
167
    fprintf([repmat('%10.5f\t', 1, size(A, 2)) '\n'], A');
168
end
169
% muestra el n�mero de iteraciones
170
fprintf('\n Y el numero de iteraciones necesarias fue: %d ', l) 
171
end % termina el codigo
172

173
% Ejecuci�n del ejemplo 2 , de la secci�n 12.1 del libro de Burden
174
% f = @(x,y)x*exp(y);
175
% a=0; b=2; c=0; d=1; m=5; n=6; tol=1e-10; N=100; 
176
% PEFD(f,a,b,c,d,m,n,tol,N)
177

178

179