NX = 200;
NY = 100;

// result is 2D array with bounary value (size is [NY + 2][NX + 2])
var result;
result = solveLaplaceEquationByJacobiMethod(NX, NY, 1.0 / NX, 1.0 / NY, makeExampleBoundaryFunc2(NX, NY), {eps: 0.01});
result = solveLaplaceEquationByGaussSeidelMethod(NX, NY, 1.0 / NX, 1.0 / NY, makeExampleBoundaryFunc2(NX, NY), {eps: 0.01});
result = solvePoissonEquationByJacobiMethod(NX, NY, 1.0 / NX, 1.0 / NY, makeExamplePoissonFunc(NX, NY), makeExampleBoundaryFunc2(NX, NY), {eps: 0.01});
result = solvePoissonEquationByGaussSeidelMethod(NX, NY, 1.0 / NX, 1.0 / NY, makeExamplePoissonFunc(NX, NY), makeExampleBoundaryFunc2(NX, NY), {eps: 0.01

// example of boundary codition
function makeExampleBoundaryFunc1(nx, ny) {
  return function(g, x, y) {
    if (x === 0) {
      return -1.0 + 2.0 * y / (ny + 1.0);
    } else if (x === nx + 1) {
      return 1.0 - 2.0 * y / (ny + 1.0);
    } else if (y === 0) {
      return -1.0 + 2.0 * x / (nx + 1.0);
    } else if (y === ny + 1) {
      return 1.0 - 2.0 * x / (nx + 1.0);
    }
  }
}

// example of boundary condition
function makeExampleBoundaryFunc2(nx, ny) {
  return function(g, x, y) {
    if (y === 0) {
      return Math.sin(x * 0.1);
    } else if (y === ny + 1) {
      return Math.cos(x * 0.1);
    } else {
      return 0.0;
    }
  }
}

// example of poisson function
function makeExamplePoissonFunc(nx, ny) {
  return function(x, y) {
    var d = Math.sqrt(Math.pow(x / nx - 0.5, 2.0) + Math.pow(y / ny - 0.5, 2.0));
    var s = 0.2;
    return 5.0 * exp(- d * d / (2.0 * s * s)) / sqrt(2.0 * Math.PI * s * s);
  }
}

/**
 * Creates 2D array
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {function} defaultFunc - function to get value
 * @return {array} - 2d array
 */
function createArray2D(nx, ny, defaultFunc) {
  var defaultFunc = defaultFunc !== undefined ? defaultFunc : function(x, y) {return 0.0;};
  var array = new Array(ny);
  for (var y = 0; y < ny; y++) {
    array[y] = new Array(nx);
    for (var x = 0; x < nx; x++) {
      array[y][x] = defaultFunc(x, y);
    }
  }
  return array;
}

/**
 * Solves 2D Laplace equation with boundary condition by Gauess SeidelMethod
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {number} dx - size of x grid
 * @param {number} dy - size of y grid
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solveLaplaceEquationByJacobiMethod(nx, ny, dx, dy, boundaryFunc, opts) {
  var dx2 = dx * dx;
  var dy2 = dy * dy;
  var func = function(f, x, y) {
    return (dy2 * (f[y][x - 1] + f[y][x + 1]) + dx2 * (f[y - 1][x] + f[y + 1][x])) / (2.0 * (dx2 + dy2));
  }

  return solveByJacobiMethodFor2D(nx, ny, func, boundaryFunc, opts);
}

/**
 * Solves 2D Laplace equation with boundary condition by Gauess SeidelMethod
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {number} dx - size of x grid
 * @param {number} dy - size of y grid
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solveLaplaceEquationByGaussSeidelMethod(nx, ny, dx, dy, boundaryFunc, opts) {
  var dx2 = dx * dx;
  var dy2 = dy * dy;
  var func = function(f, x, y) {
    return (dy2 * (f[y][x - 1] + f[y][x + 1]) + dx2 * (f[y - 1][x] + f[y + 1][x])) / (2.0 * (dx2 + dy2));
  }

  return solveByGaussSeidelMethodFor2D(nx, ny, func, boundaryFunc, opts);
}

/**
 * Solves 2D Poisson equation with boundary condition by Jacobi method
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {number} dx - size of x grid
 * @param {number} dy - size of y grid
 * @apram {function} poissonFunc - function for poisson equation
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solvePoissonEquationByJacobiMethod(nx, ny, dx, dy, poissonFunc, boundaryFunc, opts) {
  var dx2 = dx * dx;
  var dy2 = dy * dy;
  var func = function(f, x, y) {
    return (dy2 * (f[y][x - 1] + f[y][x + 1]) + dx2 * (f[y - 1][x] + f[y + 1][x]) - dx2 * dy2 * poissonFunc(x, y)) / (2.0 * (dx2 + dy2));
  }

  return solveByJacobiMethodFor2D(nx, ny, func, boundaryFunc, opts);
}


/**
 * Solves 2D Poisson equation with boundary condition by Gauess Seidel method
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {number} dx - size of x grid
 * @param {number} dy - size of y grid
 * @apram {function} poissonFunc - function for poisson equation
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solvePoissonEquationByGaussSeidelMethod(nx, ny, dx, dy, poissonFunc, boundaryFunc, opts) {
  var dx2 = dx * dx;
  var dy2 = dy * dy;
  var func = function(f, x, y) {
    return (dy2 * (f[y][x - 1] + f[y][x + 1]) + dx2 * (f[y - 1][x] + f[y + 1][x]) - dx2 * dy2 * poissonFunc(x, y)) / (2.0 * (dx2 + dy2));
  }

  return solveByGaussSeidelMethodFor2D(nx, ny, func, boundaryFunc, opts);
}

/**
 * Solves 2D equation with boundary condition by Jacobi method
 * @param {number} nr - number of repetiation
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {function} func - function which represents equation
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solveByJacobiMethodFor2D(nx, ny, func, boundaryFunc, opts) {
  var nr = opts && opts.nr !== undefined ? opts.nr : null;
  var eps = opts && opts.eps !== undefined ? opts.eps : null;
  if (nr === null && eps === null) {
    nr = 100;
    eps = 0.01;
  }
  
  var pg = createArray2D(nx + 2, ny + 2);
  var ng = createArray2D(nx + 2, ny + 2);
  
  var setBoundary = function(g) {
    for (var x = 1; x < nx + 1; x++) {
      g[0][x] = boundaryFunc(g, x, 0);
      g[ny + 1][x] = boundaryFunc(g, x, ny + 1);
    }
    for (var y = 1; y < ny + 1; y++) {
      g[y][0] = boundaryFunc(g, 0, y);
      g[y][nx + 1] = boundaryFunc(g, nx + 1, y);
    }
    g[0][0] = boundaryFunc(g, 0, 0);
    g[0][nx + 1] = boundaryFunc(g, nx + 1, 0);
    g[ny + 1][0] = boundaryFunc(g, 0, ny + 1);
    g[ny + 1][nx + 1] = boundaryFunc(g, nx + 1, ny + 1);
  }
  setBoundary(pg);
  setBoundary(ng);

  var step = function() {
    var error = 0.0;
    var temp = pg;
    pg = ng;
    ng = temp;
    for (var y = 1; y < ny + 1; y++) {
      for (var x = 1; x < nx + 1; x++) {
        var v = func(pg, x, y);
        error += Math.abs(pg[y][x] - v);
        ng[y][x] = v;
      }
    }
    return error;
  }
  
  if (nr) {
    for (var r = 0; r < nr; r++) {
      var e = step();
      if (eps && e < eps) {break;}
    }
  } else {
    while(true) {
      var e = step();
      if (e < eps) {break;}
    }
  };

  return ng;
}

/**
 * Solves 2D equation with boundary condition by Gauss-Seidel method
 * @param {number} nr - number of repetiation
 * @param {number} nx - number of rows
 * @param {number} ny - number of columns
 * @param {function} func - function which represents equation
 * @param {function} boundaryFunc - function to set boundary
 * @param {Object} opts - options
 * @param {number} opts.nr - number of repetiation 
 * @param {number} opts.eps - minimum error to stop repetiation
 * @return {array} - 2D array with boundary (size = array[ny + 2][nx + 2])
 */
function solveByGaussSeidelMethodFor2D(nx, ny, func, boundaryFunc, opts) {
  var nr = opts && opts.nr !== undefined ? opts.nr : null;
  var eps = opts && opts.eps !== undefined ? opts.eps : null;
  if (nr === null && eps === null) {
    nr = 100;
    eps = 0.01;
  }
  
  var g = createArray2D(nx + 2, ny + 2);
  
  // set up boundary
  for (var x = 1; x < nx + 1; x++) {
    g[0][x] = boundaryFunc(g, x, 0);
    g[ny + 1][x] = boundaryFunc(g, x, ny + 1);
  }
  for (var y = 1; y < ny + 1; y++) {
    g[y][0] = boundaryFunc(g, 0, y);
    g[y][nx + 1] = boundaryFunc(g, nx + 1, y);
  }
  g[0][0] = boundaryFunc(g, 0, 0);
  g[0][nx + 1] = boundaryFunc(g, nx + 1, 0);
  g[ny + 1][0] = boundaryFunc(g, 0, ny + 1);
  g[ny + 1][nx + 1] = boundaryFunc(g, nx + 1, ny + 1); 

  var step = function() {
    var error = 0.0;
    for (var y = 1; y < ny + 1; y++) {
      for (var x = 1; x < nx + 1; x++) {
        var v = func(g, x, y);
        error += Math.abs(g[y][x] - v);
        g[y][x] = v;
      }
    }
    return error;
  }
  
  if (nr) {
    for (var r = 0; r < nr; r++) {
      var e = step();
      if (eps && e < eps) {break;}
    }
  } else {
    while(true) {
      var e = step();
      if (e < eps) {break;}
    }
  };

  return g;
}