General Exact Schemes for Second-Order Linear Differential Equations Using the Concept of Local Green Functions

: In this paper, we introduce a special system of linear equations with a symmetric, tridi-agonal matrix, whose solution vector contains the values of the analytical solution of the original ordinary differential equation (ODE) in grid points. Further, we present the derivation of an exact scheme for an arbitrary mesh grid and prove that its application can completely avoid other errors in discretization and numerical methods. The presented method is constructed on the basis of special local green functions, whose special properties provide the possibility to invert the differential operator of the ODE. Thus, the newly obtained results provide a general, exact solution method for the second-order ODE, which is also effective for obtaining the arbitrary grid, Dirichlet, and/or Neumann boundary conditions. Both the results obtained and the short case study conﬁrm that the use of the exact scheme is efﬁcient and straightforward even for ODEs with discontinuity functions.


Introduction
Consider the linear second order inhomogeneous ordinary differential equation (ODE): where κ(x) ≥ κ 0 > 0 is a positive, arbitrary function, q(x) and f (x) are continuous functions and the function u(x) is the solution of ODE on the finite closed interval [0, ]. The use of boundary value problem (BVP) with Equation (1) for modeling physical phenomena is very rare in the relevant literature. A large part of the possible implicit application areas are the initial value problems (IVPs) prescribed by Equation (1) with constant coefficients. BVPs including the coefficients of the resistance (R), capacity (C) and inductance (L) can be used to model RLC circuits [1]. It is also well known that the IVP for Equation (1) is also used in harmonic analysis, where the coefficients are the spring constant, damping and mass [2]. Another group of applications consist of reaction diffusion equations, which, however, even in the simplest case, contain some nonlinear terms [3]. The third major field of science where equations such as Equation (1) can be found is the Sturm-Liouville theory [4]. In this case, although BVPs are solved for non-constant coefficient ODEs, equations can be interpreted as eigenvalue problems, and solutions are also known in special cases [4]. It follows that the BVPs defined by Equation (1) are a special group of ODEs. Perhaps the most well-known of these is the so-called Cauchy-Euler equation [5], where κ(x) = x, q(x) = 1/x and f (x) = 0.
Although the solutions of the mathematical problems mentioned above are well known, since several mathematical methods have been developed for their production, solutions for non-constant coefficient equations are known only in special cases (for example, the Cauchy-Euler equation [5]). Common methods for solving non-constant coefficient IVPs are the constant variation method, the method of undetermined coefficients, etc. [5].
The BVPs with Equation (1) are usually solved using well-known numerical methods, such as the finite difference method (FDM) [6], the finite volume method (FVM) [7] and the boundary element method (BEM) [8], by approximating the analytical solution of the ODE. Among the classical numerical methods for solving second-order differential equations, the finite element method (FEM) is presently the most widely used [9]. Nothing better illustrates the challenges of BVP defined by Equation (1) than the fact that, even in recent years, we can find publications where the authors researched methods for solving the non-constant coefficient problem [10][11][12].
This research is motivated by two main research directions. One of these directions is the study of symbolic approaches related to non-linear ODEs and systems, the so-called reaction-diffusion, cross-diffusion equations [13,14]. Further, the research group is ambitious in studying other types of non-linear ODEs, which are still intensively researched today [15][16][17]. The second research area is the inverse problems involving second-order differential equations, which are also non-linear mathematical problems [18]. Still, compared to the classical differential equation solvers, they require a completely different handling, since the functions to be determined are the coefficients of the differential equation [19][20][21]. All these research activities are supported by a solid foundation in exact schemes, which, as effective and robust mathematical methods, provide new perspectives for basic research.
Finding and then applying an analytical solution to solve the physical problems minimizes the additional errors induced by mathematical techniques. These efforts have already motivated many researchers who, as a result of their work, have presented some higher-order schemes for first-and second-order differential equations, as were presented in [22]. Other authors have published papers about using the beneficial properties of the Bessel equation [5] to solve homogenous equations with a second-order ODE in order to obtain the analytic solution [23]. In the case of singular secondary ODEs, solvability is still the subject of mathematical studies. Kiguradze et al. [24][25][26] published several related studies that thoroughly examine the difficulties of these mathematical problems to be solved.
This research also aims to develop a new analytical solving techniques to solve the BVPs in case of Equation (1). During our research work, we developed an exact solution technique that uses the inverse of the differential operator to solve the ODE. The methods we have deployed propose schemes that provide the values for the exact solution (We define the exact solution as the value of the analytical solution at any given location.) of Equation (1) at arbitrarily chosen grid points. The specific properties of the exact schemes are the main novelty of the presented research, with the introduced strength of their application compared to other alternative methods. In fact, there are two possible implementations. In the first case, without a mesh, we can use them to define symbolically the analytical solution itself. In the second case, using a mesh, we produce the values of the analytical solution in terms of grid points, without needing to determine the analytical solution itself [27]. The main advantage of this characteristic is its robustness and efficiency, which makes it irrelevant for the method how many grid points are chosen inside the domain, since the calculated values will always be the values of the analytical solution in grid points [27]. Compared with alternative solution methods, the use of exact schemes also leads to the solution of a linear system of equations with outstanding properties, still the matrix and the right-hand side of the system of equations are constructed with special functions, and the integration leads to the stability required to solve the linear system of equations. All these beneficial properties are guaranteed by the so-called local green functions, whose concept was presented by Vizvari et al. [27]. Local Green functions differ from global green functions in the following. While global green functions can be used to symbolically represent the analytical solution of the ODE, on the other hand, the local green functions are defined on the mesh intervals with the same properties as global green functions and they are used to form the system of linear equations that are generated with the exact scheme.
The application concept of the local Green functions [27] opens the possibility to develop a mathematical method deriving a recursive formula which allows us to compute the values of the analytic solution in grid points using simple linear algebraic methods. In the self-adjoint case, when q(x) = 0 for Equation (1), it has already been shown that integral functions of 1 κ(x) can be used to derive exact schemes using the apparatus of the FEM [27]. Using the electrical analogy, these local Green functions are integral functions of the material's resistance, more precisely the resultant resistance of the material under investigation. In this paper, we also intend to use strictly positive test functions that generate the exact schemes based on the FEM apparatus for Equation (1) [9]. However, compared to the FEM, local green functions have the ability to invert the differential operator of the ODE as a result of performing the necessary partial integration. Based on the previous results [27], the local green functions have the same properties as the integrals of 1 κ(x) ; however, in this case, the local green functions are the results of solving a special initial value problem, since this guarantees the existence of the necessary and important properties that the exact scheme exploits. These facts complicate the results published in the previous work [27], since the solution of the initial value problem for a homogeneous equation can only be solved by using symbolic techniques, since here the goal is to generate functions. After all, the implementation of the exact scheme is possible even by fully numerical methods, since the use of local green functions and the necessary mathematical manipulations lead to discrete values of the analytical solution. The robustness and efficiency of the application of local green functions and exact schemes are best demonstrated by the solvability of ODEs with discontinuity functions and/or singular points, which will be presented in detail in the case studies.
Thus, the research of exact schemes is an important aspect in the domain of mathematical physics, and increasing the accuracy of such schemes is a continuous challenge for researchers. Furthermore, increasing the accuracy also increases the complexity and computational burden of the schemes, resulting in difficulties in applicability. However, our proposed schemes can overcome these limitations, as they provide the exact solution at an arbitrary location independent from the space discretization. Based on this fact, our research team have developed exact schemes for the symmetric, self-adjoint case (q(x) = 0), that also work robustly for Dirichlet and mixed boundary conditions [27]. In this paper, as a continuation of our research work [27] and as a generalization of the self-adjoint case (q(x) =0), we present our latest results in the field of exact schema research work. In the following sections, we will introduce the application of the exact scheme to Equation (1). This paper is structured as follows. The first section is the introduction, and the second section presents the exact schemes developed for two cases. The first subsection introduces the concept of local green functions, the second subsection defines the exact scheme with the Dirichlet boundary condition [5] prescribed in both endpoints, and the third subsection shows the handling of mixed boundary conditions, where Dirichlet at one and Neumann boundary condition [5] at the other endpoint are prescribed. The third section introduces a brief case study with discontinuous function κ(x). From these results, the fourth section draws the conclusions and establishes the future work.

