Diferencia entre revisiones de «Contaminación del lago (AMLS)»
De MateWiki
| Línea 5: | Línea 5: | ||
*Sara Berti.}} | *Sara Berti.}} | ||
{|style="margin: 0 auto;" | {|style="margin: 0 auto;" | ||
| − | | [[Archivo: | + | | [[Archivo:Lago1.jpg|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre Contaminación en el lago en formato póster realizado por el grupo AMLS.]] |
|} | |} | ||
{|style="margin: 0 auto;" | {|style="margin: 0 auto;" | ||
| − | | [[Archivo: | + | | [[Archivo:Lago2.jpg|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre Contaminación en el lago en formato póster realizado por el grupo AMLS.]] |
|} | |} | ||
{|style="margin: 0 auto;" | {|style="margin: 0 auto;" | ||
| − | | [[Archivo: | + | | [[Archivo:Lago3.jpg|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre Contaminación en el lago en formato póster realizado por el grupo AMLS.]] |
|} | |} | ||
| − | + | ||
| − | + | ||
| − | + | ||
Revisión actual del 14:10 18 nov 2025
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Contaminación del lago (AMLS) |
| Asignatura | MNEDP |
| Curso | 2025-26 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Códigos de Matlab
close all
clear all
%-------------------------------MALLA 1------------------------------------
%%%----------PASO 1: VARIABLES-----------
% Término fuente f(x,y)
% --- Fuente física: vertido puntual normalizado (gaussiana) ---
Q = 1.0; % Caudal total de contaminante
x0 = 0.5; y0 = 0.5; % Punto de vertido del contaminante
sigma = 0.05; % Anchura de la fuente
f = @(x,y) (Q/(2*pi*sigma^2)) .* exp(-((x-x0).^2 + (y-y0).^2) / (2*sigma^2));
% Parámetros físicos
D = 0.01; % Coeficiente de difusión
k = 0.1; % Coeficiente de reacción
%%%----------PASO 2: CREACIÓN DE LAS MATRICES DE NODOS Y TRIANGULOS-----------
Vertices = load('vertices1.txt');
Triangulos = load('triangulos1.txt'); %Cargamos los ficheros
Ne = length(Triangulos); % número de Triangulos
M = length(Vertices); % número de nodos (vertices)
% Coordenadas de los nodos
coord1 = zeros(M,2);
coord1(:,1) = Vertices(:,2);
coord1(:,2) = Vertices(:,3);
% Conectividad (qué nodos forman cada triángulo)
nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos1=nodos;
%%%----------PASO 3: MATRICES LOCALES DE REFERENCIA-----------
% Matriz de masa local (para el término de reacción)
L0 = [1/12 1/24 1/24;
1/24 1/12 1/24;
1/24 1/24 1/12];
% Matrices para los términos de derivadas parciales
Lxx = [1/2 -1/2 0;
-1/2 1/2 0;
0 0 0];
Lyy = [1/2 0 -1/2;
0 0 0;
-1/2 0 1/2];
Lxy = [1/2 0 -1/2;
-1/2 0 1/2;
0 0 0];
% Vector de integración de referencia (para f)
l0 = [1/6; 1/6; 1/6];
%%%----------PASO 4: ENSAMBLADO DE LA MATRIZ GLOBAL Y EL VECTOR DE CARGAS-----------
A = zeros(M); % matriz global (rigidez + reacción)
b = zeros(M,1); % vector global de cargas
for e = 1:Ne
% Nodos del triangulo e
nodosK = nodos(e,:);
% Matriz del elemento (de coordenadas locales a reales)
Bk = [coord1(nodosK(2),1)-coord1(nodosK(1),1), coord1(nodosK(3),1)-coord1(nodosK(1),1);
coord1(nodosK(2),2)-coord1(nodosK(1),2), coord1(nodosK(3),2)-coord1(nodosK(1),2)];
detBk = det(Bk); % área*2 del triángulo
Ck = inv(Bk) * (inv(Bk))'; % matriz Ck = (Bk^-1)(Bk^-T)
% Matriz elemental (difusión + reacción)
Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
% Ensamblado en la matriz global
A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
A1=A;
% Cálculo del término fuente medio en el triángulo
fk = mean([ ...
f(coord1(nodosK(1),1), coord1(nodosK(1),2)), ...
f(coord1(nodosK(2),1), coord1(nodosK(2),2)), ...
f(coord1(nodosK(3),1), coord1(nodosK(3),2)) ]);
% Vector local de cargas
lk = fk * abs(detBk) * l0;
% Ensamblado en el vector global
b(nodosK) = b(nodosK) + lk;
end
%%%----------PASO 5: APLICAR CONDICIONES DE CONTORNO DIRICHLET HOMOGÉNEAS (c=0)-----------
% Detectamos los nodos de frontera:
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));
% Impongo c=0 en los nodos de la frontera
for n = boundary_nodes'
A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
b(n) = 0;
end
%%%----------PASO 6: RESOLUCIÓN DEL SISTEMA-----------
u1 = A\b; % solución aproximada en los vertices
%%%----------PASO 7: GRAFICAMOS-----------
figure(1)
trisurf(nodos, coord1(:,1), coord1(:,2), u1)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
%--------------------------------------MALLA 2--------------------------------------------------------
Vertices = load('vertices2.txt');
Triangulos = load('triangulos2.txt');
Ne = length(Triangulos);
M = length(Vertices);
coord2 = zeros(M,2);
coord2(:,1) = Vertices(:,2);
coord2(:,2) = Vertices(:,3);
nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos2=nodos;
L0 = [1/12 1/24 1/24;
1/24 1/12 1/24;
1/24 1/24 1/12];
Lxx = [1/2 -1/2 0;
-1/2 1/2 0;
0 0 0];
Lyy = [1/2 0 -1/2;
0 0 0;
-1/2 0 1/2];
Lxy = [1/2 0 -1/2;
-1/2 0 1/2;
0 0 0];
l0 = [1/6; 1/6; 1/6];
A = zeros(M);
b = zeros(M,1);
for e = 1:Ne
nodosK = nodos(e,:);
Bk = [coord2(nodosK(2),1)-coord2(nodosK(1),1), coord2(nodosK(3),1)-coord2(nodosK(1),1);
coord2(nodosK(2),2)-coord2(nodosK(1),2), coord2(nodosK(3),2)-coord2(nodosK(1),2)];
detBk = det(Bk);
Ck = inv(Bk) * (inv(Bk))';
Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
A2=A;
fk = mean([ ...
f(coord2(nodosK(1),1), coord2(nodosK(1),2)), ...
f(coord2(nodosK(2),1), coord2(nodosK(2),2)), ...
f(coord2(nodosK(3),1), coord2(nodosK(3),2)) ]);
lk = fk * abs(detBk) * l0;
b(nodosK) = b(nodosK) + lk;
end
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));
for n = boundary_nodes'
A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
b(n) = 0;
end
u2 = A\b;
figure(2)
trisurf(nodos, coord2(:,1), coord2(:,2), u2)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
%---------------------------------------MALA 3-----------------------------------------------------------
Vertices = load('vertices3.txt');
Triangulos = load('triangulos3.txt');
Ne = length(Triangulos);
M = length(Vertices);
coord3 = zeros(M,2);
coord3(:,1) = Vertices(:,2);
coord3(:,2) = Vertices(:,3);
nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos3=nodos;
L0 = [1/12 1/24 1/24;
1/24 1/12 1/24;
1/24 1/24 1/12];
Lxx = [1/2 -1/2 0;
-1/2 1/2 0;
0 0 0];
Lyy = [1/2 0 -1/2;
0 0 0;
-1/2 0 1/2];
Lxy = [1/2 0 -1/2;
-1/2 0 1/2;
0 0 0];
l0 = [1/6; 1/6; 1/6];
A = zeros(M);
b = zeros(M,1);
for e = 1:Ne
nodosK = nodos(e,:);
Bk = [coord3(nodosK(2),1)-coord3(nodosK(1),1), coord3(nodosK(3),1)-coord3(nodosK(1),1);
coord3(nodosK(2),2)-coord3(nodosK(1),2), coord3(nodosK(3),2)-coord3(nodosK(1),2)];
detBk = det(Bk);
Ck = inv(Bk) * (inv(Bk))';
Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
A3=A;
fk = mean([ ...
f(coord3(nodosK(1),1), coord3(nodosK(1),2)), ...
f(coord3(nodosK(2),1), coord3(nodosK(2),2)), ...
f(coord3(nodosK(3),1), coord3(nodosK(3),2)) ]);
lk = fk * abs(detBk) * l0;
b(nodosK) = b(nodosK) + lk;
end
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));
for n = boundary_nodes'
A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
b(n) = 0;
end
u3 = A\b;
figure(3)
trisurf(nodos, coord3(:,1), coord3(:,2), u3)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
%---------------------------------------MALLA 4-------------------------------------------------------------------
Vertices = load('vertices4.txt');
Triangulos = load('triangulos4.txt');
Ne = length(Triangulos);
M = length(Vertices);
coord4 = zeros(M,2);
coord4(:,1) = Vertices(:,2);
coord4(:,2) = Vertices(:,3);
nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos4=nodos;
L0 = [1/12 1/24 1/24;
1/24 1/12 1/24;
1/24 1/24 1/12];
Lxx = [1/2 -1/2 0;
-1/2 1/2 0;
0 0 0];
Lyy = [1/2 0 -1/2;
0 0 0;
-1/2 0 1/2];
Lxy = [1/2 0 -1/2;
-1/2 0 1/2;
0 0 0];
l0 = [1/6; 1/6; 1/6];
A = zeros(M);
b = zeros(M,1);
for e = 1:Ne
nodosK = nodos(e,:);
Bk = [coord4(nodosK(2),1)-coord4(nodosK(1),1), coord4(nodosK(3),1)-coord4(nodosK(1),1);
coord4(nodosK(2),2)-coord4(nodosK(1),2), coord4(nodosK(3),2)-coord4(nodosK(1),2)];
detBk = det(Bk);
Ck = inv(Bk) * (inv(Bk))';
Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
A4=A;
fk = mean([ ...
f(coord4(nodosK(1),1), coord4(nodosK(1),2)), ...
f(coord4(nodosK(2),1), coord4(nodosK(2),2)), ...
f(coord4(nodosK(3),1), coord4(nodosK(3),2)) ]);
lk = fk * abs(detBk) * l0;
b(nodosK) = b(nodosK) + lk;
end
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));
for n = boundary_nodes'
A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
b(n) = 0;
end
u4 = A\b;
figure(4)
trisurf(nodos, coord4(:,1), coord4(:,2), u4)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
