2D-уравнение Лапласа с JavaScript

Применение метода конечных разностей. Решение с помощью Math.js и построение графика с помощью Plotly.js

Двумерное уравнение Лапласа и метод центрированных разностей

Вот двумерное уравнение Лапласа

Идея состоит в том, чтобы выполнить метод центрированной разности как в направлениях x, так и в направлении y, и сложить их вместе.

Рассмотрим точку сетки, обозначенную i, j в координатах x и y, метод центрированной разности дает результат в направлении x (при условии, что Δx = Δy),

метод центрированной разности дает в направлении Y дает

Складывая два уравнения рядом друг с другом, мы будем иметь

Это уравнение говорит о том, что в точке сетки i и j решение определяется соседними значениями сверху, слева, справа и снизу точки сетки.

Граничное условие

Так же, как и в одномерном случае, мы устанавливаем граничные условия. Представьте себе стены на севере, востоке, юге и западе.

Я считаю, когда i = 0, как Запад, j = 0 как Юг, i = N как Восток, j = N как Север. Таким образом, левый нижний угол — это координата (0, 0).

Построить матрицу

Например, 4 x 4 точки сетки. Всего точек сетки 16. Внутренние 4 из 16 точек можно решить с помощью уравнения. Все точки внешней сетки определяются граничными условиями. Следовательно, нужно решить 4 уравнения.

Эту систему линейных уравнений можно преобразовать в следующую

Тогда его можно записать в матричной форме

или просто скажем

Точно так же, как мы делали это в одномерном случае, мы можем решить матричное уравнение с помощью решателя.

Кстати, я изначально использовал итерационный метод, такой как метод Гаусса-Зейделя, и он тоже был хорош. Изменение формы раствора не требуется. Здесь я использую решатель, используемый в одномерном случае, так как мне кажется, что он более согласован.

Получив решение u, нам нужно изменить форму u для графика, поскольку u задается в векторном формате в матричном уравнении. Также обратите внимание, что это u не включает границы. Поэтому я называю это inner_u в разделе кодирования. Включая границу, я называю это u. Просто для уточнения.

Кодирование

Основной поток:

Установите размер сетки и граничное условие
Построить матрицу A и вектор b
Найди решение
Изменить форму
Сюжет

1. Установите размер сетки и граничное условие

Квадратная сетка, ширина сетки 4
Северное граничное значение равно 2, восток 2, юг 1 и запад 1.

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

2. Построить матрицу A и вектор b

Я вызываю функцию create2DLaplaceMatrices, и она принимает указанные выше условия.

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

Чтобы заставить это работать, я считаю следующее.

Номер решения определяется как (Ширина сетки -2) x (Ширина сетки -2) из-за граничного условия
Имя (Ширина сетки — 2) как «Inner_Grid_Width»
Вектор решения u или inner_u в порядке от северо-западного угла к юго-восточному.
Установите индекс решения, чтобы отразить вышеуказанный порядок (индексы 2d для индекса 1d)
В точке (i, j) элемент Matrix равен -4, соседи (сверху, справа, снизу, слева) равны 1.
Если сосед находится на границе, установите граничное условие
Обратите внимание, что соседи указаны с точки зрения индекса решения, поскольку решение находится в векторной форме. Следовательно, верх и низ — это либо плюс, либо минус внутренней ширины сетки.

Так должно выглядеть:

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 };
}

Это даст желаемый результат.

3. Решите это!

Чтобы решить матричное уравнение, используйте тот же решатель, что и в одномерном случае, функцию math.lusolve. Поэтому перед запуском убедитесь, что математическая библиотека загружена. Для решения просто

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

Просто помните, что решение находится в виде вектора, а значения хранятся в массиве в массиве. Я называю решение inner_u, потому что inner_u не содержит граничных значений.

4. Изменить форму

Прежде чем построить решение, измените форму одномерного массива inner_u на двумерный массив u. Я называю эту функцию «изменить форму».

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

Чтобы создать двумерный массив, создайте массив u и сначала вставьте массив Grid Width с северным граничным условием. Затем вставьте решения из inner_u, затем вставьте западное и восточное граничные значения и добавьте южные граничные значения. Наконец, углы брали средние значения двух соседей.

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;
}

Сюжет

Наконец, Plotly.js упрощает построение трехмерных поверхностей. Давайте создадим функцию графика поверхности «plotSurface», взяв решение u.

plotSurface(u);

Убедитесь, что тег div с идентификатором «myDiv» подготовлен в html-файле.

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

Если все готово, загрузите подготовленный index.html в одномерном случае, вы должны увидеть

После примечания

Используя решатель, вы должны почувствовать задержку, как только вы увеличите ширину сетки до 20 или около того, в зависимости от вашей машины. Итерационный метод с Gauss Seidel и SOR значительно ускорил процесс. Даже должна быть лучшая библиотека для решения линейных уравнений. Что-то рассмотреть.

Что попробовать:

Итерационный метод решения
Решите уравнение Пуассона