Deep Theory of Functional Connections: A New Method for Estimating the Solutions of PDEs

TL;DR

This paper introduces deep Theory of Functional Connections (TFC), a novel neural network-based method that transforms PDE boundary value problems into unconstrained optimization tasks, enabling grid-free, analytical solutions.

Contribution

The work combines neural networks with TFC to solve PDEs without domain discretization, providing a closed-form, differentiable solution across the entire domain.

Findings

Successfully applied to four PDE problems with various boundary conditions.

Produces a continuous, analytical approximation of solutions after training.

Avoids grid-based discretization and interpolation used in traditional methods.

Abstract

This article presents a new methodology called deep Theory of Functional Connections (TFC) that estimates the solutions of partial differential equations (PDEs) by combining neural networks with TFC. TFC is used to transform PDEs with boundary conditions into unconstrained optimization problems by embedding the boundary conditions into a "constrained expression." In this work, a neural network is chosen as the free function, and used to solve the now unconstrained optimization problem. The loss function is taken as the square of the residual of the PDE. Then, the neural network is trained in an unsupervised manner to solve the unconstrained optimization problem. This methodology has two major differences when compared with popular methods used to estimate the solutions of PDEs. First, this methodology does not need to discretize the domain into a grid, rather, this methodology randomly samples points from the domain during the training phase. Second, after training, this methodology represents a closed form, analytical, differentiable approximation of the solution throughout the entire training domain. In contrast, other popular methods require interpolation if the estimated solution is desired at points that do not lie on the discretized grid. The deep TFC method for estimating the solution of PDEs is demonstrated on four problems with a variety of boundary conditions.

Figures

Equations

