Diferencia entre revisiones de «Ecuación del calor MMA»

De MateWiki
Saltar a: navegación, buscar
 
(No se muestra una edición intermedia de otro usuario)
Línea 7: Línea 7:
 
[[Archivo: foto MMA.jpeg||800px]]
 
[[Archivo: foto MMA.jpeg||800px]]
  
−
=== Serie de Fourier Función Exponencial ===
+
=== Dispersión de partículas ===
  
  
−
{{matlab|codigo=
+
<syntaxhighlight lang="python">
  
−
clear; clc; close all;
+
import numpy as np
 +
import matplotlib.pyplot as plt
  
−
% Dominio
+
# Parámetros
−
x = linspace(-pi, pi, 2000);
+
np.random.seed(42)
 +
n_particles = 400
 +
times = [0.05, 0.4, 1.0]
  
−
% Función original
+
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
−
f = exp(x);
+
  
−
% Número máximo de modos
+
for ax, t in zip(axes, times):
−
Nmax = 100;
+
    # Movimiento aleatorio (normal)
 +
    x = np.random.normal(0, np.sqrt(t), n_particles)
 +
    y = np.random.normal(0, np.sqrt(t), n_particles)
  
−
% ---- Coeficientes trigonométricos ----
+
    ax.scatter(x, y, s=10, alpha=0.7)
−
a0 = (1/pi) * trapz(x, f);
+
    ax.set_title(f"t = {t}")
 +
    ax.set_aspect("equal")
  
−
an = zeros(1, Nmax);
+
    # Límites visuales
−
bn = zeros(1, Nmax);
+
    ax.set_xlim(-4, 4)
 +
    ax.set_ylim(-4, 4)
 +
    ax.grid(True)
  
−
for n = 1:Nmax
+
fig.suptitle("Dispersión de partículas")
−
    an(n) = (1/pi) * trapz(x, f .* cos(n*x));
+
plt.tight_layout()
−
    bn(n) = (1/pi) * trapz(x, f .* sin(n*x));
+
plt.show()
−
end
+
  
−
% Valores de N que queremos visualizar
+
</syntaxhighlight>
−
N_values = [5 10 15 50];
+
  
−
figure; hold on;
 
  
−
% Dibujar función original
 
−
plot(x, f, 'k', 'LineWidth', 1);
 
  
−
% Colores automáticos
 
−
colors = lines(length(N_values));
 
  
−
for j = 1:length(N_values)
 
−
    N = N_values(j);
 
  
−
    % Suma parcial trigonométrica
 
−
    SN = (a0/2) * ones(size(x));
 
−
    for n = 1:N
 
−
        SN = SN + an(n)*cos(n*x) + bn(n)*sin(n*x);
 
−
    end
 
  
−
    plot(x, SN, 'Color', colors(j,:), 'LineWidth', 0.8);
+
=== Trayectorias de movimiento browniano ===
−
end
+
  
−
legend('e^x', 'S_5', 'S_{10}', 'S_{15}', 'S_{50}', 'Location', 'best');
+
<syntaxhighlight lang="python">
−
title('Serie de Fourier de e^x');
+
import numpy as np
−
xlabel('x');
+
import matplotlib.pyplot as plt
−
ylabel('Valor');
+
−
grid on;
+
−
}}
+
  
 +
# Parámetros
 +
np.random.seed(42)
 +
T = 1
 +
n_steps = 300
 +
dt = T / n_steps
  
 +
t = np.linspace(0, T, n_steps + 1)
  
 +
plt.figure(figsize=(8, 5))
  
 +
# Varias trayectorias
 +
for _ in range(12):
 +
    increments = np.random.normal(0, np.sqrt(dt), n_steps)
 +
    B = np.concatenate([[0], np.cumsum(increments)])
 +
    plt.plot(t, B)
  
 +
plt.title("Trayectorias de movimiento browniano")
 +
plt.xlabel("tiempo")
 +
plt.ylabel("posición")
 +
plt.grid(True)
  
