Theory of Functional Connections Applied to Linear ODEs Subject to Integral Constraints and Linear Ordinary Integro-Differential Equations

Abstract

This study shows how the Theory of Functional Connections (TFC) allows us to obtain fast and highly accurate solutions to linear ODEs involving integrals. Integrals can be constraints and/or terms of the differential equations (e.g., ordinary integro-differential equations). This study first summarizes TFC, a mathematical procedure to obtain constrained expressions. These are functionals representing all functions satisfying a set of linear constraints. These functionals contain a free function, $g\left(x\right)$, representing the unknown function to optimize. Two numerical approaches are shown to numerically estimate $g\left(x\right)$. The first models $g\left(x\right)$ as a linear combination of a set of basis functions, such as Chebyshev or Legendre orthogonal polynomials, while the second models $g\left(x\right)$ as a neural network. Meaningful problems are provided. In all numerical problems, the proposed method produces very fast and accurate solutions.

1. Introduction

This paper shows how to solve linear Ordinary Differential Equations (ODEs) and linear Integro-Differential Equations (IDEs) using a new mathematical framework to perform functional interpolation, called Theory of Functional Connections (TFC). TFC derives functionals, called constrained expressions, containing a free function and representing all possible functions satisfying a set of linear constraints [1,2,3,4]. The most important feature of the constrained expressions is: they always satisfy all the constraints no matter what the free function is.

Although it was recently developed (2017), TFC has already found several applications, especially in solving differential equations [5,6,7,8]. The free function can be expressed by any set of linearly independent functions, such as an expansion of orthogonal polynomials (e.g., Chebyshev, Legendre, etc.) or Neural Networks (NN), such as shallow NN with random features, or Deep NNs (DNNs). When the free function is expanded by a set of orthogonal polynomials the method is here identified as “classic TFC”. When shallow NNs with random features are used, the method has been identified as Extreme-TFC (X-TFC) [9], and when DNNs are used the method is called Deep-TFC [10].

The best expansion of the free function depends on the differential equation considered. Usually, for linear/nonlinear ODEs and simple bivariate PDEs, orthogonal polynomials represent the best choice in terms of accuracy. Indeed, orthogonal polynomials are generally the best mathematical tools for approximating and for convergence properties [11,12]. Nevertheless, for multivariate and more complex PDEs, orthogonal polynomials suffer the “curse of dimensionality”. For these problems, NNs represent a better choice as they better tolerate the curse of dimensionality.

In all previous applications, the free function has been expressed as a linear combination of known basis functions and unknown constant coefficients. These coefficients are then estimated by least-squares for linear Differential Equations (DEs) [5] and by nonlinear least-squares for nonlinear DEs [6]. The problems considered in this work will be solved with classic TFC using Chebyshev orthogonal polynomials and X-TFC.

Note that X-TFC and Deep-TFC are also identified as Physics-Informed Neural Networks (PINN) frameworks [9]. PINNs are recently developed machine-learning methods that employ NNs for data-physics-driven regression, data-physics-driven solution of DEs, and data-physics-driven discovery of parameters governing DEs [13]. Since X-TFC and Deep-TFC use NNs as free function, they are considered part of the PINNs family. Thanks to the constrained expression, developed by TFC, both X-TFC and Deep-TFC are more robust and accurate than the standard PINN frameworks as introduced by Raissi et al. [13].

TFC has been developed for univariate and multivariate scenarios [2,7,8] to solve a variety of mathematical problems: a homotopy continuation algorithm for dynamics and control problems [14], domain mapping [15], data-driven parameters discovery applied to epidemiological compartmental models [16], transport theory problems such as radiative transfer [17] and rarefied-gas dynamics [18], nonlinear programming under equality constraints [19], Timoshenko–Ehrenfest beam [20], boundary-value problems in hybrid systems [21], eighth-order boundary value problems [22], and in Support Vector Machine [23]. TFC has been widely used for solving optimal control problems for space application, solved via indirect methods [24]: orbit transfer and propagation [25,26,27,28], energy-optimal in relative motion [29], energy-optimal and fuel-efficient landing on small and large planetary bodies [30,31], the minimum time-energy optimal intercept problem [32].

The aim of this paper is to show how TFC can accommodate integral constraints in linear differential equations and solve linear integro-differential equations, is organized as follows. In Section 2, a summary of the TFC framework is provided, with the explanation of how to derive constrained expressions. In Section 3, the application of TFC to solve ODEs with integral constraints is shown, and some case problems are reported as examples. In Section 4, the TFC framework is used to solve linear integro-differential equations.

In all test problems considered in this article, the results have been obtained using MATLAB R2020a software on an Intel Core i7-9700 CPU PC with 64 GB of RAM. The results’ accuracy is provided in terms of the absolute error. That is,

$$err=|y_{\mathrm{TFC}}-y_{\mathrm{true}}|,$$

where $y_{\mathrm{TFC}}$ is the TFC approximated solution, and $y_{\mathrm{true}}$ is the true analytical solution.

Since the classic TFC uses Chebyshev or Legendre orthogonal polynomials as a basis set, the final Appendix A provides a definition, orthogonality, derivatives, and integral expressions and properties for both Chebyshev and Legendre orthogonal polynomials.

2. Theory of Functional Connections Summary

A mathematical generalization of interpolation, called Theory of Functional Connections, has recently been developed and successfully applied to solve, by least-squares, initial, boundary, and multi-value problems of linear [5] and nonlinear [6] ODEs and PDEs [8,10]. The theory has been developed for univariate [1] and multivariate rectangular domains [2,8].

