Approximating Partial Differential Equations with Physics-Informed Legendre Multiwavelets CNN

1. Introduction

PDEs, as a powerful mathematical tool in classical physics, enable the comprehensive depiction, simulation, and modeling of various dynamic phenomena, including heat flow, diffusion processes, wave propagation, and the motion of fluid substances [1], owing to their inherent capability to capture the fundamental characteristics of intricate systems. However, obtaining the analytical solutions of PDEs is challenging. To address this issue, conventional numerical methods like the finite volume method (FVM) [2], finite difference method (FDM) [3], and finite element method (FEM) [4], are utilized to approximate the solutions. Despite the strength and rigor of these approaches, they may experience an exponential rise in computational complexity with grid refinement and “dimension explosion” as the independent variable’s dimension increases.

With the rapid development of deep learning (DL), which is a form of machine learning (ML) that utilizes NNs with multiple hidden layers [5], has emerged as a new paradigm in scientific computing due to its universal approximation [6] and expressive capabilities. Solving PDEs through DL has gained significant momentum. Recent researches have shown that DL holds promise in building meta-models for efficiently predicting dynamic systems. By training deep neural networks (DNNs), especially with Artificial Intelligence (AI)-based approaches, systems or even families of PDEs can be expressed, leading to exponential advancements in computing efficiency during practical usage. Unlike traditional numerical methods, DNN-based approaches can simultaneously approximate target solutions in multiple dimensions, including the time dimension. They can learn direct mappings between physical states and spatial/temporal coordinates without the need for repeated iterations at each time step [7].

Currently, an increasing number of researchers are utilizing DL methods to investigate PDEs [6,8]. Among the various studies, an important contribution that cannot be overlooked is the PINNs proposed by Raissi et al. [9]. This NN model enhances the performance of the NN model by taking into account the physical laws found in PDEs and encoding them into the NN as regularization terms. Researchers have been paying more attention to PINNs lately, and they are gradually applying them to a wider range of research fields [10]. Although PINNs can be trained with little to no ground truth data, they often fail to effectively generalize to domains that were not encountered during training [11]. However, CNN can learn the inherent laws of the data and thus have better generalization ability for new data. Furthermore, the structure and characteristics of CNN facilitate the solution of PDE. For example, CNN can simplify the solution of complex PDE by reducing the dimension, leverages powerful nonlinear approximation capabilities to handle intricate initial and boundaries conditions (ICs/BCs), and, when integrated with parallel computing and GPU acceleration, significantly improves the computational efficiency of PDE solutions. Moreover, CNN can employ multi-input methods to manage multi-variable problems, and, through end-to-end learning, it minimizes the need for manual intervention and the complexity of manual settings [9,12]. Nevertheless, CNN relies on large amounts of training data.

In this paper, we propose PiLMWs-CNN for deriving continuous solutions that can be trained with a physics-informed loss only. This method merges the benefits of i) utilizing PINNs to eliminate the necessity for extensive training data [9,13] and ii) capitalizing on the enhanced inference speed and improved generalization abilities provided by CNN [14]. The main contributions of our work are as follows:

• Construct a set of new standard orthogonal compact supported LMWs.
• Prove the compact support property of LMWs integrated twice.
• Propose a novel network called PiLMWs-CNN to obtain LMWs coefficients (LMWCs) for approximating PDEs.

The rest of the paper is organized as follows. In Section 2, we conduct a survey of related works on solving PDEs using DNNs. We analyze and discuss these works from three different perspectives: physics-informed methods, neural operators (NO) learning, and advantages and applications of combining wavelets with CNNs. In Section 3, we provide the research methods of the paper. First, we conduct a feasibility analysis on combining the Legendre wavelet bases with CNNs to solve PDEs, including the role of compactness and orthogonality. Next, we provide relevant mathematical foundations, including the construction of the LMWs. Then, we give the network architectures and pipeline used in this study. In Section 4, we apply this method to approximate the wave equation, fluid equation, and heat conduction equation, with a particular focus on the loss of physical information and the stability of the solution. In Section 5, we give the conclusions.

2. Related Work

In this section, we will discuss different approaches to solve PDEs using NNs based on various mathematical foundations.

2.1. Physics-Informed Methods

This type of methods generally involves using PDEs as the constraints on training DNNs directly [9,15]. In numerous practical applications, constraints are often imposed by generalized prior knowledge, which take the form of differential equations (DEs) and are incorporated into the model training as the components of the loss function [16]. More recently, Raissi et al. [9] introduced the concept of PINNs, which have demonstrated their potential in various scenarios [17]. There are also specialized DL libraries, like DeepXDE [16], developed specifically for solving DEs.

Treating PDEs as training constraints was initially referred to as "soft constraints" [18], as the knowledge in the form of PDEs only enforces the correctness of the predicted results without guiding the optimization process. Additionally, the modeling process takes into account additional constraints in the form of ICs/BCs. Unlike conventional numerical approaches that rely on discretizing grids [19], utilizing PDE-based constraints to train DNNs ensures continuous solutions over the entire research domain, leading to highly precise results. PINNs are particularly suitable when the underlying physics of the problem is known beforehand, as they exploit the governing laws of the physical processes during training.