$$\begin{gathered} m{y}^{\prime \prime }+ky=0\text{subject to:} \\ \{\begin{array}{l}y(0)=y_0\\ y^{\prime }(0)=y_0^{\prime }\end{array} \end{gathered}$$

$$f(\text{x})={\mathcal{M}}_{i_1,i_2,\dots ,i_n}(c(\text{x}))v_{i_1}(x_1)v_{i_2}(x_2)\dots v_{i_n}(x_n)+$$

$$g(\text{x})-{\mathcal{M}}_{i_1,i_2,\dots ,i_n}(g(\text{x}))v_{i_1}(x_1)v_{i_2}(x_2)\dots v_{i_n}(x_n)$$

$$\mathcal{M}(c(x_1))=\left[\begin{array}{cccc}c(0)& c_{x_1}(0)& c(1)& c_{x_1}(1)\end{array}\right].$$

$$\mathcal{M}(c(x_1))$$

$$v$$

$$f(x_{1})=p^{1}(x_{1})\Big{(}b(0)-g(0)\Big{)}+p^{2}(x_{1})\Big{(}b_{x_{1}}(0)-g_{x_{1}}(0)\Big{)}+g(x_{1})$$

$$b(0)=p^{1}(0)\Big{(}b(0)-g(0)\Big{)}+p^{2}(0)\Big{(}b_{x_{1}}(0)-g_{x_{1}}(0)\Big{)}+g(0)$$

$$b_{t}(0)=p_{x_{1}}^{1}(0)\Big{(}b(0)-g(0)\Big{)}+p^{2}_{x_{1}}(0)\Big{(}b_{x_{1}}(0)-g_{x_{1}}(0)\Big{)}+g_{x_{1}}(0)$$

$$f(x_{1})=g(x_{1})+y_{0}-g(0)+t\Big{(}y^{\prime}_{0}-g_{x_{1}}(0)\Big{)}$$

$$\mathcal{L}(x_{1})=\Big{(}mf^{\prime\prime}(x_{1})+kf(x_{1})\Big{)}^{2}.$$

$${}^{k}b_{p}^{d}[f(\textbf{x})]=\frac{\partial^{d}f(\textbf{x})}{\partial x_{k}^{d}}\bigg{|}_{x_{k}=p}$$

$${\mathcal{M}}_{1,\dots ,1,i_k,1,\dots ,1}=\left\{\begin{array}{c}0{,}^k{c}_{{\text{p}}_k}^{{\text{d}}_k}\end{array}\right\}.$$

$$\displaystyle\mathcal{M}_{i_{1}11}=\big{[}0,c(0,x_{2},x_{3}),c(1,x_{2},x_{3})\big{]}^{\tiny{T}}$$

$$\displaystyle\mathcal{M}_{1i_{2}1}=\big{[}0,c(x_{1},0,x_{3}),c_{x_{2}}(x_{1},0,x_{3})\big{]}^{\tiny{T}}$$

$$\displaystyle\mathcal{M}_{11i_{3}}=\big{[}0,c(x_{1},x_{2},0),c_{x_{3}}(x_{1},x_{2},0)\big{]}^{\tiny{T}}.$$

$${\cal M}_{i_{1}i_{2}\dots i_{n}}=\,^{1}b^{\textbf{d}^{1}_{i_{1}-1}}_{\textbf{p}^{1}_{i_{1}-1}}\bigg{[}\,^{2}b^{\textbf{d}^{2}_{i_{2}-1}}_{\textbf{p}^{2}_{i_{2}-1}}\bigg{[}\dots\bigg{[}\,^{n}b^{\textbf{d}^{n}_{i_{n}-1}}_{\textbf{p}^{n}_{i_{n}-1}}[c(\textbf{x})]\bigg{]}\dots\bigg{]}\bigg{]}(-1)^{m+1},$$

$$M_{133}=-c_{x_2{x}_3}(x_1,0,0)$$

$$M_{221}=-c(0,0,x_3)$$

$$M_{332}=c_{x_2}(1,0,0)$$

$${\text{v}}_k=\left\{\begin{array}{ccccc}1,& \sum _{i=1}^{{\ell }_k}{\alpha }_{i1}{h}_i(x_k),& \sum _{i=1}^{{\ell }_k}{\alpha }_{i2}{h}_i(x_k),& \dots ,& \sum _{i=1}^{{\ell }_k}{\alpha }_{i{\ell }_k}{h}_i(x_k)\end{array}\right\},$$

$$\left[\begin{array}{cccc}{}^k{b}_{p_1}^{d_1}[h_1]& {}^k{b}_{p_1}^{d_1}[h_2]& \dots & {}^k{b}_{p_1}^{d_1}[h_{{\ell }_k}]\\ {}^k{b}_{p_2}^{d_2}[h_1]& {}^k{b}_{p_2}^{d_2}[h_2]& \dots & {}^k{b}_{p_2}^{d_2}[h_{{\ell }_k}]\\ ⋮& ⋮& \ddots & ⋮\\ {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}[h_1]& {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}[h_2]& \dots & {}^k{b}_{p_{{\ell }_k}}^{d_{{\ell }_k}}[h_{{\ell }_k}]\end{array}\right]\left[\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& \dots & {\alpha }_{1{\ell }_k}\\ {\alpha }_{21}& {\alpha }_{22}& \dots & {\alpha }_{2{\ell }_k}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{{\ell }_{k}1}& {\alpha }_{{\ell }_{k}2}& \dots & {\alpha }_{{\ell }_k{\ell }_k}\end{array}\right]=\left[\begin{array}{cccc}1& 0& \dots & 0\\ 0& 1& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1\end{array}\right].$$

$$\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]⟹\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]$$

$$v_1=\left\{\begin{array}{ccc}1& 1-x_1& x_1\end{array}\right\}.$$

$${\nabla }^{2}z(x,y)=e^{-x}(x-2+y^3+6y)\text{subject to:}$$

$$\{\begin{array}{l}z(x,0)=c(x,0)=x{e}^{-x}\\ z(0,y)=c(0,y)=y^3\\ z(x,1)=c(x,1)=e^{-x}(x+1)\\ z(1,y)=c(1,y)=(1+y^3)e^{-1}\end{array}$$

$$\text{where}(x,y)\in [0,1]\times [0,1]$$

$${\mathcal{M}}_{i_{1}1}$$

$${\mathcal{M}}_{1i_2}$$

$${\mathcal{M}}_{22}$$

$${\mathcal{M}}_{32}$$

$${\mathcal{M}}_{i_1{i}_2}=\left[\begin{array}{ccc}0& c(0,y)& c(1,y)\\ c(x,0)& -c(0,0)& -c(1,0)\\ c(x,1)& -c(0,1)& -c(1,1)\end{array}\right]=\left[\begin{array}{ccc}0& y^3& (1+y^3)e^{-1}\\ x{e}^{-x}& 0& -e^{-1}\\ e^{-x}(x+1)& -1& -2e^{-1}\end{array}\right]$$

$$\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]⟹\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]$$

$$v_{i_1}=\left\{\begin{array}{ccc}1& 1-x& x\end{array}\right\}.$$