The generalization of interpolation consists of a mathematical procedure to obtain analytical expressions describing all possible functions subject to n linear constraints. These expressions are functionals that are called constrained expressions. Two formally equivalent approaches [1,8] can be used to derive them. These are,

$$\begin{array}{cc}\hfill y(x,g(x\left)\right)& =g\left(x\right)+\sum _{k=1}^n{\eta }_k(x,g\left(x\right))s_k\left(x\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill y(x,g(x\left)\right)& =g\left(x\right)+\sum _{k=1}^n{\varphi }_k(x,\mathbf{s}\left(x\right)){\rho }_k(x,g\left(x\right)),\hfill \end{array}$$

(2)

where n is the number of the linear constraints, $g\left(x\right)$ is a function that can be freely chosen, ${\eta }_k(x,g\left(x\right))$ are functional coefficients, $\mathbf{s}\left(x\right)=\{s_1\left(x\right),\cdots ,s_n\left(x\right)\}$ is a set of n user-defined support functions that must be linearly independent (necessary conditions), ${\varphi }_k$ are switching functions, and ${\rho }_k(x,g\left(x\right))$ are projection functionals.

By imposing the n constraints to Equation (1), the values of the n functional coefficients, ${\eta }_k(x,g\left(x\right))$, are obtained for some set of n user-assigned functions, $s_k\left(x\right)$. Once the n functional coefficients, ${\eta }_k(x,g\left(x\right))$, are computed, by substituting them back into Equation (1) we obtain the constrained expression. To give an example, using $s_1\left(x\right)=x$ and $s_2\left(x\right)=x^2$, the constrained expression,

$$f(x,g\left(x\right))=g\left(x\right)+\underset{{\varphi }_1\left(x\right)}{\underbrace{{\displaystyle \frac{x(2x_2-x)}{2(x_2-x_1)}}}}\underset{{\rho }_1(x,g\left(x\right))}{\underbrace{(f_{x1}-g_x\left(x_1\right))}}+\underset{{\varphi }_2\left(x\right)}{\underbrace{{\displaystyle \frac{x(x-2x_1)}{2(x_2-x_1)}}}}\underset{{\rho }_2(x,g\left(x\right))}{\underbrace{(f_{x2}-g_x\left(x_2\right))}},$$

where the switching functions, ${\varphi }_k\left(x\right)$, and projection functionals, ${\rho }_k(x,g\left(x\right))$, are identified, always satisfies the two derivative constraints, $f_x\left(x_1\right)=f_{x1}$ and $f_x\left(x_2\right)=f_{x2}$, no matter what $g\left(x\right)$ is. In this article, the use of both formulations, provided by Equations (1) and (2), will be shown.

The projection functional, ${\rho }_k(x,g\left(x\right))$, projects the free function, $g\left(x\right)$, on the k-th constraint. Here is an example of projection functionals associated to a set of constraints,

$$\left\{\begin{array}{c}f\left(3\right)=9\hfill \\ 3f\left(0\right)-f\left(2\right)=0\hfill \\ \ddot{f}\left(1\right)=-1\hfill \end{array}\right.\to \left\{\begin{array}{c}{\rho }_1=9-g\left(3\right)\hfill \\ {\rho }_2=-3g\left(0\right)+g\left(2\right)\hfill \\ {\rho }_3=-1-\ddot{g}\left(1\right)\hfill \end{array}\right..$$

As for the switching functions, ${\varphi }_k\left(x\right)$, they are expressed as a linear combination of a set of support functions, $\mathit{s}\left(x\right)$. These functions satisfy ${\varphi }_k\left(x\right)=1$ when the k-th constraint is verified and ${\varphi }_k\left(x\right)=0$ when any other constraint is verified. For instance, given the support functions,

$$s_1\left(x\right)=x,s_2\left(x\right)=e^x,\mathrm{and}{s}_3\left(x\right)=sinx,$$

the switching functions are computed as ${\displaystyle {\varphi }_j\left(x\right)=\sum _{i=1}^3{\alpha }_{ji}{s}_i\left(x\right)}$, where the ${\alpha }_{ji}$ coefficients are computed by inverting the support matrix. For this example,

$$\left[\begin{array}{ccc}{\alpha }_{11}& {\alpha }_{12}& {\alpha }_{13}\\ {\alpha }_{21}& {\alpha }_{22}& {\alpha }_{23}\\ {\alpha }_{31}& {\alpha }_{32}& {\alpha }_{33}\end{array}\right]\left[\begin{array}{ccc}{s}_1\left(3\right)& 3s_1\left(0\right)-s_1\left(2\right)& {\ddot{s}}_1\left(1\right)\\ s_2\left(3\right)& 3s_2\left(0\right)-s_2\left(2\right)& {\ddot{s}}_2\left(1\right)\\ s_3\left(3\right)& 3s_3\left(0\right)-s_3\left(2\right)& {\ddot{s}}_3\left(1\right)\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$$

Reference [4] presents, in detail, how to derive the TFC constrained expressions using the formalism defined by Equation (2).

Constrained expressions have been used to solve differential equations. This is done by expressing the solution using the TFC constrained expressions from Equation (1) or, equivalently, from Equation (2). These functionals allow us to reduce the whole functions space to only the functions subspace satisfying the constraints. This is particularly important when solving differential equations. In fact, when substituting these functionals into the DE, a new differential equation is obtained in terms of $g\left(x\right)$. This new DE is subject to no constraints because the constrained expression fully satisfies the constraints. The unknown free function, $g\left(x\right)$, is then expressed as a linear combination of basis functions. In the classic TFC, $g\left(x\right)$ is a linear combination of orthogonal polynomials. That is,

$$g(x,\mathbf{\xi })=\sum _{j=0}^m{\xi }_j{h}_j\left(x\right)={\mathbf{\xi }}^{\mathrm{T}}\mathit{h}\left(x\right),$$

where m is the number of the basis functions $h_j\left(x\right)$ (orthogonal polynomials). In X-TFC, $g(x,\mathbf{\xi })$ is a shallow NN trained with Extreme Learning Machine (ELM) [33],

$$g(x,\mathbf{\xi })=\sum _{j=1}^L{\xi }_j{\sigma }_j\left(w_{j}x+b_j\right)={\mathbf{\xi }}^{\mathrm{T}}\mathbf{\sigma }\left(x\right)={\mathbf{\xi }}^{\mathrm{T}}\left\{\begin{array}{c}{\sigma }_1\left(x\right)\\ ⋮\\ {\sigma }_L\left(x\right),\end{array}\right\}$$

where L is the number of hidden neurons, $w_j\in \mathbb{R}$ is the input weights vector connecting the j-th hidden neuron and the input nodes, ${\xi }_j\in \mathbb{R}$ with $j=1,\cdots ,L$ is the j-th output weight connecting the j-th hidden neuron and the output node, and $b_j$ is the bias of the j-th hidden neuron, ${\sigma }_j(·)$ are activation functions, and $\mathbf{\sigma }={\{{\sigma }_1,\cdots ,{\sigma }_L\}}^{\mathrm{T}}$. According to the ELM algorithm [33], biases and input weights are randomly selected and not tuned during the training, thus they are known hyper-parameters. The activation functions, ${\sigma }_j(·)$, are also known as they are user selected. Thus, the only unknown NN hyper-parameters to compute are the output weights $\mathbf{\xi }={\left\{{\xi }_1,\cdots ,{\xi }_L\right\}}^{\mathrm{T}}$. Therefore, for both classic TFC and X-TFC, the unknown vector, $\mathbf{\xi }$, which actually constitutes the only unknown of our problem, is then estimated numerically, as for instance by least-squares for linear [5] and nonlinear [6] differential equation s.

In general, the independent variable of the DE (for instance, time) is defined in the $t\in [t_0,t_f]$ range, while the selected basis functions may be defined as a different range, $x\in [x_0,x_f]$ (for instance, Chebyshev and Legendre orthogonal polynomials are defined in the $x\in [-1,+1]$). Thus, a change of variable is needed. The most simple mapping between these two variables is linear,

$$x=x_0+{\displaystyle \frac{x_f-x_0}{t_f-t_0}}(t-t_0)⟷t=t_0+{\displaystyle \frac{t_f-t_0}{x_f-x_0}}(x-x_0).$$

(3)

By setting the range ratio, $c={\displaystyle \frac{x_f-x_0}{t_f-t_0}}$, the derivatives in terms of the new variable are

$${\displaystyle \frac{{\mathrm{d}}^{k}f}{\mathrm{d}{t}^k}}=c^k{\displaystyle \frac{{\mathrm{d}}^{k}f}{\mathrm{d}{x}^k}},$$

(4)

and the derivative constraints can be written as:

$${\left.{\displaystyle \frac{{\mathrm{d}}^{k}f}{\mathrm{d}{t}^k}}\right|}_{t_i}=f_{t_i}^{\left(k\right)}=c^k{\left.{\displaystyle \frac{{\mathrm{d}}^{k}f}{\mathrm{d}{x}^k}}\right|}_{x_i}=c^k{f}_{x_i}^{\left(k\right)}.$$

(5)

The change of variable in integrals takes advantage from the fact that the mean value is independent from the independent variable. Therefore,

$${\displaystyle {\displaystyle \frac{1}{t_f-t_0}}{\int }_{t_0}^{t_f}f\left(t\right)\mathrm{d}t={\displaystyle \frac{1}{x_f-x_0}}{\int }_{x_0}^{x_f}f\left(x\right)\mathrm{d}x\to c{\int }_{t_0}^{t_f}f\left(t\right)\mathrm{d}t={\int }_{x_0}^{x_f}f\left(x\right)\mathrm{d}x.}$$

When expressing the free function as $g(x,\mathbf{\xi })={\mathbf{\xi }}^{\mathrm{T}}\mathit{h}\left(x\right)$, then, the derivatives of $g(x,\mathbf{\xi })$ can be written as:

$${\displaystyle \frac{{\mathrm{d}}^{k}g}{\mathrm{d}{x}^k}}={\mathbf{\xi }}^{\mathrm{T}}{\displaystyle \frac{{\mathrm{d}}^k\mathit{h}}{\mathrm{d}{x}^k}}\mathrm{and}{\left.{\displaystyle \frac{{\mathrm{d}}^{k}g}{\mathrm{d}{x}^k}}\right|}_{x_i}={\mathbf{\xi }}^{\mathrm{T}}{\left.{\displaystyle \frac{{\mathrm{d}}^k\mathit{h}}{\mathrm{d}{x}^k}}\right|}_{x_i}.$$

This procedure can be applied to linear differential equations with non-constant coefficients. The final expression obtained is linear in terms of the unknown vector, $\mathbf{\xi }$, and can be written as:

$${\mathit{a}}^{\mathrm{T}}\left(x\right)\mathbf{\xi }=b\left(x\right).$$

(6)

To solve the problem numerically, this new equation is discretized for a set of N distinct values of x. A different discretization scheme can be used. Usually, the discretization points are either randomly selected or linearly uniformly spaced. When using Chebyshev polynomials the best discretization is to use the zeros of the Chebyshev polynomials, also called Chebyshev points or nodes or, more formally, Chebyshev–Gauss points [12]. They are defined by the cosine distribution,

$$x_k=cos\left({\displaystyle \frac{(2k-1)\pi }{2N}}\right)\mathrm{for}k=1,\cdots ,N.$$

By specifying Equation (6) for these $x_k$ values a system of N linear equations is obtained in m (or L) unknowns that is then solved for $\mathbf{\xi }$ by least-squares. Several least-squares methods can be used to solve (6). The optimal least-squares method to use for each problem depends on the problem itself and on what TFC technique is used to solve the problem (e.g., in this paper, classic TFC or X-TFC). For instance, for the classic TFC and linear DEs, our analysis identifies the QR decomposition on a scaled coefficient matrix as the best approach minimizing the condition number of the matrix to invert.

3. TFC for ODEs with Integral Constraints

In this section, we show how TFC is applied to solve linear ODEs with integral constraints. After showing how to derive the constrained expression when dealing with a general integral constraint, we will solve a couple of linear ODEs subjects to boundary conditions and integral constraints.

For the first two examples in this section, we will derive the constrained expression by following the formulation of Equation (1). For subsequent examples and problems, however, the second formulation, Equation (2), will be adopted.

3.1. Definite Integral Constraint

Let us consider the integral constraint,

$${\displaystyle {\int }_a^{b}f\left(x\right)\mathrm{d}x=I.}$$

(7)

The constrained expression has the form,

$$f(x,g\left(x\right))=g\left(x\right)+{\eta }_1{s}_1\left(x\right),$$

(8)

where $g\left(x\right)$ is a free function and $s_1\left(x\right)$ is a user-defined support function, and ${\eta }_1$ is a coefficient that is derived by imposing the constraint. By integrating Equation (8) the following equation:

$${\displaystyle I={\int }_a^{b}g\left(x\right)\mathrm{d}x+{\eta }_1{\int }_a^b{s}_1\left(x\right)\mathrm{d}x}$$

is obtained, which can be rearranged to obtain the expression for ${\eta }_1$,

$${\eta }_1={\displaystyle \frac{{\displaystyle I-{\int }_a^{b}g\left(x\right)\mathrm{d}x}}{{\displaystyle {\int }_a^b{s}_1\left(x\right)\mathrm{d}x}}}.$$

Substituting this value into Equation (8), we obtain a constrained expression for $f\left(x\right)$, that is, an expression that always satisfies the integral constraint for any expression of $g\left(x\right)$,

$$f(x,g\left(x\right))=g\left(x\right)+{\displaystyle \frac{s_1\left(x\right)}{{\displaystyle {\int }_a^b{s}_1\left(x\right)\mathrm{d}x}}}\left({\displaystyle I-{\int }_a^{b}g\left(x\right)\mathrm{d}x}\right).$$

(9)

This function becomes undefined when ${\displaystyle {\int }_a^b{s}_1\left(x\right)\mathrm{d}x=0}$, and thus $s_1\left(x\right)$ must be selected to avoid this condition. By simply selecting $s_1\left(x\right)=1$ Equation (9) becomes,

$${\displaystyle f(x,g\left(x\right))=g\left(x\right)+\frac{1}{b-a}\left({\displaystyle I-{\int }_a^{b}g\left(x\right)\mathrm{d}x}\right).}$$

3.2. Integral and Linear Constraints

In this second example, let us consider the more complex case where, in addition to the integral constraint given in Equation (7), we also consider the additional linear constraints,

$$\alpha f\left(x_0\right)+\beta f\left(x_f\right)=1,$$

where $\alpha$, $\beta$, $x_0$, and $x_f$ are all assigned. A sketch of this example is shown in Figure 1.

The constrained expression has the form,

$$f(x,g\left(x\right))=g\left(x\right)+{\eta }_1{s}_1\left(x\right)+{\eta }_2{s}_2\left(x\right).$$

(10)

Then, by applying the constraints the system of linear equations,

$$\left\{\begin{array}{c}1-\alpha g_0-\beta g_f\\ {\displaystyle I-{\int }_a^{b}g\left(x\right)\mathrm{d}x}\end{array}\right\}=\left[\begin{array}{cc}\alpha s_1\left(x_0\right)+\beta s_1\left(x_f\right)& \alpha s_2\left(x_0\right)+\beta s_2\left(x_f\right)\\ {\displaystyle {\int }_a^b{s}_1\left(x\right)\mathrm{d}x}& {\displaystyle {\int }_a^b{s}_2\left(x\right)\mathrm{d}x}\end{array}\right]\left\{\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right\}$$

(11)

is obtained. This system tells us that $s_1\left(x\right)$ and $s_2\left(x\right)$ functions can be any functions with the exception of those making the matrix singular. The matrix singularity occurs when

$${\displaystyle \left(\alpha s_1\left(x_0\right)+\beta s_1\left(x_f\right)\right){\int }_a^b{s}_2\left(x\right)\mathrm{d}x=\left(\alpha s_2\left(x_0\right)+\beta s_2\left(x_f\right)\right){\int }_a^b{s}_1\left(x\right)\mathrm{d}x.}$$

For instance, by selecting $s_1\left(x\right)=1$ and $s_2\left(x\right)=x$, the previous condition becomes,

$$(\alpha +\beta )(b+a)=2(\alpha x_0+\beta x_f)$$

which imply that, in order to avoid singularity, the following relationship must be satisfied,

$$\left\{\begin{array}{cc}\alpha ,& \beta \end{array}\right\}\left\{\begin{array}{c}b+a-2\alpha x_0\\ b+a-2\beta x_f\end{array}\right\}\ne 0.$$

When this condition is satisfied, the matrix of Equation (11) can be inverted and the coefficients, ${\eta }_1$ and $\eta$, can be computed. Then, the constrained expression for this problem is obtained by substituting the expressions found for ${\eta }_1$ and ${\eta }_2$ in Equation (10).

Problem #1

Now, let us consider the following DE to be solved on the $t\in [0,\pi ]$ range, with a constraint in the form of an integral over the integration range,

$$\ddot{f}+f=0\mathrm{subject}\mathrm{to}:\left\{\begin{array}{c}f\left(0\right)=1\hfill \\ {\displaystyle {\int }_0^{\pi }f\left(t\right)\mathrm{d}t=\pi }\hfill \end{array}\right.,$$

whose analytical solution is

$$f={\displaystyle \frac{\pi }{2}}sint+cost.$$

In this problem, we build the constrained expression according to the formulation of Equation (2), where the projection functionals ${\rho }_k(x,g\left(x\right))$ are given by

$$\left\{\begin{array}{c}{\rho }_1=1-g_0\hfill \\ {\displaystyle {\rho }_2=\pi -{\int }_0^{\pi }g\left(t\right)\mathrm{d}t}\hfill \end{array}\right..$$

Given the support functions,

$$s_1\left(t\right)=1,\mathrm{and}{s}_2\left(t\right)=t,$$

the switching functions are computed as ${\displaystyle {\varphi }_j\left(t\right)=\sum _{i=1}^2{\alpha }_{ji}{s}_i\left(t\right)}$, where the ${\alpha }_{ji}$ coefficients are computed just by inverting a matrix. For this problem,

$$\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]\left[\begin{array}{cc}{s}_1\left(0\right)& {\displaystyle {\int }_0^{\pi }{s}_1\left(\tau \right)\mathrm{d}\tau }\\ s_2\left(0\right)& {\displaystyle {\int }_0^{\pi }{s}_2\left(\tau \right)\mathrm{d}\tau }\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right].$$

Then, the constrained expression can be built as:

$$f(t,g\left(t\right))=g\left(t\right)+{\varphi }_1\left(t\right)(1-g_0)+{\varphi }_2\left(t\right)\left({\displaystyle \pi -{\int }_0^{+\pi }g\left(\tau \right)\mathrm{d}\tau }\right),$$

where the switching functions are

$${\varphi }_1\left(t\right)={\displaystyle \frac{\pi -2t}{\pi }}\mathrm{and}{\varphi }_2\left(t\right)=\frac{2t}{{\pi }^2}.$$

Now, expressing the free function according to $g(t,\mathbf{\xi })={\mathbf{\xi }}^{\mathrm{T}}\mathit{h}\left(x\left(t\right)\right)$ the terms of the DE become:

$$\begin{array}{cc}\hfill f(x,\mathbf{\xi })& ={\left\{{\displaystyle \mathit{h}\left(x\right)-{\varphi }_1\left(x\right){\mathit{h}}_0-{\varphi }_2\left(x\right)\frac{1}{c}{\int }_{-1}^1\mathit{h}\left(x\right)\mathrm{d}x}\right\}}^{\mathrm{T}}\mathbf{\xi }+{\varphi }_1\left(x\right)+\pi {\varphi }_2\left(x\right)\hfill \\ \hfill \ddot{f}(x,\mathbf{\xi })& =c^2\ddot{\mathit{h}}{\left(x\right)}^{\mathrm{T}}\mathbf{\xi },\hfill \end{array}$$

where the boundary conditions have been mapped to $x\in [x_0,x_f]$. For classic TFC $x\in [-1,+1]$, and for X-TFC $x\in [0,1]$ according to Equations (3)–(5). The mapping coefficient is $c={\displaystyle \frac{x_f-x_0}{t_f-t_0}}$.

With the function now expressed in terms of $\mathbf{\xi }$ the differential equation can now be transformed into a function of $\mathbf{\xi }$,

$${\left\{{\displaystyle c^2\ddot{\mathit{h}}\left(x\right)+\mathit{h}\left(x\right)-{\varphi }_1\left(x\right){\mathit{h}}_0-{\varphi }_2\left(x\right)\frac{1}{c}{\int }_{-1}^1\mathit{h}\left(x\right)\mathrm{d}x}\right\}}^{\mathrm{T}}\mathbf{\xi }=-{\varphi }_1\left(x\right)-\pi {\varphi }_2\left(x\right).$$

This equation can then be specified for a discrete values of x. This discretization yields to an over-determined linear system that can be solved by least-squares.

$$\mathbb{A}\mathbf{\xi }=\mathit{b}\to \mathbf{\xi }={\left({\mathbb{A}}^{\mathrm{T}}\mathbb{A}\right)}^{-1}{\mathbb{A}}^{\mathrm{T}}\mathit{b}.$$

(12)

In this example, using $n=100$ and $m=20$, the classic TFC was executed with a computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-16}\right)$, the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-16}\right)$. For X-TFC, we set the following hyper-parameters $n=100$, $L=100$, Gaussian activation function, input weight and bias were sampled from $\mathcal{U}[-10,+10]$. The computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-16}\right)$, the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-16}\right)$.

