Exact schemes for second-order linear differential equations in self-adjoint cases

When working with mathematical models, to keep the model errors as small as possible, a special system of linear equations is constructed whose solution vector yields accurate discretized values for the exact solution of the second-order linear inhomogeneous ordinary differential equation (ODE). This case involves a 1D spatial variable x with an arbitrary coefficient function κ (x) and an arbitrary source function f (x) at each grid point under Dirichlet or/and Neumann boundary conditions. This novel exact scheme is developed considering the recurrence relations between the variables. Consequently, this scheme is similar to those obtained using the finite difference, finite element, or finite volume methods; however, the proposed scheme provides the exact solution without any error. In particular, the adequate test functions that provide accurate values for the solution of the ODE at arbitrarily located grid points are determined, thereby eliminating the errors originating from discretization and numerical approximation.


Introduction
Consider the following linear second-order inhomogeneous ordinary differential equation (ODE) in self-adjoint form: where κ(x) ≥ κ 0 > 0 is a positive function ensuring the existence of the integrals applied for the solution procedure. Equation (1.1) is a fundamental formula in the mathematics of physical laws. The practical impact of the solution of Eq. (1.1) is profound since many important physical phenomena used in various industrial applications are described by this equation. Problems in this form arise in various fields, including thermal energy transport, diffusion, electrostatics, and electrodynamics. In electrical applications, κ is the conductivity distribution, u is the electric potential, and f is the source term; all these parameters depend on the 1D spatial variable x. It is assumed that the classic solution u(x) exists on [0, ] with the appropriate boundary conditions.
Finding and applying an analytical solution to Eq. (1.1) ensures that the error in the mathematical algorithm used to solve the physical problem is solely numerical. Therefore, our mathematical problem in the case of Eq. (1.1) is finding and applying the inverse of the differential operator to construct an analytical ODE solver. On the basis of our research work, this paper proposes a scheme that provides the values for the exact solution a of Eq. (1.1) in arbitrarily chosen grid points. Although some researchers have attempted to develop exact schemes, practically applicable solutions have been reported only for homogeneous second-order ODEs [1]. Furthermore, some researchers [2] used the nonstandard finite difference method to develop a scheme by using the Green function corresponding to the operator. Based on the fundamental solutions corresponding to differential operators, several schemes applicable to first-order ODEs have been reported [3].
In addition, some higher-order schemes for first-and second-order differential equations were presented in [4]. Furthermore, the authors of [4] reported an exact scheme originating from Samarskii applicable to a special case, namely an equidistant distribution of space. Delkhosh et al. proposed analytical methods for solving a second-order ODE, the effectiveness of which was demonstrated in their published work [5,6]. The method they described uses the beneficial properties of the Bessel equation to solve homogenous equations with a second-order self-adjoint differential operator [5,6]. Based on the state-ofthe-art literature, it is generally concluded that analytical methods based on the current state of science provide solutions for some specific version of Eq. (1.1) (f (x) = 0 or κ(x) = 1).
Well-known numerical methods, such as the finite difference method (FDM) [7], finite volume method [8], and boundary element method [9], can be used to resolve equations in form (1.1) by approximating the analytical solution of the ODE. Among the classic numerical methods for second-order differential equations, at present, the finite element method (FEM) is the most widely used [10]. This method has also been used to develop approaches involving a higher order of convergence, particularly with the application of higher-order elements [10].
Thus, research on exact schemes is an important aspect in the domain of mathematical physics, and increasing the accuracy of such schemes represents a continuous challenge for researchers. It has been reported [11] that there is no unified theoretical foundation for the construction of exact schemes. Furthermore, increasing the accuracy also increases the complexity and computational burden of these schemes, resulting in difficulties during their application. However, the proposed scheme can overcome these limitations as it provides the exact solution at an arbitrary location independent of the spatial discretization.
The subsequent sections present the details pertaining to the mathematical construction of the proposed method. Furthermore, we demonstrate the advantages of the method through several examples.