$$\left[\begin{array}{cc}1& 0\\ 1& 1\end{array}\right]\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]⟹\left[\begin{array}{cc}{\alpha }_{11}& {\alpha }_{12}\\ {\alpha }_{21}& {\alpha }_{22}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right]$$

Full text

1 Introduction

Partial differential equations are a powerful mathematical tool that is used to model physical phenomena, and their solutions are used in the design and verification processes of a variety of systems. For example, the solution of PDEs may be used to verify that the loads experienced by a beam throughout its lifetime will not cause the beam to fail. Many methods exist to approximate the solutions of PDEs. The most famous of these methods is the finite element method (FEM) [1, 2, 3]. FEM has been incredibly successful in approximating the solution to PDEs in a variety of fields including structures, fluids, and acoustics. However, FEM does have some drawbacks.

FEM discretizes the domain into elements. This works well for low-dimensional cases, but the number of elements grows exponentially with the number of dimensions. Therefore, the discretization becomes prohibitive as the number of dimensions increases. Another issue is that FEM solves the PDE at discrete nodes, but if the solution is needed at locations different than these nodes, an interpolation scheme must be used.

Reference [4] explored the use of neural networks to solve PDEs, and showed that the use of neural networks avoids these problems. Rather than discretizing the entire domain into a number of elements that grows exponentially with the dimension, neural networks can sample points randomly from the domain. Moreover, once the neural network is trained, it is a closed form, analytical, differentiable approximation of the PDE. Therefore, no interpolation scheme is needed when estimating the solution at points that did not appear during training, and further analytical manipulation of the solution can be done with ease. Furthermore, Ref. [4] compared the neural network method with FEM on a set of test points that did not appear during training (i.e. points that were not the nodes in the FEM solution), and showed that the solution obtained by the neural network generalized well to points outside of the training data. In fact, the maximum error on the test set of data was never more than the maximum error on the training set of data. In contrast, the FEM had more error on the test set than on the training set. In one case, the test set had approximately 3 orders of magnitude more error than the training set.

Reference [4] is not the first article to explore the use of neural networks to solve PDEs. One of the early papers on the topic was Ref. [5]. Although Ref. [4] improved what was presented in [5] in almost every way, one item that was missing was exact satisfaction of the boundary constraints. The solution for PDEs presented by Ref. [5] used a method similar to the Coon’s Patch [6] to satisfy the boundary constraints exactly.

Exact boundary constraint satisfaction is of interest for a variety of problems, especially when the boundary value information is known to a high degree of precision. This is especially useful for physics informed problems. Take for example, using the heat equation to describe the temperature within a rod. If the temperature on the edges of the rod is know to a high degree of precision, and the heat equation is being used to estimate the temperature inside the rod, then the only information that is known for sure a priori is the boundary conditions. Hence, it is desired to have the boundary conditions of the rod met exactly. Moreover, embedding the boundary conditions means that the neural network only needs to sample points from the interior of the domain, not the domain and the boundary. Thus, a method for embedding boundary conditions for PDEs of arbitrary dimension is desired. Luckily, a framework for embedding constraints into machine learning algorithms has already been invented: the Theory of Functional Connections (TFC) [7]. In Ref. [7], TFC was used to embed constraints into support vector machines, but left embedding constraints into neural networks to future work. This research shows how to embed constraints into neural networks, and leverages this technique to numerically estimate the solutions of PDEs.

TFC is a framework that is able to satisfy many types of boundary conditions while maintaining a function that can be freely chosen. This free function can be chosen, for example, to minimize the residual of a differential equation. TFC has already been used to solve differential equations with initial value constraints, boundary value constraints, relative constraints, integral constraints, and infinite constraints [8, 9, 10, 11]. Recently, the framework was extended to $n$ dimensions [12] for constraints on the value and flux of $n-1$ dimensional manifolds. This means the framework can now generate constrained expressions that satisfy the boundary constraints of multidimensional PDEs [13].

2 Theory of Functional Connections

The Theory of Functional Connections (TFC) is a mathematical framework designed to turn constrained problems into unconstrained problems. This is accomplished through the use of the constrained expression, which is a mathematical expression that represents the family of all possible functions that satisfy the problem’s constraints. This technique is especially useful when solving PDEs, as it transforms the PDE into an unconstrained optimization problem. TFC has two major steps: 1) embed the boundary conditions of the problem into the constrained expression 2) solve the now unconstrained optimization problem. The paragraphs that follow will explain these steps in more detail.