The Concept of Local Green Functions
As a first step in the concept of local green functions and by obtaining the solution of the initial value problems in Equations (3) and (4), we define (2n) non-negative functions on different subintervals. The exact scheme is generated based on a special symmetric system of the local Green functions. Now, let us consider an arbitrary discretization of the interval [0, ] into (n + 1) subintervals: . Using these discretization, the local green functions are defined as the solutions of the initial value problems described in Equations (3) and (4): Based on this, in order to produce the exact scheme, we need a positive monotonically increasing function (ψ i−1 (x)) and a positive, but monotonically decreasing function (ϕ i (x)), which we have defined as a solution to the initial value problem corresponding to a homogeneous equation and the appropriate initial conditions. From the derivatives listed in Equations (3) and (4), it is clear that the test functions ψ i increase monotonically and ϕ i decreases symmetrically on the interval of their domain [x i , x i+1 ] for all i = 1, 2, · · · , n − 1. On the first interval [x 0 , x 1 ], only the test function ψ 0 is considered; and on the last interval [x n , x n+1 ], we define only the test function ϕ n similar to the function in the FEM [9]. Based on the equalities in (3) and (4), these test functions can be considered as local green functions [27], which have the following useful property: (3) and (4). The test functions ψ i (x) and ϕ i (x) are equal at the appropriate endpoints of the [x i , x i+1 ] interval, i.e.,: where i = 1, 2, · · · , (n − 1) and a i is the i th coefficient of the constructed system of linear equations.
Proof. To prove the Theorem 1, in the first step, multiply Equation (1) by ψ i (x) and in the second step by ϕ i (x), subtract them from each other and integrate the difference on interval [x i , x i+1 ] with respect to x (using the properties presented in the Equations (3a) and (4a)): where the prime ( ) denotes the derivative with respect to x. After simplification and applying integration by parts on the right-hand side, the following form is obtained: Taking advantage of the special properties of the test functions (Equations (3b), (4b), (3c) and (4c)), the formula is simplified to: This proves Theorem 1.
Having introduced the symmetric concept of local green functions, the properties presented in Equations (3) and (4) and Theorem 1 allow us to construct the exact schemes as follows.

