# Chebyshev Spectral Neural Networks for Solving Partial Differential
Equations

Pengsong Yin  School of Mathematical Sciences, Shanghai Jiao Tong University,
Shanghai 200240, P.R. China  Shuo Ling  School of Mathematical Sciences,
Shanghai Jiao Tong University, Shanghai 200240, P.R. China  Wenjun Ying
Corresponding author.
Email addresses<EMAIL_ADDRESS>School of Mathematical Sciences, MOE-LSC
and Institute of Natural Sciences, Shanghai Jiao Tong University, Minhang,
Shanghai 200240, P.R. China.

###### Abstract

The purpose of this study is to utilize the Chebyshev spectral method neural
network (CSNN) model to solve differential equations. This approach employs a
single-layer neural network wherein Chebyshev spectral methods are used to
construct neurons satisfying boundary conditions. The study uses a feedforward
neural network model and error backpropagation principles, utilizing automatic
differentiation (AD) to compute the loss function. This method avoids the need
to solve non-sparse linear systems, making it convenient for algorithm
implementation and solving high-dimensional problems. The unique sampling
method and neuron architecture significantly enhance the training efficiency
and accuracy of the neural network. Furthermore, multiple networks enables the
Chebyshev spectral method to handle equations on more complex domains. The
numerical efficiency and accuracy of the CSNN model are investigated through
testing on elliptic partial differential equations, and it is compared with
the well-known Physics-Informed Neural Network(PINN) method.

Keywords: Chebyshev spectral method neural network, differential equations,
single-layer neural network, automatic differentiation, complex domains

##  1 Introduction

Partial differential equations (PDEs) play a crucial role in various fields of
engineering and science. Many problems in mathematics, physics, engineering,
economics, and other domains can be modeled using differential equations [19].
In most cases, obtaining analytical solutions for partial differential
equations is challenging. Therefore, various numerical methods such as finite
difference method (FDM) [22], finite element method (FEM) [18], finite volume
method (FVM) [14], etc., are employed to solve these equations. While the
aforementioned numerical methods provide good approximate solutions, they can
be challenging for high-dimensional and complex domain problems.

Computational methods based on machine learning and neural networks have been
widely researched in recent years. In fact, the increasing computational power
and abundant data have made it possible to apply neural networks to
increasingly complex problems. Furthermore, the versatility and simplicity of
neural networks greatly solidify their powerful tool status in various
domains. Meanwhile, the structure of neural networks is often reformulated to
adapt to the diverse environmental requirements necessitated by different
problems.

###  1.1 Artificial Neural Networks for PDEs

Recently, a class of Artificial Neural Networks(ANN) has been introduced,
known as Physics-Informed Neural Networks (PINN) [17, 20]. The training of
PINN is achieved by minimizing the loss function, which, in addition to
observational data, incorporates the governing equations, initial conditions,
and boundary conditions. As is well-known, neural networks possess universal
approximation capabilities [6, 26]. Therefore, following this process,
approximate solutions to initial/boundary value problems may be effective.

Compared to traditional numerical methods, ANNs offer several advantages in
providing approximate solutions. For example, such methods generally exhibit
strong versatility, capable of solving both linear and nonlinear ordinary
differential equations as well as partial differential equations.
Additionally, many neural network approaches demonstrate certain superiority
in solving high-dimensional and general domain problems.

While traditional neural networks have achieved remarkable results in solving
differential equations, there are some drawbacks associated with using
conventional neural network structures. Issues such as low solution accuracy
and poor generality need to be overcome, and training may be challenging in
certain situations, especially when the solution of the equation exhibits
high-frequency characteristics [24]. Therefore, many scientists are exploring
approaches to address these challenges by incorporating mathematical prior
knowledge, modifying network structures, training methods, and other
techniques.

###  1.2 Single Layer Artificial Neural Network Model

Susmita Mall and S. Chakraverty proposed a fast learning single-layer neural
network called Functional Link Artificial Neural Network (FLANN) to solve
differential equations [8, 9]. In FLANN the hidden layer is replaced by a
functional expansion block for enhancement of the input patterns using
Chebyshev polynomials that is computationally more efficient than the multi-
layer perceptron network [16, 26]. While this single-layer Chebyshev network
exhibits extremely fast optimization speed and, to some extent, has a mature
approximation theory [10], using such a network structure alone for any PDE is
not ideal. An important factor is that it does not take into account the need
to satisfy the boundary conditions or initial conditions required for solving
the equation.

To align the network structure more closely with the inherent properties of
the equation, this paper proposes a new neural network solution for PDE based
on the Chebyshev spectral method [4]. In the subsequent sections, we will
refer to it as the Chebyshev spectral neural network(CSNN). Its advantage lies
in its naturally satisfying the boundary conditions of the equation, as it
avoids sampling at the boundaries, thereby reducing errors near the
boundaries. Additionally, due to the excellent properties of spectral methods,
such as exponential convergence for smooth functions, achieving high accuracy
with fewer parameters of networks, and absence of phase errors [21], this
method exhibits significant improvements in solution accuracy. In some cases,
it even surpasses Physics-Informed Neural Networks (PINN) [17].

###  1.3 Spectral Methods for Complex Domain Problems

For the traditional spectral methods, addressing complex domain problems is
crucial to expanding the application scope. Current mainstream solutions can
be broadly categorized into two types: coordinate transformation and
partitioning methods. Orszag introduced the concept of linear transformations
[12]. Wang introduced general coordinates to spectral methods[25], making it
possible to compute for arbitrary shapes. Partitioning methods distribute the
global nature of spectral methods to different regions to adapt to more
complex geometries. Morchoishne first proposed spectral multi-domain methods
[11]. Patera introduced spectral element methods [15], where spectral
functions replace polynomials on each element, utilizing variational
principles for finite element approximations. Although there are various
methods available for solving complex domain problems, solving non-sparse
matrices or intricate domain partitions remains a computationally expensive
and challenging issue. CSNN employs variable substitution, transforming the
polar coordinate set region that the equation needs to solve into a
rectangular region for sampling and training. It uses an additional neural
network to learn the transformed boundary values to achieve the solution, thus
somewhat simplifying the computational complexity. While the form of the
equation may undergo transformation, the powerful learning and generalization
capabilities of neural networks enable effective learning even for nonlinear
equations with low error. The subsequent chapters will provide a detailed
explanation of this method.

The outline of this paper is as follows. Section 2 introduces some prior
knowledge about the Chebyshev spectral method. Section 3 \- 4 present the
specific structure of the CSNN proposed by us and the algorithm for solving
PDEs using it. Following that, section 5 provides a certain error analysis.
Section 6 is the numerical experiments section, validating the feasibility and
efficiency of our method. The conclusion and discussion part is in section 7.

##  2 Preliminaries

This section introduces some knowledge about the traditional Chebyshev
spectral method, serving as the foundation for CSNN.

###  2.1 Chebyshev Polynomials

The Chebyshev differential equation yields a set of orthogonal polynomials
under the integral inner product on $[-1,1]$ with the specific weight function
$w(x)=(1-x^{2})^{-\frac{1}{2}}$, which are known as the Chebyshev polynomials.
The first two Chebyshev polynomials are known as