The TFC framework is easiest to understand when explained via an example, like a simple harmonic oscillator. Equation (1) gives an example of a simple harmonic oscillator problem.

[TABLE]

Step one of the TFC framework is to build the constrained expression. The constrained expression consists of two parts; the first part is a function that satisfies the boundary constraints, and the second part projects the free-function, $g(\textbf{x})$, onto the hyper-surface of functions that are equal to zero at the boundaries. For problems with boundary and derivative constraints, the constrained expression can be written compactly using Eq. (2).

[TABLE]

where x is a vector of the independent variables, $\begin{Bmatrix}x_{1}&x_{2}&\dots x_{n}\end{Bmatrix}^{\tiny{T}}$, $\mathcal{M}$ is an $n$-th order tensor containing the boundary conditions $c(\textbf{x})$, $v_{i_{1}}\dots v_{i_{n}}$ are vectors whose elements are functions of the independent variables, $g(\textbf{x})$ is the free-function that can be chosen to be anything111More specifically, the free function $g(\textbf{x})$ can be chosen to be any function as long as that function is defined at the constraints., and $f(\textbf{x})$ is the constrained expression. The first term in $f(\textbf{x})$ satisfies the boundary conditions, and the final two terms project the free-function, $g(\textbf{x})$, onto the space of functions that vanish at the constraints.

For one dimension with independent variable $x_{1}$ on the domain $[0,1]$ and for boundary and first-derivative constraints only, $\mathcal{M}$ has the following form,

[TABLE]

The values of $\mathcal{M}$ that are unused are eliminated, and the vector, $v$, is created afterwards with the appropriate size. For the one dimensional simple harmonic oscillator, $\mathcal{M}$ and $v$ are given by Eq. (3).

[TABLE]

where $p^{1}({x_{1}})$ and $p^{2}(x_{1})$ are functions that need to be solved for. Substituting everything into the constrained expression, $f(x_{1})$, and simplifying yields Eq. (4).

[TABLE]

The functions $p_{1}(x_{1})$ and $p_{2}(x_{1})$ are solved for by setting the constrained expression, $f(x_{1})$, equal to the boundary constraints at the boundaries. For the simple harmonic oscillator this yields the set of simultaneous equations given by Eqs. (5) and (6).

[TABLE]

[TABLE]

Equation (5) shows that $p^{1}(0)=1$ and $p^{2}(0)=0$, and Eq. (6) shows that $p^{1}_{x_{1}}(0)=0$ and $p^{2}_{x_{1}}(0)=1$. One solution to this set of simultaneous equations is $p^{1}(x_{1})=1$ and $p^{2}(x_{1})=x_{1}$. Substituting these results back into the constrained expression and substituting $b(0)=y_{0}$ and $b_{x_{1}}(0)=\dot{y}_{0}$ yields Eq. (7).

[TABLE]

Notice that for any function $g(x_{1})$, the boundary conditions will always be satisfied exactly. Therefore, solving the differential equation has now become an unconstrained optimization problem. This unconstrained optimization problem could be cast in the following way. Let the function to be minimized, $\mathcal{L}$, be equal to the square of the residual of the differential equation,

[TABLE]

This function is to be minimized by varying the function $g(x_{1})$. One way to do this is to make $g(x_{1})$ a linear combination of a series of basis functions, and minimize the coefficients that multiply those basis functions via least-squares or some other optimization technique. For an example of this implementation via Chebyshev Orthogonal Polynomials see [9]. Also, [10] shows how this can be done using non-linear least squares for non-linear ODEs.

2.1 $n$-Dimensional Constrained Expressions

The previous example derived the constrained expression by creating and solving a series of simultaneous algebraic equations. This technique works well for constrained expressions in one dimension; however, a different, more mechanized formalism exists that is useful for deriving these expression in $n$ dimensions for constraints on the value and flux of $n-1$ dimensional manifolds [12]. This formalism shows how to construct the two portions of the constrained expression shown in Eq. (2): 1) $n$-th order tensor $\mathcal{M}$ 2) the $n$ vectors, $v$.