Table 1. Absolute errors on uniform test points for problem #1.

$\mathit{t}/\mathit{\pi }$	TFC	X-TFC
0.0	0.0	0.0
0.1	0.0	$2.22\times {10}^{-16}$
0.2	0.0	$2.22\times {10}^{-16}$
0.3	0.0	$0.0$
0.4	0.0	$0.0$
0.5	0.0	$0.0$
0.6	$4.44\times {10}^{-16}$	$4.44\times {10}^{-16}$
0.7	$1.11\times {10}^{-16}$	$2.22\times {10}^{-16}$
0.8	$6.66\times {10}^{-16}$	$0.0$
0.9	$1.66\times {10}^{-16}$	$2.77\times {10}^{-16}$
1.0	$2.22\times {10}^{-16}$	$6.66\times {10}^{-16}$

reports the absolute error with respect to the analytical solutions obtained with classic TFC and X-TFC on 11 points uniformly distributed in the $[0,1]$ range.

3.3. Mixed Constraints

Consider a function subject to a value-level, relative, and integral constraint,

$${\displaystyle f\left(t_1\right)=f_1,f\left(t_0\right)=f\left(t_f\right),\mathrm{and}{\displaystyle \frac{1}{t_f-t_0}}{\int }_{t_0}^{t_f}f\left(\tau \right)\mathrm{d}\tau =I.}$$