−
{{matlab|codigo=
+
plt.show()
  
−
clear; clc; close all;
+
</syntaxhighlight>
  
−
% Intervalo
 
−
a = -pi;
 
−
b = pi;
 
  
−
% Mallado fino para integración
 
−
x = linspace(a,b,5000);
 
−
dx = x(2)-x(1);
 
  
−
% Función
 
−
f = sign(t);
 
  
−
% Número máximo de términos
 
−
Nmax = 100;
 
  
−
error_L2 = zeros(1,Nmax);
 
  
−
for N = 1:Nmax
 
−
   
 
−
    % Coeficientes
 
−
    a0 = (1/pi)*trapz(x,f);
 
−
   
 
−
    an = zeros(1,N);
 
−
    bn = zeros(1,N);
 
−
   
 
−
    for n = 1:N
 
−
        an(n) = (1/pi)*trapz(x, f.*cos(n*x));
 
−
        bn(n) = (1/pi)*trapz(x, f.*sin(n*x));
 
−
    end
 
−
   
 
−
    % Suma parcial
 
−
    SN = a0/2*ones(size(x));
 
−
   
 
−
    for n = 1:N
 
−
        SN = SN + an(n)*cos(n*x) + bn(n)*sin(n*x);
 
−
    end
 
−
   
 
−
    % Error L2
 
−
    error_L2(N) = sqrt(trapz(x, (f-SN).^2));
 
−
   
 
−
end
 
  
−
% Gráfica del error
 
−
figure;
 
−
plot(1:Nmax, error_L2,'LineWidth',2)
 
−
ylim([0 max(error_L2)])
 
−
xlabel('N'); ylabel('Error L^2');
 
−
title('Error L^2 (e^x)'); grid on
 
  
−
fprintf('Error L2 en N=40: %.16e\n', error_L2(100));
+
=== Evolución de la densidad ===
−
}}
+
  
 +
<syntaxhighlight lang="python">
  
 +
import numpy as np
 +
import matplotlib.pyplot as plt
  
 +
# Dominio
 +
x = np.linspace(-6, 6, 1000)
  
 +
# Distintos tiempos
 +
times = [0.2, 0.8, 1.6, 5,50]
  
−
=== Serie de Fourier Función Signo ===
+
plt.figure(figsize=(8, 5))
  
 +
for t in times:
 +
    # Solución de la ecuación del calor (gaussiana)
 +
    p = 1 / np.sqrt(2 * np.pi * t) * np.exp(-x**2 / (2 * t))
 +
    plt.plot(x, p, label=f"t = {t}")
  
−
{{matlab|codigo=
+
plt.title("Evolución de la densidad")
 +
plt.xlabel("x")
 +
plt.ylabel("densidad")
 +
plt.legend()
 +
plt.grid(True)
  
−
clear; clc; close all;
+
plt.show()
 +
</syntaxhighlight>
  
−
% Dominio
 
−
x = linspace(-pi, pi, 4000);
 
  
−
% Función discontinua: signo (escalón)
 
−
f = ones(size(x));
 
−
f(x < 0) = -1;
 
  
−
% Fourier: máximo modo que necesitamos
 
−
Nmax = 100;
 
−
nvals = -Nmax:Nmax;
 
−
c = zeros(size(nvals));
 
  
−
% Coeficientes c_n = (1/2pi) int f(x) e^{-inx} dx
 
−
for k = 1:length(nvals)
 
−
    n = nvals(k);
 
−
    integrando = f .* exp(-1i*n*x);
 
−
    c(k) = trapz(x, integrando) / (2*pi);
 
−
end
 
  
−
% N a visualizar
 
−
N_values = [5 10 15 50 100];
 
  
−
figure; hold on;
 
  
−
% Original
 
−
plot(x, f, 'k', 'LineWidth', 1);
 
  
−
% Colores automáticos
+
=== Evolución de la densidad con ruido blanco ===
−
colors = lines(length(N_values));
+
  
−
for j = 1:length(N_values)
+
<syntaxhighlight lang="python">
−
    N = N_values(j);
+
import numpy as np
−
    SN = zeros(size(x));
+
import matplotlib.pyplot as plt
−
   
+
−
    for k = 1:length(nvals)
+
−
        if abs(nvals(k)) <= N
+
−
            SN = SN + c(k)*exp(1i*nvals(k)*x);
+
−
        end
+
−
    end
+
−
   
+
−
    plot(x, real(SN), 'Color', colors(j,:), 'LineWidth', 0.8);
+
−
end
+
  
−
legend('sign(x)', 'S_5', 'S_{10}', 'S_{15}', 'S_{50}','S_{100}', 'Location', 'best');
+
# Dominio espacial
−
title('Serie de Fourier de sign(x)');
+
x = np.linspace(-6, 6, 500)
−
xlabel('x'); ylabel('Valor');
+
−
grid on;
+
−
}}
+
  
 +
# Parámetro de difusión
 +
D = 1
  
 +
# Tiempos
 +
times = [0.2, 0.8, 1.6]
  
 +
# Solución del calor
 +
