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) containing 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.

Introduction

This paper has been motivated by several articles dedicated to estimating the solutions of high-order boundary-value problems (BVPs) including, fourth-order [1], sixth-order [2], eighth-order [3-19], 2m-order [20], and higher-order [21-24] BVPs. This paper focuses specifically on eighth-order BVPs because of the volume of research done on them, which is covered in Refs. [3-19]. These references list many existing scientific problems requiring solutions of high-degree BVPs. For example, eighth-order BVPs appear in the physics of specific hydrodynamic stability problems (infinite horizontal layer of fluid heated from below and under rotation) when instability sets in as overstability [25], and in orthotropic cylindrical shells under line load [26]. From a theoretical point of view, the study of the existence and uniqueness for the solutions of high-order boundary value problems is presented in Ref. [22] and studied further in Ref. [27].

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, where g(x) is the free function, η are unknown coefficients to be solved by imposing the n constraint conditions, and s(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(x) support function set has been selected as the monomial set.

The η coefficients are computed by imposing the constraints using Equation (1). Then, once the expressions of the η 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 energy-optimal 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.

Derivation of the Constrained Expression for Eighth-Order Boundary-Value Problems

In this paper, we consider the solution of an eighth-order BVP via the TFC. In general, the problem can be posed in its implicit form as, where the notation is used to denote the k 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, where and, using monomial support functions,

Now, according to the theory developed in Ref. [28], a system of equations can be constructed by evaluating the candidate function defined by Equation (3) at the constraint locations and setting the function equal to the specified constraint value. For example, the constraint on the function at the initial value (i.e., y(x) = y) is applied as such,

This can be done for the remaining seven constraint conditions, and the resulting system of equations can be expressed in a compact form,

This system of equations can be solved for the unknown coefficients and organized in the form, where the β (x) terms, called switching functions, only depend on the independent variable. This technique is general for any domain x ∈ [x, x], for example the general expression for β1(x) is,

However, for ease of presentation, since all problems presented in this paper, except for Problem #5, are defined on the domain x ∈ [0, 1], we will express these switching functions in terms of this integration range. For completeness, the support functions for two general points x and x (i.e., the switching functions for Equation (4)) are provided in Appendix A. The β terms for x ∈ [0, 1] are summarized below in Equations (5)–(12).

With the solution of the β (x) terms, the constrained expression is fully solved and represents all possible functions satisfying the boundary-value constraints. More specifically, by substituting the constrained expression and its derivatives into the original differential equation a new differential equation in terms of g(x) and its derivatives is obtained. This new differential equation, which has no constraints, can be written in the compact form,

To solve this differential equation, prior work [29,30,32,33,37] has expanded g(x) as a linear combination of m basis functions, where z = z(x), is an m × 1 vector of unknown coefficients, and (z) is an m × 1 vector containing the m basis functions (in this paper Chebyshev orthogonal polynomials are used). Particular attention must be paid when using this expansion and least-squares. For example, the basis functions in (z) must be linearly independent of the support functions (x) in order to solve the system via least-squares. If any of the terms in (z) are not linearly independent of the support functions, then the matrix to be inverted in the least squares step will be ill-conditioned. 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 x0 to x7; therefore, the Chebyshev polynomial expansion must start from x8. 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,

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, 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 (x). This can be seen by expanding Equation (17),

The subsequent derivatives follow by simply taking the derivatives of the (z) and β (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 is the n-th derivative of the (x) function; the (x) function also includes the derivative of (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 , 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 (z) as Chebyshev orthogonal polynomials, the optimal distribution of N + 1 points is provided by Chebyshev-Gauss-Lobatto collocation points [40,41], defined as 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 ) 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 is composed of a linear combination of (x) and its derivatives discretized over x where . Note, by our definition the domain is x ∈ [x, x], where x (k = 0) = x0 = x is the initial value and x(k = N) = x = x is the final value, and k is defined in the description of Chebyshev-Gauss-Lobatto collocation points. Furthermore, it follows that () 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 represents the current step’s estimated parameter. The parameter update is provided by, where the ∆ 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 L2 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 least-squares 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.

Figure 1.

Initialization error of the solution of Problem #4 by imposing 0 = 0.

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.

Although the TFC approach to solve differential equations has typically used (N ∈ [100, 200] and m ∈ [20, 80]) [29,30,34,37], in order to make a commensurable comparison with the techniques developed in Refs. [3-19], a grid of N = 11 points was used and m = 10 basis terms were selected: one less than the number of points selected. The following sections introduce the six most commonly solved eighth-order differential equations in Refs. [3-19], and compare the TFC solution accuracy with that of the most accurate solution from these references. Furthermore, to supplement the theory given in the previous sections, a step-by-step procedure for TFC is laid out for the first linear and non-linear problem.

Linear Eighth-Order Problems

Problem #1

Consider the linear eighth-order differential equation solved in Refs. [3,13,19] subject to which has the exact solution y(x) = (1 − x)e.

From Equations (16) and (18), we are shown that the estimated solution and eighth-order derivative take on the forms, and

Thus, using TFC, the differential equation can be re-written as,

Discretizing the problem into points, x where k ∈ [0, N], and collecting terms yields, where,

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.

Problem #1: Absolute solution error.

x	TFC Absolute Error	Ref. [19] Absolute Error
0	0	0
0.1	2.2204 × 10−16	6.3 × 10−11
0.2	1.1102 × 10−16	6.5 × 10−10
0.3	1.1102 × 10−16	2.0 × 10−09
0.4	1.1102 × 10−16	3.3 × 10−09
0.5	1.1102 × 10−16	3.9 × 10−09
0.6	6.6613 × 10−16	3.4 × 10−09
0.7	2.7756 × 10−15	2.0 × 10−09
0.8	3.8858 × 10−15	6.9 × 10−10
0.9	8.4932 × 10−15	7.6 × 10−11
1	0	0

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-15,18,19] subject to which has the exact solution y(x) = x(1 − x)e.

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.

Problem #2: Absolute solution error.

x	TFC Absolute Error	Ref. [13] Absolute Error
0	0	0
0.1	0	1.63 × 10−10
0.2	8.3267 × 10−17	1.63 × 10−09
0.3	0	4.90 × 10−09
0.4	1.1102 × 10−16	8.46 × 10−09
0.5	5.5511 × 10−17	1.01 × 10−08
0.6	3.8858 × 10−16	8.68 × 10−09
0.7	3.3307 × 10−16	5.15 × 10−09
0.8	3.3307 × 10−16	1.76 × 10−09
0.9	8.0769 × 10−15	Not reported
1	0	0

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] subject to which has the exact solution y(x) = (x2 − 1) sin(x).