Dirichlet Boundaries
In the case of Dirichlet boundaries [5], the values of the solution at the endpoints are obtained as follows: where B D is the set of Dirichlet boundary conditions [5]. The scheme that is used to obtain exact solutions to the Dirichlet problem in Equation (1) is described by the following theorem: Theorem 2 (Exact Dirichlet scheme). If ψ i−1 (x) and ϕ i (x) are the test functions [9] satisfying the Equation (3), the Equation (4), respectively, and Theorem 1, the solution vector from following the exact scheme yields to the values of the analytical solution of Equation (1) at the grid points with boundary conditions in Equation (9). Equations (3) and (4) are special IVPs defined based on the homogeneous equation corresponding to the original ODE (Equation (1)).
Let us consider a system of n ≥ 3 linear equations: −a n−1 u n−1 + b n u n = a n β + a n−1 G n−1 + a n H n with indexes i = 2, 3, · · · , (n − 1). The solution vector of Equation (10) U = (u 1 , u 2 , . . . , u n ) T leads to the same values as the solution u(x) of the second-order ODE (1) with boundary conditions (9) at the interior grid points (2) without any error; that is, u(x i ) = u i . The matrix of the system of linear equations constructed by recursion (Equation (10)) is a tridiagonal, symmetric, diagonally dominant (if q(x) = 0), n × n matrix, whose inverse exists in all cases. The coefficients in Equation (10) are defined as where i = 1, 2, · · · , n.
Proof. In the first step, we multiply the ODE (1) by ψ i−1 (x) and integrate it over I i−1 : x i After applying partial integration, the following expression is derived using the properties of ψ i−1 and multiplying with a i−1 (based on Equation (3)): Similarly, we multiply ODE (1) by ϕ i (x) and integrate it over I i to obtain By applying partial integration, the following relation is obtained using the properties of ϕ i and multiplying with a i (based on Equation (4)): The fundamental recursive formula is obtained by adding Equation (16) and Equation (18), and using the equalities and obtaining from Lψ i−1 (x) = 0 and Lϕ i (x) = 0, respectively, we obtain: If we substitute i = 1, 2, 3, · · · , n and u(0) = α, u( ) = β, we obtain Equation (10) in Theorem 2. It follows that: where u(x i ) is the value of the analytical solution of the Equation (1) in the i th mesh point and u i is the i th element of the solution vector U.