def u(x, t, D=1):
 +
    return (1 / np.sqrt(4 * np.pi * D * t)) * np.exp(-x**2 / (4 * D * t))
  
 +
# Semilla
 +
np.random.seed(4)
  
 +
plt.figure(figsize=(8,5))
  
−
{{matlab|codigo=
+
for t in times:
−
 
+
    u_det = u(x, t, D)
−
clear; clc; close all;
+
−
 
+
−
% Intervalo
+
−
a = -pi;
+
−
b = pi;
+
−
 
+
−
% Mallado fino para integración
+
−
x = linspace(a,b,5000);
+
−
dx = x(2)-x(1);
+
−
 
+
−
% Función signo
+
−
f = sign(x);
+
−
 
+
−
% Número máximo de términos
+
−
Nmax = 100;
+
−
 
+
−
error_L2 = zeros(1,Nmax);
+
−
 
+
−
for N = 1:Nmax
+
 
      
 
      
−
     % Coeficiente a0
+
     # Ruido proporcional
−
     a0 = (1/pi)*trapz(x,f);
+
     noise = 0.02 * np.random.randn(len(x))
 +
    u_noise = u_det + noise
 
      
 
      
−
     an = zeros(1,N);
+
     plt.plot(x, u_noise, label=f"t = {t}")
−
    bn = zeros(1,N);
+
−
   
+
−
    for n = 1:N
+
−
        an(n) = (1/pi)*trapz(x, f.*cos(n*x));
+
−
        bn(n) = (1/pi)*trapz(x, f.*sin(n*x));
+
−
    end
+
−
   
+
−
    % Suma parcial
+
−
    SN = a0/2*ones(size(x));
+
−
   
+
−
    for n = 1:N
+
−
        SN = SN + an(n)*cos(n*x) + bn(n)*sin(n*x);
+
−
    end
+
−
   
+
−
    % Error L2
+
−
    error_L2(N) = sqrt(trapz(x, (f-SN).^2));
+
−
   
+
−
end
+
−
 
+
−
% Gráfica del error
+
−
figure;
+
−
plot(1:Nmax, error_L2,'LineWidth',2)
+
−
xlabel('N'); ylabel('Error L^2');
+
−
title('Error L^2');
+
−
grid on
+
−
 
+
−
fprintf('Error L2 en N=100: %.16e\n', error_L2(100));
+
−
}}
+
−
 
+
−
 
+
−
 
+
−
 
+
−
 
+
−
 
+
−
{{matlab|codigo=
+
−
 
+
−
clear; clc; close all;
+
−
 
+
−
x = linspace(-pi, pi, 5000);
+
−
 
+
−
% Onda cuadrada
+
−
f = ones(size(x));
+
−
f(x < 0) = -1;
+
−
 
+
−
Nmax = 100;  % más grande para que Gibbs sea muy claro
+
−
nvals = -Nmax:Nmax;
+
−
c = zeros(size(nvals));
+
−
 
+
−
% Coeficientes de Fourier
+
−
for k = 1:length(nvals)
+
−
    n = nvals(k);
+
−
    integrando = f .* exp(-1i*n*x);
+
−
    c(k) = trapz(x, integrando)/(2*pi);
+
−
end
+
−
 
+
−
N_values = [5 10 15 50 100];
+
−
 
+
−
figure; hold on;
+
−
 
+
−
plot(x, f, 'k', 'LineWidth', 2);
+
−
 
+
−
colors = lines(length(N_values));
+
−
 
+
−
for j = 1:length(N_values)
+
−
    N = N_values(j);
+
−
    SN = zeros(size(x));
+
−
   
+
−
    for k = 1:length(nvals)
+
−
        if abs(nvals(k)) <= N
+
−
            SN = SN + c(k)*exp(1i*nvals(k)*x);
+
−
        end
+
−
    end
+
−
   
+
−
    plot(x, real(SN), 'Color', colors(j,:), 'LineWidth', 1.5);
+
−
end
+
−
 
+
−
legend('sign(x)', 'S_5', 'S_{10}', 'S_{15}', 'S_{50}', 'S_{100}', 'Location', 'best');
+
−
title('Fenómeno de Gibbs');
+
−
xlabel('x'); ylabel('Valor');
+
−
grid on;
+
−
 
+
−
xlim([-1 1])
+
−
ylim([-1.5 1.5])
+
−
}}
+
−
 
+
−
 
+
−
 
+
−
=== Sumas de Cesàro y de Abel ===
+
−
 