By using Equation (2), the constrained expression has the form,

$$f(t,g\left(t\right))=g\left(t\right)+{\varphi }_1\left(t\right){\rho }_1(t,g\left(t\right))+{\varphi }_2\left(t\right){\rho }_2(t,g\left(t\right))+{\varphi }_3\left(t\right){\rho }_3(t,g\left(t\right))$$

where the projection functionals, ${\rho }_k(t,g\left(t\right))$, are given by

$$\left\{\begin{array}{c}{\rho }_1=f_1-g\left(t_1\right)\hfill \\ {\rho }_2=g\left(t_f\right)-f\left(t_0\right)\hfill \\ {\displaystyle {\rho }_3=I-{\displaystyle \frac{1}{t_f-t_0}}{\int }_{t_0}^{t_f}g\left(\tau \right)\mathrm{d}\tau }\hfill \end{array}\right..$$

Given the support functions,

$$s_1\left(t\right)=1,s_2\left(t\right)=t,\mathrm{and}{s}_3\left(t\right)=t^2$$

the ${\alpha }_{ji}$ coefficients of the switching functions, ${\displaystyle {\varphi }_j\left(t\right)=\sum _{i=1}^3{\alpha }_{ji}{s}_i\left(t\right)}$, are computed by inverting the support matrix. Specifically,

$$\left[\begin{array}{ccc}{\alpha }_{11}& {\alpha }_{12}& {\alpha }_{13}\\ {\alpha }_{21}& {\alpha }_{22}& {\alpha }_{23}\\ {\alpha }_{31}& {\alpha }_{32}& {\alpha }_{33}\end{array}\right]\left[\begin{array}{ccc}{s}_1\left(t_1\right)& s_1\left(t_0\right)-s_1\left(t_f\right)& \frac{1}{t_f-t_0}{\int }_{-\pi }^{+\pi }{s}_1\left(\tau \right)\mathrm{d}\tau \\ s_2\left(t_1\right)& s_2\left(t_0\right)-s_2\left(t_f\right)& \frac{1}{t_f-t_0}{\int }_{-\pi }^{+\pi }{s}_2\left(\tau \right)\mathrm{d}\tau \\ s_3\left(t_1\right)& s_3\left(t_0\right)-s_3\left(t_f\right)& \frac{1}{t_f-t_0}{\int }_{-\pi }^{+\pi }{s}_3\left(\tau \right)\mathrm{d}\tau \end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$$

Once the ${\varphi }_j\left(t\right)$ are known, the constrained expression can be obtained as:

$$\begin{array}{cc}\hfill f(t,g(t\left)\right)=g\left(t\right)& +{\varphi }_1\left(t\right)\left(f_1-g\left(t_1\right)\right)+{\varphi }_2\left(t\right)\left(f\left(t_f\right)-g\left(t_0\right)\right)+\hfill \\ \hfill & +{\varphi }_3\left(t\right)\left({\displaystyle I-{\displaystyle \frac{1}{t_f-t_0}}{\int }_{t_0}^{t_f}g\left(\tau \right)\mathrm{d}\tau }\right).\hfill \end{array}$$

Problem #2

As an ultimate “stress-test” of the TFC method, a mixed-constraint case is considered where point, relative, and integral constraints are used in the solution of a differential equation, on the range $t\in [-\pi ,\pi ]$

$$\stackrel{⃛}{y}+sint\ddot{y}+(1-t)\dot{y}+ty=f\left(t\right)\mathrm{subject}\mathrm{to}:\left\{\begin{array}{c}y\left(t_0\right)=0\hfill \\ y\left(t_f\right)=y\left(t_0\right)\hfill \\ {\displaystyle {\int }_{-\pi }^{+\pi }y\left(\tau \right)\mathrm{d}\tau =-2\pi }\hfill \end{array}\right.$$

where the forcing term is

$$f\left(t\right)=(t-1){sin}^{2}t+(2+2t-t^2-2cost)sint+t(t-1)cost.$$

This problem admits the analytical solution,

$$y=(1-t)sint.$$

Following the procedure explained, the constrained expression for this problem is,

$$y(t,g\left(t\right))=g\left(t\right)+{\varphi }_1\left(t\right)(-g_0)+{\varphi }_2\left(t\right)(-g_f)+{\varphi }_3\left(t\right)\left({\displaystyle -2\pi -{\int }_{-\pi }^{+\pi }g\left(\tau \right)\mathrm{d}\tau }\right),$$

where the ${\alpha }_{ji}$ coefficients of the switching functions ${\displaystyle {\varphi }_j\left(t\right)=\sum _{i=1}^3{\alpha }_{ji}{s}_i\left(t\right)}$ are computed by solving the following system:

$$\left[\begin{array}{ccc}{\alpha }_{11}& {\alpha }_{12}& {\alpha }_{13}\\ {\alpha }_{21}& {\alpha }_{22}& {\alpha }_{23}\\ {\alpha }_{31}& {\alpha }_{32}& {\alpha }_{33}\end{array}\right]\left[\begin{array}{ccc}{s}_1\left(t_0\right)& s_1\left(t_f\right)& {\displaystyle {\int }_{-\pi }^{+\pi }{s}_1\left(\tau \right)\mathrm{d}\tau }\\ s_2\left(t_0\right)& s_2\left(t_f\right)& {\displaystyle {\int }_{-\pi }^{+\pi }{s}_2\left(\tau \right)\mathrm{d}\tau }\\ s_3\left(t_0\right)& s_3\left(t_f\right)& {\displaystyle {\int }_{-\pi }^{+\pi }{s}_3\left(\tau \right)\mathrm{d}\tau }\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$$

Now, expressing the free function according to $g\left(t\right)={\mathbf{\xi }}^{\mathrm{T}}\mathit{h}\left(x\left(t\right)\right)$ and rearranging terms leads us to the final constrained expression,

$$y(x,\mathbf{\xi })={\left\{{\displaystyle \mathit{h}\left(x\right)-{\varphi }_1\left(x\right){\mathit{h}}_0-{\varphi }_2\left(x\right){\mathit{h}}_f-{\varphi }_3\left(x\right)\frac{1}{c}{\int }_{x_0}^{x_f}\mathit{h}\left(x\right)\mathrm{d}x}\right\}}^{\mathrm{T}}\mathbf{\xi }-2\pi {\varphi }_3\left(x\right),$$

and by substituting this constrained expression into the differential equation (by evaluating its derivatives accordingly), we obtain a linear system that can be solved by least-squares using Equation (12).

Even in this case, the differential equation and boundary conditions must be mapped onto $x\in [x_0,x_f]$. For classic TFC $x\in [-1,+1]$, and for X-TFC $x\in [0,1]$ according to Equation (3), where the mapping coefficient is $c={\displaystyle \frac{x_f-x_0}{t_f-t_0}}$.

For classic TFC we set $n=35$ and $m=30$. The computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-15}\right)$, the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-30}\right)$. For X-TFC we set the following hyper-parameters $n=40$, $L=40$, sinusoidal activation function, input weight and bias were sampled from $\mathcal{U}[-20,+20]$. The computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-15}\right)$, the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-30}\right)$.