$T_{0}(x)=1,$ $T_{1}(x)=x.$

The well-known recursive formula can be employed to generate higher-order
Chebyshev polynomials [1]

$T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x),\qquad n\geq 1,$ (2.1)

where $T_{n}(x)$ denotes n-th order Chebyshev polynomial and $-1<x<1$ is the
argument of the polynomials.

An important feature of the Chebyshev polynomials is that Chebyshev-Gauss-
Lobatto quadrature points and wights can be expressed as follows

$x_{0}=1,\quad x_{N}=-1,\quad x_{j}=\cos{\frac{\pi j}{N}},\quad 1\leq j\leq
N-1,$ (2.2)
$\omega_{0}=\omega_{N}=\frac{\pi}{2N},\qquad\omega_{j}=\frac{\pi}{N},\quad
1\leq j\leq N-1.$ (2.3)

For Gauss quadrature defined with $\\{x_{j},\omega_{j}\\}_{j=0}^{N}$, define a
discrete inner product in $C[a,b]$ and its associated norms:

$(u,v)_{N,\omega}=\sum_{j=0}^{N}u(x_{j})v(x_{j})\omega_{j},\qquad\|u\|_{N,\omega}=(u,u)_{N,\omega}^{\frac{1}{2}}$
(2.4)

###  2.2 Basis Functions for Spectral Method

To introduce the construction method of basis functions, we exemplify with a
(one-dimensional) two-point boundary-value problem. Without loss of
generality, we consider the equation with homogeneous boundary conditions

$\displaystyle a_{-}u(-1)+b_{-}u^{\prime}(-1)$ $\displaystyle=0,$ (2.5)
$\displaystyle\qquad a_{+}u(1)+b_{+}u^{\prime}(1)$ $\displaystyle=0.$

For if the right-hand side is not zero, it can be always homogenized(see in
Appendix A). Our objective is to identify a basis function of the following
structure that satisfies (2.5)

$\phi_{k}(x)=T_{k}(x)+a_{k}T_{k+1}(x)+b_{k}T_{k+2}(x).$ (2.6)

Since $T_{k}(\pm 1)=(\pm 1)^{k}$ and $T_{k}^{{}^{\prime}}(\pm
1)=(\pm)^{k-1}k^{2}$, substitute (2.6) into (2.5) we have

