% steepest descent algorithm for the FD discretization of Poisson
% eqn
%
% -u'' = f
%
% on a 1D uniform grid.

m = 80; a = 0; b = 1;
x = linspace(a,b,m+2); x = x(:);
h = (b-a)/(m+1);
A = spdiags(repmat([-1,2,-1]/h^2,m,1),...
            [-1:1],m,m);
b = sin(pi*x(2:end-1)) + 8*sin(4*pi*x(2:end-1));
uu= A\b;

u=randn(m,1);
r=b-A*u; p=A*r;

i = 0;

while norm(r)/norm(b) > 1.0E-06,
    i = i+1;
    alpha = dot(r,r)/dot(p,r);
    u = u+alpha*r;
    r = r-alpha*p;
    p = A*r;
    
    % visualisation
    subplot(1,2,1);
    plot(x(2:end-1),uu,'k',x(2:end-1),u,'r');
    title(sprintf('Iteration: %5d',i));
    xlabel('x'); ylabel('u');
    subplot(1,2,2);
    plot(x(2:end-1),uu-u);
    title(sprintf('Iteration: %5d',i));
    xlabel('x'); ylabel('e');
    drawnow; pause(0.0001);
end


