% Matlab code makemgdemo1D.m
% For "Applied Numerical Linear Algebra",  Question 6.16
% Written by James Demmel, Apr 22, 1995
%                Modified, Jun  2, 1997
%
% Initialize b, x, jac1, jac2, iter to demonstrate testfmg1D
%
k=7;
kk=2^k+1;
[x,y]=meshgrid(1:kk,1:kk); %
% b=sin((0:kk-1)*3*pi/kk)'; 
t=2*eye(kk-2)-diag(ones(kk-3,1),1)-diag(ones(kk-3,1),-1);
[v,d]=eig(t); %
[ds,is]=sort(diag(d)); %
vs=v(:,is); %
% b=vs*(1.02.^(0:-1:-kk+1)'.*randn(kk,1)); %
sol=zeros(kk-2,1);
sol((kk-1)/8:3*(kk-1)/8) = (1:(kk-1)/4 + 1);
sol(5*(kk-1)/8:7*(kk-1)/8) = -2*((kk-1)/4 + 1:-1:1);
xsol=[0;sol;0];
b=[0;t*sol;0];
nb = norm(b,'inf');
b=b/nb;
sol = sol/nb;
xsol = xsol/nb;
[nb,mb]=size(b);
hold off; clf
% x=b;b=zeros(size(b));
figure(1)
hold off, clf
subplot(121)
% axes('position',[.1,.1,.35,.5])
plot(b,'k')
axis([1,nb,-1,1])
title('right hand side b')
subplot(122)
% axes('position',[.55,.1,.35,.5])
plot(xsol,'k');
axis([1,nb,-1,1])
title('solution x')
figure(2)
hold off, clf
plot(vs'*t*sol,'k');
axis([1,nb,-1,1])
title('frequency components of b')
x=zeros(size(b));
jac1=1;jac2=1;iter=8;