$\left\\{\begin{aligned}
&\left(a_{+}+b_{+}(k+1)^{2}\right)a_{k}+\left(a_{+}+b_{+}(k+2)^{2}\right)b_{k}=-a_{+}-b_{+}k^{2}\\\
&-\left(a_{-}-b_{-}(k+1)^{2}\right)a_{k}+\left(a_{-}-b_{-}(k+2)^{2}\right)b_{k}=-a_{-}+b_{-}k^{2}\end{aligned}\right.$

Thus, we obtain the expressions for the computation of $a_{k}$ and $b_{k}$.

$\displaystyle a_{k}=$
$\displaystyle-\left\\{\left(a_{+}+b_{+}(k+2)^{2}\right)\left(-a_{-}+b_{-}k^{2}\right)\right.$
(2.7)
$\displaystyle\left.-\left(a_{-}-b_{-}(k+2)^{2}\right)\left(-a_{+}-b_{+}k^{2}\right)\right\\}/\text{
DET }_{k},$ $\displaystyle b_{k}=$
$\displaystyle\left\\{\left(a_{+}+b_{+}(k+1)^{2}\right)\left(-a_{-}+b_{-}k^{2}\right)\right.$
$\displaystyle\left.+\left(a_{-}-b_{-}(k+1)^{2}\right)\left(-a_{+}-b_{+}k^{2}\right)\right\\}/\text{
DET }_{k},$

with

$\mathrm{DET}_{k}=2a_{+}a_{-}-2b_{+}b_{-}(k+1)^{2}(k+2)^{2}+(a_{-}b_{+}-a_{+}b_{-})[(k+1)^{2}+(k+2)^{2}].$

Specifically, for the homogeneous Dirichlet conditions, substituting
$a_{k}=1$, $b_{k}=0$, $c_{k}=-1$ into (2.7), the basis functions can be easily
obtained as

$\phi_{i}(x)=T_{i}(x)-T_{i+2}(x)$ (2.8)

For one-dimensional problems, we can directly take $\phi_{k}(x)$ as the basis
for the approximation function space

$\mathcal{X}_{N}=span\\{\phi_{i}(x):i=0,1,\dots,N-2\\}.$ (2.9)

For a higher-dimensional cube (or hypercube) domain, we can compute polynomial
basis functions for each coordinate corresponding to the homogeneous boundary
conditions required on each dimension. Then, we can take the product of these
basis functions to form new high-dimensional basis functions. For instance, an
approximation function space for two-dimensional rectangle problem with
homogeneous could be

$\mathcal{X}_{N}=span\\{\phi_{i}(x)\phi_{j}(y):i,j=0,1,\dots,N-2\\}.$ (2.10)

###  2.3 PDEs in Complex Geometries

In this subsection, let us discuss how to solve equations in non-rectangular
domains. Taking 2D problem as examples, assuming that we already have the
polar coordinate representation for a complex domain in $\mathbb{R}^{2}$:

$\displaystyle\mathrm{f}_{1}(\theta)\leq\mathrm{r}\leq\mathrm{f}_{2}(\theta),$
(2.11) $\displaystyle 0\leq\theta<2\pi,$

where $(r,\theta)$ are polar coordinates which can be transformed into the
rectangle

$\displaystyle-1\leq z\leq 1,$ (2.12) $\displaystyle 0\leq\theta<2\pi,$

by the simple stretching transformation

$z=2\frac{r-f_{1}(\theta)}{f_{2}(\theta)-f_{1}(\theta)}-1.$ (2.13)

The complexity of employing a coordinate transformation arises in the
coefficients of the differential equation within the transformed domain[12].
Specifically, derivatives undergo transformation as per

$\displaystyle\left.\frac{\partial u}{\partial
r}\right|_{\theta}=\left.\frac{2}{f_{2}(\theta)-f_{1}(\theta)}\frac{\partial
u}{\partial z}\right|_{\theta}$ $\displaystyle\left.\frac{\partial
u}{\partial\theta}\right|_{r}=\left.\frac{\partial
u}{\partial\theta}\right|_{z}-\left.\frac{\left[z\left(f_{2}^{\prime}-f_{1}^{\prime}\right)-\left(f_{2}^{\prime}+f_{1}^{\prime}\right)\right.}{f_{2}-f_{1}}\frac{\partial
u}{\partial z}\right|_{\theta}.$

It is worth noting that nonlinear variable substitutions can sometimes
introduce unexpected nonlinearities and singularities into the resulting
differential equations. However, thanks to the powerful generalization and
learning capabilities of neural networks, smaller errors can still be achieved
when solving these differential equations. This will be discussed in detail in
the Numerical Experimental sections.

##  3 Chebyshev Spectral Method Neural Network(CSNN)

Figure 1 illustrates the structure of CSNN, which consists of an input layer
with d inputs (input $\mathbf{x}\in\mathbb{R}^{d}$, in this figure and later
discussion, we illustrate with the case where $d=2$), followed by neurons
formed by basis functions based on Chebyshev polynomials, and a single output.

xy Chebyshev Expansion
$\phi_{0}(x)$$\phi_{N-2}(x)$$\psi_{0}(y)$$\psi_{N-2}(y)$$\vdots$$\vdots$$\phi_{0}(x)\cdot\psi_{0}(y)$$\phi_{0}(x)\cdot\psi_{N-2}(y)$$\phi_{N-2}(x)\cdot\psi_{0}(y)$$\phi_{N-2}(x)\cdot\psi_{N-2}(y)$$\vdots$$\vdots$output$\omega_{0,0}$$\omega_{0,N-2}$$\omega_{N-2,0}$$\omega_{N-2,N-2}$
Figure 1: The structure of CSNN model(only the blue lines contain the neural
network weight parameters that need to be trained).

The CSNN model comprises three parts: The first part transforms the non-
homogeneous problem into the corresponding homogeneous problem; the second
part computes the Chebyshev basis functions that satisfy the conditionsas
neurons; the third part discusses the learning and optimization algorithms of
the neural network.

The first two parts have been elaborated in detail in Section 2.2. It should
be noted that if the input elements satisfy the two-point boundary
conditions(e.g., $z$ in (2.12)), the basis functions can directly take the
form of equation (2.6). If the input elements are periodic elements introduced
by variable substitution (e.g., $\theta$ in (2.12)), due to their inherent
periodicity, the spectral expansion is more suitable. We can choose the basis
by $\phi_{k}(\theta)=\cos{k\theta}+\sin{k\theta}$, namely

$\mathcal{\theta}_{N}=span\\{\cos(i\theta):i=0,1,\dots,N-2\\}\cup
span\\{\sin(i\theta):i=0,1,\dots,N-2\\}.$ (3.1)

Experiments show that this choice indeed reduces the error more effectively
than using Chebyshev basis functions for all cases directly.

The learning algorithm is employed to update the network parameters and
minimize the loss function. In this case, the error backpropagation algorithm
is utilized to update the weights of the CSNN. To accomplish this, the
gradient of an loss function with respect to the network parameters
$\mathbf{p}$ is computed.

The network output with input data $\mathbf{x}$ and parameters (weights)
$\mathbf{p}$ can be computed as:

$\mathcal{N}(\mathbf{x};\mathbf{p})=\sum_{j=0}^{N-2}\sum_{i=0}^{N-2}w_{ij}\phi_{i}(x)\psi_{j}(y).$

Where $\mathbf{x}=\left(x,y\right)$ represents the input data,
$\phi_{i}(x)\psi_{j}(y)$ and $w_{ij}$ denote the basis function in (2.10) and
the weight vectors of CSNN, respectively.

Upon observation, it is evident that the function approximation space obtained
using a single-layer CSNN neural network resembles that of spectral methods.
In the subsequent numerical experiments, we also employ a single-layer network
for experimentation. In fact, numerical results indicate that the single-layer
CSNN exhibits high accuracy.

##  4 Solving PDEs by CSNN

The general form of representing ordinary differential equations (ODEs) or
partial differential equations (PDEs) for the solution u($\mathbf{x}$) defined
on a domain $\Omega\subset\mathbb{R}^{d}$ is as follows:

$f(\mathbf{x},u(\mathbf{x}),\nabla
u(\mathbf{x}),\nabla^{2}u(\mathbf{x}),\dots,\nabla^{n}u(\mathbf{x}))=0\qquad\mathbf{x}\in\Omega$
(4.1)

where $f$ is a function that defines the structure of differential equation,
$u(X)$ and $\nabla$ denote the solution and differential operator
respectively. With suitable boundary conditions

$\mathcal{B}(u,\mathbf{x})=0\qquad on\quad\partial\Omega$ (4.2)

where $\mathcal{B}(u,\mathbf{x})$ could be Dirichlet, Neumann, Robin boundary
conditions.

In the following discussion in section 4, we will take $d=2$ as an example
while the case for higher dimensions is similar.

###  4.1 PDEs in Rectangular Domains

We begin by considering a simple rectangular domain. We first create a CSNN
$\hat{u}(\mathbf{x};\mathbf{\theta})$ as an approximation to the solution
$u({\mathbf{x}})$. Then, we can take the derivatives of
$\hat{u}(\mathbf{x};\mathbf{\theta})$ with respect to the input $\mathbf{x}$
via automatic differentiation (AD) [13].

In fact, due to the specific structure of the CSNN, the network’s output
automatically satisfies the boundary conditions of the equation. Therefore, we
only need to calculate the residual of the equation within the domain, without
further considering the boundary conditions. This not only reduces the
computational load of the loss function but also avoids the need to adjust the
weights of the boundary and internal residuals as in traditional PINN methods.

As for sampling, if the training basis functions are Chebyshev
polynomials(2.6), then $N-2$ Chebyshev-Gauss-Lobatto quadrature points (2.2)
are taken in each dimension. If the training basis functions are Fourier
polynomials, uniform sampling is directly used. This is a significant
difference between CSNN and PINN. This sampling method greatly reduces the
training workload.

We denote the sampling set as
$\mathcal{T}=\\{\mathbf{x_{1}},\mathbf{x_{2}},\dots,\mathbf{x_{\mathcal{|T|}}}\\}$.
To quantify the difference between the neural network $u({\mathbf{x}})$, we
define the loss function as the weighted summation of the $L^{2}$ norm of
residuals for the equation

$\mathcal{L}(\mathbf{\theta};\mathcal{T})=\frac{1}{|\mathcal{T}|}\sum_{\mathbf{x}\in\mathcal{T}}\|f(X,\hat{u}(X;\mathbf{\theta}),\nabla\hat{u}(X;\mathbf{\theta}),\nabla^{2}\hat{u}(X;\mathbf{\theta}),\dots,\nabla^{n}\hat{u}(X;\mathbf{\theta}))\|_{2}^{2}$
(4.3)

Finally, by minimizing the loss function, we find the most suitable parameters
$\mathbf{\theta}$. Considering that the loss with respect to $\mathbf{\theta}$
is highly nonlinear and non-convex [2], we typically utilize gradient-based
optimizers such as gradient descent, Adam [7], and L-BFGS [5] to minimize the
loss function.

In conclusion, the process of solving PDEs on rectangular domains using the
CSNN method is illustrated in Algorithm 1.

Algorithm 1 CSNN on rectangular domains

1: For each dimension whose boundary condition is not of the type of periodic,
homogenize boundary conditions into (2.5);

2: For each dimension whose boundary condition is not of the type of periodic,
compute the values of $a_{k},b_{k},c_{k}$ using (2.7);

3: Construct the Chebyshev spectral method neural Network using the computed
$\phi_{i}$ (2.6);

4: Compute the residual points $\mathcal{T}$ in each dimension;

5: Train the neural network for a limited number of iterations with the loss
function $\mathcal{L}(\mathbf{\theta};\mathcal{T})$ given by (4.3).

###  4.2 PDEs in Complex Domains

In this section, we will introduce how to solve PDEs on general parameterized
domains, which is one of the major advantages of the CSNN method over
traditional spectral methods and also one of our highlights.

####  4.2.1 PDEs After Coordinate Transformation on Variables

For more complex polar coordinate domains, when the inner boundary of the
original region has prescribed values (such as in the case of an annular
region (2.11)), the transformation poses no extra issues. This is because the
boundary values on the rectangular region obtained from the transformation are
determined by the inner and outer boundary conditions of the original region.
Therefore, apart from the change in the form of the original equation due to
the variable substitution and a new periodic boundary condition arise by
$\theta$ direction, all other aspects are identical to solving PDEs on a
rectangular domain(e.g. (2.12)).

However, when the original domain only has an outer boundary, such as

$\displaystyle\mathrm{r}\leq\mathrm{f}_{2}(\theta),$ (4.4) $\displaystyle
0\leq\theta<2\pi.$

The values on one side of the transformed rectangular region are determined by
the original region’s boundary conditions, while the values on the other side
become undetermined unknowns (Figure 3 illustrates the original domain of the
equation. Figure 3 illustrates the parameter domain in polar coordinates. The
rectangular orange border corresponds to the boundary values of the original
equation, while the green border corresponds to the unknown solution at the
center of the original equation.)

Figure 2: the original domain

Figure 3: the parameter domain(with $r$ on the horizontal axis and $\theta$ on
the vertical axis)

Traditional spectral methods may find it challenging to solve in such cases.
In CSNN, we can treat the boundary unknowns $u_{0}$ as an additional parameter
$\hat{u_{0}}$ which can be trained during the training process together with
$\hat{u}(\mathbf{x};\theta)$. In addition, in the next subsection, we will
introduce a correction technique for $u_{0}$ to assist in adjusting its value
during the training process.

####  4.2.2 Correction for $\hat{u_{0}}$ Network

Through multiple numerical experiments, we have found that after several
iterations, the primary error in CSNN’s solution to the differential equation
arises from the discrepancy between the output values of the $\hat{u_{0}}$
network and the true solution at the origin. This error stems from never
sampling the boundary points during training. We need to make certain
correction to the $\hat{u_{0}}$ network to ensure that the equation is as
close as possible to satisfying our original equation at the origin.

Predefine an iteration threshold value $M$. When the number of iterations
exceeds $M$, we perform a correction to $u$ every 100 iterations. Take 2D
Possion equation as an example, we employ a five-point finite difference
scheme for the correction

$-\frac{u(h,0)-2u(0,0)+u(-h,0)}{h^{2}}-\frac{u(0,h)-2u(0,0)+u(0,-h)}{h^{2}}\approx
f(0,0),$ (4.5)

where $h$ is a very small positive number. Thus, we obtain the correction form
for the unknown value $u_{0}$ at the origin as

$\bar{u_{0}}=\frac{\hat{u}(h,0)+\hat{u}(-h,0)+\hat{u}(0,h)+\hat{u}(0,-h)}{4}+h^{2}f(0,0).$
(4.6)

Here, $\hat{u}$ is the CSNN at the current stage, and $\bar{u}$ is the
corrected prediction. In fact, we only need to update the value of $\bar{u}$
to the parameter of the single-parameter $\hat{u_{0}}$ network.

In conclusion, the process of solving PDEs on general parametric domains using
the CSNN method is illustrated in Algorithm 2.

Algorithm 2 CSNN on general parametric domains

1: By variable substitution, the PDEs on general domains are transformed into
new PDEs on rectangular domains;

2: Determining the boundary conditions for the new PDEs and creating
additional trainable parameters $\hat{u_{0}}$ caused by a new periodic
boundary condition;

3: Run Algorithm 1 on the new PDEs defined on the rectangular domain, and
correct $\hat{u_{0}}$ every certain number of training steps.

##  5 Approximation Theory and Error Analysis for CSNN

Let $\mathcal{F}$ denote the family of all the functions that can be
represented by the architecture of CSNN, and u represent the true solution of
the equation. We define
$u_{\mathcal{T}}=arg\min_{f\in\mathcal{F}}\mathcal{L}(f;\mathcal{T})$ as the
neural network whose loss is at global minimum on the trainning set
$\mathcal{T}$. Since finding $u_{\mathcal{T}}$ through minimizing the loss is
typically computationally challenging, we denote that our optimizer returns an
approximate solution $\tilde{u_{\mathcal{T}}}$. Then, the total error
$\mathcal{E}$ can be represented as [3]

$\mathcal{E}:=\|\tilde{u_{\mathcal{T}}}-u\|\leq\|\tilde{u_{\mathcal{T}}}-u_{\mathcal{T}}\|+\|u_{\mathcal{T}}-u\|.$
(5.1)

To simplify the problem, we assume that when the training samples are large
enough, the global minimum point can be reached using the optimization
algorithm of the neural network.i.e.

$\|\tilde{u_{\mathcal{T}}}-u_{\mathcal{T}}\|\leq\epsilon$ (5.2)

holds for any arbitrarily small $\epsilon$.

Here, for the approximation of $\|u_{\mathcal{T}}-u\|$ we only consider the
one-dimensional case. Suppose we need to solve the following model problem

$\displaystyle-u^{{}^{\prime\prime}}+\alpha u=f,\qquad in\quad I=(-1,1)$ (5.3)
$\displaystyle u(1)=u(-1)=0$

In this case, $\mathcal{F}=\\{v\in P_{N}:v(-1)=v(1)=0\\}$. According to the
assumptions, we have

$\displaystyle\alpha
u_{\mathcal{T}}(x_{j})-u_{\mathcal{T}}^{{}^{\prime\prime}}(x_{j})=f(x_{j}),\qquad
j=1,2,\dots,N-1$ (5.4) $\displaystyle
u_{\mathcal{T}}(x_{0})=u_{\mathcal{T}}(x_{N})=0$

where $x_{j}=\cos(\frac{j\pi}{N})$. We can rewrite (5.4) in a variational
formulation with discrete inner product (2.4) [23]:

$\alpha(u_{\mathcal{T}},v_{\mathcal{T}})_{N,\omega}-(u_{\mathcal{T}}^{{}^{\prime\prime}},v_{\mathcal{T}})_{N,\omega}=(I_{N}f,v_{\mathcal{T}})_{N,\omega}\qquad\forall
v_{\mathcal{T}}\in\mathcal{F}$ (5.5)

where $I_{N}:C(-1,1)\to P_{N}$ is the interpolating operator associated with
Gauss-Lobatto points.

According to [23], for $u\in H_{\omega^{-\frac{3}{2},-\frac{3}{2},*}}^{m}(I)$
and $f\in H_{\omega^{-\frac{1}{2},-\frac{1}{2},*}}^{k}(I)$, we have

$\displaystyle\left\|u-u_{\mathcal{T}}\right\|_{1,\omega}\lesssim
N^{1-m}\left\|\partial_{x}^{m}u\right\|_{\omega^{m-3/2,m-3/2}}+N^{-k}\left\|\partial_{x}^{k}f\right\|_{\omega^{k-\frac{1}{2},k-\frac{1}{2}}},\quad
m,k\geqslant 1,$ (5.6)
$\displaystyle\left\|u-u_{\mathcal{T}}\right\|_{\omega}\lesssim
N^{-m}\left\|\partial_{x}^{m}u\right\|_{\omega^{m-1/2,m-1/2}}+N^{-k}\left\|\partial_{x}^{k}f\right\|_{\omega^{k-\frac{1}{2},k-\frac{1}{2}}},\quad
m,k\geqslant 1.$

where $\omega$ represents the Chebyshev weight $(1-x^{2})^{-\frac{1}{2}}$, and
the definition of the Jacobi-weighted Sobolev space can be found in Appendix
B.

It indicates that if we ignore the error introduced by the neural network
optimization algorithm, as long as the solution of the equation and the right-
hand side $f$ are sufficiently smooth and the width of the neural network is
wide enough, the approximate solution can converge to the true solution to any
order.

##  6 Numerical Experiments

In this section, based on the aforementioned CSNN method, we validate the
effectiveness and generality of the method by solving differential equations
from one to four dimensions. All numerical experiments are conducted using the
NVIDIA GeForce RTX 3060(Laptop) GPU, and the Adam optimizer is used to train
the neural network. We begin with several simple examples.

We evaluate the accuracy by using the absolute $L^{\infty}$ error, absolute
$L^{2}$ error and relative $L^{2}$ error defined as follows:

$\|e\|_{absolute-L^{\infty}}=\max_{i=1,\cdots,K}|\hat{u}(x_{i})-u^{*}(x_{i})|,$
$\|e\|_{absolute-L^{2}}=\sqrt{\frac{\sum_{i=1}^{K}|\hat{u}(x_{i})-u^{*}(x_{i})|^{2}}{K}},$
$\|e\|_{relative-L^{2}}=\frac{\sqrt{\sum_{i=1}^{K}|\hat{u}(x_{i})-u^{*}(x_{i})|^{2}}}{\sqrt{\sum_{i=1}^{K}|u^{*}(x_{i})|^{2}}},$

where $\hat{u}$ is the neural network prediction, $u^{*}$ is the true
solution, and $K$ is the number of test points.

###  6.1 The 1D Robin boundary condition differential equation

Consider the following problem:

$u{{}^{\prime\prime}}(x)+xu^{{}^{\prime}}(x)-u(x)=(24+5x)e^{5x}+(2+2x^{2})\cos(x^{2})-(4x^{2}+1)\sin(x^{2})$
(6.1)

with the Robin boundary conditions:

$u(-1)-u^{{}^{\prime}}(-1)=-4e^{-5}+\sin(1)+2\cos(1),\qquad
u(1)+u^{{}^{\prime}}(1)=6e^{5}+\sin(1)+2\cos(1).$ (6.2)

The exact solution is $u(x)=e^{5x}+\sin(x^{2})$.

We take $N=6,8,10,12,14$ to construct CSNN networks for computation. The epoch
size for the training process of each network is set as $300\times N$. The
initial learning rate is set to $0.1$. A learning rate scheduler with a
reduction factor of $0.7$ is utilized. This scheduler dynamically adjusts the
learning rate if the loss fails to decrease for 600 consecutive iterations.

The test points are uniformly sampled with 501 points in the interval from -1
to 1. The results are shown in Table 1.

$N$ | $6$ | $8$ | $10$ | $12$ | $14$
---|---|---|---|---|---
Training Time (seconds) | 10.44 | 17.38 | 23.51 | 31.30 | 41.21
Absolute $L^{\infty}$ Error | 7.38E-2 | 5.46E-3 | 2.85E-4 | 1.10E-5 | 1.01E-6
Table 1: Result for 1D Differential Equation

It is observed that as the number of CSNN sampling points increases, the error
of the solution decreases at a very high rate. In contrast, the slight
increase in sampling points for PINNs has minimal impact on the equation’s
solution.

###  6.2 The 2D Poisson equation on a rectangular domain

In this part, we solve the two-dimensional Poisson equation with Dirichlet
boundary conditions:

$-\Delta u(x,y)=f(x,y)\qquad(x,y)\in\Omega=[-1,1]^{2}.$ (6.3)

The right-hand side term $f$ of the equation are determined by the true
solution $u=e^{-x}\sin{(\pi y)}$. Let the boundary conditions of the equation
be

$\left\\{\begin{array}[]{l}\alpha_{0}(y)=u(-1,y),\\\ \alpha_{1}(y)=u(1,y),\\\
\beta_{0}(x)=u(x,-1),\\\ \beta_{1}(x)=u(x,1),\end{array}\right.$ (6.4)

Figure 4: the problem on the boundary

which Satisfy the continuity property

$\left\\{\begin{aligned} \alpha_{0}(-1)&=u(-1,-1)=\beta_{0}(-1),\\\
\alpha_{0}(1)&=u(1,1)=\beta_{1}(1),\\\
\alpha_{1}(-1)&=u(1,-1)=\beta_{0}(1),\\\
\alpha_{1}(1)&=u(1,1)=\beta_{1}(1).\end{aligned}\right.$ (6.5)

The trial solution is given by the expression,

$\psi(x,y)=A(x,y)+\tilde{u}(x,y;\mathbf{\theta}).$ (6.6)

Here, $\tilde{u}(x,y;\mathbf{\theta})$ is the output obtained by the CSNN
neural network, and $A(x,y)$ represents the construct function which satisfies
the non-homogeneous boundary conditions, where

$\displaystyle A(x,y)=$
$\displaystyle\left(\frac{1-x}{2}\alpha_{0}(y)+\frac{1+x}{2}\alpha_{1}(y)\right)$
(6.7)
$\displaystyle+\left(\frac{1-y}{2}\left(\beta_{0}(x)-\left(\frac{1-x}{2}\beta_{0}(-1)+\frac{1+x}{2}\beta_{0}(1)\right)\right)\right)$
$\displaystyle+\left(\frac{1+y}{2}\left(\beta_{1}(x)-\left(\frac{1-x}{2}\beta_{1}(-1)+\frac{1+x}{2}\beta_{1}(1)\right)\right)\right)$

The optimization algorithm chosen is the Adam, with a learning rate 0.1, and
the number of iterations set to 2000. The test points are on a uniformly
spaced grid of 300 cells on each dimension in the x and y space. The
$L^{\infty}$ error and training duration for different network complexities N
are shown in Table 2.

$N$ | $6$ | $8$ | $10$
---|---|---|---
Training Time (seconds) | $19.16$ | 26.18 | 38.67
Absolute $L^{\infty}$ Error | 4.466E-4 | 1.019E-5 | 8.345E-7
Table 2: Result for 2D Poisson Equations Figure 5: predict for 2D Poisson
equation(the left image shows the predicted values by the CSNN network, the
middle image displays the accurate solution of the equation, and the right
image represents the computed error).

A major advantage of the CSNN network is the selection of a large learning
rate and the rapid convergence speed. The computational results for N = 10 are
shown in Figure 5. The error distribution is relatively uniform across the
entire domain, with the error gradually decreasing near the boundary until it
approaches zero.

###  6.3 The 2D Poisson equation on an annular domain

In this part, we consider solving Poisson equation on the annular region
$\Omega=\\{(r,\theta)|0.5\leq r\leq 1$, $0\leq\theta\leq 2\pi\\}$. The
original equation is given as follows:

$-\Delta u(x,y)=f(x,y)\qquad(x,y)\in\Omega.$ (6.8)

The right-hand side term $f$ and the Dirichlet boundary conditions on the
inner and outer boundaries are both determined by the true solution
$u=e^{x}\cos{(\pi y)}$.

Take $x=\frac{s+3}{4}\cos((t+1)\pi)$, $y=\frac{s+3}{4}\sin((t+1)\pi)$, the
above equation (6.8) can be transformed into

$-\frac{16}{(s+3)}u_{s}-16u_{ss}-\frac{16}{{\left(\pi\cdot(s+3)\right)}^{2}}u_{tt}=f(\frac{{3+s}}{4}\cos\left((t+1)\pi\right),\frac{3+s}{4}\cos\left((t+1)\pi\right)),$
(6.9)

where $t,s\in[-1,1]$.

In the s-direction, we employ Chebyshev basis functions (as in (2.9)), and the
sampling points are still chosen as the Lobatto-Gauss-Legendre points (2.2)
(with N=10). For the t-direction, considering its periodicity, we use
trigonometric basis functions (as in (3.1)), and the sampling points range
uniformly from -1 to 1 (with 500 sampling points).

We denote the sampling set in the parameter space as
$\mathcal{T}=\\{\mathbf{x_{1}},\mathbf{x_{2}},\dots,\mathbf{x_{\mathcal{|T|}}}\\}$.
Constructing the loss function

$\mathcal{L}(\mathbf{\theta};\mathcal{T})=\frac{1}{|\mathcal{T}|}\sqrt{\sum_{\mathbf{x}\in\mathcal{T}}|\frac{16}{(s+3)}u_{s}+16u_{ss}+\frac{16}{{\left(\pi\cdot(s+3)\right)}^{2}}u_{tt}+f|^{2}}$
(6.10)

The initial learning rate is set to $0.01$, and a learning rate scheduler with
a reduction factor of $0.6$ is utilized. This scheduler dynamically adjusts
the learning rate if the loss fails to decrease for 500 consecutive
iterations. As the training progresses, the loss function steadily
decreases(see in Figure 6). It can be observed that even with a relatively
small initial learning rate, the CSNN network converges quickly.

Figure 6: loss for 2D Poisson Equations on annular Domain

The test points are on a uniformly spaced grid of 300 cells on each dimension
in the s and t space. The solution of the equation is depicted in Figure 7 and
the absolute $L^{\infty}$ error reaches $9.07\times 10^{-4}$ after 2000
training iterations.

Figure 7: The 2D Poisson equation on an annular domain(the left image shows
the predicted values by the CSNN network, the middle image displays the
accurate solution of the equation, and the right image represents the computed
error.)

###  6.4 The 2D Helmholtz equation on an annular domain

We consider the 2D Helmholtz equation on the $\Omega=\\{(r,\theta)|0.5\leq
r\leq 1$, $0\leq\theta\leq 2\pi\\}$

$\Delta u(x,y)+u(x,y)=f(x,y).$ (6.11)

We test our method with $f$ and the Dirichlet boundary condition given by the
exact solution $u(x,y)=e^{x}\sin(\pi y)$. Traditional methods for solving the
2D Helmholtz equation may encounter some difficulties. In finite difference,
finite element, higher wave numbers may lead to numerical dissipation or
numerical waveguide issues.

In this example, we take $N=13$, and the training parameters and test samples
are the same as the previous subsection. The final numerical results are shown
in Figure 8 whose absolute $L_{\infty}$ error is $1.0192\times 10^{-5}$.

Figure 8: The 2D Helmholtz equations on an annular domain(the left image shows
the predicted values by the CSNN network, the middle image displays the
accurate solution of the equation, and the right image represents the computed
error.)

To compare the efficiency and accuracy differences between CSNN and PINNs, we
applied the PINNs method to solve the examples in this section. The network of
the PINNs method consists of 4 hidden layers, each with 128 neurons, using the
$\tanh$ activation function. After 2000 iterations with the same learning rate
$lr=0.001$, the $L_{\infty}$ error of the solution given by the PINNs method
is $0.014$, which is over 1000 times higher than CSNN. It can be observed that
CSNN method has a higher convergence speed and can achieve higher accuracy as
well. Simultaneously, during training, CSNN does not exhibit irregular
oscillations in the solution error like PINNs; instead, the error steadily
decreases (see in Fig 9).

Figure 9: Comparison with PINNs

###  6.5 The 2D Poisson equation on a circular domain

In this example, we consider solving the Poisson equation (6.8) on the
circular region $\Omega=\\{(r,\theta)|r\leq 1$, $0\leq\theta\leq 2\pi\\}$ with
exact solution $u(x,y)=e^{x}\cos(\pi y)$.

The initial learning rate is set to $lr=8\times 10^{-3}$, and a learning rate
scheduler with a reduction factor of $0.7$ is utilized. This scheduler
dynamically adjusts the learning rate if the loss fails to decrease for 800
consecutive iterations. Additionally, we require an extra single-parameter
neural network $\hat{u_{0}}$ to learn the value $u_{0}$ at the center of the
circle. We set the initial learning rate for $\hat{u_{0}}$ as $6\times
10^{-2}$ and decrease the step size by multiplying it by 0.8 if the loss fails
to decrease for 600 consecutive iterations.

The test points are on a uniformly spaced grid of 300 cells on each dimension
in the s and t space. After 10000 iterations of training, the $L^{2}$ error of
the solution can reach $9\times 10^{-4}$(see in Fig 10).

Figure 10: the 2D Poisson equation on a circular domain(the left image shows
the predicted values by the CSNN network, the middle image displays the
accurate solution of the equation, and the right image represents the computed
error.)

However, it must be acknowledged that the feasibility of the CSNN method is
strongly associated with the sizes of the learning rates for both networks.
Incorrect step size settings can lead to slow convergence of the networks or
even failure to reach the global minimum. To address this issue, we propose a
network correction for $\hat{u}$ (as discussed in Section 4.2.2), where we
apply adjustments to the network parameters to satisfy the equation conditions
within a certain number of convergence steps. This measure significantly
improves the resilience of the CSNN method to the learning rate settings.

###  6.6 The 3D variable coefficients elliptic PDE

For the three-dimensional case, we are interested in the numerical solution of
the following form of the elliptic problem

$\displaystyle\nabla\cdot(\sigma(\mathbf{x})\nabla
u)-\kappa(x)u=f(\mathbf{\mathbf{x}})\qquad$
$\displaystyle\text{in}\quad\Omega$ (6.12) $\displaystyle u=g^{D}$
$\displaystyle\text{on}\quad\partial\Omega$

with $\sigma=\cos(\pi(x+y+z))+2$ and $\kappa=\sin(\pi(x+y+z))+2$. The
Dirichlet boundary data $g(x,y,z)$ and source data $f(x,y)$ are chosen such
that the exact solution to the problem is given by
$u(x,y,z)=e^{\frac{2}{3}x+\frac{2}{3}y+\frac{1}{3}z}\sin(\pi(x+y+z))$. The
computational domain $\Omega$ is a rolled-up cylindrical shell with an inner
radius of 0.5, an outer radius of 1, and a height ranging from -1 to 1(see Fig
11).

Figure 11: The computational domain of Example 6.6

If we take $x=r\cos\theta$, $y=r\sin\theta$, $z=z$, the aforementioned domain
can be transformed into a standard rectangular region
$\Omega^{{}^{\prime}}=\\{(r,\theta,z)\in[0.5,1]\times[0,2\pi)\times[-1,1]\\}$.
And the differential operator is transformed into

$\nabla\cdot\left(\sigma\nabla
u(x,y,z)\right)=\frac{1}{r}\frac{\partial}{\partial r}(r\sigma\frac{\partial
u}{\partial
r})+\frac{1}{r^{2}}\frac{\partial}{\partial\theta}(\sigma\frac{\partial
u}{\partial\theta})+\frac{\partial}{\partial z}(\sigma\frac{\partial
u}{\partial z})$ (6.13)

The initial learning rate is set to $0.01$, and a learning rate scheduler with
a reduction factor of $0.7$ is utilized. This scheduler dynamically adjusts
the learning rate if the loss fails to decrease for 100 consecutive
iterations. $r$ and $z$ directions are sampled at $N$ CGL points, and the
$\theta$ direction is uniformly sampled with 100 points. The test points are
on a uniformly spaced grid of 128 cells on each dimension in the $r,\theta,z$
space. In this experiment, $N=8$, after training 78 seconds for 500
iterations, the $L_{\infty}$ error and $L_{2}$ error are $9.5367\times
10^{-7}$ and $1.9601\times 10^{-7}$, respectively. This demonstrates that even
when solving nonlinear PDEs in irregular domains, CSNN still maintains a high
level of accuracy.

As an example, for the cross-section at $r=0.7$, the numerical results and
$L^{\infty}$ errors on the cylindrical surface are shown in Figure 12.

Figure 12: The 3D variable coefficients elliptic(the left image shows the
predicted values by the CSNN network, the middle image displays the accurate
solution of the equation, and the right image represents the computed error.)

###  6.7 The 4D Poisson Equation

To illustrate the applicability of CSNN in higher dimensions, we solve the
Poisson equation on a four-dimensional hypercube, testing the speed and
accuracy of the method and comparing it with the traditional PINNs approach.
The Poisson equation is given as following

$\begin{cases}-\Delta u(\mathbf{x})=4\pi^{2}sin(\pi x_{1})sin(\pi
x_{2})sin(\pi x_{3})sin(\pi x_{4})&\mathbf{x}\in\Omega=[-1,1]^{4}\\\
u(\mathbf{x})=0&\mathbf{x}\in\partial\Omega\end{cases}.$ (6.14)

The exact solution is $u(\mathbf{x}=(x_{1},x_{2},x_{3},x_{4}))=sin(\pi
x_{1})sin(\pi x_{2})sin(\pi x_{3})sin(\pi x_{4})$ which can be verified
easily.

The result for our method on different values of $N$ after training 1000
epochs is given in Table 3. The test points are on a uniformly spaced grid of
$11$ cells on each dimension in the 4D hypercube.

$N$ | $6$ | $8$ | $10$
---|---|---|---
Training Time (seconds) | $79.20$ | 90.71 | 142.3
Absolute $L^{2}$ Error | 3.324E-2 | 7.419E-4 | 3.284E-5
Relative $L^{2}$ Error | 4.063E-2 | 9.069E-4 | 4.014E-5
Table 3: Result for 4D Poisson Equations

From the table, it can be observed that our method performs excellently in
solving this PDE, especially for $N=8$ and $10$. It demonstrates impressive
speed and accuracy in numerical solutions for PDEs in high-dimensional spaces.

In addition, it is worth mentioning that the performance of the PINNs method
in this case is catastrophic. Specifically, we employe a fully connected
neural network with 3 hidden layers, each containing 128 neurons. The hyper-
parameter $\beta$ for balancing PDE residual and boundary residual is set to
1, and the optimization scheme involves the widely used Adam optimizer. We
randomly sample 65,536 points within $\Omega$ every 100 training steps to
compute the corresponding $L^{2}$ norm inside the domain. Additionally, 1,024
points are randomly sampled along the boundaries to compute the corresponding
$L^{2}$ norm on the boundary. Despite such substantial computational expenses,
there is no corresponding improvement. Even after training for more than
10,000 iterations over 3,000 seconds, the loss function remains above 0.8.
This results in the neural network failing to adequately fit the desired
solution, with an absolute $L^{2}$ error reaching 0.88. This highlights the
superiority of our method.

##  7 Conclusion and Discussion

We propose a novel neural network framework for solving PDEs. Our experience
with this method suggests that it offers the following advantages:

  * 1)

High convergence rates and capability to fit high-frequency solutions.

  * 2)

Reduced sampling point requirements, eliminating the need for boundary
sampling points.

  * 3)

Training using only a single-layer neural network greatly reduces the number
of parameters, making convergence analysis feasible.

  * 4)

