% numerical relaxation animation to find equilibrium for charges on a line
% Barnett 10/9/08
n = 20;                                    % # nodes - 1
eps = 1e-4; maxF = 30;                     % movement speed, max force
x = -2:.05:2; y = -1.5:0.05:1.5; figure;
[xx yy] = meshgrid(x, y); zz = xx + 1i*yy; % zz is rect array of complex #s
xj = [-1, rand(1,n-1), 1];                 % initial node positions
F = zeros(1,n+1);                          % forces on charges
for i=1:1e2
  a = poly(xj);                            % coeffs of monic poly with roots xj
  for m=1:50                               % inner loop
    for j=2:n                              % don't move the end guys
      F(j) = sum(1./(xj(j) - [xj(1:j-1) xj(j+1:end)]));
    end
    xj = xj + eps*max(min(F,maxF),-maxF);  % truncate to max force
  end
  phi = log(polyval(a,zz))/(n+1);
  contour(x,y,phi,[-1:.1:1]);
  hold on; plot(xj, 0, 'k.', 'markersize', 20); hold off;
  text(0,1,sprintf('max force on any node = %.3g', max(abs(F))))
  xlabel('Re z'); ylabel('Im z'); title('phi(z) contours'); drawnow
end