Partial di!erential equations (PDEs) of elliptic and parabolic type are ubiquitous in the mathematical modeling of many physical phenomena, and thus see wide application to many real-world problems. Classical numerical methods (finite di!erence methods, finite element methods, etc.) introduce a computational mesh over which matrix representations of di!erential operators are defined. This results in large linear systems whose e"cient solution presents a number of computational challenges. Recently, neural networks have gained popularity in solving PDEs, partially thanks to the universal approximation property of neural networks [1,2]. Popular approaches include Physics Informed Neural Networks (PINNs) [3], the Deep Galerkin method [4], and the Deep Ritz method [5,6]. These methods incorporate the underlying physics of the problem by introducing the di!erential operators of the PDE or an equivalent variational principle into the loss function. The common feature of these methods is the requirement of training data as these methods belong to the class of supervised learning methods. This requires su"cient training data, which is often generated using more classical numerical algorithms. An important drawback is the poor generalization in neural networks for supervised learning.
Our proposed method di!ers from the methods described above because it does not require training data. This helps mitigate the generalization problem of existing approaches. Our unsupervised learning method is an iterative method that computes successive approximations of the solution of the PDE given a source function and boundary data. The iterations are constructed to minimize a finite di!erence based loss function. To reduce the computational cost, the PocketNet paradigm is applied to the neural network architecture [7]. The method is presented for both steady-state and time-dependent problems. Finite di!erence methods have inspired the construction of various neural networks [8,9], but again these methods need training data. The recent work [10] solves elliptic PDEs with constant diffusion and it shares the following similarities with our proposed numerical method: no training data is needed and the loss function is derived from the five-point stencil of the finite di!erence approximation. However, our loss function is di!erent as it incorporates the Dirichlet boundary condition in a weighted fashion. We also define the method for elliptic PDEs with nonconstant di!usion coe"cient and we extend it to time-dependent problems.
Another key feature of our method is the use of small neural networks, which makes it computationally e"cient.
An outline of the paper is as follows. Section 2 introduces the model problem, the network architecture and the loss function. Section 3 describes the algorithms for both steady-state and time-dependent problems. Results and discussion are presented in Section 4. Conclusions follow.
We propose a fully unsupervised method for estimating finite di!erences solutions to partial di!erential equations via convolutional neural networks. In contrast to existing deep learning-based methods, our approach is fully unsupervised. In other words, our approach does not require training data and estimates the solution to a given PDE directly via the optimization process, also called the training process.
Let # → R 2 be an open set. Given a function f : # ↑ R and di!usion coe"cients ω, we consider the 2D Poisson problem with solution u :
We begin with the case when ω is constant and assume that ω = 1. For readibility, we assume that # is the square domain (0, L) 2 , that is partitioned into a uniform grid of N ↗ N squares with size h = L/N . The methodology presented in the paper can be easily extended to rectangular domains. We first recall the standard finite di!erence method based on the five-point stencil, applied to (1). Let V 0 (resp. V ω ) be the set of indices (i, j) such that the point (x i , y j ) = (ih, jh) belongs to the interior (resp. boundary) of #.
The finite di!erence solution u h is a vector with entries u i,j satisfying (2)-( 3) and it is known that for smooth enough solutions u and small enough h, the value u i,j is a good approximation of u(x i , y j ). One can rewrite (2) in terms of a discrete convolutional operator ϑ and convolutional kernel K ! defined as:
We use the shorthand notation f i,j = f (x i , y j ). The equivalent form for ( 2) is
For a given kernel K ≃ R 3↑3 , the convolution operator ϑ is defined by
With this convolution-based discretization, we can now reformulate the finite di!erence approximation of (1) as a convex optimization problem given by ( 7) arg min
.
We obtain an approximate solution û to (7) by training a neural network N ε : R N ↑N ↑ R N ↑N with trainable parameters ω ≃ R M . We propose to train N ε using the unsupervised loss function L ϑ : R N ↑N ↑ R given by
where ϖ = h 2 /4 is a weighting term. The idea behind this loss function is to directly estimate the finite di!erence approximation to (1) without using any training data. Specifically, we use f as the input to the CNN and we perform M updates of the parameters in the neural network N ω , using the loss function defined above. Algorithm 1 in Section 3 outlines this training procedure in more detail.
In the case where ω is piecewise constant, ω = (ω i,j ), we use the following finite di!erence discretization
where the interface values ω i+ 1 2 ,j , ω i→ 1 2 ,j , ω i,j+ 1 2 , and ω i,j→ 1 2 are the harmonic averages of ω i,j given by
.
The interface values necessitate the imposition of a half-grid (also called dual grid) onto our original grid. As a result, we modify our neural network prediction via a dilation operator D : R N ↑N ↑ R (2N →1)↑(2N →1) such that
Additionally, we map ω from the original grid to the dual grid by applying the following steps:
Step 1: apply the dilation operator to obtain ω ↔ ≃ R (2N →1)↑(2N →1) :
Step 2: update the values of ω ↔ on the half-grid:
, if i odd and j even 2
, if i even and j odd (13) Step 3: set to zero the values of ω ↔ that would correspond to the values on the original grid:
In the equation above, 1 N ↑N is an N ↗ N matrix with all entries equal to one, and ⇐ is the Hadamard or pointwise product.
We now introduce the following convolution kernels
0 0 1 0 0 0 0 0 0 0 0 0 ↓1 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
which will allow us to rewrite (9) as follows:
which is then used in the first term of the loss function (8) to train our neural network. The same training process outlined in Algorithm 1 applies here. Note that in (16), the convolution operator ϑ uses a stride of two because of the dual grid and the use of larger convolution kernels.
We generalize the method described in the previous section to timedependent problems. We propose a modified loss function that incorporates the backward Euler method to address the time derivative. Consider the following parabolic problem on a square domain # and for the time interval [t 0 , T ].
Let ϱ > 0 be the time step size. At each time step, t n = nϱ , the finite di!erence solution denoted by u n h solves the following optimization problem arg min
Above we use the short-hand notation f n i,j = f (x i , y j , t n ) and g n i,j = g(x i , y j , t n ). Here, the unknowns u n i,j are the approximations of the solution u evaluated at the point (x i , y j ) and time t n . We then use the following loss function to train our CNN at each time step
Again, we will choose the weight ϖ = h 2 /4. The key di!erence between ( 8) and ( 19) is that we minimize at each time step. The input to the CNN at each time step is the previous solution ûn→1 and the output is the solution at the time step t n . Algorithm 2 in Section 3 outlines this training procedure in more detail.
We select as our architecture the widely used U-Net [11]. The U-Net was initially designed for image segmentation tasks but is successful in several scientific machine learning tasks [12,13,14]. The architecture consists of convolution block operators (blocks), a downsampling (or encoding) path, and an upsampling (or decoding) path. A channel-wise concatenation operation links each layer in the downsampling and upsampling paths. Each block consists of two 2D convolution layers. Unlike most implementations of the U-Net architecture [11,15,16], we omit normalization (i.e., batch or instance normalization) and do not apply a non-linear activation function like ReLU and use the identity function instead (see Section 4). Additionally, we use the PocketNet approach proposed in [7] and leave the number of feature maps (channels) at each resolution within our architecture constant, namely 32. Figure 1 sketches the U-Net architecture for di!erent network depths. We define the depth, d, of each network to be the number of downsampling operations present in the network. The size of the feature tensor is halved from depth d to depth d + 1.
We test our proposed approach for elliptic problems on four di!erent test cases. Unless otherwise specified, the domain is the unit square.
The Bubble Function. In the first test case, (1) is solved with data chosen such that the exact solution is the bubble function
This function is symmetric, and we apply homogeneous zero Dirichlet boundary conditions. Note that for this function, one can show that the five-point stencil scheme given by (2) applied to the bubble function gives exactly ↓$u evaluated at the grid point (x i , y j ), which implies that the finite di!erence approximation is exact.
The "Peak" Function. In the second case, we want to test the e!ectiveness of our strategy when the solution has high gradients and is non-symmetric. Hence, we select data such that the exact solution is the "peak" function
The Exponential-Trigonometric Function. We test our approach on a function with non-homogeneous Dirichlet boundary conditions in the third case.
We select data such that the exact solution is the exponential-trigonometric function given by u(x, y) = e →x 2 →y 2 sin(3ςx) sin(3ςy) + x.
Non-Constant and Discontinuous Di!usion. To test how our approach fairs on problems with low regularity, we consider an example with a non-constant, discontinuous di!usion coe"cient taken from [17]. The problem is posed on the square domain # = (↓1, 1) 2 which we divide into four subdomains # i corresponding to the quadrants of the Cartesian plane. On each subdomain, ω is constant and takes the values ω 1 = ω 3 = 5 and ω 2 = ω 4 = 1. The exact solution takes the form
where (r, ↼) are the polar coordinates of a given point in #, a i , b i are constants that depend on the subdomains (see [17] for exact values). This discontinuity in the di!usion coe"cient introduces a singularity in the solution at the origin, namely the function belongs to the Sobolev space H 1+ς (#), with φ ⇑ 0.5354.
We test our proposed approach for parabolic problems on three di!erent test cases.
The Trigonometric Function. In the first two cases, ( 17) is solved with data such that the exact solution is given by the trigonometric function u(x, y, t) = cos(t) sin(nςx) sin(nςy).
(
In the first case, we select n = 1, which gives a symmetric function with a single peak. In the second case, we set n = 4, which results in 16 peaks and troughs in our domain.
The Gaussian Function. In the third case, we test our approach on a Gaussian function centered on the point 1 2 , 1 2 u(x, y, t) = cos(t)e →50((2x→1
We use the Adam optimizer [18] and 5↗5 convolutional kernels in each layer. The initial learning rate is 0.001 for the steady-state elliptic problems and 0.0001 for parabolic problems. During training, we use L2 regularization with a penalty of 10 →7 and the norm of the network's gradient is clipped so that it is no greater than 10 →2 . Our models are implemented in Python using TensorFlow (v2.12.0) and trained on an NVIDIA Quadro RTX 6000 GPU [19]. All network weights are initialized using the default TensorFlow initializers. All other hyperparameters are left at their default values. The code for our network architecture is available at https://github.com/aecelaya/pde-nets.
To assess the accuracy of our predicted solutions, we use the following norms of the error between the exact solution u and its approximate û for steady-state problems:
For time-dependent problems, similar errors are computed at the final time T .
We use Algorithm 1 to approximate a solution û to (1). In this algorithm, we start by initializing the neural network N ε with trainable parameters ω. In this case, we use the U-Net architecture described in Section 2.3. We set the maximum number of optimization steps M as the stopping criterion for the algorithm. We then begin the minimization (or training) process by generating a prediction from our neural network with the source term f as the input. Given the current prediction, we compute the loss value and update the network weights via backpropagation. If the loss value from the current prediction is less than the previous best value, we update our best loss and save the current prediction. Algorithm 1 fully outlines this procedure.
We use Algorithm 2 to approximate a solution ûN T to (17), where N T is an integer such that T = N T ϱ . Like with Algorithm 1, we randomly initialize a U-Net and set a maximum number of optimization steps as a stopping
Update ω using backpropagation on ⇀
9: if ⇀ < ⇀ ↔ then 10: û↔ ⇓ û ↽ Save best prediction 11: ⇀ ↔ ⇓ ⇀ ↽ Update smallest loss value 12:
k ⇓ k + 1 criterion. However, in this case, the maximum number of steps applies to each time step. We then set the initial condition as the previous solution.
We apply nearly the same steps for each time step as with Algorithm 1. Namely, we produce a candidate prediction from our neural network with the previous solution as the input. We compute the loss value and update the network weights via backpropogation on the loss. If the loss value is lower than the previously observed lowest loss value, then we update our lowest loss value and save the current prediction. At the end of this process, we set the previous solution to the best current solution. This process continues for every time step. Algorithm 2 fully outlines this procedure. Note that in this algorithm, we do not reinitialize our neural network weights at each time step. Instead, we use the previous network configuration as the initial state for the next time step.
Using the methods described in Section 2 and Algorithm 1, we approximate the finite di!erence solution to the elliptic problems defined in Section 2.4.1 with constant di!usion coe"cients (i.e., ω = 1). Table 1 shows the accuracy of our unsupervised predictions for a varying grid size and number of optimization steps. We fix the depth of our U-Net architecture to three in while k < M do 9:
nϱ) ↽ Compute loss 11: Update ω using backpropagation on ⇀ 12: if ⇀ < ⇀ ↔ then 13: ûn ⇓ w ↽ Save best prediction 14: ⇀ ↔ ⇓ ⇀ ↽ Update smallest loss value 15: k ⇓ k + 1 16:
n ⇓ n + 1 this case. For reference, Table 3 shows the accuracy of the finite di!erence method for varying grid sizes. For the peak and exponential-trigonometric functions, our approach provides a good approximation to the finite di!erence solution. We do not, however, see this for the bubble function example.
In that case, the solution of the finite di!erence method is exact. Under our selected settings (i.e., depth and optimization steps), our method stops at an error of approximately 10 →6 in the ||•|| 2,h norm for the bubble function example. The cause of this discrepancy between our method applied to the bubble function and the finite di!erence solution may be that current neural network optimizers like Adam are stochastic, gradient-based, and do not have optimal step sizes. We hypothesize that our neural network is getting stuck in a neighborhood around the global minimum because of these previously listed limitations of optimizers like Adam. Hence, using secondorder, Hessian-based optimizers or better step-size computations may allow for faster convergence to more accurate solutions for our neural network-based approach.
Figure 2 shows the true solution, predicted solution, and absolute di!erence for the bubble, peak, and exponential-trigonometric cases on a 128↗128 grid. This figure shows that for every case, our unsupervised algorithm produces visually indistinguishable solutions from the true solution in each case.
We also want to study the network depth's e!ect on our predictions' accuracy. Table 2 shows the accuracy of our unsupervised neural network predictions for varying grid sizes and U-Net depths. We set the number of optimization steps to 4,000. Here, we see that the depth of the U-Net architectures does not appear to have a significant e!ect on the accuracy of our predictions. The only notable exceptions are for the bubble and exponentialtrigonometric functions on a 128↗128 grid with a depth equal to two. This discrepancy may be caused by the network not capturing su"ciently rich features on the finer grid. The fact that we see a decrease in the errors between our predictions and the true solution to those comparable to the finite di!erence method as we increase the depth to three or greater supports this hypothesis. Note that for coarser grids (i.e., fewer grid points), we do not test higher network depths since the size of the features (the output of each layer) at the coarser grids in such networks would be smaller than the size of the convolutional kernels in the network.
The results of Algorithm 1 for the non-constant di!usion problem defined in Section 2.4.1 are shown in Table 4. Here, we see more significant errors than with the constant di!usion problems. Figure 3 shows the true solution, predicted solution, and absolute di!erence for the non-constant diffusion problem on a 128↗128 grid. This figure shows that the errors are mostly concentrated around the discontinuities in ω. We do not compare this to the finite di!erence method because the exact solution exhibits a singularity at the origin. Indeed, the gradient of the exact solution blows up at the point (0, 0).
Because the finite di!erence solution solves the convex minimization problem (7), the finite di!erence approximation's accuracy limits our approach's accuracy. This explains the relatively poor performance in the case of discontinuous di!usion. It is well known that, due to the singularity at the origin, the exact solution to the discontinuous di!usion problem belongs to the Sobolev space H 1+ς (#), where φ ⇑ 0.5354 [17]. Hence, the finite di!erence method, which requires more regularity, performs poorly.
Using the methods described in Section 2 and Algorithm 2, we approximate the finite di!erence solution to the parabolic problems defined in Section Problem Table 2: Accuracy of unsupervised predictions for varying grid sizes N and U-Net depths d. We set the number of optimization steps to 4,000.
2.4.2. The time step is ϱ = 0.1. Except for the first time step, we set the number of optimization steps per time step to 250. Because we start with randomly initialized weights in the first time step, we set the number of optimization steps to 1,000 in the first time step. Table 6 shows the accuracy of our unsupervised predictions for varying grid sizes. For reference, Table
Problem
1.3582e-16 2.8449e-16 2.4524e-16 6.8695e-16 8.1134e-16 1.9776e-15 Peak 1.7597e-02 1.1151e-01 5.3441e-03 3.1282e-02 1.3988e-03 7.9696e-03 Exp-Trig 2.2986e-03 7.3266e-03 5.7262e-04 1.8251e-03 1.4303e-04 4.5588e-04 Table 3: Finite di!erence errors for selected elliptic problems for comparison.
2 283,713 3.9709e-02 1.5798e-01 3.1294e-02 1.2265e-01 3.3859e-02 1.4377e-01 3 412,225 3.9713e-02 1.5794e-01 3.1224e-02 1.2258e-01 4.9868e-02 1.1281e-01 4 541,889 --3.1110e-02 1.2282e-01 3.8937e-02 1.0907e-01 5 669,249 ----3.6566e-02 1.0511e-01
Table 5: Accuracy of unsupervised predictions on non-constant di!usion problem for varying grid sizes N and U-Net depths d. We set the number of optimization steps to 4,000.
7 shows the accuracy of the finite di!erence method with backward Euler for varying grid sizes. Here, we see that Algorithm 2 achieves comparable accuracy to the finite di!erence method with backward Euler. However, we do not see the same convergence as with the finite di!erence approach for the trigonometric functions. The small number of optimization steps can explain this lack of convergence and, again, the weakness of first-order optimizers like Adam. Figures 4 through 6 show the true solution, predicted solution, and absolute di!erence for the first ten time steps of each problem. These figures show that our unsupervised approach produces visually accurate solutions. Finally, Table 8 displays the errors for several choices of activation functions at time t = 0.5. The errors are similar for ReLU, Tanh, Swish [20], and identity activation functions. These results indicate that linear activation functions (i.e., the identity) are su"cient for learning solutions to time-dependent problems.
Problem t
0.5 7.6825e-04 1.5388e-03 1.0661e-03 2.1703e-03 1.1562e-03 2.4109e-03 1.0 5.2059e-04 1.0459e-03 7.0987e-04 1.4189e-03 7.7541e-04 1.5364e-03 2.5 6.1710e-04 1.2354e-03 8.6825e-04 1.7483e-03 9.4440e-04 1.9382e-03 5.0 1.6629e-04 3.3075e-04 2.5692e-04 5.2649e-04 3.0442e-04 6.2141e-04 Trig. (n = 4) 0.5 6.0011e-03 1.1977e-02 1.3953e-03 2.9254e-03 3.2371e-04 9.0033e-04 1.0 3.6984e-03 7.3841e-03 8.5877e-04 1.8117e-03 1.9062e-04 5.3841e-04 2.5 5.4592e-03 1.0897e-02 1.3062e-03 3.1186e-03 8.4891e-04 2.8031e-03 5.0 1.9352e-03 3.8744e-03 4.8639e-04 1.3043e-03 9.8322e-04 2.5854e-03 Gaussian 0.5 3.0932e-03 3.8826e-02 7.1997e-04 1.0341e-02 1.8902e-04 2.5298e-03 1.0 1.9069e-03 2.3921e-02 4.4389e-04 6.3892e-03 1.3421e-04 1.4701e-03 2.5 2.8197e-03 3.5385e-02 6.5362e-04 9.5078e-03 2.7103e-04 2.3168e-03 5.0 1.0004e-03 1.2579e-02 2.8684e-04 2.6686e-03 3.1119e-04 9.0614e-04 Table 6: Accuracy of unsupervised predictions for varying grid sizes and time steps for parabolic problems. The depth of the network is set to three, the number of optimization iterations at each time step to 250, and the time step to 0.1. Table 8: Accuracy of unsupervised predictions for varying grid sizes and activation functions at time 0.5. The depth of the network is set to three, the number of optimization steps to 250, and the time step is 0.1.
labeled dataset. We randomly initialize the neural network weights in Algorithms 1 and 2 (first time step only). This random initialization may result in predictions with high errors at the beginning of training. Without optimal step sizes in our optimizer, these poor solutions may result in slow convergence towards an optimal solution. Using transfer learning with a pre-trained network that is trained in a supervised setting (i.e., with labeled training data), the initial predictions from our algorithms may be already close to optimal, resulting in faster convergence. One could view the use of these pre-trained networks as a sort of preconditioner for our neural network-based approach. Indeed, we see the benefits of transfer learning with Algorithm 2 after the first time step. Instead of reinitializing the weights for the subsequent time steps, we resuse the weights from the previous solution. That allows us to use a small number of optimization steps (i.e., M = 250) at each time for parabolic problems.
It is important to note that our methods for estimating the finite di!erence solutions do not explicitly construct the finite di!erence matrix. Instead, we use our neural networks to map the source term f to the approximate finite di!erence solution. In this sense, we are learning the inverse mapping of the finite di!erence matrix. From a numerical point of view, this provides the advantage of not having to construct or store a finite di!erence matrix. Instead, we implement the required stencil(s) for our problem in the loss function and apply them appropriately. Only having to implement the stencils and not construct a finite di!erence matrix is also an advantage from an implementation point-of-view.
The justification for the lack of non-linear activation functions in our network is that the relationship between the source term f and the finite di!erences solution u h is given by a system of linear equations of the form Au h = f . In other words, we know a priori that our goal is to learn a linear relationship between f and our estimated solution. In many machine learning applications, the exact nature of the input and output relationship is unknown and assumed to be highly non-linear. Hence, including non-linear activation functions results in a network that learns a non-linear relationship between inputs and outputs. The use of the identity function as an activation puts us outside of the scope of the universal approximation theorem [21,22,23]. However, in our case, we do not need our neural network to be a universal approximator. The role of the network is to learn a linear relationship for a single example.
We utilize the PocketNet approach proposed by [7] in our proposed algorithms. This approach takes advantage of the similarity between the U-Net architecture and geometric multigrid methods to drastically reduce the number of parameters, while maintaining the same accuracy as conventional CNNs for medical imaging and scientific machine learning tasks [7,24,25]. Additionally, we replace transposed convolution with bilinear upsampling. We find that these changes save time and memory and yield the same accuracy that we see using conventional CNNs (i.e., doubling the number of channels at every depth). These results indicate that smaller neural networks (in terms of parameters) can achieve high accuracy for scientific machine learning tasks.
We see in Tables 2 and 5 that the depth of the U-Net architecture does not generally have a significant e!ect on our results. This indicates that, regardless of depth, the U-Net architecture is su"ciently expressive to learn a mapping from the right hand side f to an approximation of the finite di!erence solution u h . However, non-U-shaped architectures like the HRNet may also produce similar or improved results [26]. Additionally, the use of residual or dense connections within the convolutions of our architecture may also be beneficial [27,28]. Such block designs have been shown to speed up convergence to lower loss values for neural networks [29]. Finally, modifying our existing architecture by adding deep supervision could also speed up convergence to lower loss values in fewer iterations [30,31].
The use of the weighting parameter ϖ in ( 8) and ( 19) is necessary to enforce the given boundary conditions. We arrived at our proposed values of ϖ via a grid search over a range of possible values. However, finding optimal values of weighting parameters like ϖ is an open question [32]. Future work will explore this topic further.
The results presented above show the e!ectiveness of our proposed unsupervised approaches for estimating the finite di!erence solution to elliptic and parabolic problems. Unlike current neural network-based methods for estimating the solutions of PDEs, our algorithms do not rely on training data. Instead, they leverage the machinery of neural networks to estimate the finite di!erence solution to a given problem directly. As a result, our approach does not su!er from the generalization gap seen in PINNs or other neural network-based approaches, which rely on training data.
Our approach could benefit finite di!erence solvers by producing better initial guesses and/or acting as a preconditioner. With a few iterations, Algorithm 1 can produce good initial guesses for iterative solvers, thereby reducing the number of iterations required to solve the linear system resulting from (2) or (9). This same idea can also apply to time-dependent problems, but with initial guesses being produced at each time step. Additionally, because we employ identity for activation functions, our neural networks are linear. Hence, it may be possible to represent a pretrained architecture as a matrix. This resulting matrix could then be used as a preconditioner for finite di!erence solvers. This use of preconditioners and further testing on other kinds of PDEs, like convection-di!usion and nonlinear problems, will be the object of future work.