+
−
{{matlab|codigo=
+
−
 
+
−
clear; clc; close all;
+
−
 
+
−
% Intervalo base [-pi, pi]
+
−
a = -pi; b = pi;
+
−
x = linspace(a,b,2000);
+
−
 
+
−
% Parámetros de la animación
+
−
N_list = 1:2:80;            % N creciente
+
−
r_list = linspace(0.60,0.98, length(N_list));  % r creciente (mismo nº de frames)
+
−
 
+
−
% Funciones a animar
+
−
cases = {
+
−
    @(t) sign(t), 'sign';
+
−
    @(t) exp(t),  'exp'
+
−
};
+
−
 
+
−
for c = 1:size(cases,1)
+
−
    f = cases{c,1};
+
−
    name = cases{c,2};
+
−
 
+
−
    outdir = fullfile(pwd, ['frames_' name]);
+
−
    if ~exist(outdir, 'dir'); mkdir(outdir); end
+
−
 
+
−
    % Precomputo coeficientes hasta Nmax para eficiencia
+
−
    Nmax = max(N_list);
+
−
    [a0, an, bn] = fourier_coef(f, a, b, Nmax);
+
−
 
+
−
    for k = 1:length(N_list)
+
−
        N = N_list(k);
+
−
        r = r_list(k);
+
−
 
+
−
        % Aproximaciones
+
−
        SN    = fourier_partial(a0, an, bn, x, N);
+
−
        CesN  = fourier_cesaro(a0, an, bn, x, N); % Fejér
+
−
        AbelN = fourier_abel(a0, an, bn, x, N, r);
+
−
 
+
−
        % Plot
+
−
        figure(1); clf;
+
−
        plot(x, f(x), 'LineWidth', 2); hold on;
+
−
        plot(x, SN,    'LineWidth', 1.2);
+
−
        plot(x, CesN,  'LineWidth', 1.2);
+
−
        plot(x, AbelN, 'LineWidth', 1.2);
+
−
        grid on;
+
−
 
+
−
        xlabel('x'); ylabel('y');
+
−
        title(sprintf('%s: N=%d, r=%.2f', name, N, r));
+
−
        legend('f(x)', 'S_N', 'Cesàro (Fejér)', 'Abel', 'Location','best');
+
−
 
+
−
        % Guardar frame
+
−
        fname = fullfile(outdir, sprintf('frame_%04d.png', k));
+
−
        exportgraphics(gcf, fname, 'Resolution', 200);
+
−
    end
+
−
end
+
−
 
+
−
disp('Frames generados.');
+
−
 
+
−
%% ================== FUNCIONES ==================
+
−
 
+
−
function [a0, an, bn] = fourier_coef(f, a, b, N)
+
−
% Coeficientes trigonométricos en [a,b], periodo T=b-a
+
−
    T = b - a;
+
−
    xx = linspace(a,b,12000);    % mallado fino para integrar
+
−
    fx = f(xx);
+
−
 
+
−
    a0 = (2/T) * trapz(xx, fx);
+
−
 
+
−
    an = zeros(1,N);
+
−
    bn = zeros(1,N);
+
−
 
+
−
    for n = 1:N
+
−
        an(n) = (2/T) * trapz(xx, fx .* cos(2*pi*n*xx/T));
+
−
        bn(n) = (2/T) * trapz(xx, fx .* sin(2*pi*n*xx/T));
+
−
    end
+
−
end
+
−
 
+
−
function SN = fourier_partial(a0, an, bn, x, N)
+
−
% S_N(x) usando primeros N coeficientes en [-pi,pi] => cos(nx), sin(nx)
+
−
    SN = (a0/2) * ones(size(x));
+
−
    for n = 1:N
+
−
        SN = SN + an(n)*cos(n*x) + bn(n)*sin(n*x);
+
−
    end
+
−
end
+
−
 
+
−
function sigmaN = fourier_cesaro(a0, an, bn, x, N)
+
−
% Cesàro/Fejér: sigma_N = (1/(N+1)) sum_{k=0}^N S_k
+
−
    sigmaN = zeros(size(x));
+
−
    Sk = (a0/2) * ones(size(x));  % S_0
+
−
    sigmaN = sigmaN + Sk;
+
−
 
+
−
    for k = 1:N
+
−
        Sk = Sk + an(k)*cos(k*x) + bn(k)*sin(k*x);
+
−
        sigmaN = sigmaN + Sk;
+
−
    end
+
−
 
+
−
    sigmaN = sigmaN/(N+1);
+
−
end
+
−
 
+
−
function Ar = fourier_abel(a0, an, bn, x, N, r)
+
−
% Abel truncado: A_r(x) = a0/2 + sum_{n=1}^N r^n(...)
+
−
    Ar = (a0/2) * ones(size(x));
+
−
    for n = 1:N
+
−
        Ar = Ar + (r^n) * ( an(n)*cos(n*x) + bn(n)*sin(n*x) );
+
−
    end
+
−
end
+
−
}}
+
−
 