3.4. Discussions

The results from these two problems emphasize the utility and convenience of using TFC. Once the constrained expression is defined, the method to solving the differential equation does not change regardless of the use of different constraints. Moreover, this makes the solution accuracy only dependent on the problem complexity. Finally, once the solution is computed (on the discretization points), we have an analytical representation of it. That is, no further manipulations are needed (e.g., interpolation) if we want to evaluate the solutions in points that are different from the discretization points. Indeed, as shown in the results of both problems we do not lose any accuracy when evaluating the solution on test points, as it would happen with some other state-of-the-art methods such as the Finite Element Method (FEM) [34].

4. TFC for Linear Ordinary Integro-Differential Equation

In this section we show how TFC is applied to solve linear ordinary Integro-Differential Equations (IDEs) using all constrained expressions derived from the formalism of Equation (1).

As first test problem, the linear Fredholm integro-differential equations is considered. Comparisons of the numerical results obtained with TFC, X-TFC, and the method published in Ref. [35] are presented.

4.1. Problem #1

Consider the following integro-differential equation:

$${\displaystyle \dot{y}\left(x\right)=1-\frac{1}{3}x+x{\int }_0^1\tau y\left(\tau \right)\mathrm{d}\tau }$$

with one constraint $y\left(0\right)=0$, for $x\in [0,1]$. The analytical solution is $y\left(x\right)=x$. The constrained expression for this problem is simply

