1 % laplace.m
   2 % Author: Lorenzo Bolla <lbolla (at) gmail dot com>
   3 % Last modified: Apr. 11, 2007
   4  
   5 xmin = 0;
   6 xmax = 1;
   7 ymin = 0;
   8 ymax = 1;
   9  
  10 nx = 200;
  11 ny = nx;
  12  
  13 dx = (xmax - xmin) / (nx - 1);
  14 dy = (ymax - ymin) / (ny - 1);
  15  
  16 x = xmin:dx:xmax;
  17 y = ymin:dy:ymax;
  18  
  19 n_iter = 100;
  20 eps = 1e-16;
  21  
  22 u = zeros(nx,ny);
  23 u(1,:)   = xmin^2 - y.^2;
  24 u(end,:) = xmax^2 - y.^2;
  25 u(:,1)   = x.'.^2 - ymin^2;
  26 u(:,end) = x.'.^2 - ymax^2;
  27 err = Inf;
  28 count = 0;
  29 while err > eps && count < n_iter
  30     count = count + 1;
  31     dx2 = dx^2;
  32     dy2 = dy^2;
  33     dnr_inv = 0.5 / (dx2 + dy2);    
  34     u_old = u;
  35     u(2:end-1, 2:end-1) = ((u(1:end-2, 2:end-1) + u(3:end, 2:end-1))*dy2 + (u(2:end-1,1:end-2) + u(2:end-1, 3:end))*dx2)*dnr_inv;
  36     v = (u - u_old);
  37     err = sqrt(v(:).'*v(:));
  38 end