Before discussing how to build the $\mathcal{M}$ tensor and $v$ vectors, let’s introduce some mathematical notation. Let $k\in[1,n]$ be an index used to denote the $k$-th dimension. Let {}^{k}c_{p}^{d}:=\frac{\partial^{d}c(\textbf{x})}{\partial x_{k}^{d}}\bigg{|}_{x_{k}=p} be the constraint specified by taking the $d$-th derivative of the constraint function, $c(\textbf{x})$, evaluated at the $x_{k}=p$ hyperplane. Further, let ${}^{k}c_{\textbf{p}_{k}}^{\textbf{d}_{k}}$ be the vector of $\ell_{k}$ constraints defined at the $x_{k}=\textbf{p}_{k}$ hyperplanes with derivative orders of $\textbf{d}_{k}$, where $\textbf{p}_{k}$ and $\textbf{d}_{k}\in\mathcal{R}^{\ell_{k}}$. In addition, let’s define a boundary condition operator ${}^{k}b_{p}^{d}$ that takes the $d$-th derivative with respect to $x_{k}$ of a function, and then evaluates that function at the $x_{k}=p$ hyperplane. Mathematically,

[TABLE]

This mathematical notation will be used to introduce a step-by-step method for building the $\mathcal{M}$ tensor. This step-by-step process will be be illustrated via a 3-dimensional example that has Dirichlet boundary conditions in $x_{1}$ and initial conditions in $x_{2}$ and $x_{3}$ on the domain $x_{1},x_{2},x_{3}\in[0,1]\times[0,1]\times[0,1]$. The $\mathcal{M}$ tensor is constructed using the following three rules.

The element $\mathcal{M}_{111}=0$. 2. 2.

The first order sub-tensor of $\mathcal{M}$ specified by keeping one dimension’s index free and setting all other dimension’s indices to 1 consists of the value [math] and the boundary conditions for that dimension. Mathematically,

[TABLE]

Using the example boundary conditions,

[TABLE] 3. 3.

The remaining elements of the $\mathcal{M}$ tensor are those with at least two indices different from one. These elements are the geometric intersection of the boundary condition elements of the first order tensors given in Eq. (2), plus a sign ($+$ or $-$) that is determined by the number of elements being intersected. Mathematically this can be written as,

[TABLE]

where $m$ is the number of indices different than one. Using the example boundary conditions we give three examples:

[TABLE]

A simple procedure also exists for constructing the $v$ vectors. The $v$ vectors have a standard form:

[TABLE]

where $h_{i}(x_{k})$ are $\ell_{k}$ linearly independent functions. The simplest set of linearly independent functions, and those most often used in the TFC constrained expressions, are monomials, $h_{i}(x_{k})=x_{k}^{i-1}$. The $\ell_{k}\times\ell_{k}$ coefficients, $\alpha_{ij}$, can be computed by matrix inversion,

[TABLE]

Using the example boundary conditions, let’s derive the $v_{1}$ vector using the linearly independent functions $h_{1}=1$ and $h_{2}=x_{1}$.

[TABLE]

[TABLE]

For more examples and a proof that these procedures for generating the $\mathcal{M}$ tensor and $v$ vectors form a constrained expression see Ref. [12].

2.2 Two-Dimensional Example

This subsection will give an in depth example for a two-dimensional TFC case. The example is originally from [5] problem 5, and is one of the PDE problems analyzed in this article. The problem is shown in Eq. (9).

[TABLE]

Following the step-by-step procedure given in the previous section we will construct the $\mathcal{M}$ tensor:

The first element is $\mathcal{M}_{11}=0$. 2. 2.

The first order sub-tensors of $\mathcal{M}$ are:

[TABLE] 3. 3.

The remaining elements of $\mathcal{M}$ are the geometric intersection of elements from the first order sub-tensors.

[TABLE]

Hence, the $\mathcal{M}$ tensor is,

[TABLE]

Following the step-by-step procedure given in the previous section we will construct the $v$ vectors. For $v_{i_{1}}$ let’s choose the linearly independent functions $h_{1}=1$ and $h_{2}=x$.

[TABLE]

[TABLE]

For $v_{i_{2}}$ let’s choose the linearly indpendent functions $h_{1}=1$ and $h_{2}=y$.

[TABLE]

[TABLE]

Now, we use the constrained expression form given in Eq. (2) to finish building the constrained expression.

[TABLE]

Notice, that Eq. (10) will always satisfy the boundary conditions of the problem regardless of the value of $g(x,y)$. Thus, the problem has been transformed into an unconstrained optimization problem where the cost function, $\mathcal{L}$, is the square of the residual of the PDE,

[TABLE]

