%% (1) Random PDE 1 - Random Source Term
clear all; clc;

% Set domain
L = 1; H = 1;

% Deterministic diffusion term
a = 'exp(0*x)';

% Random source term
f = 'exp(randn(1) * sin(pi*(x+y)) + randn(1) * cos(pi*(x+y)))';

% ... and its expectation
Ef = 'exp(0.5 * sin(pi*(x+y)).^2 + 0.5 * cos(pi*(x+y)).^2)';
Ef = 'exp(0.5) * exp(0*x)';
    
% Random boundary data 
s = '(x-x.^2) .* sin(rand(1) * (1+x+y))';
g{1} = s;     % Southern boundary data
g{2} = s;     % Eastern boundary data
g{3} = s;     % Northern boundary data
g{4} = s;     % Western boundary data

% ... and its expectation 
Es = '(x-x.^2) .* 1./(1+x+y) .* (1-cos(1+x+y))';
Eg{1} = Es;     % Southern boundary data
Eg{2} = Es;     % Eastern boundary data
Eg{3} = Es;     % Northern boundary data
Eg{4} = Es;     % Western boundary data

% Set number of elements and generate mesh
nx = 64; ny=64;
[pt, el] = fem2d_lin_mesh(L,H,1+nx,1+ny);
x = pt(1,:)'; y = pt(2,:)';

% Estimate mean of FE-solution u
M = 1e3;
meanS = zeros(length(x),1);
meanF = zeros(length(x),1);
meanU = zeros(length(x),1);
for m = 1:M,
    f_sample = eval(f);
    s_sample = eval(s);
    u_sample = fem2d_lin(L,H,1+nx,1+ny,a,f_sample,g);
    
    % online computation of mean
    meanS = (m-1)/m * meanS + 1/m * s_sample;
    meanF = (m-1)/m * meanF + 1/m * f_sample;
    meanU = (m-1)/m * meanU + 1/m * u_sample;
end
figure(); 
fem2d_lin_plot(u_sample,pt,el); xlabel('x'); ylabel('y');

% Plot both
figure(); 
subplot(1,2,1); 
fem2d_lin_plot(meanF,pt,el); xlabel('x'); ylabel('y');
title('Empirical mean of f')
subplot(1,2,2); 
fem2d_lin_plot(eval(Ef),pt,el); xlabel('x'); ylabel('y');
title('Theoretical mean of f')

% Plot both
figure(); 
subplot(1,2,1); 
fem2d_lin_plot(meanS,pt,el); xlabel('x'); ylabel('y');
title('Empirical mean of s')
subplot(1,2,2); 
fem2d_lin_plot(eval(Es),pt,el); xlabel('x'); ylabel('y');
title('Theoretical mean of s')

% Compute FE-solution when f is replaced by its mean
Eu = fem2d_lin(L,H,1+nx,1+ny,a,Ef,Eg);

% Plot both: empirical and theoretical mean of u
figure(); 
subplot(1,2,1); 
fem2d_lin_plot(meanU,pt,el); xlabel('x'); ylabel('y');
title('Empirical mean of u')
subplot(1,2,2); 
fem2d_lin_plot(Eu,pt,el); xlabel('x'); ylabel('y');
title('Theoretical mean of u')

%% (2): Random PDE 1 - Random Coefficient
clear all; clc;

% Set domain
L = 1; H = 1;

% Random diffusion term
a = 'exp(randn(1) * sin(pi*x) .* cos(pi*y))';

% ... and its expectation
Ea = 'exp(0.5 * sin(pi*x).^2 .* cos(pi*y).^2)';

% Deterministic source term
f = 'exp(0*x)';
    
% Deterministic boundary data (homogeneous)
g{1} = '0*x';     % Southern boundary data
g{2} = '0*y';     % Eastern boundary data
g{3} = '0*x';     % Northern boundary data
g{4} = '0*y';     % Western boundary data

% Set number of elements and generate mesh
nx = 64; ny=64;
[pt, el] = fem2d_lin_mesh(L,H,1+nx,1+ny);
x = pt(1,:)'; y = pt(2,:)';

% Estimate mean of FE-solution u
M = 1e3;
meanA = zeros(length(x),1);
meanU = zeros(length(x),1);
for m = 1:M,
    a_sample = eval(a);
    u_sample = fem2d_lin(L,H,1+nx,1+ny,a_sample,f,g);
    
    % online computation of mean
    meanA = (m-1)/m * meanA + 1/m * a_sample;
    meanU = (m-1)/m * meanU + 1/m * u_sample;
end

% Compute FE-solution when a is replaced by its mean
Eu = fem2d_lin(L,H,1+nx,1+ny,Ea,f,g);

% Plot both: mean of u, u for mean of f
figure(); 
subplot(1,2,1); 
fem2d_lin_plot(meanU,pt,el); xlabel('x'); ylabel('y');
title('Empirical mean of u')
subplot(1,2,2); 
fem2d_lin_plot(Eu,pt,el); xlabel('x'); ylabel('y');
title('u for mean of a')
