%% Example 5.9
% One-Dimensional V Cycle Multigrid
%

%% Setup
close all;
clear N n x0 xN h X b u num_vcycles v u r R R1 nv1;

N = 64;
n = N+1;
x0 = 0;
xN = 1;
h = 1/N;

X = linspace(x0,xN,N+1)';
b = (sin(pi*X) + sin(16*pi*X))/2;

%% V-cycles
u = zeros(N+1,1);
num_vcycles = 16;

for v = 1:num_vcycles
u = mgv1n(u,b);
r = zeros(n,1);
r(2:n-1) = b(2:n-1) + (2*u(2:n-1)-u(1:n-2)-u(3:n)) / h^2;
R(v) = max(abs(r));
end
R1 = R;
nv1 = num_vcycles;

%% V-cycles - simple restriction
u = zeros(N+1,1);
num_vcycles = 28;

for v = 1:num_vcycles
u = mgv1r(u,b);
r = zeros(n,1);
r(2:n-1) = b(2:n-1) + (2*u(2:n-1)-u(1:n-2)-u(3:n)) / h^2;
R(v) = max(abs(r));
end

%% Plot
figure(1)
semilogy(1:nv1,R1,'o-')
hold on
semilogy(1:num_vcycles,R,'ro-')

xlabel('Iteration')
ylabel('max |r_i|')
legend('Average Restriction','Simple Restriction')