These methods leverage the idea that multilayer feedforward NNs (FFNNs) can serve as universal approximators for high-dimensional functions [20]. Both DNNs and solving DEs share a common objective: to achieve accurate approximations. Whether it is the FEM using basis functions to approximate solutions within a unit or DNNs with complex structures fitting high-dimensional functions, the computed or predicted values at each unit center or sampling point should be extremely close to the actual values. The goal for DNNs is to minimize the loss function value, or residual, which should ideally approach zero. In the numerical approach, this residual represents the error in the computed solution.

However, data-driven algorithms are often the only practical choice for natural processes because the underlying physics governing them is unknown. These algorithms do not generalize beyond the distribution of training data, even though they do not ensure the preservation of the problem’s governing physics [21].

2.2. Neural Operators Learning

In response to the challenges mentioned above, the concept of NO was introduced by Lu et al. [22]. NOs are designed to learn the relationship between two infinite-dimensional function spaces. After training, these NOs can accurately predict the solution for any input function presented to them. In comparison to conventional PDE solvers, NOs offer significant computational advantages. These NOs are based on the concept of universal operator approximation, as proposed by Chen et al. [23], which is similar to the theory of universal function approximation. Kovachki et al. [24] have conducted rigorous mathematical analysis and provided theoretical guarantees for the effectiveness of NOs.

By employing this methodology, DL can successfully obtain mesh-free, infinite-dimensional linear and nonlinear operators. Chen et al. [23] developed an initial prototype of NO methods based on operator theory, drawing inspiration from linear algebra and functional analysis. Currently, there are three popular instances in this area, which are Deep Operator Network (DeepONet) [25], Graph Neural Operator (GNO) [26], and Fourier Neural Operator (FNO) [27]. The DeepONet architecture comprises of two NNs: the branch NN and the trunk NN. The branch NN handles the input function, while the trunk NN is responsible for generating the output function at a specific sensor point. In DeepONet, the output is derived by taking the inner product of these NNs. The GNO focuses on learning mappings in infinite-dimensional spaces by combining nonlinear activation functions with a specific category of integral operators. Nevertheless, GNO may become unstable as the number of hidden layers increases. These advanced networks possess exceptional feature representation capabilities, and when combined with the idea of approximating functionals using the prototype network, the NOs become infused with boundless vitality. The FNO is an innovative approach inspired by the well-known Fourier transform (FT). Building upon the GNO, the FNO was introduced as an enhanced method to learn network parameters in Fourier space [27]. The FNO employs the fast FT (FFT) to perform spectral decomposition on the input signal, subsequent to calculating the convolution integral kernel in the Fourier domain. One notable limitation of FNO stems from the frequency localization of basis functions in FFTs, which hinders its ability to provide spatial resolution [28]. Consequently, FNO’s performance may be adversely affected when dealing with complex BCs.

An alternative method to tackle this issue is to employ wavelets, which exhibit both spatial and frequency localization properties [28]. By incorporating spatial information, wavelets demonstrate enhanced ability in handling signals with discontinuities and spikes, thus surpassing the pattern learning capabilities of FFTs when specifically dealing with image data. Wavelets have found applications in various domains, such as compression of fingerprints, iris recognition, denoising of signals, analysis of motion, and detection of faults, among other areas. The literature also highlights the utilization of wavelets in NNs [29]. In their work, Gupta et al. [30] proposed an operator that utilizes MW decomposition (MWD). This approach involves the use of four distinct NNs to calculate wavelet kernels. Specifically, the network architecture includes a fully connected NN (FNN), a CNN, and a Fourier integral layer similar to FNO. As the FNN represents the coarse scale, the combination of the latter three networks computes the detailed wavelet coefficients. As a result, MWD can be considered an improvement over FNO. Given that the solution of a PDE is obtained by calculating the inverse operator map between the input and the solution space, Gupta et al. aim to transform the problem of learning PDEs into a domain where a more compact representation of the operator can be obtained. In order to achieve this, they proposed an approach for learning the operator, based on MWs, which allows for compression of the associated kernel. By explicitly incorporating the inverse MW filters, they learn the projection of the kernel onto fixed MW polynomial bases. The use of MWs offers a natural basis for the MW subspace, using OPs with respect to a specific measure, and employing suitable scale/shifts to capture locality at different resolutions. Tripura et al. [31] introduced a NO, termed Wavelet Neural Operator (WNO), which utilizes spectral decomposition to learn mappings between infinite-dimensional function spaces. The proposed WNO demonstrates its effectiveness in dealing with highly nonlinear families of PDEs that involve discontinuities and abrupt changes in both the solution domain and boundary.

2.3. Wavelets and CNN

It is essential to highlight that the wavelet transform (WT) is a signal processing technique that exhibits the properties of localization in both the time and frequency domain. This attribute makes it a powerful tool for analyzing local features [32]. Consequently, Li et al. [33] developed a deep wavelet NN (DWNN) based on the PINNs approach. The employment of wavelets enables the extraction of multi-scale and detailed features, resulting in enhanced performance in solving PDEs.

However, it is well known that wavelet functions are not simultaneously symmetric, orthogonal, and compact supported. These challenges are overcome by LMWs, also referred to as more general, vector-valued polynomial types of wavelets [34], which simultaneously have all three of these properties. This leads to high smoothness and high approximation order and suggests that LMWs can outperform scalar wavelets in a variety of applications [35].