Local Green functions of ODE for arbitrary partitioning
Consider an arbitrary discretization of the interval [0, ] into (n + 1) subintervals using arbitrarily distributed node points: where n may be 1, 2, 3, . . . . We represent the subintervals as I i-1 = [x i-1 , x i ] for the indexes i = 1, 2, . . . , n, n + 1; thus, We define (2n) nonnegative integral functions of 1 κ(x) on the different subintervals represented as follows: for each i = 1, 2, . . . , n, assuming that the integrals are existing and finite. In the subsequent sections, these functions are used as test functions in a similar sense to the case of the FEM [10].
The following properties are fulfilled for the test functions based on the fundamental theorems of calculus: x 0 = 0, Then the local Green functions can be defined as follows: (2.7) In Fig. 1, the graphs of functions ψ 0 , ψ 1 , and ψ 2 are represented as solid lines, while those of functions ϕ 1 , ϕ 2 , and ϕ 3 are represented as dotted lines on the domains of their definitions.
From the derivatives listed in (2.4) and (2.5), it is clear that test functions ψ i increase monotonously and that test functions ϕ i decrease on the interval of their domain I i for all i = 1, 2, . . . , n -1. On the first interval I 0 , only the test function ψ 0 is considered, while on the last interval I n , we define only the test function ϕ n similar to the FEM [10]. Initially, it appears as though the local Green functions are not all independent because Furthermore, as indicated by (2.8) and demonstrated in the following subsections, the test functions are handled in a similar manner as in the FEM [10].
Remark 1 If κ(x) ≡ 1 is selected as the constant function and the uniform mesh x i = ih (i = 0, 1, . . . , n + 1) is used on the interval [0, 1] with a step size of h = 1 n+1 , the test functions form linear functions resembling a saw-tooth pattern, as shown in Fig. 2, and the shape functions can be calculated as follows:

Fundamental recursive relation: flux elimination process
In this section, we demonstrate that the test functions defined in (2.3a)-(2.3b) are adequate to create a recursive relation among consecutive values of the solution u(x). This approach enables the elimination of the derivative u (x) from the equations. The objective The coefficients of the linear recursive relation are defined using the following integral formula: Using the definition of the test functions in (2.3a)-(2.3b), equivalent formulas for the following coefficients can be defined: The recursive relation can be developed essentially in two steps. An arbitrary pair of consecutive subintervals In the first step, we multiply ODE (1.1) by the test function ψ i-1 (x) and then integrate with respect to x over the interval I i-1 : where the prime symbol ( ) denotes the derivative with respect to x and i = 1, 2, . . . , n.
Applying integration by parts to the left-hand side of (3.3), the derivative of ψ i-1 (x) and the anti-derivative of the remainder can be obtained as follows: The first term becomes zero by substituting the lower limit according to the third identity reported in (2.4); however, by substituting the upper limit, we can obtain the coefficient (3.2). The integrand of the second term can be simplified using the factor κ(x), and based on the fundamental theorems of calculus, the following expression is obtained: By multiplying this equation with the coefficient a i-1 , we obtain the following relation: which is satisfied for all i = 1, 2, . . . , n.
It is desirable to eliminate the first term from (3.5), the so-called flux, because the derivative u (x) is unknown. This elimination is a necessary second step to obtain the recursive formula.
Therefore, we multiply ODE (1.1) by the test function ϕ i (x) and integrate over the interval I i : Applying integration by parts to the left-hand side of (3.6), the derivative of factor ϕ i (x) and the anti-derivative of the remainder can be obtained as follows: (3.7) The first term becomes zero by substituting the upper limit according to the third identity reported in (2.5); furthermore, by substituting the lower limit, we obtain the coefficient 1 a i = ϕ i (x i ), as discussed in (3.2). The integrand of the second term can be simplified using the factor κ(x), and from the fundamental theorems of calculus, we can obtain the following: Multiplying this equation by the coefficient a i , the following relation is obtained: which is satisfied for all i = 1, 2, . . . , n. Now, we can eliminate the flux from Eqs. (3.5) and (3.9). By adding these two equations, the basic recursive relation can be obtained: for all indexes i = 1, 2, . . . , n.

Exact scheme for Dirichlet boundary conditions
We use the basic recurrence relation of Eq. (3.10) to obtain the exact values of solution (1.1) at all (n + 2) node points: However, only n equations are available in (3.10). Therefore, two values must be prescribed arbitrarily to obtain a unique solution. One possible approach to provide two independent values is to apply Dirichlet boundary conditions, in which the values of the solution at the endpoints are obtained as follows: In this case, we substitute u(x 0 ) = α into the first equation and u(x n+1 ) = β into the last equation in (3.10) to obtain the key result as follows.

Theorem 1 Consider a system of n ≥ 3 linear equations:
-a n-1 u n-1 + (a n-1 + a n )u n = a n β + a n-1 G n-1 + a n H n where the coefficients a i can be obtained using integrals

. , n) and the coefficients G i-1 and H i can be obtained using double integrals
(4.5) The matrix-vector format LU = F D of system (4.3) can be defined, in which the discrete Laplacian matrix (L) is as follows: L = a 0 + a 1 -a 1 0 0 -a 1 a 1 + a 2 -a 2 0 -a 2 a 2 + a 3 -a 3 0 -a n-2 a n-2 + a n-1 -a n-1 0 0 a n-1 a n-1 + a n This matrix is symmetric and positive definite and has a tridiagonal shape with dimensions of n × n. The vector on the right-hand side can be defined for the Dirichlet boundary conditions as follows: . . a n-2 G n-2 + a n-1 H n-1 βa n + a n-1 G n-1 + a n H n  We can assemble the entries of the Laplacian matrix based on the tridiagonal form (4.6) by using accurate values of the exponential function e x in (4.9) without approximating with rounded values: By calculating the integrals G i and H i in (4.5) and using the local Green functions from Example 2.1, the following expressions are obtained:   Now, we can collect the entries of the right-hand side vector based on the formula and because the boundary conditions are homogeneous, i.e., α = β = 0. After simplification, we obtain  The solution u(x) is represented in Fig. 3 where the kernel function G(x, t) is the (global) Green function (4.17) We should note that this derivation from local Green functions towards the global Green function is reversible; i.e., we could have started with a solution of the form (4.16) with the Green function (4.17) applied to the partitions I i-1 and I i , and we would still have arrived at the basic recursion relations (3.10).