Exact Scheme for Dirichlet and Neumann Boundaries
In this case, we obtain the solution of the ODE (Equation (1)) using the following boundary conditions: Owing to the definition of the boundary conditions in Equation (23), the Dirichlet scheme cannot be used directly to obtain the values of the exact solution. However, the value of u( ) can be calculated using β by defining the Neumann-to-Dirichlet transformation.
Consider that the same manipulations lead to the identity (16), albeit on the entire interval [0, ]. In this case, we obtain: After the substitution κ( )u ( ) = δ, u(0) = α and isolating u( ) from Equation (24), the Neumann boundary condition κ( )u ( ) = δ is transformed into the Dirichlet boundary value u( ), which can be substituted directly in the result presented in Theorem 2. The result of the isolation is

Implementation of the Exact Scheme to Solve Case Studies
In this section, we demonstrate the effectiveness and robustness of the proposed exact scheme in the case of a discontinuous κ(x) and a regular singular ODE, that appear in several physical and mathematical problems based on Equation (1). However, before solving the case studies, consider generally the application of the exact scheme in a stepwise sequence:

1.
After defining the ODE to be solved and the necessary boundary conditions (Dirichlet problem, mixed problem), based on Equation (2), we perform an arbitrary partition of the interval [0, ]. The indexes of the nodes are: i = 0, 1, . . . , n − 1, n, n + 1, where 0 and n + 1 are the indices of the boundary points.

4.
Boundary conditions are also easily managed: • For Dirichlet boundary conditions, we use the substitution defined in Equation (9), • For mixed boundary conditions, we apply the Neumann to Dirichlet transformation based on the Equation (25) in order to calculate the potential value at the boundary point where the Neumann boundary condition is defined. (10), using the coefficient values and boundary conditions calculated before, and solve it by inverting the tridiagonal matrix. 6.

Construct the linear system of equations defined in Equation
The elements of the vector obtained in this way are equal to the values of the analytic solution of the ODE defined in the first step taken in the grid points, according to the boundary conditions (according to Theorem 2).
The case studies presented in this work have been evaluated symbolically and numerically using the MAPLE software package in order to be able to show the local green functions in a symbolic way.