For one-dimensional ODEs, the minimization of the cost function was accomplished by making $g$ the summation of orthogonal polynomials, and performing least-squares or some other optimization technique to find the coefficients that multiply those orthogonal polynomials. For two dimensions, one could make $g(x,y)$ the product of two sums of these orthogonal polynomials, calculate all of the cross-terms, and then solve for the coefficients that multiply all terms and cross-terms using least-squares or non-linear least-squares. However, this will become prohibitive as the dimension increases. An alternative solution, and the one explored in this article, is to make the free function, $g(x,y)$, a neural network.

3 PDE Solution Methodology

As with the previous section, the easiest way to describe the methodology is with an example. The example used throughout this section will be the PDE given in Eq. (9).

As mentioned in the previous section, deep TFC approximates solutions to PDEs by finding the constrained expression for the PDE and choosing a neural network as the free function. For all of the problems analyzed in this article, a simple, fully connected neural network was used. These networks consist of non-linear activation functions composed with affine transformations of the form $\mathcal{A}=W\cdot x+b$, where $W$ are neuron weights, $b$ are neuron biases, and $x$ is a vector of inputs from the previous layer (or the inputs to the neural network if it is the first layer). The weights and biases of the entire neural network make up the tunable parameters, $\theta$. In this article, neural networks will be represented with the symbol $\mathcal{N}$. For example, a neural network with inputs $x$ and $y$ would be given as $\mathcal{N}(x,y;\theta)$. Thus, the constrained expression, given originally in Eq. (10), now has the form given in Eq. (11).

[TABLE]

In order to estimate the solution to the PDE, the parameters of the neural network have to be optimized to minimize the loss function, which is taken to be the square of the residual of the PDE. For this example,

[TABLE]

The attentive reader will notice that training the neural network will require, for this example, taking two second order partial derivatives of $f(x,y)$ to calculate $\mathcal{L}_{i}$, and then taking gradients of $\mathcal{L}$ with respect to the neural network parameters, $\theta$, in order to train the neural network.

To take these higher order derivatives, TensorFlow’s™ gradients function was used [14]. This function uses automatic differentiation [15] to compute these derivatives. However, one must be conscientious when using the gradients function to ensure they get the correct gradients.

When taking the gradient of a vector, $y_{j}$, with respect to another vector, $x_{i}$, TensorFlow™ computes,

[TABLE]

where $z_{i}$ is a vector of the same size as $x_{i}$. The only place where this may be an issue in the example used in this section is when computing $\nabla^{2}f_{i}$. The desired output of this calculation is the following vector,

[TABLE]

where $z_{i}$ has the same size as $f_{i}$ and $(x_{i},y_{i})$ is the pair used to generate $f_{i}$. TensorFlow’s™ gradients function will compute the following vector,

[TABLE]

However, because $f_{i}$ only depends on $(x_{i},y_{i})$ and the derivative operator commutes with the sum operator, TensorFlow’s™ gradients function will compute the desired vector. Moreover, the size of the output vector will be correct because the input vectors, $x_{i}$ and $y_{i}$, have the same size as $f_{i}$.

3.1 Training the Neural Network

Three methods were tried when optimizing the parameters of the neural networks:

Adam optimizer [16] 2. 2.

Broyden-–Fletcher–-Goldfarb-–Shanno (BFGS) optimizaion algorithm [17] 3. 3.

Hybrid method

The first method, Adam, is a variant of stochastic gradient descent (SGD) that combines the advantages of two other popular SGD variants: AdaGrad [18] and RMSProp [19].

The second method, BFGS, is a quasi-Newton method designed for solving unconstrained, non-linear optimization problems. This method was chosen based on its performance when optimizing neural network parameters to estimate PDE solutions in Ref. [5]. The hybrid method uses both methods in series. For all four problems shown in this article, the BFGS optimizer gave the best results.

4 Results

This section compares the estimated solution found using deep TFC with with the analytical solution. Four PDE problems are analyzed. The first is the example PDE given in Eq. (9), and the second is the wave equation. The third and fourth PDEs are simple solutions to the incompressible Navier-Stokes equations.

4.1 Problem 1

The first problem analyzed was the PDE given by Eq. (9), copied below for the readers’ convenience.

[TABLE]

The neural network used to estimate the solution to this PDE was a fully connected neural network with five hidden layers and 30 neurons per layer. The non-linear activation function used was the hyperbolic tangent. Other neural network sizes and non-linear activation functions were tried, but this size and activation function combination did the best. The biases of the neural network were all initialized as zeros, and the weights were initialized using TensorFlow’s™ implementation of the Xavier initialization with uniform random initialization [20]. Training pairs, $(x,y)$, were created for this problem by sampling $x$ and $y$ as independent and identically distributed (IID) random variables on the uniform distribution on $[0,1]$. The network was trained using the BFGS method on a batch size of 10,000 training pairs.

