Least-squares Solutions of Eighth-order Boundary Value Problems using the Theory of Functional Connections

This paper shows how to obtain highly accurate solutions of eighth-order boundary-value problems of linear and nonlinear ordinary differential equations. The presented method is based on the Theory of Functional Connections, and is solved in two steps. First, the Theory of Functional Connections analytically embeds the differential equation constraints into a candidate function (called a $constrained \, expression$) that contains a function that the user is free to choose. This expression always satisfies the constraints, no matter what the free function is. Second, the free-function is expanded as a linear combination of orthogonal basis functions with unknown coefficients. The constrained expression (and its derivatives) are then substituted into the eighth-order differential equation, transforming the problem into an unconstrained optimization problem where the coefficients in the linear combination of orthogonal basis functions are the optimization parameters. These parameters are then found by linear/nonlinear least-squares. The solution obtained from this method is a highly accurate analytical approximation of the true solution. Comparisons with alternative methods appearing in literature validate the proposed approach.

The technique presented in this paper is rooted in functional interpolation expressions. These expressions are particularly well suited to solve differential equations. This has been shown in Ref. [28], the seminal paper on the Theory of Functional Connections (TFC), and in subsequent articles, showing its application to solve linear [29] and nonlinear ODEs [30]. The TFC formalized the method of analytical constraint embedding (a.k.a. functional interpolation), since it provides expressions representing all functions satisfying a set of specified constraints.
The general equation to derive these interpolating expressions, named constrained expressions, follows as, y x, g x = g x + ∑ k = 1 n η k s k x (1) where g(x) is the free function, η k are unknown coefficients to be solved by imposing the n constraint conditions, and s k (x) are "support functions," which are a set of n linearly independent functions. In prior work [29,30] as well as in this paper, the s k (x) support function set has been selected as the monomial set.
The η k coefficients are computed by imposing the constraints using Equation (1). Then, once the expressions of the η k coefficients are obtained, they are back substituted into Equation (1) to produce the constrained expression, a functional representing all possible functions satisfying the specified set of constraints. The use of this constrained expression has already been applied to many areas of study, including solving low-order differential equations [29,30], hybrid systems [31], and optimal control problems, including energyoptimal landing, energy-optimal intercept [32], and fuel-optimal landing [33]. Furthermore, this technique has been successfully used to embed constraints into machine learning frameworks [34,35], in quadratic and nonlinear programming [36], and in a variety of other applications [37]. In addition, this technique has been generalized to n-dimensions [38,39], providing functionals representing all possible n-dimensional manifolds subject to constraints on the value and arbitrary order derivative of n − 1 dimensional manifolds. F x, y, y′, …, y 8 = 0 subject to: y k x i = y i k y k x f = y f k for k = 0, 1, 2, 3 (2) where the notation y k : = d k y k dx k is used to denote the k th derivative of y(x) with respect to x. Now, in order to embed the eight constraints, we can set n = 8 in Equation (1) leading to the expression, y x, g x = g x + η T s x (3) where in terms of this integration range. For completeness, the support functions for two general points x i and x f (i.e., the switching functions for Equation (4)) are provided in Appendix A. The β k terms for x ∈ [0, 1] are summarized below in Equations (5)- (12).
Particular attention must be paid when using this expansion and least-squares. For example, the basis functions in h(z) must be linearly independent of the support functions s(x) in order to solve the system via least-squares. If any of the terms in h(z) are not linearly independent of the support functions, then the matrix to be inverted in the least squares step will be illconditioned. Thus, the terms that are not linearly independent of the support functions must be skipped in the expansion of g(x). In this problem, our support functions span from the monomial term x 0 to x 7 ; therefore, the Chebyshev polynomial expansion must start from x 8 .
Furthermore, in general, the basis functions may not be defined on the same range as the problem domain (i.e., Chebyshev and Legendre polynomials are defined on z ∈ [−1, +1], Fourier series is defined on z ∈ [−π, +π], etc.). Therefore, the basis domain (z) must be mapped into the problem domain (x), which can be done using the simple linear equations, Furthermore, all subsequent derivatives of the free-function g(x) are defined as, where by defining, the expression can be simplified to, d n g dx n = c n ξ T d n h z dz n = c n ξ T h n z . (15) This defines all mappings from the basis domain into the problem domain. With the expression of g(x) in terms of a known basis, we can rewrite the constrained expression given in Equation (4) in the form, y x, ξ = a x, ξ + b x (16) where a(x, ξ) is a function that is zero where the constraints are defined and b(x) is a function that equals the constraints where they are defined. In fact, if g(x) is selected such that g(x) = 0 (meaning ξ = 0) the expression would simplify to an interpolating function y(x, 0) = b(x). Furthermore, since the ξ vector shows up linearly in the expression of a(x, ξ), Equation (16) can also be written as, where a(x, ξ) now becomes a vector equation a(x). This can be seen by expanding Equation (17), The subsequent derivatives follow by simply taking the derivatives of the h(z) and β k (x) terms. That is, the form of subsequent derivatives of the constrained expression remains the same and we can generally write the constrained expression up to the eighth-order derivative as shown in Equation (16), where a n x : = d n a x dx n is the n-th derivative of the a(x) function; the a(x) function also includes the derivative of h(z), which follows Equation (15) such that where c is defined in Equation (14). With this adjustment to the constrained expression, the transformed differential equation defined in Equation (13) can be reduced to a function of only x and the unknown vector ξ, F x, ξ = 0, (19) which may be linear or nonlinear in the unknown parameter ξ. Lastly, in order to solve this equation numerically, we must discretize the domain into N + 1 points. Since in this paper we consider the linear basis h(z) as Chebyshev orthogonal polynomials, the optimal distribution of N + 1 points is provided by Chebyshev-Gauss-Lobatto collocation points [40,41], defined as z k = − cos kπ N for k = 0, 1, …, N, and the map from z → x has been previously defined. As compared to the uniform point distribution, the collocation point distribution allows a much slower increase of the condition number as the number of basis functions, m, increases. In general, we can define the residual of our differential equation in Equation (19) for each discretized point, For a linear differential equation F (and therefore a linear differential equation F ) the constrained expression and its derivatives will show up linearly and therefore will remain linear in the unknown ξ term. This leads to the form where the matrix A is composed of a linear combination of a(x) and its derivatives discretized over x k where x = x 0 , ⋯, x k , ⋯, x N T . Note, by our definition the domain is x ∈ [x i , x f ], where x (k = 0) = x 0 = x i is the initial value and x(k = N) = x N = x f is the final value, and k is defined in the description of Chebyshev-Gauss-Lobatto collocation points. Furthermore, it follows that b(x) is composed of a linear combination of b(x), its derivatives, and a potential forcing term f (x) for the discrete values of x. This system can now be easily solved with any available least-squares technique. All numerical solutions in this paper utilize a scaled QR method to perform the least-squares.
In the case of a nonlinear differential equation, Equation (20) can be expressed as a loss function at each discretization point, and the system can be solved by an iterative least-squares method where the Jacobian is defined as, where ξ i represents the current step's estimated ξ parameter. The parameter update is provided by, where the ∆ξ i can be defined using classic least-squares, or in this paper through a QR decomposition method. This process is repeated until either the absolute value of the loss function is below some tolerance ϵ, or until the L 2 norm of the loss function continues to increase, which is specified by the following conditions, In this paper, this tolerance was set as twice the value of machine-level precision for double point precision, ϵ = 4.4409 × 10 −16 .