Effective computation of mixed boundary conditions.

  * 5)

Capable of solving equations in non-rectangular domains.

There are several disadvantages that need to be addressed in the future work:

  * 1)

In solving problems in complex domains, it is currently challenging to
effectively determine the learning rates for each neural network.

  * 2)

Inability to solve problems in non-parametrizable complex domains.

In future research, we will focus on solving more complex equations, such as
the Stokes equations and the elasticity equations. We will also attempt to
apply our methods to interface problems, hoping to achieve good numerical
results. Additionally, more flexible and complex boundaries are also of
interest to us. We hope to apply our method to any form of complex boundaries,
not just parametric boundaries.

## Acknowledgements

This work is financially supported by the Strategic Priority Research Program
of Chinese Academy of Sciences(Grant No. XDA25010405). It is also partially
supported by the National Key R & D Program of China, Project Number
2020YFA0712000, the National Natural Science Foundation of China (Grant No.
DMS-11771290) and the Science Challenge Project of China (Grant No.
TZ2016002).

## Appendix A

Assume that U is the solution to a two-point boundary-value problem satisfying
the following boundary conditions

$\displaystyle a_{-}U(-1)+b_{-}U^{\prime}(-1)$ $\displaystyle=c_{-},$ (1)
$\displaystyle\quad a_{+}U(1)+b_{+}U^{\prime}(1)$ $\displaystyle=c_{+}.$

