Implementation of a finite-difference scheme for solving the Laplace equation

In this example, we will solve the following Laplace equation

using a 5-point finite difference solution scheme.

Creating a grid

Suppose we are interested in solving a problem on a specific grid defined in the form of an L-shaped group of points on a plane.

Pkg.add(["LinearAlgebra"])

gr()

Plots.GRBackend()

First, we will create a simple matrix of zeros, then fill in the necessary areas with ones, also adding a 1-element "border" to set Dirichlet boundary conditions. Function spy from a set of graphs gr() It is well suited for visualizing areas where the values of matrix elements are nonzero.

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 )

If we need to number all these points, how does the function do it? numgrid for the MATLAB platform, use the following code:

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

You can see which element numbers of the input matrix should be involved in calculating which results. With a 5-element scheme, elements 1, 2, and 5 (adjacent elements) participate in the calculations of element 1.

The Laplace operator

Now let's load a function that sets the Laplace operator on an arbitrary grid, using a 4-connected domain for calculations.

using SparseArrays, LinearAlgebra
include( "delsq.jl" ); # Создание лапласиана

Let's use the function spy to better see the individual markers.

L = lshape(9);
G = delsq(L);
spy( G .!= 0, colorbar=false, markersize=3, markercolor=1 )

Although the standard output for sparse matrices is already quite good at showing areas of different values.

delsq( lshape(9) )

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

And for a 32*32 matrix (we will output a part of this matrix):

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

How many points are there in total in this function (roughly speaking, how many calculations it sets):

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

675

Solving the PDE with a boundary condition

Now we are ready to solve the above differential equation with the Dirichlet boundary condition (1 if the point is inside the domain, 0 if it is on the boundary). Julia easily allows you to solve linear systems of equations with coefficients in the form of sparse matrices (using the [UMFPack] library (http://faculty.cse.tamu.edu/davis/suitesparse.html ) – the same as MATLAB).

N = 60                   # Зададим правую часть уравнения
L = lshape(N)
A = -delsq(L)
b = ones( size(A,1) ) / N^2

u = A \ b;          # Решаем систему линейных уравнений

Let's draw a contour graph of this solution.:

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 )

Now let's output the solution as a three-dimensional grid.:

wireframe( U[1:3:end,1:3:end],
           zflip=true, camera=(240,30),
           fmt=:png, aspect_ratio=:equal )

Conclusion

We solved a simple equation on an L-shaped calculation grid, which was set using a function we wrote. In addition, since we used sparse matrices, our task should be solved faster and take up less memory than when using dense (conventional) matrices, thus, theoretically, we get the opportunity to solve many times larger tasks.