Parameter Initialization for Nonlinear Problems
For nonlinear problems the ξ must be initialized at the beginning of the iterative leastsquares process. In this paper, the initialization was chosen to be ξ 0 = 0 for all nonlinear problems. Setting the coefficient vector equal to zero is synonymous with selecting g(x) = 0, or in other words, choosing the constrained expression with the simplest interpolating polynomial satisfying all of the problem constraints. In the case of BVPs, the solution lies somewhere around this initial guess. Although introduced and solved in a later section, consider Problem #4 which involves solving the differential equation, This problem is highlighted in this section because it had the largest initialization error of the three nonlinear differential equations presented. Figure 1 displays the error of the solution due to the initialization technique of ξ = 0. As it can be seen in this figure, the error is on the order of 10 −7 for this specific case. The iterative least-squares will then reduce this error to close to machine-level precision. An astute reader will notice that the TFC method at initialization produces a more accurate solution than the techniques developed in [3,18]. A more detailed explanation of this result is discussed in the conclusion.

Numerical Solution
This section compares the TFC method with competing methods on a variety of problems, both linear and nonlinear. For each problem, the differential equation and the boundary conditions are presented, followed by a table that compares the absolute error of the two methods on a grid of 11 equidistantly spaced points that span the domain. In addition, each table includes the reference where the competing method's solution error was found.

