2D Laplace Equation with JavaScript

Application of Finite Difference Method. Solving with Math.js, and plot with Plotly.js

2D Laplace Equation and Centered Difference Method

Here is the two dimensional Laplace equation

The idea is to perform centered difference method both in x and y directions, and put them together.

Consider the grid point denoted by i, j in x and y coordinate, the centered difference method yields in x direction yields (assuming Δx=Δy),

the centered difference method yields in y direction yields

Adding the two equations side by side together, we will have

This equation is telling that at grid point i and j, the solution is given by its neighbor values, top, left, right and bottom of the grid point.

Boundary Condition

Just like one dimensional case, we set the boundary conditions. Imagine there are walls at north, east, south, and west side.

I consider when i = 0 as West, j = 0 as South, i = N as East, j = N as North. So the left lower corner is the (0, 0) coordinate.

Build Matrix

For example, 4 x 4 grid points. There are total of 16 grid points. The inner 4 of 16 points can be solved using the equation. All the outer grid points are determined by boundary conditions. Hence there are 4 equations to solve.

This system of linear equations can be rearranged to the following

Then it can be written in Matrix form

or simply put it

Just like we did it in one dimensional case, we can solve the Matrix equation using the solver.

By the way, I initially used iterative method such as Gauss Seidel method, and it was good as well. Reshaping of the solution is not necessary. Here, I use the solver used in one dimensional case as I felt like it is more consistent.

Once get the solution of u, we need to reshape u for the plot because u is given in vector format in the Matrix equation. Also note that this u does not include boundaries. So I call it inner_u below coding section. Including boundary is I call it u. Just to clarify.

Coding

Basic flow:

1. Set Grid size and boundary condition
2. Build Matrix A and Vector b
3. Solve it
4. Reshape
5. Plot

1. Set Grid size and boundary condition

• Square grid, and Grid width is 4
• North Boundary Value is 2, East 2, South 1 and West 1

const GRID_WIDTH = 4
const DIRICHLET = {
  NORTH: 2,
  EAST: 2,
  SOUTH: 1,
  WEST: 1
}

2. Build Matrix A and Vector b

I call the function ‘create2DLaplaceMatrices’ and it takes the conditions set above.

let { A, b } = create2DLaplaceMatrices(GRID_WIDTH, DIRICHLET);

To make this to work, I consider the following.

• Number of solution is given by (Grid Width -2) x (Grid Width -2) because of boundary condition
• Name(Grid Width — 2) as “Inner_Grid_Width”
• Solution vector u, or inner_u, in the order from North West corner to South East
• Set the solutionIndex to reflect the above order (2d indices to 1d index)
• At point (i, j), the Matrix element is -4, neighbors (top, right, bottom, left) is 1.
• If neighbor is at the boundary, set boundary condition
• Note neighbors are indicated in terms of solutionIndex because solution is in vector form. Hence, top and bottom are either plus or minus of Inner Grid Width.

So it should look:

function create2DLaplaceMatrices(GRID_WIDTH, DIRICHLET) {
  const INNER_GRID_WIDTH = GRID_WIDTH - 2;
  const NUMBER_OF_SOLUTIONS = (INNER_GRID_WIDTH)*(INNER_GRID_WIDTH);
  let solutionIndex = 0;
  let A = [];
  let b = new Array(NUMBER_OF_SOLUTIONS).fill(0);
  for (j = 0; j < INNER_GRID_WIDTH; j++) {
    for (i = 0; i < INNER_GRID_WIDTH; i++) {
      A.push(new Array(NUMBER_OF_SOLUTIONS).fill(0));
      solutionIndex = i + (j) * (INNER_GRID_WIDTH);
      A[solutionIndex][solutionIndex] = -4;
      if (i == 0) {
        b[solutionIndex] -= DIRICHLET.WEST;
      } else {
        A[solutionIndex][solutionIndex - 1] = 1
      }
      if (j == 0) {
        b[solutionIndex] -= DIRICHLET.NORTH;
      } else {
        A[solutionIndex][solutionIndex - (INNER_GRID_WIDTH)] = 1;
      }
      if (i == (INNER_GRID_WIDTH - 1)) {
        b[solutionIndex] -= DIRICHLET.EAST;
      } else {
        A[solutionIndex][solutionIndex + 1] = 1;
      }
      if (j == (INNER_GRID_WIDTH - 1)) {
        b[solutionIndex] -= DIRICHLET.SOUTH;
      } else {
        A[solutionIndex][solutionIndex + (INNER_GRID_WIDTH)] = 1;
      }
    }
  }
  return { A, b };
}

This will make the desired result.

3. Solve it!

To solve the Matrix equation, use the same solver used in one dimensional case, the function math.lusolve. So make sure the math library is loaded before run. In order to solve, simply

let inner_u = math.lusolve(A, b);

Just to remember, the solution is in the form of vector, and the values are stored in array in the array. I call the solution inner_u because inner_u does not contain boundary values.

4. Reshape

Before plot the solution, reshape the one dimensional array inner_u to two dimensional array u. I call the function ‘reshape’.

let u = reshape(inner_u, GRID_WIDTH, DIRICHLET);

To make two dimensional array, create an array u, and first push a Grid Width array with North boundary condition. Then insert the solutions from the inner_u, then insert West and East boundary value, and add the South boundary values. Finally, corners took the average values of two neighbors.

function reshape(inner_u, GRID_WIDTH, DIRICHLET) {
  const INNER_GRID_WIDTH = GRID_WIDTH — 2;
  let u = [];
  u.push(new Array(GRID_WIDTH).fill(DIRICHLET.NORTH));
  while (inner_u.length) u.push((inner_u.splice(0, INNER_GRID_WIDTH)).flat());
  for (j = 1; j < GRID_WIDTH — 1; j++) {
    u[j].splice(0, 0, DIRICHLET.WEST);
    u[j].splice(u[j].length, 0, DIRICHLET.EAST);
  }
  u.push(new Array(GRID_WIDTH).fill(DIRICHLET.SOUTH));
  u[0][0] = (u[0][1] + u[1][0]) / 2;
  u[0][GRID_WIDTH — 1] = (u[0][GRID_WIDTH — 2] + u[1][GRID_WIDTH — 1]) / 2;
  u[GRID_WIDTH — 1][0] = (u[GRID_WIDTH — 2][0] + u[GRID_WIDTH — 1][1]) / 2;
  u[GRID_WIDTH — 1][GRID_WIDTH — 1] = (u[GRID_WIDTH — 2][GRID_WIDTH — 1] + u[GRID_WIDTH — 1][GRID_WIDTH — 2]) / 2;
  return u;
}

Plot

Finally, Plotly.js makes three dimensional surface plot easy. Let’s make surface plot function ‘plotSurface’ taking the solution u.

plotSurface(u);

Make sure div tag with id named “myDiv” was prepared in html file.

function plotSurface(u) {
  let data = { z: u, type: ‘surface’ };
  let layout = {
    font: { size: 16 }
  };
  let config = { responsive: true };
  Plotly.newPlot(‘myDiv’, [data], layout, config);
}

If everything is ready, load the index.html prepared in one dimensional case, you should see

After Note

Using the solver, you should feel delay once you increase the Grid Width to 20 or so, depending on your machine. Iterative method with Gauss Seidel and SOR significantly sped up the process. Even there should be better library for linear equation solver. Something to consider.

Things to try:

• Iterative Method to solve
• Solve Poisson Equation

Reference

• Numerical methods for Laplace’s equation
• Finite Difference Method for Laplace Equation