%
% This tests the use of multilevel Monte Carlo for a simple
% elliptic PDE with random forcing:
%
% u(0)=0, u(1)=0, u''(x)= 50 Z sin(pi x) where Z = N(0,1)
%
% the output is P = \int_0^1 u^2 dx
%

function ellip

close all; clear all;

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

% global      % commented out because not needed in this example

%
% call mlmc_test to perform MLMC tests
%

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

N     = 100000; % samples for convergence tests
L     = 8;      % levels for convergence tests 

N0    = 200;   % 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(@ellip_l, M, N,L, N0,Eps, nvert);


%-------------------------------------------------------
%
% level l estimator for elliptic solver
%

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

% global

nf = 2^(l+1);
hf = 1/nf;
Af = hf^(-2)*spdiags(repmat([-1 2 -1],nf-1,1),-1:1,nf-1,nf-1);
xf = 50*sin(pi*(1:nf-1)*hf)';

if l>0
  nc = nf/2;
  hc = 1/nc;
  Ac = hc^(-2)*spdiags(repmat([-1 2 -1],nc-1,1),-1:1,nc-1,nc-1);
  xc = 50*sin(pi*(1:nc-1)*hc)';
end

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

for N1 = 1:100:N         % compute samples 100 at a time
  N2 = min(100,N-N1+1);

  Z  = rand(1,N2);

  uf = Af \ (xf*Z);
  Pf = sum(hf*uf.^2,1);

  if l==0
    Pc = zeros(1,N2);
  else
    uc = Ac \ (xc*Z);
    Pc = sum(hc*uc.^2,1);
  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