Linear Eighth-Order
Discretizing the problem into points, x k where k ∈ [0, N], and collecting terms yields, This system can be solved by least-squares to yield the unknown coefficients, ξ, which can then be substituted back into the constrained expression to give the TFC estimate of the solution. Table 1 shows the absolute error of the TFC solution and the solution from Ref. [19] at each of the 11 points. Table 1 shows that the TFC solution error is orders of magnitude lower than the solution from Ref. [19] at all of the points in the domain, except the boundaries. At the boundaries, each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Problem #2-Consider the linear eighth-order differential equation solved in
Refs. [3,[13][14][15]18,19] Table 2 shows the absolute error of the TFC solution and the solution from Ref. [13] at each of the 11 points. Reference [13] did not report the solution at x = 0.9, so that entry in the table is labeled "not reported." Table 2 shows that the TFC solution error is orders of magnitude lower than the solution from Ref. [13] at all of the points in the domain, except the boundaries. At the boundaries, each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Problem #3-Consider the linear eighth-order differential equation solved in
Refs. [3,6,13,19]    each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Problem #4-Consider the nonlinear eighth-order differential equation solved in
Refs. [3,18]  Equations (16) and (18) show that the estimated solution, together with its third-order and eighth-order derivative, take the form Thus, using TFC, the differential equation can be re-written as, This system can be solved by an iterative least-squares method as shown in Section 2 to yield the unknown coefficients, ξ, which can then be substituted back into the constrained expression to give the TFC estimate of the solution. Table 4 shows the absolute error of the TFC solution, which was obtained in three iterations, and the solution from Ref. [18] at each of the 11 points. Table 4 shows that the TFC solution error is orders of magnitude lower than the solution from Ref. [18] at all of the points in the domain, except the boundaries. At the boundaries, each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Problem #5-Consider the nonlinear eighth-order differential equation solved in
Refs. [3,[14][15][16]18] Table 5 shows the absolute error of the TFC solution, which converged in 2 iterations, and the solution from Ref. [15] at each of the 11 points. Table 5 shows that the TFC solution error is orders of magnitude lower than the solution from Ref. [15] at all of the points in the domain, except the boundaries. At the boundaries, each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Problem #6-Consider the nonlinear eighth-order differential equation solved in
Refs. [3,15,18,19] Table 6 shows the absolute error of the TFC solution, which was obtained in 2 iterations, and the solution from Ref. [19] at each of the 11 points. Table 6 shows that the TFC solution error is orders of magnitude lower than the solution from Ref. [19] at all of the points in the domain, except the boundaries. At the boundaries, each of the methods has zero error, because each of the methods satisfies the boundary conditions exactly.

Accuracy of the Derivatives
The previous section compared the absolute error of TFC at a number of points along the domain with the absolute error of previous methods. This section discusses the accuracy of the derivatives when using TFC. One of the major advantages of TFC compared to other methods is that the estimated solution is analytical. As a result, further manipulation of the estimated solution is easily achieved, such as taking derivatives. As an example, consider problem #5. Table 7 shows the mean absolute error of y and its derivatives up to order eight. The second column of Table 7 used 10 basis functions to compute the solution and 11 equidistant points to compute the error, while the third column used 30 basis functions to compute the solution and 100 equidistant points to compute the error.
If enough Chebyshev orthogonal polynomials are used in the free function to estimate the solution, the error in subsequent derivatives should increase by an order of magnitude or less. Table 7 shows that when 10 basis functions are used, the error increases as the order of the derivative increases. In this case, there were not enough Chebyshev orthogonal polynomials used, as indicated by the large mean error in the eighth derivative. In other words, the number of basis functions used was not nearly enough to accurately estimate the solution of the eighth derivative.
When 30 basis functions are used, the mean error increases as the order of the derivative increases, until the eighth derivative is reached. In derivatives one through seven, the mean error increases approximately an order of magnitude or less when compared to the previous derivative. However, the eighth derivative has less error than the seventh derivative, because the eighth derivative shows up in the differential equation, and thus in the residual. Hence, the eighth derivative is directly affected when computing the solution, whereas the other derivatives are not.
In problem #5, the differential equation only contains the function and the eighth derivative. As a different example, consider the problem solved in [17], which has the exact solution y(x) = (x 2 − 1) sin(x). From hereon, we shall refer to this problem as problem #7. Table 8 shows the mean absolute error of y and its derivatives up to order eight for problem #7. The second column of Table 8 used 10 basis functions to compute the solution and 11 points to compute the error, while the third column used 30 basis functions to compute the solution and 100 points to compute the error. Table 8 shows that when all derivatives are included in the differential equation, the anomalous decrease in mean error as subsequent derivatives are taken disappears (i.e., the mean solution error from derivative seven to derivative eight increases as expected).
These results are similar to the results obtained in earlier studies on first-and second-order linear [29] and nonlinear [30] differential equations. This application to higher-order systems further highlights the power and robustness of this technique.
In Section 3, a discussion of the initialization of the TFC approach for nonlinear differential equations was provided. In the initialization, the coefficient vector ξ is set to zero, which implies that g(x) = 0. It was found that this still solved the differential equation with an accuracy on the order of O 10 −7 . This can be explained by an equation first presented in the seminal paper on TFC [28]. In this paper an equation (Equation (9) in Ref. [28]) is presented which describes the general expression for the interpolating expression for the function and its first n derivatives. This equation simplifies to the expression of a Taylor series when g(x) = 0. While this is not the exact case for the boundary-value constraints in this paper, the nature of the aforementioned equation points to the fact that when g(x) = 0, the constrained expression derived in this paper acts as a Taylor series approximation for two points. The relationship between the TFC method and Taylor series has yet to be explored, and the extension of this Taylor series-like expansion about n points will be the focus of future work.
Section 5 discussed the accuracy of the derivatives of the estimated solution. The solution accuracy was reduced with each subsequent derivative, but overall the accuracy of the final derivatives only lost a few orders of magnitude, resulting in an overall error on the order of O 10 −10 − 10 −12 , provided that enough basis functions were used when estimating the solution. In addition, it was shown that the accuracy of a given derivative depends, in part, on whether it explicitly shows up in the differential equation.

Abbreviations
The following abbreviations are used in this manuscript: BVP boundary value problem ODE ordinary differential equation