$$y(x,\mathbf{\xi })={\left(\mathit{h}\left(x\right)-{\mathit{h}}_0\right)}^{\mathrm{T}}\mathbf{\xi },$$

where ${\mathit{h}}_0$ is the vector of the basis function computed at $x=0$.

For classic TFC $n=100$ discrete points and $m=29$ basis functions were used. The computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-5}\right)$, and the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-9}\right)$. For X-TFC the hyper-parameters were set to $n=20$, $L=20$, the activation function was sinusoidal, and the input weight and bias were sampled from $\mathcal{U}[-1,+1]$. The computational time obtained was of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points of $\mathcal{O}\left({10}^{-4}\right)$, the variance of the absolute error on the discretization points of $\mathcal{O}\left({10}^{-8}\right)$. The absolute errors obtained with classic TFC and X-TFC on 10 test points are reported in Figure 2 and compared with those reported by Ref. [35].

4.2. Problem #2

The second test problem is on the following integro-differential equation:

$${\displaystyle \dot{y}\left(x\right)=x{e}^x+e^x-x+x{\int }_0^{1}y\left(\tau \right)\mathrm{d}\tau ,}$$

with one constraint $y\left(0\right)=0$, for $x\in [0,1]$. The analytical solution is $y\left(x\right)=x{e}^x$. Also for this this problem the constrained expression is