A Case Study for Discontinuous κ(x) with Increasing Mesh Resolution
First, we solve a Dirichlet and then a Dirichlet-Neumann problem using a mesh grid, wherein one of the grid points is provided at the breakpoint of the κ(x) function. Let the functions in Equation (1) be a q(x) = 2 constant, and discontinuous functions. The solution of a problem prescribed in this manner is boundary conditions. In this case study, our goal is to show the efficiency and robustness of the exact scheme by solving an ODE with κ(x), f (x) (presented in the Equations (26) and (27)) and constant q(x) functions for several mesh partitions with increasing resolution. Consider the mesh point coordinates in the interval [0, 2] as follows: The values of u(x) are determined using Theorem 2 in the grid points of the mesh. For the partitioning described in Equation (30), we divide the domain into 4 subintervals, i.e., n + 1 = 4, hence we define the values of the analytic solution at n = 3 points. The point of interest is x 2 = 1, which is exactly the breakpoint of κ(x), where the analytic solution can also be determined using the exact scheme. Consider first the mesh in Equation (30), which is defined based on instructions in Section 3.1 (point 1). According to point 2 in Section 3.1, Equations (3) and (4), the obtained local Green functions are defined in the Equation (31).   Figure 1 illustrates the properties of the local green functions ψ i−1 and ϕ i . On each interval, ψ i−1 starts from zero and monotonically increases, according to the method used to generate it, and ϕ i starts from the maximum value and monotonically decreases until it reaches zero. In Figure 1, the symmetry property of the local green functions introduced in Theorem 1 is also evident. These properties of the local green functions, explained and exploited in Theorems 1 and 2, provide their physical meaning. To construct the exact schema defined in Theorem 2, the local green functions guarantee continuity boundary conditions in the connection of the subintervals. In fact, given the pair of test functions ψ i−1 (x) and ϕ i (x) with finite support on the interval [x i−1 , x i+1 ], their zeros and special properties provide the flux elimination necessary to produce the recursive formula shown in Equation (10). Moreover, the continuity boundary conditions in the connection points of the subintervals, for the equality of the values of u(x), arise from the special properties of the derivatives of the local green functions (Equations (3c) and (4c)), which will be exploited in the proof of Theorem 2.
Once the local green functions are determined, the coefficient values a i−1 , b i , G i−1 and H i can be calculated from Equations (11)- (14) to construct the linear system of equations described in Equation (10) (i = 2, 3). Based on the Theorem 2, the system in Equation (10) is expressed as follows: The values of the solution vector (u 1 , u 2 , u 3 ) T = (45/64, 0, 5/8) T (obtained from Equation (32)) are equal to the values of u(x) that are obtained in points x 1 = 1/4, x 2 = 1 and x 3 = 3/2. The results of the case study are illustrated in Figure 2.
As it can be seen in Figure 2, the values of u i clearly coincide with the function u(x). This illustrates, that the results obtained by applying the exact scheme described in Theorem 2 are equal to the values of the solution of ODE at the same grid points.
In the next step, we increase the resolution of the mesh by using the MAPLE random number generator to choose points within the interval [0, 2], under the condition that we retain x = 1 where is the discontinuity of κ(x). Considering all this, by increasing the mesh generation resolution step by step, we solve the ODE defined (in terms of Equations (26) and (27) and q(x) = 2) under the boundary conditions in Equation (29) in three different cases:
n + 1 = 101 subintervals and the analytic solution is calculated in n = 100 points. In the next step, we solve the problem by using the considerations defined in Theorem 2. The solutions are shown in Figure 3. Figure 3 clearly shows, that although the handling of the exact scheme is very similar to the classical numerical methods, the behavior of the solution is quite different. Comparing the graphs in Figure 3, it can be seen that the calculated potential values highlighted in red are visually identical to the grid point values of the analytical solution, i.e., the points highlighted in red coincide with the corresponding values of the function u(x) indicated in blue. Moreover, since the aforementioned facts can be observed in all three figures, we can conclude that the outstanding property proved in Theorem 2 is stable regardless of the mesh resolution.
In the last step, we can also solve the same example problem above with the functions κ(x) defined in the Equation (26), q(x) = 2, f (x) defined in the Equation (27) and u(0) = 1, and κ(2)u (2) = 14 mixed boundary conditions, since these can be defined by u(x). Subsequently, using the Dirichlet-to-Neumann transformation (Equation (24)), we can easily verify that u(2) = 3. During the application of the Equation (24), the continuity of the solution and flux must be ensured to determine ψ 0 (x). Our statement can easily be verified by substituting into the formula of the Dirichlet-to-Neumann transformation and then isolating u(2) from the formula.

