% jacobi - Program to solve the Laplace equation using 

% Jacobi's method on a square grid

clear; help jacobi;  % Clear memory and print header

N = input('input number of grid points on a side - ');

L = 1;        % System size (length)

h = L/(N-1);  % Grid spacing

phi0 = 1;     % Potential at y=L

fprintf('Potential at y=L equals %g \n',phi0);

fprintf('Potential is zero on all other boundaries\n');

change_want = 1e-4;   % Stop when the change is given fraction

%%%%% Set initial conditions and boundary conditions %%%%%

coeff = pi/L; 

x = (0:N-1)*h;     % x coordinate

y = (0:N-1)*h;     % y coordinate

% Initial guess is first term in separation of variables soln.

phi = phi0 * 4/(pi*sinh(pi)) * sin(coeff*x')*sinh(coeff*y); 

phi(:,N) = phi0*ones(N,1);    % Set the boundary condition

%%%%% MAIN LOOP %%%%%

flops(0);  % Reset the flops counter to zero;

newphi = phi;              %%!!%% Jacobi only

%omega = 2/(1+sin(pi/N));  %%!!%% SOR only

max_iter = N^2;  % Set max to avoid excessively long runs

for iter=1:max_iter

  temp = 0;

  for i=2:(N-1)      % Loop over interior points only

   for j=2:(N-1)

     %% Jacobi method %%

     newphi(i,j) = .25*(phi(i+1,j)+phi(i-1,j)+ ...

                               phi(i,j-1)+phi(i,j+1));

     temp = temp + abs(1-phi(i,j)/newphi(i,j));

    %% G-S method %%

%    newphi = .25*(phi(i+1,j)+phi(i-1,j)+ ...

%                               phi(i,j-1)+phi(i,j+1));

%    temp = temp + abs(1-phi(i,j)/newphi);

%    phi(i,j) = newphi;

    %% SOR method %%

%    newphi = .25*omega*(phi(i+1,j)+phi(i-1,j)+ ...

%               phi(i,j-1)+phi(i,j+1))  +  (1-omega)*phi(i,j);

%    temp = temp + abs(1-phi(i,j)/newphi);

%    phi(i,j) = newphi;

   end

  end

  phi = newphi;   %%!!%% Reset phi (Jacobi only)

  % Break out of the main loop if the average fractional change

  % in an iteration is less than change_want

  change(iter) = temp/(N-2)^2;

  fprintf('After %g iterations, fractional change = %g\n',...

                            iter,change(iter));

  if( change(iter) < change_want ) 

    disp('Desired accuracy achieved; breaking out of main loop');

    break;

  end

end

fprintf('Number of flops = %g\n',flops);

subplot(121)

  cs = contour(rot90(phi),9,x,y);  % Contour plot with labels

  xlabel('x'); ylabel('y'); clabel(cs,[.2 .4 .6 .8])

subplot(122)

  mesh(phi);   % Wire-mesh plot of potential

subplot(111)

title('Potential');

pause;  % Pause between plots; strike any key to continue

semilogy(change);

xlabel('iteration');  ylabel('fractional change');