Exact scheme for Dirichlet and Neumann boundary conditions
Consider ODE (1.1) with a Dirichlet boundary condition at x = 0 and a Neumann boundary condition at x = : Because the flux value β is given at x = , we specify a function mapping from the flux β to the solution value u( ). This function is termed the Neumann-to-Dirichlet map at the boundary x = .
Considering the same manipulations leading to identity (3.5) (albeit on the whole interval [0, ]), in this case, we have Substituting the Neumann boundary condition and isolating u( ) from this equation results in the following: Therefore, the Neumann boundary condition (5.1) at x = is transformed into a Dirichlet boundary condition. Applying the scheme from the Dirichlet boundary conditions and using the Dirichlet boundary condition (5.2) instead of u( ) = β leads to the following theorem.

Theorem 2
The solution vector U = (u 1 , u 2 , . . . , u n ) T of the system of linear equations, that is, -a n-1 u n-1 + (a n-1 + a n )u n = a n u( ) + a n-1 G n-1 + a n H n ,  (5.5) and u( ) is defined in (5.2).
The Laplacian matrix (4.6) for the Dirichlet-Neumann boundary conditions (5.1) is the same as in the case of the Dirichlet boundary conditions. Only the last element of the vector on the right-hand side must be changed as follows: . . a n-2 G n-2 + a n-1 H n-1 a n u( ) + a n-1 G n-1 + a n H n  Therefore, the right-hand side vector F DN in 5.6 is the same as the vector F D in Example 4.1.
Consequently, the solution values of the system of linear equations are the same as well.

Case study for piecewise constant κ(x)
In this section, we demonstrate the effectiveness and robustness of the proposed exact scheme in the case of an ODE with a discontinuous conductivity function, which appears frequently in practical problems. Consider a physical problem based on (1.1) with homogeneous Dirichlet boundary conditions, in which the source term is f (x) = x and the is not continuous at x = 1, as shown in Fig. 4. Because the proposed method eliminates the numerical difficulties arising from the discontinuity, no smoothing (i.e., a sigmoid approximation of the Heaviside function) is required to obtain the exact solution. In particular, the discontinuity of κ(x) does not create issues for the implementation of the proposed method. Although the κ(x) function is discontinuous, its integrals exist, and thus the necessary calculations can be performed. The integral of the conductivity function on the whole domain is Thus, the shape functions are expressed as follows: as shown in Fig. 5. The different parts of the Green function are defined as follows: (6.5) Thus, based on the shape functions (6.3), (6.4) and (6.5), we can construct the Green functions corresponding to conductivity function (6.1) as follows: (6.6) Figure 6 shows the surface plot of the Green function G(x, t).
The solution of the problem u(x) can be derived by using (4.16); however, in general, the solution is defined using

Conclusion
This paper presents a practically useful and easily applicable robust scheme for calculating the exact solutions of a second-order self-adjoint ODE on any grid point corresponding to an arbitrary discretization. Through detailed derivations, we proved our theorems and assertions, and the effectiveness of the method was demonstrated through several relevant examples involving continuous and discontinuous κ(x). By rearranging the scheme, an implicit form of the solution of the differential equation can be obtained, which can be used to derive the analytical solution if it can be evaluated symbolically; in addition, by numerically evaluating the integrals, the proposed approach can provide the exact numerical values of the solution at any given point. Because the proposed scheme is easy to implement in various mathematical software environments (Mathematica, Maple, MAT-LAB, etc.), for problems possessing the given ODE structure, the exact solution can be determined, thereby avoiding the use of and problems associated with various numerical approximation methods. For all these reasons, the proposed method can be used not only to effectively solve mathematical problems with the differential operators in Eq. (1.1) but also to solve physical processes that can be described by the ODE in Eq. (1.1). In this case, if the integrals resulting from the applied scheme can be evaluated symbolically, the physical problem can be solved analytically, but if the integrals are numerically evaluated, we can obtain the exact solution to the physical problem at arbitrarily located grid points. At the moment, the proposed scheme works for the specific ODE structure discussed; a future direction can be the extension of this method to a broader range of differential operators.