Therefore, some researchers have focused on investigating how to solve PDEs using the Legendre wavelets method. It is worth mentioning that Legendre wavelets have both spectral accuracy and orthogonality, as well as other properties mentioned about wavelets [36]. For example, the wavelet technique can transform complex problems into a system of algebraic equations. Heydari et al. [37] presented a numerical technique utilizing two-dimensional Legendre wavelets to solve fractional PDEs with Dirichlet BCs. The method involved using operational matrices for fractional integration and derivatives to obtain numerical solutions to the problems under consideration. The approach presented by Abbas et al. [38] primarily focuses on converting the fundamental DEs into integral equations through integration. This is accomplished by approximate representation of the diverse signals embedded in the equations using truncated orthogonal series $\phi \left(x\right)$. Subsequently, the operation matrix of integration, denoted as P, is employed to remove the integral operations.

Meanwhile, we found in the research that there are many applications for the combination of wavelets and CNN. It was discovered that a CNN can be considered a simplified version of multiresolution analysis (MRA). This revelation highlights the fact that conventional CNNs fail to capture a substantial portion of the spectral information accessible through MRA. Therefore, Fujieda et al. [39] introduced the idea of wavelet CNNs, which combine MRA and CNNs within a single model. This strategy enables the exploitation of spectral information that is frequently neglected in conventional CNNs but is essential in various image processing tasks. Zhao et al. [40] suggested a Wavelet-Attention CNN (WA-CNN) for image classification purposes. They explored Discrete WT (DWT) in the frequency domain and devised a novel WA block, which solely implements attention in the high-frequency domain. DWT exhibits superior down-sampling information quality in the image processing domain, thereby significantly minimizing the loss of down-sampling information in CNN. Guo et al. [41] developed a method for intelligent fault diagnosis of rolling bearings, utilizing WT and deformable CNN. Their research highlights the significance of selecting appropriate wavelet bases in the WT. Different wavelet bases have a substantial influence on the resulting time-frequency map. In the context of CNNs, WT is developed and employed for breast cancer detection [42]. The adaptive wavelet pooling layers, as proposed by the Wolter et al. [43], make use of fast WT (FWT) to lower the feature resolution. Through FWT, the input features are decomposed into various scales, resulting in a reduction of feature dimensions by eliminating the fine-scale subbands. This approach introduces additional flexibility by optimizing the wavelet basis functions and weighting coefficients at different scales, thus avoiding high repetition.

The findings above indicate that CNNs exhibit strong performance when handling multidimensional data, such as images and speech signals. In fact, these networks have been employed in a PINN, primarily within the domain of image processing. Inspired by the above researches, we consider combining LMWs with CNN to approach PDEs, and conduct experiments on three different types of equations to validate the effectiveness and accuracy of this approach.

3. Method

3.1. Feasibility Analysis

By integrating LMWs with CNN for solving PDEs, this approach capitalizes on four key advantages, thereby enhancing the modeling and analysis capabilities in PDE solving and delivering more precise predictions and solutions.

• Multiscale Representation: PDEs often involve spatial features at different scales. By using the LMWs, it provides the capability of multiscale representation, which can better capture the features of the input data at different scales. This is crucial for analyzing and solving PDEs involving multiscale spatial features.
• Local Correlation Modeling: The mathematical models in PDEs often assume that the system’s behavior is influenced by its local neighborhood. The LMWs have good properties in capturing the local correlation of input data effectively. By using convolution operations in CNN, it can leverage this local correlation and perform local feature extraction on the input data, which helps in describing the local behavior in PDEs more accurately. Bedford et al. [44] prove rigorously that good approximation "locally" guarantees good approximation globally.
• Translation Invariance: Analysis of spatial translation invariance is often required for input data in PDEs. The translation invariance property of LMWs ensures that the translation of input signals does not affect the representation of their MW coefficients (MWCs). In CNN, convolution operations can preserve the translation invariance of input data, better meeting the requirement of translation invariance in PDEs.
• Explanation and Interpretability: LMWs have good explanatory and interpretable properties, providing insights into the features and structure of PDEs. By combining LMWs with CNN, we can leverage the explanatory power of MWs and the learning capability of CNN to better understand and interpret the feature representations of PDEs. In fact, Restricted Boltzmann machines (RBMs) [45] and other variants of DL models, such as deep belief networks and autoencoders, have a clear advantage in providing explanations and interpretability.

Furthermore, the compact support property of LMWs offers several advantages [46,47] when applied to solving PDEs with NNs.

• Reduced computational cost: The compact support of LMWs enables a smaller set of wavelet coefficients to be considered, thereby reducing the computational burden associated with representing the solution as a wavelet series.
• Localized feature extraction: LMWs are well-suited for extracting localized features within the solution, facilitating the NN’s ability to comprehend the underlying structure of the PDE.
• Regularization effect: LMWs’ compact support property inherently creates a regularization effect in the NN, which reduces overfitting and improves generalization performance.
• Adaptability to different scales: LMWs readily adapt to various scales by adjusting the wavelet coefficients, permitting the NN to capture features at different resolutions.
• Robustness to noise: LMWs’ compact support enhances their resilience to noise and other imperfections in the input data, enabling more accurate solutions to be obtained when employed in conjunction with NNs.

