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

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

采用5点有限差分解方案。

创建网格

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

Pkg.add(["LinearAlgebra"])

gr()

Plots.GRBackend()

首先，我们将创建一个简单的零矩阵，然后用一个填充必要的区域，还添加一个1元素的"边界"来设置Dirichlet边界条件。功能 spy 从一组图 gr() 它非常适合可视化矩阵元素值非零的区域。

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 )

如果我们需要对所有这些点进行编号，函数是如何做到的？ numgrid 对于MATLAB平台，使用以下代码:

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]库（http://faculty.cse.tamu.edu/davis/suitesparse.html ）-与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形计算网格上求解了一个简单的方程，该方程是使用我们编写的函数设置的。此外，由于我们使用稀疏矩阵，我们的任务应该比使用密集（常规）矩阵时更快地解决并且占用更少的内存，因此，理论上，我们有机会解决许多倍大的任务。