CC BY
21
Discrete & Continuous Models
#& Applied Computational Science 2022, 30 (1) 62-78
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Review article
UDC 535:535.3:681.7
DOI: 10.22363/2658-4670-2022-30-1-62-78
Finite-difference methods for solving 1D Poisson problem
Serge Ndayisenga1, Leonid A. Sevastianov1,2, Konstantin P. Lovetskiy1
1 Peoples' Friendship University of Russia (RUDN University) 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation 2 Bogoliubov Laboratory of Theoretical Physics Joint Institute for Nuclear Research 6, Joliot-Curie St., Dubna, Moscow Region, 141980, Russian Federation
(received: January 21, 2022; revised: February 2, 2022; accepted: February 18, 2022)
Abstract. The paper discusses the formulation and analysis of methods for solving the one-dimensional Poisson equation based on finite-difference approximations — an important and very useful tool for the numerical study of differential equations. In fact, this is a classical approximation method based on the expansion of the solution in a Taylor series, based on which the recent progress of theoretical and practical studies allowed increasing the accuracy, stability, and convergence of methods for solving differential equations. Some of the features of this analysis include interesting extensions to classical numerical analysis of initial and boundary value problems. In the first part, a numerical method for solving the one-dimensional Poisson equation is presented, which reduces to solving a system of linear algebraic equations (SLAE) with a banded symmetric positive definite matrix. The well-known tridiagonal matrix algorithm, also known as the Thomas algorithm, is used to solve the SLAEs. The second part presents a solution method based on an analytical representation of the exact inverse matrix of a discretized version of the Poisson equation. Expressions for inverse matrices essentially depend on the types of boundary conditions in the original setting. Variants of inverse matrices for the Poisson equation with different boundary conditions at the ends of the interval under study are presented — the Dirichlet conditions at both ends of the interval, the Dirichlet conditions at one of the ends and Neumann conditions at the other. In all three cases, the coefficients of the inverse matrices are easily found and the algorithm for solving the problem is practically reduced to multiplying the matrix by the vector of the right-hand side.
Key words and phrases: 1D Poisson equation, finite difference method, tridiagonal matrix inversion, Thomas algorithm, Gaussian elimination
1. Introduction
Applied mathematical models are mainly based on the use of partial differential equations [1]. The solution must satisfy a given equation of mathematical physics and some additional relations, which are, first, boundary
© NdayisengaS., Sevastianov L.A., Lovetskiy K. P., 2022
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
and initial conditions. The most important for applications [2] are second-order equations — elliptic, parabolic, and hyperbolic. Currently for equations of mathematical physics, methods of numerical solution and the appropriate software [3], [4], as well as computer algebra systems (CASs) such as Sage, Mathematica, Maxima and Maple are actively developed to implement these methods. Many features of stationary problems of mathematical physics described by elliptic equations of the second order can be illustrated by considering the simplest boundary value problems for an ordinary differential equation of the second order. Perhaps the simplest second-order elliptic equation is the Poisson equation.
Let us consider some methods for the numerical solution of this equation and compare the investigated methods.
The Poisson equation [1] is a special case of the heat conduction equation describing the dependence of the temperature of a medium on spatial coordinates and time, and the heat capacity and thermal conductivity of the medium (in the general case, inhomogeneous) are considered to be given. We will consider the problem of finding the steady-state distribution of density or temperature (e.g., when the distribution of sources does not depend on time). In this case, terms with time derivatives are eliminated from the non-stationary equation and a stationary heat equation is obtained, which belongs to the class of elliptic equations. A two-point boundary value problem is the problem of finding a solution to an ordinary differential equation or second-order systems in the interval a < x < b. Additional conditions are imposed on the solution at any two points of the interval, e.g., a and b — the 'boundaries' of the segment (hence the name of the problem).
Consider a second-order differential equation
It is called the one-dimensional stationary heat conduction equation and arises in the mathematical modeling of many important processes. For example, this equation describes the steady-state temperature distribution u (x) in a heat-conducting rod of length I = b — a. In this case, k (x) is the
du
thermal conductivity coefficient; w (x) = —k (x) — is the heat flux density,
p (x) is the heat transfer coefficient (pu is the heat sink power proportional to the temperature u); f (x) is the density of heat sources (at f ^0 it is the density of heat sinks).
The boundary value problem is much harder to solve than the Cauchy problem, and various approaches are used for this purpose. The most common are various sampling methods that allow replacing the original problem with a certain discrete analog. The resulting discrete boundary value problem is a system of equations (possibly nonlinear) with a finite number of unknowns and can be numerically solved using special direct or iterative methods. One of the simplest discretization algorithms often used in applied scientific and technical calculations is the method of finite differences [5].
The most commonly used method for solving difference equations arising in the approximation of boundary value problems for equations of mathematical physics is the sweep method [6], [7], or the Thomas method [8].
(1)
Below we will show how the difference method is applied to solve the boundary value problem (1), restricting ourselves, for simplicity, to an equation with a constant coefficient k(x) = 1. In this case, the boundary value problem with Dirichlet boundary conditions takes the form
u" (x) — p (x) u (x) — f (x), a < x < b, (2)
u (a) — a, u (b) — fi. (3)
Introduce on [a, b] a grid b, which for
simplicity is assumed uniform. Let us approximately express the second derivative of the solution in terms of the values of the future solution at the grid nodes un — u(xn). We use the simplest symmetric difference approximation
u" (xn) « {un_! — 2un + un+1 ), h — xn+1 — xn — const. (4)
Using such an approximation at each internal grid node xn, 1 < n < N — 1 and substituting it into the differential equation (2), we transform the differential equation (1) into a system of finite-difference equations, i.e., into a system of approximate linear algebraic equations, the solution of which will be an approximate solution yn « u (xn). Finite-difference equations cannot be written at the boundary nodes n — 0, n — N, otherwise the indices of the nodes will go beyond the permissible limits [5]. Denoting pn — p(xn) and fn — f(xn), we get a system of (N — 1) linear equations with respect to the approximate values of the solution at grid nodes
yn-1 — (2 + h2 Pn)Vn + Vn+1 — h2 la, 1^n^N—1. (5)
The number of unknowns yn, 0 < n < N equals (N + 1), i.e., it is greater than the number of equations (5). The lacking two equations are to be obtained from the boundary conditions (3)
Vo —a, Vn — (6)
Solving the algebraic system (5), (6) we get an approximate solution of the boundary value problem (2), (3).
Further analysis of the described algorithm of solving the boundary value problem is to answer three important questions.
— What are the conditions for the existence of a solution to the system of algebraic equations?
— Does the solution of the system of algebraic equations tend to the exact solution of the boundary value problem upon reducing the grid step?
— Is it possible to develop an algorithm (procedure) for finding the solution with given accuracy by reducing the grid step?
It is known [5, P. 66], that for a rather wide class of the boundary value problem coefficients it is possible to prove the existence of a finite-difference solution and its convergence to the exact solution. The following theorem takes place.
Theorem 1. Let p(x), f(x) are twice continuously differentiate on [a,b], p (x) > m, where the constant m^ 0. Also let the step h be small enough, so that h < 2. Then the finite-difference solution exists, its difference from the exact solution by the norm c being of the order of O (h2).
Remark 1. The matrix of the system (5), (6) is tridiagonal. It is not difficult to solve the system by the Gaussian method for a strip matrix or by sweep method. These are direct methods. They allow finding a solution, executing about nine arithmetic operations for each node. By virtue of the conditions of the theorem, the solution of the system of equations by the sweep method exists, is unique and found without accumulating round-off errors.
Remark 2. The conditions of the theorem are sufficient, but not necessary. Even if the conditions are not met, in most cases the finite-difference solution exists and converges to the exact one. Under additional assumptions, it is possible to construct an asymptotically accurate estimate of the error. Then it is possible to apply the grid refinement and Richardson's method to find the posterior estimate of the error and calculations with control of the accuracy.
2. Finite-difference scheme
The problem in matrix form can be represented as
2 + Pi 1 0 - 0 \ (Ul ) (hh2 fi -u
1 -2 + P2 1 - 0 U2 h2 f2
0 1 -2+p3 - 0 u3 h2 f3
0 0 1 - 0 u4 = h2 h
0
0 0 0 • 1
V 0
0
0
1 -2+Pn J
\uN J
(7)
\h2fN -ub J
When applying the sweep method to systems of the form (7), during a forward sweep, both the coefficients of the matrix and the elements of the vector on the right-hand side are recalculated. The matrix is thus reduced to two-diagonal form. During the backward sweep, the components of the solution are calculated at the second stage. Tridiagonal matrices, which are inverted using the simple sweep method, often arise when solving differential equations of two independent variables by the finite-difference method, e.g., when solving a linear one-dimensional heat equation.
For such systems, the solution can be obtained in operations instead of required by the Gaussian elimination method. The first sweep of the method calculates the sweep coefficients, based on which the inverse substitution yields the solution. Examples of such matrices usually arise from discretization of the one-dimensional Poisson equation and interpolation by the natural cubic spline.
a
For the simplest one-dimensional Poisson equation in the case when p(x) = 0, the authors of Refs. [9], [10] proposed a solution based on the analytical (exact) representation of the inverse matrix coefficients.
3. The exact formulation of the inverse of the tridiagonal matrix for solving the 1D Poisson equation with the finite difference method
Consider a method for solving the one-dimensional Poisson equation using the finite difference method based on exact formulas for the inverse of the Laplacian tridiagonal matrix. In the method proposed in Ref. [11], formulas for the coefficients of the inverse matrix are directly derived. Thus, the procedure of solving the one-dimensional Poisson equation becomes very accurate and very fast. This method is a very important tool for solving many physical and technical problems, where the Poisson equation often appears when describing (modeling) various physical phenomena.
3.1. The finite difference method for solving the Poisson equation with Dirichlet—Dirichlet boundary conditions
Consider a function u (x), that satisfies the Poisson equation u" (x) = f(x) on the interval ]a,b[, where f(x) is a given function. We require that the function u (x) satisfy the Dirichlet-Dirichlet boundary conditions: u (a) = a, u (b) = ft. On the considered interval [a, b] we specify a one-dimensional grid xi = a + i ■ Ax, i = 0,..., N + 1, where the uniform step of the grid is
b — Qj
calculated as Ax = ^ ^ = h. We denote by ui = u (xi) and fi = f(xi),
i = 0,..., N + 1 the values of the approximate solution and the function in the right-hand side.
Replacing the second derivative by symmetric difference expressions, we obtain the following system for internal nodes:
u%_1 — 2ui + u%+1 = h2fi, i = 1,..., N. (8)
In matrix form, the system of linear algebraic equations (8), taking into account the boundary conditions, can be written in the form Au = F, where
F = (h2 fi — ua, h2 f2,..., h2fN_i, h2 fN — ub)T, or
(-2 1 0 0 0 ...... 0 ] f Ul ) (h2 fi -ua \
1 -2 1 0 0 ...... 0 U2 h2 Î2
0 1 -2 1 0 ...... 0 u3 h2 f3
0 0 1 -2 1 ••. ••• 0 X u4 h2 f4
0 0 0 1 -2 • • 0 u5 h2 h
0 0 0 0 1 UN-1 h2fN-i
0 0 0 0 0 0 1 -2) V UN ) \h2fN -ub)
Thus, the solution of the one-dimensional Poisson equation is reduced to the inversion of the tridiagonal symmetric negative definite matrix
A=(aZJ), i,j=l,...,N.
The inverse matrix which we denote by
B=(bZJ), i,j=l,...,N,
is also symmetric.
The elements of matrix A may be briefly written as
-2, i= j,
= Jl, l%-3l = l, x = l,...,N (10)
a
li-3l>l,
and the elements of matrix B are related by the following formulas:
-2b%1 +b%2 = 5},
bij-i - 2b%J + blJ+i = Sj, l < j < N, (11)
-biN-1 — 2^iN = $i1,
where is the Kronecker symbol.
3.2. Calculating the inverse matrix Relations (11) allow deriving the following interesting dependencies
btJ+1 =blJ +bt1, btj = Jbz1 + U-l). (12)
From relations (12) it follows that the elements of inverse matrix B are unambiguously determined by the value of the element b11. This coefficient can be determined based on the behavior of matrix B at different dimensionalities N:
b11 = -N/(N+l). (13)
From relations (12) and (13), it is easy to express the elements of the first row and the first column of the inverse matrix
bu = -(N-(J-l))/(N + l), bt1 =-(N-(i-l))/(N+l). ( )
These relations allow completing the accurate and full determination of the coefficients of the inverse matrix B = (b^), i,j = l,..., N:
h ={-j(N-(i-1))/(N+1), i>j, ^ { -l(N-U-1))/(N+1), i<j; (5)
B =
1
'(N+1)
x
f N
N-1
x
N-1 2 (N-1)
V
4 2
N-(j-1) 2[N-(j-1)]
N-(i-1) 2[N-(t-1)] - i[N-(j-1)]
2j j
2
4
2i
1 2
2 (N-1) N-1 N-1 N J
With the inverse matrix elements known, it is easy to get the solution of the one-dimensional Poisson equation by mere multiplication of the matrix by the right-hand side vector u = BF.
3.3. Classification of media
Taking into account the specific form of the inverse matrix and its persym-metry makes it easy to express the solution uN at the point xn
N
UN =
(N + 1)-i £i-Ft.
(16)
i=i
The direct search for the solution _i at the point xN-i leads to the
expression
nN_i =-(N+1)_
N_i
T,2i-Fi
i=i
+ (N-1)Fj
N
(17)
In a similar way, it is possible to derive the expressions for calculating the rest components of the solution in the form
uN_k = -(N+1) 1 x
x
(k+1)
'N_k i=1
+ (N-k)
N
£ (N-(i-1))F,
i=N_k+1
k = 0,1,...,N-1 (18)
or in the form
uk = -(N +1)_1 x
(N-k+1)
x
i=1
+ k
N
£ (N-(i-1))F,
i=k+1
k=1,...,N. (19)
From the computational point of view, it is preferable to use Eqs. (15), when programming the procedure of calculating the solution.
Let us consider the numerical solution of the problem of finding a scalar potential given on the interval [-1,1] and satisfying the Poisson equation d2 $ (x)
A$ (x) = , ; 7 = f(x) = -cos2 (-K (x - 0.5)) and the Dirichlet-Dirichlet ox2
boundary conditions: $ (-1) = -0.2, $ (1) = 0.1. The exact solution is expressed by the formula
x
^exact = — +
cos (k (x — 0.5)) 2k
x
+ - — 0.1(x + 1)+0.3. 4
(20)
The software implementation of the algorithm consists of several lines, namely, filling the vector on the right-hand side of Eq. (9) and multiplying the inverse matrix B by this vector using Eqs. (14).
Figure 1 illustrates the results of the numerical experiment.
-1,00
-0,05 -0,10 -0,15 -0720^
(a) Exact solution
(b) Calculation error
Figure 1. The maximal error at points x = ±0.5 is 0.36 at N = 30 and decreases to 0.036
at N = 300
4. Solving the 1D Poisson equation with the Neumann—Dirichlet and Dirichlet—Neumann boundary conditions
The problem is to determine the scalar potential u (x) satisfying the one-dimensional Poisson equation Au(x) = f(x) on the interval ]a,b[, where f(x) is a given function. It is necessary to find the solution satisfying the Neumann-Dirichlet boundary conditions u' (a) = u'a and u (b) = ub. Let us consider a special uniform grid for the finite difference method with the step
Ax = = h, consisting of N + 1 points. The coordinates of the grid nodes
(x,i) are determined by the expression xi = a + (i - 1) • h, i = 0,1,... N + 1. We denote by ui the approximate values of the desired solution at point
(xi), and by fi the value of the given function in the right-hand side at the same point. In addition, let us denote by u\ = u' (xi) and u'l = u" (xi) the values of the first and second derivatives of the sought solution at the grid node at the same point. Replacing the derivatives with symmetric finite-difference expressions [12], we arrive at the approximation formulas of the second order of accuracy for the first derivatives
u< = Ut+1 -hU%-1 +0(h2), i = 1,2,3,...,N (21)
and for the second derivatives
< = U%-1 ~ 2h2 + U%+1 + 0(h2), i = 1,2,...,N. (22)
The system of linear equations for the internal nodes of the interval looks as
uz-i — 2ui + ul+i = h2fi, i = 1,..., N. (23)
4.1. The Neumann—Dirichlet boundary conditions
Let us derive equations complementing the system with the boundary conditions at the left and right ends of the interval taken into account. Assuming the use of Eqs. (21) and (23) possible and combining them at i = 0, we eliminate u_1 from the system of equations.
,2 fi
—u1 + u2 = h2-^ + hu'a,
t = 1,...,N.
(24)
Thus, introducing into consideration an additional virtual point x0 = a —h allows using the central differences with the order of approximation O (h2) for the sought solution even at the boundary point of the interval.
We introduce the vector F with the components expressed as
Fi = h22 + hu'a, Fn = h2fN — ub, Fi = h2fz, i = 2,3,...,N — 1. (25)
As a result, the system of equations that determines the solution components reduces to the form
-1 1 0 0 0
0 0
1 -2 1 0 0
0 0
0 1 -2 1 0
0 0
0 0 1 -2 1
0 0
0 0 0 1 -2
0 0 1
0 0 0 0
0 1
-2J
x
Un
Un
UA
Uï
UN_1 \ UN
tjr A + hK\
h2 f: h2
h21
h2 h
h2 fN_1 \ h2 fN -ub J
(26)
where the matrix A = {airj}, i,j = 1,..., N of system (8) is symmetric tridiag-onal negative definite and possesses the property of diagonal transformation. The presence of diagonal dominance in the coefficient matrix guarantees the stability of the sweep method; however, in this case, there is a way to calculate the elements of the inverse matrix.
4.2. Calculation of the inverse matrix elements
Let us write down the properties of the inverse matrix B = {birj}, i,j = l,...,N, B = A-1, following directly from its definition. It must be symmetrical and its elements must satisfy the following relations:
f—bij + b2i = 5], 1 < j < N,
J2j
ba — 1bl2 + bi3 = 52, 1<г<N,
bi-ij — 2hij + h+ij = ÔÎ, 1 <i,j <N, biN-i — 2biN = öf, 1<x< N,
(27)
where is the Kronecker symbol.
The elements of the inverse matrix also satisfy the relations
^ =
bii + Ü—1),
bu + (i — 1), i> j.
(28)
The analysis of behavior of the system determinant allows deriving the expressions
det (B) = (—1)
N
bii =
N■(—1)
N-1
(—1)
N
= —N, and bNN =
(—1)
N-1
(—1)
N
= —1=b
(29)
1N.
Using Eqs. (27)-(29), we can exactly determine the elements of the inverse matrix, which is related to the search for the approximate solution in the case of the Neumann-Dirichlet boundary conditions. Thus, the elements of matrix B are determined by the expressions
^ =
— [N — U—1)], i^j,
-[N-(i-l)], t>j. The elements of the inverse matrix can be alternatively expressed as
'(i + j) + \i-j\
= — [N— [max (z,f) — 1]] = —
N —
—1
(30)
(31)
The expressions (30) and (31) are equivalent. However, for software implementation, the first one is preferable.
As a result of the transformations carried out, explicit expressions for the elements of the inverse matrix are obtained, and the solution of the Poisson
problem with Neumann-Dirichlet boundary conditions can be obtained using a simple multiplication of the inverse matrix by the vector of the right-hand side: U = BF, where
B = —
( N N — 1 N — 2 ■■■ •2 A
N — 1 N — 1 N — 2 • •2 1
N — 2 N — 2 N — 2 • •2 1
2 2 2 •2 1
1 1 1 •1 1/
Each solution component can be expressed directly using the formula
uk =
(N — k+1)
i=1
+
N
£ (N-(i-1))-Fi
i=k+1
k = 1,2,..., N. (32)
Formula (32) gives a simple analytical expression for the solution of the Poisson equation with Neumann-Dirichlet boundary conditions. It is very easy to program it either directly or based on Eq. (30). One double loop will be enough to compute the entire solution.
4.3. Example
Consider a numerical solution of the problem of finding a scalar potential defined on the interval [-a, b] and satisfying the Poisson equation
d2 $ (x)
A$ (x) = d^2J = f(x) = Vo cos (kx + <po),
where a, b, V0, k and are given constants, and the Neumann-Dirichlet
d$
boundary conditions —— (a) = and $ (b) = $b.
ax
The known exact solution is expressed as
$
exact
(*) =
V0
$'a — -0 sin (ka + <p0 )
(x — b) —
Vn
— -0 [cos (kx + ) — cos (kb + )] + $b. (33) kz
Let us consider the finite-difference solution at a = ——, b = —, V0 = 1,
k = ïï and ^0 = ^.
We define the computational grid with the following parameters: N=100, Ax = h = = (i-1)Ax, = &(Xi)
and fi = f(xi) = cos (kxi + ). The solution is assumed to satisfy the Neumann-Dirichlet conditions specified as: = 1/4 and = -1/2.
We calculate the solution of the Poisson problem multiplying the inverse matrix with the elements determined by expressions (30) by the right-hand side vector, corrected using Eqs. (26).
The software implementation of the algorithm consists of a few lines: filling the right-hand side vector (26) and multiplying the inverse matrix B by this vector using Eqs. (30).
Figure 2 illustrates the results of the numerical experiment.
-3^-1-L -1,571 -1,071 -0,571 -0,071 0,429 0,929
(a) Exact solution (b) Calculation error
Figure 2. The maximal error at x = -1.571 is 1.55E — 04 for N = 100 and decreases to
1.55^-06 for N = 1000
5. Dirichlet—Neumann boundary conditions 5.1. Discretization and matrix equation
By analogy with the case of the Neumann-Dirichlet boundary conditions, we consider the symmetric case with the Dirichlet-Neumann boundary conditions. Let us first define a suitable sampling grid on the interval [a, b]. Grid points [xi, i = 0,1,..., N + 1} are specified as xi = a + ih. The boundary conditions ua and u'h, complementing the Poisson equation redefine the system of finite-difference equations (23). The solution value uN+1 at the 'virtual' point xN+1 is expressed using the boundary condition for the derivative, approximating the latter by symmetric central differences. As the last equation of the system, we get
UN-1 -uN = h2 ^ -hu'b. (34)
The transformed right-hand side vector F is presented as
Fn = h2^ -hu'b, F1 = h2 f1 — ua, Fi = h2U, t = 2,...,N—1. (35)
Like in the previous case of Neumann-Dirichlet boundary conditions, the resulting matrix of the equation
(—2 1 0 0 0 ...... 0\ i Ul \ (hh2Î1 —Ua \
1—2 1 0 0 ...... 0 U2 h2 Î2
0 1—2 1 0 ...... 0 u3 h2 f3
0 0 1 —2 1 • • 0 X u4 h2 h
0 0 0 1 —2 ■■. • u5 h2 h
0
0 0 0 0 ••. • • 1 UN-1 h2 fN-1
V 0 0 0 0 0 0 1 —1J V un \h2li — hub J
symmetric three-diagonal negative definite, with the dominant main
(36)
nal.
With respect to the antidiagonal, this matrix is symmetric to the matrix used in the solution of the Poisson problem with the Neumann-Dirichlet boundary conditions. The system is definite and has a unique solution for any right-hand side.
Using the antidiagonal symmetry with respect to the Neumann-Dirichlet problem, we construct the inverse matrix for the Dirichlet-Neumann case:
B = -
1111 12 2 2 1 2
12 12 12
1 2
N-2 N — 2 N-2 N-2 N—1 N — 1 N — 2 N—1 N J
Therefore, the exact solution of the system of equations (36) can be written very simply (in a single line)
uk =
Y,'-*
i=1
+ k
n i=k+1
k = 1,2,...,N.
(37)
The software implementation of the method reduces to simple multiplication of the inverse matrix by the right-hand side vector.
Figure 3 illustrates the results of the numerical experiment.
(a) Exact solution (b) Calculation error
Figure 3. The maximal error in this case at point x = 0.785 is 1.69E — 04 for N = 100 and reduces to 1.69E — 6 for N= 1000
5.2. Example
An example of the previous section is considered, which differs only in that the boundary conditions set earlier at the left end of the interval are transferred to the right and vice versa. The software implementation of the algorithm consists of several lines: filling in the vector of the right-hand side (34) and multiplying the inverse matrix B by this vector using Eq. (37).
6. Conclusion
The paper gives examples of practical problems, in the simulation of which it is necessary to solve second-order elliptic equations with different boundary conditions. The case of the one-dimensional Poisson equation and its finite-difference solution are described in detail. Estimates of the complexity of the sweep algorithm in the case of a uniform grid are given. An approach to solving the one-dimensional Poisson equation using explicitly calculated coefficients of inverse matrices for various types of boundary conditions is also described. The Dirichlet and Neumann boundary conditions in various combinations are considered.
A comparative analysis of the computational complexity of methods for solving the one-dimensional Poisson equation, based on the use of the sweep method and methods using an explicit representation of inverse matrices is presented.
Direct calculation shows that to implement calculations by right-sweep formulas, approximately 8N arithmetic operations are required, whereas in the Gauss method for fully filled matrices this number is approximately (2/3) N3. It is also important that the tridiagonal structure of the matrix of the system makes it possible to use for its storage only an array of real variables of dimension 3N — 2.
The assertion of the author of Ref. [4] that the method he proposed using the explicit form of inverse matrices allows solving the Poisson equations with different boundary conditions faster and more accurately is, to put it mildly, incorrect. Provided that the stability conditions of the sweep method are met,
the speed of solving the problem by the sweep (Thomas) method is an order of magnitude higher due to a much smaller number of required operations.
However, unlike the sweep method [13], the practical implementation of the proposed method does not imply the allocation of additional arrays for software implementation, since the elements of the inverse matrix have a very simple form and their calculation within the loop determining the components of the solution is not difficult.
References
[1] A. N. Tikhonov and A. A. Samarskii, Equations of Mathematical Physics [Uravneniya matematicheskoy fizikij, 7th ed. Moscow: Moscow State University, Nauka, 2004, in Russian.
[2] D. A. Yakovlev, V. G. Chigrinov, and H. S. Kwok, Modeling and optimization of LCD optical performance. New York: Wiley, 2015. DOI: 10.1002/9781118706749.
[3] L. N. Trefethen, Approximation theory and approximation practice. Philadelphia: SIAM - Society for Industrial and Applied Mathematics, 2019.
[4] M. Planitz et al., Numerical Recipes: The Art of Scientific Computing, 3rd ed. New York: Cambridge University Press, 2007.
[5] N. N. Kalitkin and P. V. Koryakin, "Numerical methods [Chislennyye
metody],"in Methods of Mathematical Physics, 1st ed. Moscow: Academia, 2013, vol. 2, in Russian.
[6] A. A. Abramov and V. B. Andreyev, "On the application of the method of successive substitution to the determination of periodic solutions of differential and difference equations," USSR Computational Mathematics and Mathematical Physics, vol. 3, no. 2, pp. 498-504, 1963. DOI: 10. 1016/0041-5553(63)90034-x.
[7] A. A. Samarskiy and A. V. Gulin, Numerical Methods of Mathematical Physics [Chislennyye metody matematicheskoy fiziki]. Moscow: Scientific world, 2003, in Russian.
[8] L. H. Thomas, Elliptic problems in linear difference equations over a network. New York: Waston Sci. Comput. Lab. Rept., Columbia University, 1949.
[9] S. B. Gueye, K. Talla, and C. Mbow, "Generalization of the exact solution of 1D Poisson equation with robin boundary conditions, using the finite difference method," Journal of Electromagnetic Analysis and Applications, vol. 6, no. 12, pp. 372-381, 2014. DOI: 10.4236/jemaa.2014.612038.
[10] S. B. Gueye, K. Talla, and C. Mbow, "Solution of 1D Poisson equation with Neumann-Dirichlet and Dirichlet-Neumann boundary conditions, using the finite difference method," Journal of Electromagnetic Analysis and Applications, vol. 6, no. 10, pp. 309-318, 2014. DOI: 10.4236/jemaa. 2014.610031.
[11] S. B. Gueye, "The exact formulation of the inverse of the tridiagonal matrix for solving the 1D Poisson equation with the finite difference method," Journal of Electromagnetic Analysis and Applications, vol. 6, no. 10, pp. 303-308, 2014. DOI: 10.4236/jemaa.2014.610031.
[12] N. N. Kalitkin and E. A. Alshina, "Numerical Methods [Chislennyye metody]," in Numerical analysis. Moscow: Academia, 2013, vol. 1, in Russian.
[13] A. Amosov, Y. Dubinsky, and N. Kopchenova, Computational Methods [Vychislitel'nyye metody], 4th ed. St. Petersburg: Lan', 2021, in Russian.
For citation:
S. Ndayisenga, L. A. Sevastianov, K. P. Lovetskiy, Finite-difference methods for solving 1D Poisson problem, Discrete and Continuous Models and Applied Computational Science 30 (1) (2022) 62-78. DOI: 10.22363/2658-4670-202230-1-62-78.
Information about the authors:
Sevastianov, Leonid A. — Doctor of Physical and Mathematical Sciences, Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University), Leading Researcher of Bogoliubov Laboratory of Theoretical Physics, JINR (e-mail: [email protected], phone: +7(495)9522572, ORCID: https://orcid.org/0000-0002-1856-4643, ResearcherlD: B-8497-2016, Scopus Author ID: 8783969400) Lovetskiy, Konstantin P. — Candidate of Physical and Mathematical Sciences, Associate Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: [email protected], phone: +7(495)9522572, ORCID: https://orcid.org/0000-0002-3645-1060 , ResearcherID: A-5725-2017, Scopus Author ID: 18634692900)
Ndayisenga, Serge — Student of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia (RUDN University) (e-mail: [email protected], phone: +7(977)9086946, ORCID: https://orcid.org/0000-0002-9297-9839, ResearcherID: AAC-3303-2022)
УДК 535:535.3:681.7
DOI: 10.22363/2658-4670-2022-30-1-62-78
Конечно-разностные методы решения 1D задачи
Пуассона
С. Ндайисенга1, Л. А. Севастьянов1,2, К. П. Ловецкий1
1 Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия 2 Лаборатория теоретической физики им. Н. Н. Боголюбова Объединённый институт ядерных исследований ул. Жолио-Кюри, д. 6, Дубна, Московская область, 141980, Россия
Аннотация. В статье обсуждается постановка и анализ методов решения одномерного уравнения Пуассона на основе конечно-разностных аппроксимаций — важного и очень полезного инструмента численного исследования дифференциальных уравнений. По сути, это классический метод аппроксимации, основанный на разложении решения в ряд Тейлора. Развитие теоретических и практических результатов на базе этого метода в последние годы позволили повысить точность, стабильность и сходимость методов решения дифференциальных уравнений. Некоторые особенности этого анализа включают интересные расширения классического численного анализа начальных и граничных задач. В первой части излагается численный метод решения одномерного уравнения Пуассона, сводящийся к решению системы линейных алгебраических уравнений (СЛАУ) с ленточной симметричной положительно определённой матрицей. В качестве метода решения СЛАУ используется широко известный метод прогонки (метод Томаса). Во второй части представлен метод решения, основанный на аналитическом представлении точной обратной матрицы дискретизированного варианта уравнения Пуассона. Выражения для обратных матриц существенно зависят от типов граничных условий в исходной постановке. Представлены варианты обратных матриц для уравнения Пуассона с различными граничными условиями на концах исследуемого интервала — условиями Дирихле на обоих концах интервала, условиями Дирихле на одном из концов и Неймана на другом. Во всех трёх случаях коэффициенты обратных матриц легко вычисляются (выписываются) и алгоритм решения задачи практически сводится к умножению матрицы на вектор правой части.
Ключевые слова: Ш уравнение Пуассона, метод конечных разностей, обращение трехдиагональной матрицы, алгоритм Томаса, исключение Гаусса