$$y(x,\mathbf{\xi })={(\mathit{h}\left(x\right)-{\mathit{h}}_0)}^{\mathrm{T}}\mathbf{\xi }.$$

For classic TFC $n=25$ discretization points and $m=14$ basis functions were adopted. The computational time obtained is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-16}\right)$, and the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-32}\right)$. For X-TFC the hyper-parameters were set to $n=50$, $L=50$, the activation function was the hyperbolic tangent, and the input weight and bias were sampled from $\mathcal{U}[-1,+1]$. The computational time is of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points is of $\mathcal{O}\left({10}^{-14}\right)$, and the variance of the absolute error on the discretization points is of $\mathcal{O}\left({10}^{-27}\right)$. The absolute errors on the test points obtained with classic TFC and X-TFC, are reported in Figure 3 and compared with those reported by Ref. [35].

4.3. Problem #3

The last test problem is a second-order integro-differential equation,

$${\displaystyle \ddot{y}\left(x\right)=e^x-x+x{\int }_0^1\tau y\left(\tau \right)\mathrm{d}\tau ,}$$

subject to the initial value problem, $y\left(0\right)=1$ and $\dot{y}\left(0\right)=1$, where $x\in [0,1]$. The analytical solution is $y\left(x\right)=e^x$. The constrained expression for this problem is:

