CC BY
6
Discrete & Continuous Models & Applied Computational Science 2023, 31 (2) 150-163
ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph
Research article
UDC 519.6:004.94 PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2023-31-2-150-163 EDN: WFZCIO
Chebyshev collocation method for solving second order ODEs using integration matrices
Konstantin P. Lovetskiy1, Dmitry S. Kulyabov1,2, Leonid A. Sevastianov1,2, Stepan V. Sergeev1
1 RUDN University, 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation 2 Joint Institute for Nuclear Research, 6, Joliot-Curie St., Dubna, Moscow Region, 141980, Russian Federation
(received: April 5, 2023; revised: May 5, 2023; accepted: June 26, 2023)
Abstract. The spectral collocation method for solving two-point boundary value problems for second order differential equations is implemented, based on representing the solution as an expansion in Chebyshev polynomials. The approach allows a stable calculation of both the spectral representation of the solution and its pointwise representation on any required grid in the definition domain of the equation and additional conditions of the multipoint problem. For the effective construction of SLAE, the solution of which gives the desired coefficients, the Chebyshev matrices of spectral integration are actively used. The proposed algorithms have a high accuracy for moderate-dimension systems of linear algebraic equations. The matrix of the system remains well-conditioned and, with an increase in the number of collocation points, allows finding solutions with ever-increasing accuracy.
Key words and phrases: ordinary differential equation, spectral methods, two-point boundary value problems
1. Introduction
Ordinary differential equations (ODEs) and systems of ODEs of the second order describe most problems in classical mechanics. Most oscillatory processes are described by second order ODEs or systems of ODEs. Second order ODE systems describe a number of optical diffraction problems (see, for example, [1]). The model of adiabatic guided wave propagation of polarized light in integrated optical waveguides is also described by a system of two coupled oscillators [2-4].
There are many different methods for exact and approximate solution of initial/boundary value problems for different classes of second order ordinary
© Lovetskiy K.P., Kulyabov D.S., Sevastianov L.A., Sergeev S.V., 2023
This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4.0/legalcode
differential equations. Among them, the spectral methods of expansion in Chebyshev polynomials consistently occupy a well-deserved place.
In 1991, L. Greengard [5] formulated a method for solving a two-point boundary value problem for second order ODEs with constant coefficients, based on expanding the solution into a series of Chebyshev polynomials of the first kind. The method became stably referred to as the pseudospectral collocation method. In the same paper, mathematical constructions were introduced, which later received the names "differentiation matrix" and "integration matrix" (or "antidifferentiation matrix"). A detailed description of the properties of matrices that determine the relationship between the expansion coefficients in a series of approximated functions and the expansion coefficients of their derivatives and antiderivatives in the same set of basis functions is given in [6]. Greengard obtained estimates for the norms of these matrices and their condition numbers — large values for differentiation matrices and small values for integration (antidifferentiation) matrices.
Despite the poor conditionality of differentiation matrices, many authors used them to solve initial and boundary problems for ODEs of various orders. This is explained by the more familiar and therefore 'convenient' representation of physical models using the language of mathematical formulas.
The instability of widely used [7, 8] algorithms has been overcome by applying methods of preconditioning to the corresponding systems of linear algebraic equations. As a result of numerous studies, methods based on integration matrices in the physical space and in the spectral representation turned out to be the most preferable [9].
It is important to note that none of the applied methods for solving ODEs based on Chebyshev integration matrices [9, 10] allows obtaining systems of linear equations with sparse matrices [5]. The dense filling of matrices is a consequence of attempts to introduce boundary conditions into the system of linear algebraic equations along with differential relations [11]. The high sparseness of the matrices can be maintained by improving the algorithm by switching to the two-stage method. In this case, at the first stage, differential conditions are considered, which allow fixing the leading coefficients in the expansion of the solution into a series, thus defining the 'general solution'. The next step uses boundary/initial conditions to determine a pair (for second order equations) of missing coefficients. This makes it possible to obtain a complete set of expansion coefficients for the desired 'particular' solution.
The results of studies [5] demonstrate that the method of Chebyshev collocation that ensures the best accuracy in solving initial and boundary value problems is the method using Chebyshev integration matrices in the spectral space. This approach effectively relies on the use of operations with sparse matrices and its computational costs are quite comparable with the Fourier spectral discretization.
2. Setting of the problem
We consider an approximate solution to the two-point boundary value problem for the second-order differential equation having the form [12]
y"(x) + p(x)y'(x) + q(x)y(x) = r(x), x G (-1,1), (1)
where p(x), q(x), r(x) are sufficiently regular functions. The uniqueness of the solution for any a, ft is ensured by the boundary conditions
a0y(-1) - a-y'(-1) = a, y(1) + ft-y'(1) = ft, (2)
the constants a0 ,a1,ft0 ,ft1 being nonnegative. For example, the condition of continuous p(x) and q(x), positive q(x) >0, x G [-1,1], and nonzero a0 + ax ^ 0, a0 + ft0 ^ 0, ft0 + ft1 ^ 0 ensures the existence of the problem (1)-(2) [13].
3. Methods
The basic idea of spectral methods is to present the solution as a truncated series in known basis functions. The linear transformation (differentiation operator) that transforms the vector of coefficients a = {ak}k>0 of the function expansion f(x) = ^k>0 ak^k (x) into the vector of coefficients b = {bk}k>0 of its derivative expansion f' (x) = ^k>0 bk0k (x) into an analogous series in
the same basis functions is known as the spectral differentiation matrix. The most widespread is the use of bases of Chebyshev functions of the first kind or Lagrange functions, which is due to high interpolative properties of these functions.
Approximation by a finite series (on accuracy when discarding terms of the series with n > N). The expansion of function f G Cn [-1,1] (n times differentiate function) in Chebyshev polynomials Tk (x) • Tk (cos 9) = cos(kd), is determined by the relation
g(x) = 1a0 T(j(x) + a1T1 (x) + ... + an Tn (x) + ..., xG[-1,1], (3)
where
^k = ~ I f(x)Tk(x)(1 - x2)-1/2dx. (4)
n J-1
The residue of truncation of the series (3) to N terms
gN (x) = 2a,0 T(s(x) + ai T- (x) + ... + aNTN (x), xG[-1,1], (5)
has an order of O (N) at N ^ ^ and at f G C^ [-1,1] it tends to zero superalgebraically [14, 15].
Remark 1. According to Eq. (4), coefficients &k are the coefficients of Fourier cosine transformation, so that all N coefficients ak can be obtained by the fast Fourier cosine transformation. And using the inverse Fourier cosine transformation, it is possible to simply calculate gn(cos Qrj) on a grid uniform in 0 G [0,*-].
Most often, the approximation of continuous functions is restricted to a certain fixed number n of the Chebyshev series, as a result of discarding
the components with such Tk (x), k > n, the magnitude of which is small [16, 17]. In contrast to the approximations obtained using other power series, the approximation using the Chebyshev polynomials minimizes the number of terms necessary to approximate the function by polynomials with a given accuracy. Related to this is also the property that the approximation based on the Chebyshev series turns out to be quite close to the best uniform approximation (among polynomials of the same degree), but it is easier to calculate. In addition, it allows you to get rid of the Gibbs effect with a reasonable choice of interpolation points.
The differentiation matrices in the implicit or explicit form are presented in many publications related to the use of pseudospectral collocation methods [68]. The ODE solution using nondegenerate differentiation matrices in the (N + 1)-dimensional physical and/or spectral space quite naturally led to poor conditioned systems of the linear algebraic equations to be solved. Refs. [5,6, 18-20] formulate the specific features of the differentiation and integration matrices, considered on similar or mutually dependent grids. Using explicitly the differentiation matrices on the Chebyshev-Gauss-Lobatto to solve ODEs allows proposing stable and economic methods for solving ODEs. WE use the integration matrices n Chebyshev-Gauss-Lobatto grids in the spectral representation. For more details on the form and properties of these matrices, see [6, 18-20].
3.1. The algorithm based on using integration matrices
Note first that the derivative of Tk (x) can be explicitly written as an expansion in Chebyshev polynomials T0,T1,... ,Tk-1 of lower order [6, 21] as a sum
(]T (r) (
k ( 1 -k[-[kodd]T0 ( x) + 2 £ Tk-1-2J (x)\, xe[-1,1], (6)
3=0
where the notation x\ means the largest integer less than x, and the expression [ kodd] takes the value equal to 1 when k is odd, and 0 when k is even.
We represent the desired function ( ), the future approximate solution of equation (1), as an expansion of the form (3),(4),(5) in a finite set of Chebyshev polynomials T0, ^,..., Tn
- n • n
y( x) = Y: akTk ), ze[-M]. (7)
k=0
By differentiating (7), it is possible to present the first derivative as a series:
n
y' (X ) = £ GkT'k (X), xe[-1,1]. (8)
k=0
At the same time, the derivative y' (x) as a polynomial of degree n can be expanded is series with respect to the initial basis T0, ^,..., Tn with coefficients b = {b0, b1,..., bn}:
y' (x ) = ^ hTk (x ), xe[-1,1]
(9)
the last expansion coefficient becoming zero, bn = 0, in accordance with formula (6) of a transition to the expansion in lower-order polynomials.
Therefore, Eq. (6) describes the relation between the expansion coefficients a = {a0, a1,..., an} of a Chebyshev polynomial of the first kind and the expansion coefficients of its derivative. In matrix form, this relation can be represented using the differentiation matrix: b = Da, where the infinite matrix of Chebyshev differentiation has the form:
A similar transformation for the coefficients of the second derivative
allows using the formula c = DDa to calculate the expansion coefficients c = {c0, cy, cn} in the matrix form.
If an algorithm is needed to determine a part of the coefficients a = {a0, a1,..., an} of the expansion of function y(x) from the known coefficients b = {b0, b1,..., bn} of its derivative expansion, the appropriate matrix form for this operation is a = D+b, where the infinite tridiagonal matrix of integration (antidifferentiation) has the form [6, 18]:
D
D Chebyshev
0 1 0 3 0 5 0 7
0 4 0 8 0 12 0
0 6 0 10 0 14
0 8 0 12 0
0 10 0 14
0 12 0
0 14 0
(10)
(11)
00
10 1
4
1 2
D+ = D+
0 1 6
1
4
(12)
Chebyshev
0 1
8
1 ~6 0 1 10
For example, using the spectral integration matrices D+ to determine coefficients a = |a0 } of function y(x) for the known coefficients
c = {c0,c1 ,...,cn} of the expansion of its second derivative y"(x) allows calculating all coefficients of the function expansion by formula a = D+D+c, except the first two coefficients. This is because the first row of matrix D+ is zero.
Multiplying from the left the integration matrix D+ by vector b = {b0,b1,... ,bn} of the known coefficients of the derivative expansion allows revealing [6] the following dependence of coefficients a = {a0,a1,... ,an} on b = {b0, b1,...,bn}, which can be written in the explicit form:
D+b =
0
0 1 4
0
00 00 00 00
1 -2 0 1
6 0
0 0 0 0
-- 0
0 0
0
0
0
0
1/2
(n-2)
0 0
0 0
0
0
0
-1/2 (n-3)
0
1/2
(n-1)
0
-1/2 (n-2)
0
0 1/2
0 0
0
0
0
0
0
-1/2 (n-1)
0
" ^ " a0
b1 a1
b2 a2
b3 (13
b4 = a4
bn-3 an-3
b-n-2 an-2
bn-1 an-1
- bn . - an -
All the above, including relations (6) and (7),(9),(11) expressed through the representations a = D+D+c and b = D+c, allows us to write equation (1) in spectral representation in the following matrix form:
Tc + diag(p)TD+c + diag(q)TD+D+c = r, x e (-1,1). (13)
Here T is the Chebyshev matrix of mapping a point (vector) from the space of coefficients to the space of function values [5]
10,0 ^
1,0
Tr
2,0
T
n,0
TQ,1 T,
1,1
T,
2,1
T
n,1
T0.2 ^
1,2
Tr
2,2
T
n,2
1,n
Tr
2,n
T
C0 P0
C1 P1
C22 P2
-cn - -Pn -
(14)
so that p = Tc is the vector of values of the desired function (also in the physical space). Here, to reduce formulas, we use the notation Tkj = Tk(Xj), k,j = 0,...,n.
The system of linear algebraic equations (13) has a well-conditioned matrix [5] for any number of collocation points. We will use the Chebyshev-Gauss-Lobatto grid [7, 8], which has proven itself well in the Chebyshev
X
pseudospectral collocation method [9]. Since the matrix is not a symmetric real matrix, instead of the very convenient Cholesky method, we will use the widely used LU method to solve system (13).
The solution of the system of linear algebraic equations (13) is the vector of expansion coefficients {c0, cy,..., cn} in the (n + 1)-dimensional space of the second derivative of the desired solution of equation (1). These components determine the set of 'general' solutions to the ordinary differential equation (1). To single out some specific 'particular' solution from this set, it is required to impose additional restrictions on the components {a0, ay}, which cannot be determined from the relation a = D+D+c.
The first two components that have not yet been found will have to be additionally determined (to obtain a 'particular' solution) from the boundary conditions (2). The remaining components of the vector a remain unchanged and allow satisfying equation (1) for any first expansion coefficients in terms of basis polynomials.
The solution of equation (13) gives us the vector of coefficients {c0, cy,..., cn} of the expansion of the second derivative of the solution of Eq. (1) in Chebyshev polynomials. Thus, the main problem is reduced to solving the simplest Poisson equation:
y"(x) = f(x), -1 < x <1, (15)
where the function f(x) is calculated at any point of the interval x G based on the known vector of coefficients {c0, cy,..., cn}.
n
f( x ) = Y, ckTk (x), xe[-1,1]. (16)
k=0
The method under consideration makes it possible to solve, depending on the type of additional conditions, both the Cauchy problem with initial conditions and the problem with boundary conditions of a general form, requiring, for example, the use of the iterative shooting method [22]. The boundary conditions of the original problem (2) allow extending the definition of the spectral coefficients of the solution. Let us consider some variants, such, e.g., as the Dirichlet conditions at both ends of the interval
y(-1) = a, y(1)=f3. (17)
Neumann-Dirichlet conditions
y' (-1) = a, y(1) = ¡3 (18) or Dirichlet-Neumann condition
y' (-1) = a, y' (1)=f3. (19)
The algorithm for finding a solution to the simplest Poisson equation (15) with one of the boundary conditions (17),(18),(19) consists of three stages:
— calculation of the coefficients of polynomial interpolation of the vector f(x) in the right-hand side of Eq. (15) on the Gauss-Lobatto grid; an efficient method is presented in [23];
— calculation of those coefficients of the solution (except for the first two), which are determined by the differential conditions (15) of the problem (the solution must satisfy the differential conditions), by multiplying the transposed Chebyshev matrix by the vector of interpolation coefficients of the function f(x);
— redefinition of solution coefficients based on boundary (or other independent additional) conditions (17),(18),(19).
In the case of Dirichlet boundary conditions (boundary conditions of the first kind): p(-1) = a, p(1) = (3, the determination of the still unknown coefficients clq ,cli is reduced to solving a system of two equations, which can be, e.g., the equations, which determine the behavior of the solution at the boundary points x = ±1:
n
a0 + a1T10(-1) + Y akTk,0(-1) = a,
k=2 (20)
a0 + aiTi,n (1) +Y akTk,n(1) = ^
k=2
If we additionally consider the fact that Chebyshev polynomials of the first kind take the values Tkj(±1) = ±1, j,k = 0,1,... at the boundary of the interval, then the solution can be written explicitly
a0 = 1[a + (3- V ak L ai = ^ I ¡3- a - V ak ). (21)
Y ak ) , ai = 2lf3-a- ^ ak
k=2,k even J \ k=2,k odd J
In the case when the boundary conditions contain expressions of higher degrees of derivatives of the desired function, one can use the relation [7]
dpTn dx'P
p-i
x=±1
= (±iya+p n
k=0
n2 — k2 2 k+1
(22)
For example, in the case of mixed Neumann-Dirichlet conditions (boundary conditions of the second and first kind): p'(—1) = a, p(1) = /3, the coefficients c0, c1 are determined by the formulas:
n n
aT = a —
k=2 k=2
a — ^(—1)k+1k2ak, ao = (3 — a,1 — ^ ak (23)
and in the case of Dirichlet-Neumann conditions
n n
= —
k=2 k=2
P — Y, k2ak, a0 = a — a1 —Y( — 1)k+1k2ak. (24)
4. Solution of model examples
To illustrate the capabilities of the proposed algorithm, consider as an example the solution of the following second order ODE with Dirichlet boundary conditions y(-1) = sin 1, y(1) = sin 1:
y" + xy' = (2 + x2) cos x, x£(-1,1), y(-l) = sin 1, y(1) = sin 1.
(25)
The exact solution is y(x) = x sin x.
The problem was solved by the collocation method using the integration matrices (see figures 1,2).
Figure 1. Ten collocation points. Solution is plotted in blue, residual - in red
residual
Solution
■1.0 -0.5 0 0.5 1.0
Argument
Figure 2. Five collocation points. Solution is plotted blue, residual - red
Comparison of the exact solution of the model equation with the numerical one is given in the table 1.
Table 1
Comparison of the exact solution of the model equation with the numerical one
Number of collocation points Mean deviation abs(Vexact (X)-Vcalc (%))/N Maximum deviation of the calculated solution from the exact one
6 1.82356474757341e-06 4.63901002387395e-06
7 1.50424878363523e-06 3.0061171892859e-06
9 5.23575446557936e-09 1.05369208025419e-08
11 1.19253244714073e-11 2.39715192140721e-11
13 1.91730425714347e-14 3.86046415624311e-14
14 4.4039495190215e-17 1.11022302462516e-16
The error was estimated numerically (see figure 3). The number of accuracy control points N was taken equal to one hundred.
Figure 3. High solution accuracy with the average and maximum deviations of the numerical solution from the exact one < 10-17 is achieved with a sufficiently small number
of collocation points (n > 13)
As can be seen from the results, the accuracy of the solution depends significantly on the number of collocation points: with an increase in the number of collocation points, the algorithm, in contrast to the method using differentiation matrices, does not lose stability. Due to the inherent property of Chebyshev polynomials, when approximating smooth functions, the accuracy of the solution rapidly increases with a slight increase in the number of basis functions. In our experiment, the most accurate solution was obtained with the number of collocation points equal to 14. With a further
increase in the number of collocation points and, consequently, the number of approximation terms in the expansion series of the solution in Chebyshev polynomials, the accuracy does not increase.
5. Conclusion
In traditional algorithms, even in the most favorable cases when using differentiation matrices on arbitrary grids, the number of arithmetic operations for solving problems with acceptable accuracy turns out to be large. This fact is a consequence of the inclusion in the SLAE, obtained by passing from differential to algebraic relations, of additional equations that specify the initial and boundary conditions.
The algorithm presented in [23] uses a modified (improved) method of pseudo-spectral collocation, i.e., the solution of the problem in two stages. At the first stage, only the 'general' solution of the ODE is found, which is determined by the leading coefficients of the spectral expansion of the solution in the polynomial basis. This approach allows constructing an algorithm that uses only matrices of a simple structure to obtain the solution of the corresponding SLAE. The missing expansion coefficients are determined at the second stage based on the additional (initial or boundary) conditions, solving a simple system of two linear equations.
In this paper, we use an algorithm based on integration matrices. The matrix of the SLAE formed in this case turns out to be well conditioned even for large dimensions of the system. A high accuracy of the solution is achieved with a sufficiently small number of collocation points. The method based on integration matrices should be chosen in cases when there is a request for high and stable accuracy of solving the problem.
Acknowledgments
The publication was supported by the RUDN University Strategic Academic Leadership Program.
References
[1] V. A. Soifer, Diffraction computer optics [Difraktsionnaya komp'yuternaya optika]. M.: FIZMATLIT, 2007, in Russian.
[2] A. A. Egorov and L. A. Sevastianov, "Structure of modes of a smoothly irregular integrated-optical four-layer three-dimensional waveguide," Quantum Electronics, vol. 39, no. 6, pp. 566-574, Jun. 2009. DOI: 10. 1070/QE2009v039n06ABEH013966.
[3] A. L. Sevastianov, "Asymptotic method for constructing a model of adiabatic guided modes of smoothly irregular integrated optical waveguides," Discrete and Continuous Models and Applied Computational Science, vol. 28, no. 3, pp. 252-273, 2020. DOI: 10.22363/2658-46702020-28-3-252-273.
[4] A. L. Sevastianov, "Single-mode propagation of adiabatic guided modes in smoothly irregular integral optical waveguides," Discrete and Continuous Models and Applied Computational Science, vol. 28, no. 4, pp. 361377, 2020. DOI: 10.22363/2658-4670-2020-28-4-361-377.
[5] L. Greengard, "Spectral integration and two-point boundary value problems," SIAM Journal on Numerical Analysis, vol. 28, no. 4, pp. 10711080, 1991. DOI: 10.1137/0728057.
[6] A. Amiraslani, R. M. Corless, and M. Gunasingam, "Differentiation matrices for univariate polynomials," Numerical Algorithms, vol. 83, no. 1, pp. 1-31, 2020. DOI: 10.1007/s11075-019-00668-z.
[7] J. P. Boyd, Chebyshev and Fourier spectral methods, second. Dover Books on Mathematics, 2013.
[8] J. C. Mason and D. C. Handscomb, Chebyshev polynomials. New York: Chapman and Hall/CRC Press, 2002. DOI: 10.1201/9781420036114.
[9] S. E. El-gendi, "Chebyshev solution of differential, integral and integro-differential equations," The Computer Journal, vol. 12, no. 3, pp. 282287, Aug. 1969. DOI: 10.1093/comjnl/12.3.282.
[10] L. N. Trefethen, "Is Gauss quadrature better than Clenshaw-Curtis?" SIAM Review, vol. 50, no. 1, pp. 67-87, 2008. DOI: 10.1137/060659831.
[11] L. A. Sevastianov, K. P. Lovetskiy, and D. S. Kulyabov, "A new approach to the formation of systems of linear algebraic equations for solving ordinary differential equations by the collocation method [Novyy podkhod k formirovaniyu sistem lineynykh algebraicheskikh uravneniy dlya resheniya obyknovennykh differentsial'nykh uravneniy metodom kollokatsiy]," Izvestiya of Saratov University. Mathematics. Mechanics. Informatics, vol. 23, no. 1, pp. 36-47, 2023, in Russian. DOI: 10.18500/ 1816-9791-2023-23-1-36-47.
[12] N. Egidi and P. Maponi, "A spectral method for the solution of boundary value problems," Applied Mathematics and Computation, vol. 409, p. 125 812, 2021. DOI: 10.1016/j.amc.2020.125812.
[13] H. B. Keller, Numerical methods for two-point boundary value problems. Boston: Ginn-Blaisdell, 1968.
[14] D. Gottlieb and S. A. Orszag, Numerical analysis of spectral methods. Philadelphia, PA: Society for Industrial and Applied Mathematics, 1977.
[15] J. F. Epperson, An introduction to numerical methods and analysis, second. John Wiley & Sons, Inc, 2013.
[16] X. Zhang and J. P. Boyd, Asymptotic coefficients and errors for Chebyshev polynomial approximations with weak endpoint singularities: effects of different bases, 2022. DOI: 10.48550/arXiv.2103.11841.
[17] J. P. Boyd and D. H. Gally, "Numerical experiments on the accuracy of the Chebyshev-Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials," Journal of Computational and Applied Mathematics, vol. 205, no. 1, pp. 281-295, 2007. DOI: 10.1016/j.cam.2006.05.006.
[18] B. Fornberg, A practical guide to pseudospectral methods. New York: Cambridge University Press, 1996.
[19] F. Rezaei, M. Hadizadeh, R. Corless, and A. Amiraslani, "Structural analysis of matrix integration operators in polynomial bases," Banach Journal of Mathematical Analysis, vol. 16, no. 1, p. 5, 2022. DOI: 10. 1007/s43037-021-00156-4.
[20] L. C. Young, "Orthogonal collocation revisited," Computer methods in Applied Mechanics and Engineering, vol. 345, pp. 1033-1076, 2019. DOI: 10.1016/j.cma.2018.10.019.
[21] S. Olver and A. Townsend, "A fast and well-conditioned spectral method," SIAM Review, vol. 55, no. 3, pp. 462-489, 2013. DOI: 10 . 1137/120865458.
[22] M. Planitz et al., Numerical recipes: the art of scientific computing, 3rd Edition. New York: Cambridge University Press, 2007.
[23] L. A. Sevastianov, K. P. Lovetskiy, and D. S. Kulyabov, "Multistage collocation pseudo-spectral method for the solution of the first order linear ODE," in 2022 VIII International Conference on Information Technology and Nanotechnology (ITNT), 2022, pp. 1-6. DOI: 10.1109/ ITNT55410.2022.9848731.
For citation:
K. P. Lovetskiy, D. S. Kulyabov, L. A. Sevastianov, S. V. Sergeev, Cheby-shev collocation method for solving second order ODEs using integration matrices, Discrete and Continuous Models and Applied Computational Science 31 (2) (2023) 150-163. DOI: 10.22363/2658-4670-2023-31-2-150-163.
Information about the authors:
Lovetskiy, Konstantin P. — Candidate of Sciences in Physics and Mathematics, Associate Professor of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University) (e-mail: [email protected], phone: +7(495)952-25-72, ORCID: https://orcid.org/0000-0002-3645-1060) Kulyabov, Dmitry S. — Professor, Doctor of Sciences in Physics and Mathematics, Professor at the Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University); Senior Researcher of Laboratory of Information Technologies, Joint Institute for Nuclear Research (e-mail: [email protected], phone: +7(495)952-02-50, ORCID: https://orcid.org/0000-0002-0877-7063) Sevastianov, Leonid A. — Professor, Doctor of Sciences in Physics and Mathematics, Professor at the Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University), Leading Researcher of Bogoliubov Laboratory of Theoretical Physics, Joint Institute for Nuclear Research (e-mail: [email protected], phone: +7(495)952-25-72, ORCID: https://orcid.org/0000-0002-1856-4643)
Sergeev, Stepan V. — PhD student of Department of Applied Probability and Informatics of Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University) (e-mail: [email protected], ORCID: https://orcid.org/0009-0004-1159-4745)
УДК 519.6:004.94 PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2023-31-2-150-163 EDN: WFZCЮ
Метод коллокации Чебышева для решения ОДУ второго порядка с использованием матриц интегрирования
К. П. Ловецкий1, Д. С. Кулябов1,2, Л. А. Севастьянов1,2,
С. В. Сергеев1
1 Российский университет дружбы народов, ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия 2 Объединённый институт ядерных исследований, ул. Жолио-Кюри, д. 6, Дубна, Московская область, 1^1980, Россия
Аннотация. Реализован метод спектральной коллокации для решения двухточечных краевых задач для дифференциальных уравнений второго порядка, основанный на представлении решения в виде разложения по полиномам Чебышева. Подход позволяет устойчиво вычислять как спектральное представление решения, так и его поточечное представление на любой необходимой сетке в области определения уравнения и дополнительных условий многоточечной задачи. Для эффективного построения СЛАУ, решение которой дает искомые коэффициенты, активно используются матрицы Чебышева спектрального интегрирования. Предложенные алгоритмы обладают высокой точностью для систем линейных алгебраических уравнений средней размерности. Матрица системы остается хорошо обусловленной и с увеличением количества точек коллокации позволяет находить решения со все возрастающей точностью.
Ключевые слова: обыкновенное дифференциальное уравнение, спектральные методы, двухточечные краевые задачи