Based on the aforementioned analysis, a strong correlation between LMWs, CNNs, and PDEs becomes apparent. Notably, the compact support property of LMWs offers numerous advantages when resolving PDEs via NNs. These distinctive features render the integration of wavelets and CNNs particularly appealing for PDE solutions, providing the authors with a promising avenue to approximate PDEs. Consequently, it piques the authors’ curiosity and motivates further exploration and in-depth experimentation with this approach.

3.2. Foundations of Multiwavelet Bases Construction in Mathematics

Families of functions ${\phi }_{a,b}\left(x\right)$,

$${\phi }_{a,b}\left(x\right)={\left|a\right|}^{-1/2}\phi \left(\frac{x-b}{a}\right),a,b\in R,a\ne 0,$$

obtained from a single function $\phi$ by dilation and translation and serve as the basis for $L^2\left(R\right)$, are known as wavelets [48].

In the following, we will provide an overview of the properties of the MW bases developed in [49] and introduce the necessary notations.

3.2.1. MW bases

In this subsection, we introduce a new class of wavelet-like bases referred to as MW bases, which enable sparse representations of smooth integral operators over finite intervals. Moreover, MW bases have orthogonality, compact support, vanishing moments, and other wavelet base properties. The basis functions do not overlap on a specific scale and are arranged into small groups of multiple functions (thus, MWs) sharing the same support. One notable benefit of this construction is its simplicity.

One-dimensional construction

Our focus will initially be limited to $[0,1]\subset R$, and we now proceed to establish a basis for $L^2[0,1]$. Each basis consists of dilates and translates of a finite set of functions ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$. Specifically, these bases consist of orthonormal systems

$${\phi }_{l,i}^k=2^{i/2}{\phi }_l(2^{i}x-k),l=1,2,\cdots ,n;i,k\in Z,$$

(3.1)

where the functions ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$ are piecewise polynomials, they become zero outside the range of $[0,1]$, and possess vanishing moments, making them orthogonal to low-order polynomials,

$${\int }_0^1{\phi }_l\left(x\right)x^{i}dx=0,i=0,1,\cdots ,n-1.$$

(3.2)

Besides, we implement the MRA approach, and suppose $n\in N^{+}$, for $i=0,1,2,\cdots$, we define a space $V_i^n$ of piecewise polynomial functions,

$$\begin{array}{cc}\hfill & V_i^n=\{f:\mathrm{the}\mathrm{restriction}\mathrm{of}f\mathrm{to}\mathrm{the}\mathrm{interval}(2^{-i}k,2^{-i}(k+1))\mathrm{is}\hfill \\ \hfill & \mathrm{a}\mathrm{polynomial}\mathrm{of}\mathrm{degree}\mathrm{less}\mathrm{than}n,\mathrm{for}k=0,1,\cdots ,2^i-1,\hfill \\ \hfill & \mathrm{and}f\mathrm{vanishes}\mathrm{elsewhere}\}.\hfill \end{array}$$

(3.3)

It is clear that the space $V_i^n$ has a dimension of $2^{i}n$ and

$$V_0^n\subset V_1^n\subset \cdots \subset V_i^n\subset \cdots .$$

(3.4)

For $i=0,1,2\cdots$, we define the $2^{i}n$-dimensional space (the MW subspace) $W_i^n$ to be the orthogonal complement of $V_i^n$ in $V_{i+1}^n$,

$$V_i^n\oplus W_i^n=V_{i+1}^n,V_i^n\perp W_i^n,$$

(3.5)

therefore, we can inductively derive the decomposition

$$V_i^n=V_0^n\oplus W_0^n\oplus W_1^n\oplus \cdots \oplus W_{i-1}^n.$$

(3.6)

Suppose that functions ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$ form an orthogonal basis for $W_0^n$, that is $W_0^n=linearspan\{{\phi }_1,{\phi }_2,\cdots ,{\phi }_n\}$. More generally,

$$\begin{array}{cc}\hfill W_i^n& =\mathrm{linear}\mathrm{span}\{{\phi }_{l,i}^k:{\phi }_{l,i}^k\left(x\right)\hfill \\ \hfill & =2^{i/2}{\phi }_l(2^{i}x-k),l=1,\cdots ,n,i=0,1,2,\cdots ,k=0,1,\cdots ,2^i-1\}.\hfill \end{array}$$

(3.7)

Since $W_0^n$ is orthogonal to $V_0^n$, the first n moments of ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$ vanish,

$${\int }_0^1{\phi }_l\left(x\right)x^{i}dx=0,\hspace{1em}\hspace{1em}i=0,1,\cdots n-1.$$

(3.8)

The $2n$-dimensional space $W_1^n$ is spanned by $2n$ orthogonal functions ${\phi }_1\left(2x\right),\cdots ,{\phi }_n\left(2x\right),{\phi }_1(2x-1),\cdots ,{\phi }_n(2x-1)$, of which n are supported on the interval $[0,1/2]$ and n on $[1/2,1]$. Generally, the space $W_i^n$ is spanned by $2^{i}n$ functions obtained from ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$ by translation and dilation. There is some freedom in choosing the functions ${\phi }_1,{\phi }_2,\cdots ,{\phi }_n$, subject to the restriction of orthogonality, we can uniquely determine them by imposing normality and further vanishing moments, except for sign.

Completeness of one-dimensional construction

We define the space

$$V^n=\bigcup _{i=0}^{\infty }{V}_i^n,$$

(3.9)

