
% problem data
dirichlet_values = [0, 0];
f = @(x) (100*x.*(1-x));

lower_f = @(x) 0*x;
upper_f = @(x) 0*x+1;

% grid points
n = 100;

% iteration steps
maxnu = 2000;




% mesh size
h = 1/(n-1);
grid = linspace(0,1,n)';

% assemble rhs by trapezoidal rule
fx = f(grid);
b = fx * h;
b(1) = fx(1) * h/2;
b(n) = fx(n) * h/2;

% assemble matrix
A_diags = [-ones(n,1), 2*ones(n,1), -ones(n,1)] / h;
A = spdiags(A_diags, [-1 0 1], n, n);
A(1,1) = 1/h;
A(n,n) = 1/h;

% incorporate dirichlet bc
A(1,:) = 0;
A(n,:) = 0;
if (norm(dirichlet_values, 'inf')==0)
    A(:,1) = 0;
    A(:,n) = 0;
end
A(1,1) = 1;
b(1) = dirichlet_values(1);
A(n,n) = 1;
b(end) = dirichlet_values(end);


% interpolate obstacles
lower = lower_f(grid);
upper = upper_f(grid);

x0 = zeros(n, 1);

% Additional last argument for switching storage of all iterates on
% Second return value is a matrix containg x_0, ..., x_maxnu as
% column vectors
[x, xcache] = pgs_solver(A, b, lower, upper, x0, maxnu, 1);



% plot solution
figure(1);
clf;
plot(grid, x, 'o-r');



% plot energy
figure(2)
clf;

E = zeros(maxnu+1,1);
for nu=0:maxnu
    E(nu+1) = 1/2 * dot(A*xcache(:,nu+1),xcache(:,nu+1)) - dot(b,xcache(:,nu+1));
end
plot(0:maxnu, E, 'o-');



% plot approximate error
figure(3);
clf;

% compute approximate error
corr_norm = zeros(maxnu+1,1);
for nu = 0:maxnu
    corr = xcache(:,nu+1) - x;
    corr_norm(nu+1) = dot(A*corr, corr)^.5;
end
% plot only nonzero error values in semi log scale
nu = 0:maxnu;
nz = find(corr_norm);
semilogy(nu(nz), corr_norm(nz), 'o-');

hold on

% plot a reference triangle with slope O(1-h)
ref_nu = round([max(nu(nz))/4 max(nu(nz))*3/4]);
ref_corr(1) = corr_norm(ref_nu(1))/10;
ref_corr(2) = ref_corr(1)*(1-h)^(ref_nu(2)-ref_nu(1));
semilogy(ref_nu, ref_corr, '-r');

