%% Histograma Adaptativo según Desviación Típica en un punto x

clear; clc; close all;

% 1. Parámetros de Control
L = 1;              
N_harmonicos = 500;  
M_iter = 3000; 
x_fijo = 0.5;       
omega = 2 * pi / L;
p=2;


% 2. Definición de la desviación típica por modo
k = 1:N_harmonicos;
sigma_k = 1 ./ (k.^(p/2)); 

% 3. Simulación de Monte Carlo
valores_f = zeros(1, M_iter);
for m = 1:M_iter
    ak = randn(1, N_harmonicos) .* sigma_k;
    bk = randn(1, N_harmonicos) .* sigma_k;
    
    % Evaluamos la serie en x_fijo
    f_x = sum( ak .* cos(k * omega * x_fijo) + bk .* sin(k * omega * x_fijo) );
    valores_f(m) = f_x;
end

% 4. Cálculo de la Varianza Teórica 
% La varianza de la suma es la suma de las varianzas de los sumandos

var_teorica_adaptada = sum(sigma_k.^2); 

% 5. Representación Gráfica
figure('Color', 'w');
h = histogram(valores_f, 50, 'Normalization', 'pdf', 'FaceColor', [0.2, 0.7, 0.5], 'EdgeColor', 'none');
hold on;

% Generar Gaussiana basada en la varianza calculada
x_axis = linspace(min(valores_f), max(valores_f), 200);
y_gauss = (1/sqrt(2*pi*var_teorica_adaptada)) * exp(-x_axis.^2 / (2*var_teorica_adaptada));

plot(x_axis, y_gauss, 'r', 'LineWidth', 2.5);

title(['Distribución Adaptada: \sigma_{total}^2 = ', num2str(var_teorica_adaptada, '%.4f')]);
subtitle(['Parámetros: p = ', num2str(p)]);
xlabel(['Valor de f(x) en x = ', num2str(x_fijo)]);
ylabel('Densidad de Probabilidad');
legend('Simulación', 'Gaussiana Teórica Adaptada');
grid on;


fprintf('Varianza Teórica Calculada: %f\n', var_teorica_adaptada);
fprintf('Varianza Observada: %f\n', var(valores_f));

%% Simulación de Función Aleatoria con coeficientes normales.

% 1. Parámetros
L = 1;              % Longitud del dominio
N = 500;            % Número de armónicos
dx = 0.001;
x = 0:dx:L;
f = zeros(size(x));
p=8;

% 2. Generación de Coeficientes
% Introducimos la desviación típica (raiz de la varianza) 
k = 1:N;
sigma_k = 1./k.^(p/2); 

% Coeficientes: Normal(0, 1) multiplicada por su desviación estándar
ak = randn(1, N) .* sigma_k;
bk = randn(1, N) .* sigma_k;

% 3. Construcción de la Serie
omega = 2 * pi / L;
for i = 1:N
    f = f + ak(i)*cos(i*omega*x) + bk(i)*sin(i*omega*x);
end

% 4. Representación Gráfica
figure('Color', 'w');
plot(x, f, 'LineWidth', 1.5, 'Color', [0.8, 0.2, 0.2]);
title(['Función Aleatoria (Desviación \sigma_k^2 = 1/k^', num2str(p), ')']);
xlabel('x');
ylabel('f_{\sigma}(x)');
grid on;

% Solución de la Ecuación del Calor con Condiciones Iniciales Aleatorias
% Caso: Dominio Periódico (Anillo) con Coeficientes de Fourier Estocásticos

clear; clc; close all;

%% 1. Parámetros del Problema

L = 2*pi;              % Longitud del dominio (anillo)
alpha = 0.5;           % Difusividad térmica
T_final = 3.0;         % Tiempo total de simulación
N_x = 500;             % Puntos en el espacio
x = linspace(0, L, N_x);
k_max = 50;            % Número de modos de Fourier a sumar
p = 1;                 % Exponente de decaimiento de varianza (1/k^p)

%% 2. Generación de la Condición Inicial Aleatoria (t = 0)
% Coeficientes de Fourier aleatorios
A0 = randn(1); % Media aleatoria
Ak = randn(1, k_max) ./ ( (1:k_max).^(p/2) ); % Varianza 1/k^p
Bk = randn(1, k_max) ./ ( (1:k_max).^(p/2) );

% Construcción de la función inicial f(x)
u_initial = A0/2;
for k = 1:k_max
    u_initial = u_initial + Ak(k)*cos(2*pi*k*x/L) + Bk(k)*sin(2*pi*k*x/L);
end

%% 3. Evolución Temporal y Graficación
t_steps = [0, 0.05, 0.2, 0.5, 1.5,3]; 
figure('Color', 'w', 'Position', [100, 100, 1000, 600]);
hold on; grid on;

colors = jet(length(t_steps)); 

for i = 1:length(t_steps)
    t = t_steps(i);
    
    % Calculamos la solución u(x,t) usando la serie de Fourier
    u_t = A0/2;
    for k = 1:k_max
        decay = exp(-alpha * (2*pi*k/L)^2 * t); 
        u_t = u_t + decay * (Ak(k)*cos(2*pi*k*x/L) + Bk(k)*sin(2*pi*k*x/L));
    end
    
    plot(x, u_t, 'Color', colors(i,:), 'LineWidth', 2, ...
        'DisplayName', sprintf('t = %.2f', t));
end

title(['Evolución de u(x,t) con C.I. Aleatoria (p = ', num2str(p), ')']);
xlabel('Posición x');
ylabel('Temperatura u');
legend show;
set(gca, 'FontSize', 12);