This includes in particular the Dirichlet $\left(a_{\pm}=1\right.$ and
$\left.b_{\pm}=0\right)$, the Neumann $\left(a_{\pm}=0\right.$ and
$\left.b_{\pm}=1\right)$, and the mixed $\left(a_{-}=b_{+}=0\right.$ or
$a_{+}=b_{-}=0$ ) boundary conditions. The following explains how to
homogenize it.

  * •

Case $1\quad a_{\pm}=0$ and $b_{\pm}\neq 0$. We set $\bar{u}=\beta
x^{2}+\gamma x$, where $\beta$ and $\gamma$ are uniquely determined by asking
$\tilde{u}$ to satisfy (1)

$\displaystyle-2b_{-}\beta+b_{-}\gamma$ $\displaystyle=c_{-},$ (2)
$\displaystyle\quad 2b_{+}\beta+b_{+}\gamma$ $\displaystyle=c_{+}.$

  * •

Case $2\quad a_{-}^{2}+a_{+}^{2}\neq 0$. We set $\tilde{u}=\beta x+\gamma$,
where $\beta$ and $\gamma$ can again be uniquely determined by asking that
$\tilde{u}$ to satisfy (1)

$\displaystyle\left(-a_{-}+b_{-}\right)\beta+a_{-}\gamma$
$\displaystyle=c_{-},$ (3)
$\displaystyle\quad\left(a_{+}+b_{+}\right)\beta+a_{+}\gamma$
$\displaystyle=c_{+}.$