$$y(x,\mathbf{\xi })={\left(\mathit{h}\left(x\right)-{\mathit{h}}_0-x{\dot{\mathit{h}}}_0\right)}^{\mathrm{T}}\mathbf{\xi }+x+1.$$

To solve this problem, the classic TFC required $n=9000$ points and $m=1000$ basis functions. The computational time was of $\mathcal{O}\left(1\right)$ second, the average absolute error on the discretization points was of $\mathcal{O}\left({10}^{-5}\right)$, and the variance of the absolute error on the discretization points was of $\mathcal{O}\left({10}^{-10}\right)$. For X-TFC the setting was $n=90$ and $L=94$ hyper-parameters, hyperbolic sine activation function, and input weight and bias were sampled from $\mathcal{U}[-1,+1]$. The computational time was of $\mathcal{O}\left({10}^{-4}\right)$ s, the average absolute error on the discretization points of $\mathcal{O}\left({10}^{-4}\right)$, and the variance of the absolute error on the discretization points of $\mathcal{O}\left({10}^{-7}\right)$. The absolute errors on the test points obtained with classic TFC and X-TFC, are reported in Figure 4 and compared with those reported by Ref. [35].

4.4. Discussions

The discussions made for the linear ODE subjects to integral constraints are also valid for linear ordinary integro-differential equations. Nevertheless, the attentive reader will notice that there is a significant difference in terms of accuracy among problems 1 and 3 compared to problem 2. This is due to the fact that the kernels in the integral terms of the equations are different than one for problems 1 and 3 (e.g., for the both problems the kernel is t). This can have a negative impact on the solution accuracy as the function to be integrated can become very complex, and this can lead to numerical issues. For instance, when classic TFC is used in problems 1 and 3, where the kernel is simply t, the coefficients in front of the integrated polynomials become enormous as the number of basis functions (e.g., the degree of the Chebyshev polynomials) reaches above ten, causing numerical issues that impact negatively on the accuracy. Moreover, in many real world applications, such as in radiative transfer equation [17,36], the kernel can be very complex, making the analytical evaluation of the integral prohibitive.

The are several ways of mitigating these issues that are worth exploring. For instance, the integral could be evaluated with a Gauss-quadrature scheme or with a Monte Carlo method. The authors of this manuscript are investigating the possibility of evaluating the integral numerically using NNs. However, the investigation of different ways to overcome these limitations is not the focus of this work, and it will be explored by the authors in future papers.

5. Conclusions

This study presents an extension of the least-squares solution of linear differential equations studied in an earlier publication [5]. This paper aims to explore solutions using different boundary constraints including multi-valued, relative, and integral constraints for linear ordinary differential equations and for linear ordinary integro-differential equations. In all problems, a constrained expression is used to embed the constraints. This allows the integration range to be independent from the constraints location. In these expressions, the free function $g\left(x\right)$ was expressed as a linear combination of known basis functions and unknown constant coefficients $\mathbf{\xi }$, either using Chebyshev polynomials, or shallow NNs for the X-TFC framework. Expressing $g\left(x\right)$ as a linear combination of basis functions and constant coefficients, the coefficient vector $\mathbf{\xi }$ remains linear and is solved by using a least-squares method. While this paper only solves linear differential equations as a test case, the proposed methodology can easily be extended to nonlinear differential equations by replacing the least-squares method with a nonlinear least-squares technique as introduced in prior research [6].

For all cases explored, the classic TFC and X-TFC methods consistently produce very fast and accurate solutions. Additionally, for all problems, the constraints are analytically satisfied since they are integrated into the constrained expression. In general, this means that the TFC methods decouple the constraints from the solution of the differential equation. Due to this, the solution range of the differential equation and, where the constraints are specified, is independent. This fact makes the TFC method a unified framework to solve differential equations with no more distinctions between the initial and boundary values problems, as well as any other constraints distribution problem.