实现拉普拉斯方程求解的有限差分方案¶

在本例中，我们将求解以下拉普拉斯方程

$$\begin{align*} - \Delta u &= 1, \quad x \in \Omega \\ u &= 0, \quad x \in \partial \Omega. \end{align*}$$

使用 5 点有限差分方案。

网格创建¶

假设我们有兴趣在平面上以 L 形点群形式给出的特定网格上解决问题。

Pkg.add(["LinearAlgebra"])

gr()

Plots.GRBackend()

首先，我们创建一个简单的零点矩阵，然后用单位填充所需的区域，并添加 1 个元素的 "边界"，以设置 Dirichlet 边界条件。图形集gr() 中的函数spy 非常适合用来显示矩阵元素值不为零的区域。

10÷2

5

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 )

如果我们需要对所有这些点进行编号，就像 MATLAB 平台的函数numgrid 所做的那样，我们将使用以下代码：

g = lshape( 10 ) .> 10

10×10 BitMatrix:
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  1  0
 0  0  0  0  1  1  1  1  1  0
 0  0  0  1  1  1  1  1  1  0
 0  0  0  1  1  1  1  1  1  0
 0  0  0  0  0  1  1  1  1  0
 0  0  0  0  0  1  1  1  1  0
 0  0  0  0  0  1  1  1  1  0
 0  0  0  0  0  1  1  1  1  0
 0  0  0  0  0  0  0  0  0  0

g = lshape( 10 )

10×10 Matrix{Int64}:
 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

您可以看到输入矩阵中的哪些元素编号应参与哪些结果的计算。在 5 元素方案中，元素 1、2 和 5（相邻元素）参与计算元素 1。

拉普拉斯算子¶

现在，让我们加载一个函数，在任意网格上定义拉普拉斯算子，使用 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) )

37×37 SparseMatrixCSC{Float64, Int64} with 157 stored entries:
⎡⡻⢎⡢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠈⠪⡛⣬⡢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠪⡱⣮⡀⠀⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠻⣦⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠢⡀⠀⠻⣦⡀⠈⠢⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠢⡀⠈⠻⢆⡀⠈⠢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠈⠻⣦⡀⠈⠢⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠈⠛⣤⡀⠈⠂⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠈⠻⣦⡀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠀⠀⠈⠁⎦

对于一个 32*32 矩阵（让我们输出该矩阵的一部分），我们可以输出以下结果

G = lshape( 32 );
D = delsq( G );
spy( D .!= 0, colorbar=false, markersize=2, markershpe=:square, markercolor=1 )

该函数中总共有多少个点（粗略地说，它设置了多少次计算）：

N = sum( G[:] .> 0 )

675

求解带边界条件的 PDE¶

现在我们准备求解上述带有 Dirichlet 边界条件的微分方程（如果点在区域内则为 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 )

结论¶

我们在一个 L 形计算网格上求解了一个简单方程，我们使用自己编写的函数定义了这个网格。此外，由于我们使用了稀疏矩阵，因此与密集（规则）矩阵相比，我们解决问题的速度更快，占用的内存更少，因此理论上我们可以解决比密集矩阵大很多倍的问题。