+
−
 
+
−
=== Apéndice ===
+
−
 
+
−
[[Archivo:SumasFourierMMA.gif||]]
+
−
 
+
−
 
+
  
−
[[Categoría:EDP]]
+
plt.xlabel("x")
−
[[Categoría:EDP25/26]]
+
plt.ylabel("densidad")
 +
plt.title("Evolución de la densidad con ruido blanco")
 +
plt.legend()
 +
plt.grid(True)
  
 +
plt.show()
  
 +
</syntaxhighlight>
  
  
 
[[Categoría:EDP]]
 
[[Categoría:EDP]]
 
[[Categoría:EDP25/26]]
 
[[Categoría:EDP25/26]]

Revisión actual del 11:20 13 abr 2026

Trabajo realizado por estudiantes
Título Ecuación del calor. Grupo MMA
Asignatura EDP
Curso 2025-26
Autores Marta Tejedor

María Romojaro

Andrea Sánchez

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


Foto MMA.jpeg

1 Dispersión de partículas

import numpy as np
import matplotlib.pyplot as plt

# Parámetros
np.random.seed(42)
n_particles = 400
times = [0.05, 0.4, 1.0]

fig, axes = plt.subplots(1, 3, figsize=(12, 4))

for ax, t in zip(axes, times):
    # Movimiento aleatorio (normal)
    x = np.random.normal(0, np.sqrt(t), n_particles)
    y = np.random.normal(0, np.sqrt(t), n_particles)

    ax.scatter(x, y, s=10, alpha=0.7)
    ax.set_title(f"t = {t}")
    ax.set_aspect("equal")

    # Límites visuales
    ax.set_xlim(-4, 4)
    ax.set_ylim(-4, 4)
    ax.grid(True)

fig.suptitle("Dispersión de partículas")
plt.tight_layout()
plt.show()




2 Trayectorias de movimiento browniano

import numpy as np
import matplotlib.pyplot as plt

# Parámetros
np.random.seed(42)
T = 1
n_steps = 300
dt = T / n_steps

t = np.linspace(0, T, n_steps + 1)

plt.figure(figsize=(8, 5))

# Varias trayectorias
for _ in range(12):
    increments = np.random.normal(0, np.sqrt(dt), n_steps)
    B = np.concatenate([[0], np.cumsum(increments)])
    plt.plot(t, B)

plt.title("Trayectorias de movimiento browniano")
plt.xlabel("tiempo")
plt.ylabel("posición")
plt.grid(True)

plt.show()





3 Evolución de la densidad

import numpy as np
import matplotlib.pyplot as plt

# Dominio
x = np.linspace(-6, 6, 1000)

# Distintos tiempos
times = [0.2, 0.8, 1.6, 5,50]

plt.figure(figsize=(8, 5))

for t in times:
    # Solución de la ecuación del calor (gaussiana)
    p = 1 / np.sqrt(2 * np.pi * t) * np.exp(-x**2 / (2 * t))
    plt.plot(x, p, label=f"t = {t}")

plt.title("Evolución de la densidad")
plt.xlabel("x")
plt.ylabel("densidad")
plt.legend()
plt.grid(True)

plt.show()





4 Evolución de la densidad con ruido blanco

import numpy as np
import matplotlib.pyplot as plt

# Dominio espacial
x = np.linspace(-6, 6, 500)

# Parámetro de difusión
D = 1

# Tiempos 
times = [0.2, 0.8, 1.6]

# Solución del calor
def u(x, t, D=1):
    return (1 / np.sqrt(4 * np.pi * D * t)) * np.exp(-x**2 / (4 * D * t))

# Semilla 
np.random.seed(4)

plt.figure(figsize=(8,5))

for t in times:
    u_det = u(x, t, D)
    
    # Ruido proporcional 
    noise = 0.02 * np.random.randn(len(x))
    u_noise = u_det + noise
    
    plt.plot(x, u_noise, label=f"t = {t}")

plt.xlabel("x")
plt.ylabel("densidad")
plt.title("Evolución de la densidad con ruido blanco")
plt.legend()
plt.grid(True)

plt.show()
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Ecuación_del_calor_MMA&oldid=104639»