Лапласиан конечной разности

Этот пример показывает, как вычислить и представлять Лапласиан конечной разности на L-образной области.

Область

Числа функции numgrid указывают в L-образной области. Функция spy является полезным инструментом для визуализации шаблона ненулевых элементов в матрице. Используйте эти две функции, чтобы сгенерировать и отобразить L-образную область.

n = 32;
R = 'L';
G = numgrid(R,n);
spy(G)
title('A Finite Difference Grid')

Покажите уменьшенный вариант матрицы как выборка.

g = numgrid(R,10)
g = 10×10

     0     0     0     0     0     0     0     0     0     0
     0     1     5     9    13    17    25    33    41     0
     0     2     6    10    14    18    26    34    42     0
     0     3     7    11    15    19    27    35    43     0
     0     4     8    12    16    20    28    36    44     0
     0     0     0     0     0    21    29    37    45     0
     0     0     0     0     0    22    30    38    46     0
     0     0     0     0     0    23    31    39    47     0
     0     0     0     0     0    24    32    40    48     0
     0     0     0     0     0     0     0     0     0     0

Дискретный лапласиан

Используйте delsq, чтобы сгенерировать дискретный Лапласиан. Используйте функцию spy снова, чтобы получить графическое ощущение элементов матрицы.

D = delsq(G);
spy(D)
title('The 5-Point Laplacian')

Определите количество внутренних точек.

N = sum(G(:)>0)
N = 675

Краевая задача Дирихле

Решите краевую задачу Дирихле для разреженной линейной системы. Настройка задач:

delsq(u) = 1 во внутренней части, u = 0 на контуре.

rhs = ones(N,1);
if (R == 'N') % For nested dissection, turn off minimum degree ordering.
   spparms('autommd',0)
   u = D\rhs;
   spparms('autommd',1)
else
   u = D\rhs; % This is used for R=='L' as in this example
end

Сопоставьте решение на L-образную сетку и постройте график его как контурной карты.

U = G;
U(G>0) = full(u(G(G>0)));
clabel(contour(U));
prism
axis square ij

Теперь покажите решение как сетчатый график.

mesh(U)
axis([0 n 0 n 0 max(max(U))])
axis square ij

Смотрите также

Была ли эта тема полезной?