%
% This tests the use of multilevel Monte Carlo for a simple
% parabolic PDE with random initial data:
%
% PDE: du/dt = d^2u/dx^2
%
% u(0,t)=0, u(1,t)=0, u(x,0) = 100 Z sin(pi x) where Z = N(0,1)
%
% the output is P = \int_0^1 u(x,0.25) dx
%

function para

close all; clear all;

%
% use global variables to pass parameters into gbm_l function
%

% global 

%
% call mlmc_test to perform MLMC tests
%

nvert = 3;      % plotting option (1 for slides, 2 for papers, 3 for full set)
M     = 8;      % refinement cost factor (2^gamma in general MLMC Thm)

N     = 1000;   % samples for convergence tests
L     = 6;      % levels for convergence tests 

N0    = 100;   % initial number of samples on each MLMC level
Eps   = [ 0.005 0.01 0.02 0.05 0.1 ];  % desired accuracies for MLMC calcs

figs = mlmc_test(@para_l, M, N,L, N0,Eps, nvert);


%-------------------------------------------------------
%
% level l estimator for parabolic solver
%

function [sum1 sum2] = para_l(l,N)

% global

lam = 0.25;  % timestep/spacing^2

nf = 2^(l+1);
hf = 1/nf;

if l>0
  nc = nf/2;
  hc = 1/nc;
end

sum1(1:4) = 0;
sum2(1:2) = 0;

for N1 = 1:100:N
  N2 = min(100,N-N1+1);
  Z  = rand(1,N2);

  uf = 50*sin(pi*(0:nf)*hf)' * Z;      % random initial data
  i  = 2:nf;

  for n = 1:nf^2                       % time-marching
    uf(i,:) = uf(i,:) + lam*(uf(i+1,:)-2*uf(i,:)+uf(i-1,:) );
  end

  Pf = sum(hf*uf.^2,1);                % output

  if l==0
    Pc = zeros(1,N2);
  else

    uc = 50*sin(pi*(0:nc)*hc)' * Z;    % random initial data
    i  = 2:nc;

    for n = 1:nc^2                     % time-marching
      uc(i,:) = uc(i,:) + lam*(uc(i+1,:)-2*uc(i,:)+uc(i-1,:) );
    end

    Pc = sum(hc*uc.^2,1);              % output
  end

  sum1(1) = sum1(1) + sum(Pf-Pc);
  sum1(2) = sum1(2) + sum((Pf-Pc).^2);
  sum1(3) = sum1(3) + sum((Pf-Pc).^3);
  sum1(4) = sum1(4) + sum((Pf-Pc).^4);
  sum2(1) = sum2(1) + sum(Pf);
  sum2(2) = sum2(2) + sum(Pf.^2);
end
