SharedProyecto_final_analisis_numerico2 / PEFD.mOpen in CoCalc
Author: Rafael Calderón
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