Реализация конечно-разностной схемы для решения уравнения Лапласа¶
В этом примере мы решим следующее уравнение Лапласа
$$\begin{align*} - \Delta u &= 1, \quad x \in \Omega \\ u &= 0, \quad x \in \partial \Omega. \end{align*}$$
используя 5-точечны конечно-разностную схему решения.
Создание сетки¶
Положим, нас интересует решение задачи на некоторой конкретной сетке, заданной в форме L-образной группы точек на плоскости.
gr()
Сперва мы создадим простую матрицу из нулей, затем заполним нужные области единицами, добавив также "окантовку" размером в 1 элемент для задания граничных условий Дирихле. Функция spy из набора графиков gr() хорошо подходит для визуализации областей, где значения элементов матриц ненулевые.
10÷2
function lshape(N)
#N = round(Int, N/2) * 2 # убедимся, что входной параметр – четный
L = zeros(Int, N, N) # создаем массив из чисел типа Int со значением 0
# заполним нужную область единицами
L[2:N÷2, 2:N-1] .= 1
L[2:N-1, N÷2+1:N-1] .= 1
# ÷ обозначает целочисленное деление
# (! julia не позволяет индексировать векторы при помощи чисел типа Float)
# Вместо единиц, подставим номера значений в матрице
items = findall( !iszero, L )
L[items] = 1:length(items)
return L
end
G = lshape( 32 );
spy( G .> 0, markersize=3, markercolor=1 )
Если нам будет нужно пронумеровать все эти точки, как это делает функция numgrid платформы MATLAB, воспользуемся следующим кодом:
g = lshape( 10 ) .> 10
g = lshape( 10 )
Можно увидеть, какие номера элементов входной матрицы должны участвовать в вычислении каких результатов. При 5-элеметной схеме, в вычислениях элемента 1 участвуют элементы 1, 2 и 5 (соседние элементы).
Оператор Лапласа¶
Теперь загрузим функцию, которая задает оператор Лапласа на произвольной сетке, используя для вычислений 4-связную область.
using SparseArrays, LinearAlgebra
include( "delsq.jl" ); # Создание лапласиана
Воспользуемся функцией spy чтобы лучше увидеть отдельные маркеры.
L = lshape(9);
G = delsq(L);
spy( G .!= 0, colorbar=false, markersize=3, markercolor=1 )
Хотя и стандартный вывод для разреженных матриц уже довольно хорошо дает увидеть области разных значений.
delsq( lshape(9) )
И для матрицы 32*32 (выведем часть этой матрицы):
G = lshape( 32 );
D = delsq( G );
spy( D .!= 0, colorbar=false, markersize=2, markershpe=:square, markercolor=1 )
Сколько точек всего в этой функции (грубо говоря, сколько расчетов она задает):
N = sum( G[:] .> 0 )
Решаем PDE с граничным условием¶
Теперь мы готовы решить приведенное выше дифференциальное уравнение с граничным условием Дирихле (1 если точка внутри области, 0 если она на границе). Julia легко позволяет решать линейные системы уравнений с коэффициентами в виде разреженных матриц (при помощи библиотеки UMFPack – той же, что и MATLAB).
N = 60 # Зададим правую часть уравнения
L = lshape(N)
A = -delsq(L)
b = ones( size(A,1) ) / N^2
u = A \ b; # Решаем систему линейных уравнений
Выведем контурный график этого решения:
U = zeros( size(L) )
U[findall(!iszero, L)] = u
contour( U, levels=8, xlimits=(0,N), ylimits=(0,N),
yflip=true, c=cgrad(:turbo, rev=true),
cbar=false, aspect_ratio=:equal )
Теперь выведем решение в качестве трехмерной сетки:
wireframe( U[1:3:end,1:3:end],
zflip=true, camera=(240,30),
fmt=:png, aspect_ratio=:equal )
Заключение¶
Мы решили несложное уравнение на Г-образной расчетной сетке, которую задали при помощи написанной нами функции. Кроме того, поскольку мы использовали разреженные матрицы, наша задача должна решаться быстрее и занимать меньше памяти, чем при использовании плотных (обычных) матриц, таким образом мы, теоретически, получаем возможность решать многократно более крупные задачи.