% ch6_mg2.m --- Demo of smoothing factor  -- Damped Jacobi

help ch6_mg2
 n=16; h = 1/n; x = (h:h:1-h)'; % Internal nodes

g = h^2*(x.^3-2*x-4);  % RHS
n = n-1
A = diag(2*ones(n,1)) + diag(-ones(n-1,1),-1) + diag(-ones(n-1,1),1);
x_sol = A \ g; % Exact solution to  -u''=g

% Now Jacobi
 w = 1/3, 
 s = 'v+o*^><sd'; use=1;

rand('seed',6);x=rand(n,1); plot(1:n, x_sol-x,'r-x'); hold on
  D = diag(diag(A));  % B = L + U
  B = diag(diag(A))-A;
  g = D\g;
for iter=1:88
  x1  = w*x + (1-w)*( D\(B*x) + g );
  if find(iter==[1 2 6 7 88])>0
    plot(1:n, x_sol-x1,['-' s(use)]);
    use=use+1; if use>length(s),use=1; end
  end
  x=x1;
end
  legend('Initial error','1 iteration','2 iterations','6 iterations',...
   '7 iterations','88 iterations',0)
  axis tight
  str=sprintf('Damped Jacobi omega = %4.2f\n',w);xlabel(str)
  hold off