Case Study for a Singular ODE with a Piecewise Source Term
In this section, we demonstrate the effectiveness and robustness of the proposed exact scheme in the case of a regular singular ODE, which appears in several physical and mathematical problems based on Equation (1).
Equation (33) is a second order Cauchy-Euler type, singular-regular differential equation that has singularities at x = 0 and x → ∞. The aim of our case study is to show that the mathematical method presented in this paper is suitable for the effective treatment of singularities in addition to creating a solution for ODE Equation (33) with the boundary conditions defined in Equation (35). Furthermore, the exact scheme also handles the behavior of the source term robustly. In particular, the non-differentiable property of f (x) does not cause any problems for the implementation of the method. Although the function f (x) is non-differentiable, its integrals still exist; thus, the necessary calculations can be performed. We will examine how u(x) behaves at the singular points of ODE (Equation (33)). To perform this, we will add new mesh points a, b to exclude singularities and then obtain function u(x) on the new mesh by using the results presented in Theorem 2.
Next, we will examine how the solution can be modified into a form in which the solution also exists at the singular points of the ODE by taking the limits a → 0 + , b → ∞.
In order to obtain the solution (u(x)) to Equation (33), we choose the Dirichlet boundary conditions as follows: where 0 < a ≤ 1/2 and 1/2 ≤ b < ∞. The shape functions are and Since functions ψ 0 (t) and ϕ 1 (t) are solutions to Equations (3) and (4), respectively, they have the following properties: ψ 0 (a) = 0, ϕ 1 (b) = 0 and ψ 0 (a) = 1 a and ϕ 1 (b) = − 1 b . Figure 4 shows functions ψ 0 (t) and ϕ 1 (t) with substitutions a = 1, b = 2 and m = 2. The local green functions in Figure 4 have the same properties as the functions in Figure 1, except that here the test function pair ψ 0 and ϕ 1 are the only pairs on the entire interval [a, b]. Obviously, this difference does not affect either the mathematical properties of the functions, or their physical properties according to the continuity boundary conditions. After obtaining a 0 , a 1 , b 1 , G 0 and H 1 from ψ 0 (t) and ϕ 1 (t) using the Equations (11)- (14) in Theorem 2, the fundamental recursive formula (Equation (10)) for n = 1 inner point, and After isolating u(x) from the Equation (38) and a necessary simplification is carried out, the solution of the ODE (in Equation (33)) with the corresponding boundary conditions (in Equation (35)) is given by Equation (39). where In the next step, we examine how the solution function u(x) behaves at the singular points of the original ODE. The solution considers that the singular points are obtained by taking the limits a → 0 + , b → ∞. We also consider the properties of the f (x) to evaluate improper integrals. If a → 0 + , the limit of u(x) can be obtained from Equations (39) and (40) as follows when m > 0: Equation (41) is a significantly simplified form of Equation (39), which shows that the substitution has eliminated the u(0) = α Dirichlet boundary condition. Furthermore, it can be seen that at this point, when taking the limit b → ∞, the convergence of improper integrals is fulfilled only for certain functions f (x). In our case study, we chose a source term that fulfills the following condition: The equation for free space is obtained as: Now, we can substitute f (x) (Equation (34)) into (43) and evaluate the integrals. Since f (x) is piecewise, the solution will be (if 0 < m = 2).
It is easy to notice that u 0 (x) satisfies the original ODE (Equation (33)). It is interesting to note that in the case of m = 2, the denominator of the solution (at interval 0 ≤ x < 1 2 ) becomes zero. It is also noticeable that the numerator of the same function also becomes zero for m = 2. Therefore, instead of substitution, we can take the limit, which when evaluated with the L'Hospital rule, leads to the following result: is constant to any m. Figure 5 shows the graph of u 0 (x) for m = 1, 2, 3, 4.  (44)) of the case study for a singular ODE with a piecewise source term (m = 1, 2, 3, 4).
As shown in Figure 5 and in Equations (44) and (45), function u 0 (x) is continuous and differentiable at x = 1 2 , even in the case of a discontinuous f (x).

Conclusions
The application of the proposed scheme provides the opportunity to generate the values of the analytical solution that is obtained in grid points. This method can be easily implemented in different software environments (such as MATLAB, Maple, and Mathematica); thus, many physical problems can be modeled to avoid errors that are caused by the numerical approximations when using arbitrary grids. Currently, the most critical issue in the application of this scheme is to obtain the solution of a homogeneous equation, which may be challenging. Therefore, in future generalizations, we will consider the importance of providing a solution to this problem to complete the scope of the application areas. Furthermore, we generalize the concept presented in this paper by the application of the Sturm-Liouville theory, focusing on the symmetry of the local green functions and the gained advantageous properties.
The application of the exact scheme to Sturm-Liouville problems will represent a milestone in the research, that will complete the study of linear second-order ODEs. As a consequence, we will be able to concentrate our research activities on the symbolic treatment of nonlinear ODEs and inverse problems, as the main two motivating areas of the research group.