and observe that $\overline{V^n}=L^2[0,1]$. Here the closure $\overline{V^n}$ is defined with respect to the $L^2$-norm, $\left|\right|f\left|\right|={〈f,f〉}^{1/2}$, where the inner product $〈f,g〉={\int }_0^{1}f\left(x\right)g\left(x\right)dx$.

Let $\{{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_n\}$ be an orthonormal basis for $V_0^n$ , in accordance with $\left(\right)$-$\left(\right)$, the set $B_n=\{{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_n\}\cup \left\{{\phi }_{l,i}^k\right\}$ forms a complete orthonormal basis for $L^2[0,1]$. We refer to $B_n$ as the MW bases of order n for $L^2[0,1]$.

3.2.2. Multiple-dimensional construction

For any positive integer d, the bases construction for $L^2[0,1]$ can be extended to some other function spaces as well, such as $L^2{[0,1]}^d$. The basis for $L^2{[0,1]}^2$, which serves as an example of the construction for any finite-dimensional space, is now given in order to outline this extension. We establish the space

$$V_i^{n,2}=V_i^n\times V_i^n,i=0,1,2,\cdots ,$$

(3.10)

where $V_i^n$ is defined by $\left(\right)$. Furthermore, we define $W_i^{n,2}=W_i^n\times W_i^n$ to be the orthogonal complement of $V_i^{n,2}$ in $V_{i+1}^{n,2}$,

$$V_i^{n,2}\oplus W_i^{n,2}=V_{i+1}^{n,2},V_i^{n,2}\perp W_i^{n,2}.$$

(3.11)

$W_0^{n,2}$ is then the space generated by the orthonormal basis

$$\{{\varphi }_i\left(x\right){\phi }_j\left(y\right),{\phi }_i\left(x\right){\varphi }_j\left(y\right),{\phi }_i\left(x\right){\phi }_j\left(y\right)\text{:}i,j=1,2,\cdots ,n\}.$$

(3.12)

For each basis element $v(x,y)$ of these $3n^2$ elements, it has no projection on any low-order polynomial,

$${\int }_0^1{\int }_0^{1}v(x,y)x^i{y}^{j}dxdy=0,\hspace{1em}\hspace{1em}i,j=0,1,\cdots ,n-1.$$

(3.13)

Dilations and translations of the $v(x,y)$ span the space $W_i^{n,2}$; low-order polynomials $\{{\varphi }_i\left(x\right){\varphi }_j\left(y\right):$

$i,j=1,2,\cdots ,n\}$ are also included in the basis of $L^2{[0,1]}^2$, which is made up of these functions.

3.3. Construct LMWs

In this part, our primary focus is the development of MW bases based on Legendre polynomials, which we refer to as LMWs. We present detailed construction procedures and provide the necessary proofs to support our approach.

3.3.1. Legendre polynomials

Legendre polynomials are polynomials defined on the $[-1,1]$ and have a recursive formula [50]:

$$\begin{array}{cc}\hfill & p_0\left(x\right)=1,p_1\left(x\right)=x,\hfill \\ \hfill & p_{n+1}\left(x\right)=\frac{(2n+1)x{p}_n\left(x\right)-n{p}_{n-1}\left(x\right)}{n+1},n=1,2,\cdots .\hfill \end{array}$$

(3.15)

In order to be consistent with the previous knowledge, we consider a transformation to fix the discussion interval at $[0,1]$ and construct the standard orthogonal basis of $L^2[0,1]$. Finally, a similar transformation can be performed to return to the desired interval. The implementation details can be found in literature [51]. Therefore, according to $\left(\right)$, Legendre polynomials on $[0,1]$ can be expressed as follows through translation and expansion transformations:

$$\begin{array}{cc}\hfill & p_0\left(x\right)=1,p_1\left(x\right)=2x-1,\hfill \\ \hfill & p_{n+1}\left(x\right)=\frac{(2n+1)(2x-1)p_n\left(x\right)-n{p}_{n-1}\left(x\right)}{n+1},n=1,2,\cdots .\hfill \end{array}$$

(3.16)

Further, $q_l\left(x\right)={\left\{\sqrt{2l+1}{p}_l\left(x\right)\right\}}_{l=0}^{\infty }$ is an orthonormal basis in $L^2[0,1]$.

3.3.2. LMWs

Now, we use Legendre polynomials $q_l\left(x\right)(l=0,1,2,3)$ as LMW multiscaling function [52], to construct the Legendre mother multiwavelet [53] as

$${\phi }_l\left(x\right)=\sum _{j=0}^3\left(h_{lj}{q}_j\left(2x\right)+g_{lj}{q}_j(2x-1)\right),l=0,1,2,3.$$

(3.17)

By the properties of vanishing moments and standard orthonormality, ${\phi }_l\left(x\right)$ can be uniquely determined. Figure 1 shows the plots of these four Legendre mother multiwavelets.

Then, we construct the LMWs by translating and dilating the mother multiwavelet ${\phi }_l\left(x\right)$ and they can be expressed as

$${\phi }_{l,i}^k\left(x\right)=2^{i/2}{\phi }_l(2^{i}x-k),i\in N^{+},k=0,1,\cdots ,2^i-1,l=0,1,2,3.$$

(3.18)

In addition, ${\phi }_{l,i}^k\left(x\right)$ is compact supported on the interval $\left[k/2^i,(k+1)/2^i\right]$, and $\left\{q_l\right\}\bigcup \left\{{\phi }_{l,i}^k\right\}$ forms a complete orthonormal basis for $L^2[0,1]$.

3.3.3. Compact support

Given the common use of second-order PDEs in physics and engineering [54], we aim to solve them in the $W_2[0,1]$ [55] by integrating ${\phi }_{l,i}^k\left(x\right)$ twice. Then, we denote

$$J_0^2{\phi }_{l,i}^k\left(x\right)\triangleq {\int }_0^x{\int }_0^x{\phi }_{l,i}^k\left(t\right)d{t}^2.$$

(3.19)

Lemma 1. 

$J_0^2{\phi }_{l,i}^k\left(x\right)={\int }_0^x(x-t){\phi }_{l,i}^k\left(t\right)dt$

Theorem 1. 

$J_0^2{\phi }_{l,i}^k\left(x\right)$is compact supported on $\left[\frac{k}{2^i},\frac{k+1}{2^i}\right]$.

Proof. For $J_0^2{\phi }_{l,i}^k\left(x\right)={\int }_0^x(x-t){\phi }_{l,i}^k\left(t\right)dt$, if $x<\frac{k}{2^i}$, $J_0^2{\phi }_{l,i}^k\left(x\right)=0;$ if $x>\frac{k+1}{2^i}$,

$$J_0^2{\phi }_{l,i}^k\left(x\right)={\int }_{\frac{k}{2^i}}^{\frac{k+1}{2^i}}(x-t){\phi }_{l,i}^k\left(t\right)dt=2^{\frac{i}{2}}{\int }_{\frac{k}{2^i}}^{\frac{k+1}{2^i}}(x-t){\phi }_l(2^{i}t-k)dt$$

(3.20)

Let $2^{i}t-k=s,$ then $t=\frac{s+k}{2^i},dt=\frac{1}{2^i}ds$. So, the transformation of (3.20) is as follows,

$$\begin{array}{cc}\hfill J_0^2{\phi }_{l,i}^k\left(x\right)& =2^{-\frac{i}{2}}{\int }_0^1\left(x-\frac{s+k}{2^i}\right){\phi }_{l,0}^0\left(s\right)ds\hfill \\ \hfill & =2^{-\frac{i}{2}}{\int }_0^1\left(x-\frac{s+k}{2^i}\right){\phi }_l\left(s\right)ds\hfill \\ \hfill & =0\hfill \end{array}$$

(3.21)

The theorem holds true. $\square$

3.3.4. Function approximation

We employ the tensor product of multiple one-dimensional ($1D$) basis functions to derive basis functions in multiple dimensions [56], building upon equation $\left(\right)$ discussed earlier. To simplify the expression, we still use ${\phi }_{l,i}^p$ as a basis function in $1D$, so that we have the following equation:

$${\phi }_{lmn,ijk}^{p,q,r}(x,y,t)={\phi }_{l,i}^p\left(x\right){\phi }_{m,j}^q\left(y\right){\phi }_{n,k}^r\left(t\right).$$

(3.22)

With discrete MWCs $c_{\widehat{x},\widehat{y},\widehat{t}}^{p,q,r}$ on a grid $\widehat{x},\widehat{y},\widehat{t}\in \widehat{X}\times \widehat{Y}\times \widehat{T}$, we obtain a continuous LMW $\tilde{u}(x,y,t)$ as follows:

$$\tilde{u}(x,y,t)=\sum _{\begin{array}{c}p,q,r\in [0,2^i-1]\times [0,2^j-1]\times [0,2^k-1]\\ \widehat{x},\widehat{y},\widehat{t}\in \widehat{X}\times \widehat{Y}\times \widehat{T}\end{array}}{c}_{\widehat{x},\widehat{y},\widehat{t}}^{p,q,r}{\phi }_{lmn,ijk}^{p,q,r}(x-\widehat{x},y-\widehat{y},t-\widehat{t}).$$

(3.23)

Our goal is to find MWCs $c_{\widehat{x},\widehat{y},\widehat{t}}^{p,q,r}$ such that $\tilde{u}(x,y,t)$ as nearly as possible resembles the PDE solution. By taking the corresponding derivatives of the LMWs, one can directly obtain the partial derivatives of $\tilde{u}$ with respect to $x,y,t$.

3.4. Neural Network Architecture

Figure 2 showcases the network architectures employed for the incompressible N-S equation and the damped wave equation. The domain’s occupancy grid at a timepoint t is represented by ${\Omega }^t$ and $u_d^t$ contains Dirichlet BCs. We utilize a U-Net architecture [11] internally for the N-S equation, the damped wave equation, however, only required a basic 3-layer CNN. These networks were employed to calculate residuals $\left(\Delta c_{ij}\right)$,which are then added to $c_{ij}^t$ to derive the MWCs of the next time-step $c_{ij}^{t+dt}$. Thus, we mean-normalize the coefficients to zero, that is ${\sum }_{\widehat{x},\widehat{y}}{c}_{\widehat{x},\widehat{y},\widehat{t}}^{0,0,0}=0$.

Pipeline

As depicted in Figure 3, PiLMWs-CNN employs a CNN (PDE Model) to map discrete LMWCs and BCs from a timepoint $\widehat{t}$ to LMWCs at a subsequent timepoint $\widehat{t}+dt$. By iteratively applying the PDE Model on the LMWCs, the simulation can be propagated in time. Efficient computation of continuous LMWs is achievable using transposed convolutions (refer to (3.23)) with LMWs, as illustrated in Figure 4.

Training procedure

As demonstrated in Figure 5, our method adopts a training process similar to that in reference [11]. Initially, we create a randomized training pool consisting of domains and LMWCs. Since all LMWCs can be set to zero at first, ground truth data is not required. The PDE model then uses the random minibatch (batch size = 50) that we extracted from the training pool to predict the LMWCs for the next time step. Next, with the goal of optimizing the PDE model weights via gradient descent, we compute a physics-informed loss inside the volume spanned by the LMWCs of the minibatch and the predicted LMWCs. We employ the Adam optimizer (learning rate = ${10}^{-4}$) for this purpose. We update the training dataset with the newly predicted LMWCs during the training process to gradually replenish the pool with more realistic training data. From time to time, we reset all of the LMWCs in the training pool to zero, which enables us to learn the warm-up phase from the ground up.

Physics-informed loss

Drawing parallels with reference [11], the PiLMWs-CNN approach integrates the benefits of both physics-constrained methodologies [48] and physics-informed strategies [9]. This method shows the ability to generalize to new domain geometries and yields continuous solutions that avoid the drawbacks of a discretized loss function based on finite differences by using a CNN to manage LMWCs, which can be thought of as a discrete hidden latent representation for a continuous implicit field description based on LMWs. We aim to minimize the integrals of the squared residuals of the PDEs over the domain / domain boundaries and time steps by optimizing the LMWCs. We uniformly randomly choose points inside the designated integration domains in order to compute these integrals. Within the domain, we require two loss terms for the damped wave equation 26:

$$L_u={\int }_{\Omega }{\int }_{\widehat{t}}^{\widehat{t}+dt}{∥{\partial }_{t}u-v_u∥}^2$$

(3.24)

$$L_v={\int }_{\Omega }{\int }_{\widehat{t}}^{\widehat{t}+dt}{∥{\partial }_t{v}_u-k\Delta u+\gamma v_u∥}^2$$

(3.25)

The boundary loss term is:

$$L_b={\int }_{\partial \Omega }{\int }_{\widehat{t}}^{\widehat{t}+dt}{∥u-u_d∥}^2$$

(3.26)

For the integral to be computed, it is imperative that the residuals have bounded variation.

4. Results

In this section, we applied PiLMWs-CNN to approximate the damped wave equation, the incompressible N-S equation, and the two-dimensional heat conduction equation. Through quantitative evaluation and stability analysis, it is evident that our method achieves higher accuracy and provides better approximation for both the governing equations and BCs. The figures below are displayed in TensorFlow after training.

Example 1. The Damped Wave Equation appears as a mathematical model in biology and physics [57]. It is important in geophysics, general relativity, quantum mechanics, and plasma physics. We consider the following equation with constant positive damping:

$$\begin{array}{cccc}\hfill {\partial }_{t}u& =v_u,\hfill & \hfill & in\Omega ,\hfill \\ \hfill {\partial }_t{v}_u& =k\Delta u-\gamma v_u,\hfill & \hfill & in\Omega ,\hfill \\ \hfill u& =u_d,\hfill & \hfill & on\partial \Omega ,\hfill \\ \hfill u& =u_0,v_u=v_{u_0},\hfill & \hfill & in\Omega ,t=0.\hfill \end{array}$$

(4.1)

Here, k is the stiffness constant, $\gamma$ is a damping constant. Wave equation solutions are essential for understanding the concepts in fluid dynamics, optics, gravitational physics and electromagnetism. Moreover, for large damping constants $\gamma$, converges towards the Laplace equation solution [11].

Table 1 compares the losses of our method with that of reference [11]. We trained with LMWs $J_0^2{\phi }_{l,i}^k\left(x\right)$ (Figure 4) in the spatial dimensions, and observed a much better performance. From the result of $L_v$, it can be seen that our method follows the physical laws well.

As shown in Figure 6, our method produces stable results for the damped wave equation.

Example 2. The Incompressible N-S equations are a set of PDEs that describe the movement of viscous fluids. These equations are derived by applying Newton’s second law of motion to an incompressible Newtonian fluid. The resulting equation is known as the N-S equation written in vector notation:

$$\begin{array}{cccc}\hfill \nabla ·\overrightarrow{u}& =0,\hfill & \hfill & in\Omega ,\hfill \\ \hfill \rho ({\partial }_t\overrightarrow{u}+\overrightarrow{u}·\nabla \overrightarrow{u})& =\mu \Delta \overrightarrow{u}-\nabla p+\overrightarrow{f},\hfill & \hfill & in\Omega .\hfill \end{array}$$

(4.2)

The equations take into account factors such as the dynamic viscosity $\mu$, pressure p , fluid velocity vector $\overrightarrow{u}$, and the fluid’s density $\rho$. They can be used to model a wide range of phenomena, like the flow of water in a pipe, blood in an artery, and air over an airplane wing, weather patterns in the atmosphere, and even the flow of stars in a galaxy. In our tests, we don’t take into account the external factors $\overrightarrow{f}$. For given initial conditions ${\overrightarrow{u}}_0$, $p_0$ and BCs, these pair of equations need to be resolved. We are taking into account the Dirichlet BC, where the velocity field $\overrightarrow{u}$ is set to ${\overrightarrow{u}}_d$ at the boundaries of the domain $\partial \Omega .$

We conducted an extensive study on the stability of PiLMWs-CNN when applied to the incompressible N-S equations. The investigation involved hundreds of iterations on the DFG benchmark [58] problem at Re = 100 (Re represents the Reynolds number). Figure 7 shows the curves of $E\left[\right|\nabla ·\overrightarrow{u}\left|\right]$ and $L_p$, where $L_p$ is a momentum loss term defined as

$$L_p={\int }_{\Omega }{\int }_{\widehat{t}}^{\widehat{t}+dt}\left|\right|\rho ({\partial }_t\overrightarrow{u}+\overrightarrow{u}·\nabla \overrightarrow{u})-\mu \Delta \overrightarrow{u}+\nabla p|{|}^2$$

(4.3)

The results clearly demonstrate the superiority of our approach in terms of both accuracy and stability. This significant improvement highlights the effectiveness of our method over the existing reference [11]. Furthermore, following the same approach as implicit PINNs and Spline nets, which involves computing losses based on physics-informed losses, PiLMWs-CNN calculated $E\left[\right|\nabla ·\overrightarrow{u}\left|\right]$ and $L_p$ and achieved lower loss values, as shown in Table 2. This demonstrates the enhanced effectiveness of our method in approximating the solutions of PDEs.

Since the proposed method in this paper utilizes physics-informed loss computation during model training, which is primarily based on the residual concept, it aligns with the idea of $\epsilon$-best approximate solution [51] or $\epsilon$-approximate approach [59] that have been proposed in recent years for numerical solutions of PDEs and fractional DEs (FDEs). Therefore, as a final example, we intend to apply the PiLMWs-CNN method to the two-dimensional heat conduction equation to observe its approximation performance.

Example 3. The Two-dimensional heat conduction equation is a PDE that governs the thermal conduction and heat transfer within a medium. Mathematically, it is defined as:

$$\begin{array}{cc}\hfill {\partial }_{t}u& ={\partial }_{x^2}u+{\partial }_{y^2}u,(t,x,y)\in [0,1]\times [0,1]\times [0,1],\hfill \\ \hfill u(0,x,y)& =sin\left(\pi x\right)sin\left(\pi y\right),\hfill \\ \hfill u(t,0,y)& =u(t,1,y)=u(t,x,0)=u(t,x,1)=0.\hfill \end{array}$$

(4.4)

The exact solution of is $u(t,x,y)=e^{-2{\pi }^{2}t}sin\left(\pi x\right)sin\left(\pi y\right).$ We conducted the experiment for Ex.3 using the same network structure as the damped wave equation. $L_{b}oundary$ represents the boundary loss, calculated in a manner similar to equation $\left(\right)$, while $L_{h}eat$ signifies the equation loss, calculated similarly to equation $\left(\right)$. From Table 3, it can be observed that the absolute errors computed by our method(right) are smaller. Figure 8 demonstrates a progressive reduction in the boundary loss of the heat conduction equation (reaching ${10}^{-9}$) as the iterations progress, while the equation’s loss stabilizes (reaching ${10}^{-6}$). This indicates that our method is capable of more accurately capturing the physical laws of the equation.

Table 3. The absolute errors $|u-u_n|$ for Ex.3($t=1$). The left subtable is computed by Shifted-Legendre orthonormal method [51], and the right is computed by PiLMWs-CNN(ours).

Table 3. The absolute errors $|u-u_n|$ for Ex.3($t=1$). The left subtable is computed by Shifted-Legendre orthonormal method [51], and the right is computed by PiLMWs-CNN(ours).

5. Conclusions

This paper presents an approach to approximate PDEs by training a continuous PiLMWs-CNN utilizing solely physics-informed loss functions. A feasibility analysis of this method was conducted, properties of LMWs were analyzed, and their advantages in solving PDEs were discussed. Furthermore, the construction theory and methods of LMWs were also provided. The PiLMWs-CNN utilizes the ability of CNN in fast modeling complex systems in solving PDEs to process LMWCs through CNN, trains the network using physics-informed loss function. The experimental results show that this method is more effective, more accurate, and yields more stable solutions.

In the future, we will conduct further research on the application of the algorithm proposed in this paper for solving FDEs. This research will encompass the exploration of solution methods and network structures for constant-order, variable-order, and fractional PDEs. In addition to accuracy and stability measures, we will also investigate and experiment with other evaluation metrics. For instance, in the context of the previously mentioned damping wave equation, we will examine the doppler effect, wave reflections, as well as the drag and lift coefficients at different Reynolds numbers in the N-S equation. This assessment from a physical application perspective will provide a deeper evaluation of the algorithm’s performance and facilitate a gradual optimization of both the algorithm and the network.

Furthermore, a multigrid LMWs CNN might be contemplated for the further refinement of solutions at the boundary layers. Our method’s fully differentiable feature might also be useful in situations involving gradient-oriented shape optimization, optimal control, reinforcement learning, and sensitivity analysis. We firmly believe that the performance of upcoming ML-based PDE solvers with generalization capabilities will be positively impacted by the shift from explicit physics-constrained losses to implicit physics-informed losses on continuous fields, supported by discrete latent descriptions like LMWCs.