Table 3 shows the absolute error of the TFC solution and the solution from Ref. [19] at each of the 11 points.

Table 3.

Problem #3: Absolute solution error.

x	TFC Absolute Error	Ref. [19] Absolute Error
0	0	0
0.1	2.7756 × 10−17	6.6 × 10−12
0.2	2.7756 × 10−17	6.9 × 10−11
0.3	0	2.1 × 10−10
0.4	5.5511 × 10−17	3.5 × 10−10
0.5	0	4.1 × 10−10
0.6	7.2164 × 10−16	3.5 × 10−10
0.7	1.3323 × 10−15	2.1 × 10−10
0.8	1.1102 × 10−15	7.2 × 10−11
0.9	3.4417 × 10−15	8.0 × 10−12
1	0	0

Table 3 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.

Nonlinear Eighth-Order Problems

Problem #4

Consider the nonlinear eighth-order differential equation solved in Refs. [3,18] subject to which has the exact solution y(x) = e.

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,

Discretizing the problem into points, x where k ∈ [0, N], leads to the loss function for some values of = . The Jacobian of the loss function with respect to is,

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.

Problem #4: Absolute solution error.

x	TFC Absolute Error	Ref. [18] Absolute Error
0	0	0
0.1	2.2204 × 10−16	2.503395 × 10−06
0.2	0	8.940697 × 10−06
0.3	2.2204 × 10−16	1.561642 × 10−05
0.4	4.4409 × 10−16	1.823902 × 10−05
0.5	2.2204 × 10−16	8.821487 × 10−06
0.6	6.6613 × 10−16	7.510185 × 10−06
0.7	3.5527 × 10−15	1.883507 × 10−05
0.8	7.5495 × 10−15	1.931190 × 10−05
0.9	1.0214 × 10−14	1.168251 × 10−05
1	0	0

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-16,18] subject to which has the exact solution y(x) = ln(1 + x).

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.