We now set $u=U-\tilde{u}$. Then $u$ satisfies the new equation with the
homogeneous boundary conditions

$\displaystyle a_{-}u(-1)+b_{-}u^{\prime}(-1)$ $\displaystyle=0,$ (4)
$\displaystyle\qquad a_{+}u(1)+b_{+}u^{\prime}(1)$ $\displaystyle=0.$

## Appendix B

Let $I=(-1,1)$ and the Jacobi weight function of index $(\alpha,\beta)$ by
$\omega^{\alpha,\beta}(x)=(1-x)^{\alpha}(1+x)^{\beta}$. We define the Jacobi
weighted Sobolev spaces:

$\displaystyle
L_{\omega^{\alpha,\beta}}^{2}(I)=\left\\{u:\int_{I}u^{2}\omega^{\alpha,\beta}\mathrm{d}x<+\infty\right\\},$
$\displaystyle H_{\omega^{\alpha,\beta}}^{l}(I)=\left\\{u\in
L_{\omega^{\alpha,\beta}}^{2}(I):\partial_{x}u,\cdots,\partial_{x}^{l}u\in
L_{\omega^{\alpha,\beta}}^{2}(I)\right\\},$ $\displaystyle
H_{\omega^{\alpha,\beta},*}^{m}(I):=\left\\{u:\partial_{x}^{k}u\in
L_{\omega^{\alpha+k,\beta+k}}^{2}(I),\quad 0\leqslant k\leqslant m\right\\},$