Figure 1 shows the difference between the analytical solution and the estimated solution using deep TFC on a grid of 100 evenly distributed points (10 per independent variable). This grid represents the test set.

The maximum error on the test set was $4.806\times 10^{-6}$ meters and the average error was $8.872\times 10^{-7}$ meters. The maximum error is relatively low, five orders of magnitude lower than the solution values, which are on the order of $10^{-1}$ meters. However, the maximum solution error using deep TFC is not as low as the maximum solution error obtained in [5], which was on the order of $10^{-7}$ meters. Although there are many possible explanations for this discrepancy, for example, different implementations of the optimizer or different weight initialization functions, the author was concerned that it may be the assumed solution form, $f(x,y)$, that was causing the increase in estimated solution error. The solution form created using deep TFC is more complex both in the number of terms and the number of times the neural network appears within the assumed solution form.

To investigate, a comparison was made between the solution form posed in [5] and the solution form posed in this article, while keeping all other variables constant. For this comparison, the neural network architecture consisted of one hidden layer with 10 neurons, and the non-linear activation function used was the hyperbolic tangent function222The sigmoid function was also tried, which is what is used in [5]; however, the average error and maximum error for both solution forms, deep TFC and [5], was almost doubled compared to using the hyperbolic tangent function.. The neural network was trained using the BFGS optimizer. The training pairs were created by randomly sampling 10,000 point pairs $(x,y)$ where $x$ and $y$ are IID random variables on the uniform distribution $[0,1]$. The test set was a grid of 100 evenly distributed points (10 per independent variable).

Figure 2 was created using the solution form posed in [5]. The maximum error on the test set was $5.481\times 10^{-6}$ meters and the average error on the test set was $1.131\times 10^{-6}$ meters.

Figure 3 was created using the deep TFC solution form. The maximum error on the test set was $8.124\times 10^{-6}$ meters and the average error on the test set was $1.696\times 10^{-6}$ meters.

Comparing Figs. 2 and 3 shows that the solution form from [5] does slightly better in terms of average error and maximum error for this problem, but the difference is not significant. However, neither of the tests conducted here is able to reduce the error to the level reported in [5]. Thus, there must be some other factor, initialization, optimizer method, etc., that is causing this difference. One important take-away from this comparison is that despite the more complex solution form created using deep TFC, which can easily be applied to higher dimensions, the final solution accuracy is very similar to the simpler solution form that was designed for lower-dimensional problems.

4.2 Problem 2

The second problem analyzed was the wave equation, shown in Eq. (12).

[TABLE]

where the constant, $c$, was chosen to be 1. The constrained expression for this problem is shown in Eq. (13).

[TABLE]

The neural network used to estimate the solution to this PDE was a fully connected neural network with two hidden layers and 30 neurons per layer. The non-linear activation function used was the hyperbolic tangent. The biases and weights were initialized using the same method as problem 1. The training points, $(x,t)$, were created by sampling $x$ and $t$ IID from the uniform distribution on $[0,1]$. The network was trained using the BFGS method on a batch size of 10,000 training pairs.

Figure 4 shows the difference between the analytical solution and the estimated solution using deep TFC on a grid of 100 evenly distributed points (10 per independent variable). This grid represents the test set.

The maximum error on the test set was $2.701\times 10^{-3}$ meters and the average error on the test set was $6.978\times 10^{-4}$ meters. The error of this solution is larger than in the problem one, while the solution values are on the same order of magnitude, $10^{-1}$ meters, as in problem one. The larger relative error in problem two is most likely due to the more oscillatory nature of the solution.

4.3 Problem 3

The third problem analyzed was a known solution to the incompressible Navier-Stokes equations, called Poiseuille flow. The problem solves the flow velocity in a two-dimensional pipe in steady-state with a constant pressure gradient applied in the longitudinal axis. Equation (14) shows the associated equations and boundary conditions.

[TABLE]

where $u$ and $v$ are velocities in the $x$ and $y$ directions respectively, $H$ is the height of the channel, $P$ is the pressure, $\rho$ is the density, and $\mu$ is the viscosity. For this problem, the values $H=1$ m, $\rho=1$ kg/m3, $\mu=1$ Pa$\cdot$s, and $\frac{\partial P}{\partial x}=-5$ N/m3 were chosen. The constrained expressions for the $u$-velocity, $f^{u}(x,y,t)$, and $v$-velocity, $f^{v}(x,y,t)$, are shown in Eq. (15).