Problem #5: Absolute solution error.

x	TFC Absolute Error	Ref. [15] Absolute Error
0	0	0
0.1	1.5266 × 10−16	2.01 × 10−07
0.2	1.5821 × 10−15	4.54 × 10−07
0.3	7.0083 × 10−14	1.52 × 10−06
0.4	2.5846 × 10−13	4.07 × 10−06
0.5	3.2330 × 10−13	6.71 × 10−06
0.6	1.3139 × 10−13	9.06 × 10−06
0.7	2.1261 × 10−14	1.00 × 10−05
0.8	2.0539 × 10−14	5.45 × 10−06
0.9	3.3307 × 10−16	2.59 × 10−06
1	0	0

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] subject to which has the exact solution y(x) = e−.

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.

Problem #6: Absolute solution error.

x	TFC Absolute Error	Ref. [19] Absolute Error
0	0	0
0.1	1.1102 × 10−16	2.9 × 10−12
0.2	1.1102 × 10−16	2.7 × 10−11
0.3	0	7.6 × 10−11
0.4	0	1.3 × 10−10
0.5	1.1102 × 10−16	1.5 × 10−10
0.6	1.1102 × 10−16	1.3 × 10−10
0.7	2.2204 × 10−16	7.6 × 10−11
0.8	3.2196 × 10−15	2.5 × 10−11
0.9	9.9920 × 10−16	2.4 × 10−12
1	0	0

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.

Table 7.

Mean absolute error of all derivatives for Problem #5.

Function	Mean Absolute Error: 10 Basis Functions	Mean Absolute Error: 30 Basis Functions
y	7.5585 × 10−14	9.6866 × 10−16
y′	1.0534 × 10−12	7.5884 × 10−15
y′′	2.0202 × 10−11	5.0360 × 10−14
y(3)	4.9228 × 10−10	4.0456 × 10−13
y(4)	1.3318 × 10−08	2.8079 × 10−12
y(5)	3.8469 × 10−07	1.3927 × 10−11
y(6)	1.3150 × 10−05	5.5250 × 10−11
y(7)	3.9359 × 10−04	2.0221 × 10−10
y(8)	1.9399 × 10−02	1.4765 × 10−12

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], subject to which has the exact solution y(x) = (x2 − 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.

Mean absolute error of all derivatives for Problem #7.

Function	Mean Absolute Error: 10 Basis Functions	Mean Absolute Error: 30 Basis Functions
y	4.8255 × 10−15	5.8919 × 10−15
y′	8.3368 × 10−15	9.9755 × 10−15
y′′	7.1054 × 10−14	5.9525 × 10−14
y(3)	4.8760 × 10−13	4.8352 × 10−13
y(4)	1.5118 × 10−12	1.6443 × 10−12
y(5)	4.1244 × 10−12	4.4041 × 10−12
y(6)	3.1934 × 10−12	3.2911 × 10−12
y(7)	8.0532 × 10−12	7.3956 × 10−12
y(8)	8.1927 × 10−11	8.7722 × 10−12

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).

Conclusions

This paper explores the application of the techniques developed in [28-30] to the solution of high-order differential equations, namely eighth-order BVPs. In all the problems presented, which span the publications [3-19], the solution accuracy ranges from . 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 . 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 , 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.