The norms in $L_{\omega^{\alpha,\beta}}^{2}(I)$ and
$H_{\omega^{\alpha,\beta}}^{l}(I)$ will be denoted by
$\|\cdot\|_{\omega^{\alpha,\beta}}$ and
$\|\cdot\|_{l,{\omega^{\alpha,\beta}}}$, respectively.
$H_{\omega^{\alpha,\beta},*}^{m}(I)$ equipped with the inner product and norm

$(u,v)_{m,\omega^{\alpha,\beta},*}=\sum_{k=0}^{m}\left(\partial_{x}^{k}u,\partial_{x}^{k}v\right)_{\omega^{\alpha+k,\beta+k}},\qquad\|u\|_{m,\omega^{\alpha,\beta},*}=(u,u)_{m,\omega^{\alpha,\beta},*}^{\frac{1}{2}}.$

## References

  * [1] RB Bhat and S Chakravety.  Numerical analysis in engineering. alpha science int ltd, 2004.
  * [2] Avrim Blum and Ronald Rivest.  Training a 3-node neural network is np-complete.  Advances in neural information processing systems, 1, 1988.
  * [3] Léon Bottou and Olivier Bousquet.  The tradeoffs of large scale learning.  Advances in neural information processing systems, 20, 2007.
  * [4] John P Boyd.  Chebyshev and Fourier spectral methods.  Courier Corporation, 2001.
  * [5] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  SIAM Journal on scientific computing, 16(5):1190–1208, 1995.
  * [6] Simon Haykin.  Neural networks: a comprehensive foundation.  Prentice Hall PTR, 1998.
  * [7] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [8] Susmita Mall and Snehashish Chakraverty.  Chebyshev neural network based model for solving lane–emden type equations.  Applied Mathematics and Computation, 247:100–114, 2014.
  * [9] Susmita Mall and Snehashish Chakraverty.  Single layer chebyshev neural network model for solving elliptic partial differential equations.  Neural Processing Letters, 45:825–840, 2017.
  * [10] John C Mason and David C Handscomb.  Chebyshev polynomials.  Chapman and Hall/CRC, 2002.
  * [11] Y Morchoisne.  Inhomogeneous flow calculations by spectral methods-mono-domain and multi-domain techniques.  NASA STI/Recon Technical Report A, 1982:14527, 1982.
  * [12] Steven A Orszag.  Spectral methods for problems in complex geometrics.  In Numerical methods for partial differential equations, pages 273–305. Elsevier, 1979.
  * [13] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * [14] Suhas Patankar.  Numerical heat transfer and fluid flow.  CRC press, 2018.
  * [15] Anthony T Patera.  A spectral element method for fluid dynamics: laminar flow in a channel expansion.  Journal of computational Physics, 54(3):468–488, 1984.
  * [16] Jagdish Chandra Patra.  Chebyshev neural network-based model for dual-junction solar cells.  IEEE Transactions on Energy Conversion, 26(1):132–139, 2010.
  * [17] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational physics, 378:686–707, 2019.
  * [18] Junuthula Narasimha Reddy.  An introduction to the finite element method.  New York, 27:14, 1993.
  * [19] Henry J Ricardo.  A modern introduction to differential equations.  Academic Press, 2020.
  * [20] Samuel Rudy, Alessandro Alla, Steven L Brunton, and J Nathan Kutz.  Data-driven identification of parametric partial differential equations.  SIAM Journal on Applied Dynamical Systems, 18(2):643–660, 2019.
  * [21] Jie Shen, Tao Tang, and Li-Lian Wang.  Spectral methods: algorithms, analysis and applications, volume 41.  Springer Science & Business Media, 2011.
  * [22] Gordon D Smith.  Numerical solution of partial differential equations: finite difference methods.  Oxford university press, 1985.
  * [23] Tao Tang.  Spectral and high-order methods with applications.  Science Press Beijing, 2006.
  * [24] Sifan Wang, Xinling Yu, and Paris Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  Journal of Computational Physics, 449:110768, 2022.
  * [25] Michiru Yasuhara, Yoshiaki Nakamura, and Jian-Ping Wang.  Computations of the hypersonic flow by the spectral method.  Journal of the Japan Society for Aeronautical and Space Sciences, 36:542–549, 01 1988.
  * [26] Jacek M Zurada.  Introduction to artificial neural network.  West Publishing Company, United States of America, 1992.