[TABLE]

The neural network used to estimate the solution to this PDE was a fully connected neural network with four hidden layers and 30 neurons per layer. The non-linear activation function used was the sigmoid. The biases and weights were initialized using the same method as problem 1. The training points, $(x,y,t)$, were created by sampling $x$, $y$, and $t$ IID from the a uniform distribution that spanned the range of the associated independent variable. For $x$, the range was $[0,1]$. For $y$, the range was $[-\frac{H}{2},\frac{H}{2}]$, and for $t$, the range was $[0,1]$. The network was trained using the BFGS method on a batch size of 1,000 training pairs. The loss function used was the sum of the squares of the residuals of the three PDEs in Eq. (14).

For one particular run333All errors are given in meters per second., the maximum error in the $u$-velocity was $3.308\times 10^{-7}$, the average error in the $u$-velocity was $9.998\times 10^{-8}$, the maximum error in the $v$-velocity was $5.575\times 10^{-7}$, and the average error in the $v$-velocity was $1.542\times 10^{-7}$. The maximum error and average error for this problem are much lower than in problems 1 and 2. However, the constrained expression for this problem essentially encodes the solution, because the initial flow condition at time zero is the same as the flow condition throughout the spatial domain at any time. Thus, if the neural network outputs a value of zero for all inputs, the problem will be solved exactly. Although the neural network does output a very small value for all inputs, it is interesting to note that none of the layers have weights or biases that are at or near zero.

4.4 Problem 4

The fourth problem is another solution to the Navier-Stokes equations, and is very similar to the third. The only difference is that in this case, the fluid is not in steady state, it starts from rest. Equation (16) shows the associated equations and boundary conditions.

[TABLE]

This problem was created to avoid encoding the solution to the problem into the constrained expression, as was the case in the previous problem. The constrained expressions for the $u$-velocity, $f^{u}(x,y,t)$, and $v$-velocity, $f^{v}(x,y,t)$, are shown in Eq. (17).

[TABLE]

The neural network used in this problem is exactly the same as the neural network used in problem three. Problem 4 used 2,000 training points that were selected the same way as in problem three, except the new ranges for the independent variables were $[0,10]$ for $x$, $[0,3]$ for $t$, and $[-\frac{H}{2},\frac{H}{2}]$ for $y$.

Figures 7 through 7 show the $u$-velocity of the fluid throughout the domain at three different times. Qualitatively, the solution should look as follows. The solution should be symmetric about the line $y=0$, and the solution should develop spatially and temporally such that after a period of time and sufficiently far from the inlet, $x=0$, the $u$-velocity will be equal or very close to the steady state $u$-velocity of problem 3. Qualitatively, the $u$-velocity field looks correct at all points. Quantitatively, the $u$-velocity at $x=10$ from Fig. 7 was compared with the known steady state $u$-velocity, and had a maximum error of $5.530\times 10^{-4}$ meters per second and an average error of $2.742\times 10^{-4}$ meters per second.

5 Conclusions

This article demonstrated how to combine neural networks with the Theory of Functional Connections into a new methodology, called deep TFC, that was used to estimate the solutions of PDEs. Results on four problems were presented that display how accurately relatively simple neural networks can approximate the solutions to some well known PDEs. The difficulty of the PDEs in these problems ranged from linear, two-dimensional PDEs to coupled, non-linear, three-dimensional PDEs. Moreover, while the focus of this article was on PDEs, the capability to embed constraints into neural networks is useful for many other applications as well.

Future work should investigate the performance of different neural network architectures on the estimated solution error. For example, Ref. [4] suggests a neural network architecture where the hidden layers contain element-wise multiplications and sums of sub-layers. The sub-layers are more standard neural network layers like the fully connected layers used in the neural networks of this article.

Another topic for investigation is reducing the estimated solution error by sampling the training points based on the loss function values for the training points of the previous iteration. For example, one could create batches where half of the new batch consists of half of the points in the previous batch that had the largest loss function value and the other half are randomly sampled from the domain. This should consistently give training points that are in portions of the domain where the estimated solution is farthest from the real solution.

Acknowledgements

Carl Leake would like to acknowledge the support of the NASA Space Technology Research Fellowship (NSTRF) 2019: grant number 80NSSC19K1152.