Realisation of the finite-difference scheme for the solution of the Laplace equation¶
In this example we will solve the following Laplace equation
$$\begin{align*} - \Delta u &= 1, \quad x \in \Omega \\ u &= 0, \quad x \in \partial \Omega. \end{align*}$$
using a 5-point finite-difference scheme.
Mesh creation¶
Suppose we are interested in solving the problem on some particular grid given in the form of an L-shaped group of points in the plane.
Pkg.add(["LinearAlgebra"])
gr()
First we create a simple matrix of zeros, then we fill the required regions with units, adding also a "border" of 1 element to set the Dirichlet boundary conditions. The function spy from the graph set gr() is well suited for visualising areas where the values of matrix elements are non-zero.
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 )
If we need to number all these points, as the function numgrid of MATLAB platform does, we will use the following code:
g = lshape( 10 ) .> 10
g = lshape( 10 )
You can see which element numbers of the input matrix should participate in the computation of which results. In a 5-element scheme, elements 1, 2 and 5 (neighbouring elements) participate in the calculation of element 1.
Laplace operator¶
Now let us load a function that defines the Laplace operator on an arbitrary grid, using a 4-connected domain for the computation.
using SparseArrays, LinearAlgebra
include( "delsq.jl" ); # Создание лапласиана
Let's use the function spy to see the individual markers better.
L = lshape(9);
G = delsq(L);
spy( G .!= 0, colorbar=false, markersize=3, markercolor=1 )
Although the standard output for sparse matrices already gives quite a good view of areas of different values.
delsq( lshape(9) )
And for a 32*32 matrix (let's 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 total in this function (roughly speaking, how many calculations it sets):
N = sum( G[:] .> 0 )
Solving the PDE with boundary condition¶
Now we are ready to solve the above differential equation with Dirichlet boundary condition (1 if the point is inside the region, 0 if it is on the boundary). Julia easily allows you to solve linear systems of equations with coefficients as sparse matrices (using the library UMFPack - the same as MATLAB).
N = 60 # Зададим правую часть уравнения
L = lshape(N)
A = -delsq(L)
b = ones( size(A,1) ) / N^2
u = A \ b; # Решаем систему линейных уравнений
Let's derive the contour plot 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 us output the solution as a three-dimensional mesh:
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 computational grid, which we defined using a function we wrote. In addition, since we used sparse matrices, our problem should solve faster and take less memory than with dense (regular) matrices, so we theoretically get the ability to solve many times larger problems.