A physics-informed deep learning approach for solving strongly degenerate parabolic problems

In recent years, Scientific Machine Learning (SciML) methods for solving partial differential equations (PDEs) have gained increasing popularity. Within such a paradigm, Physics-Informed Neural Networks (PINNs) are novel deep learning frameworks for solving initial-boundary value problems involving nonlinear PDEs. Recently, PINNs have shown promising results in several application fields. Motivated by applications to gas filtration problems, here we present and evaluate a PINN-based approach to predict solutions to strongly degenerate parabolic problems with asymptotic structure of Laplacian type. To the best of our knowledge, this is one of the first papers demonstrating the efficacy of the PINN framework for solving such kind of problems. In particular, we estimate an appropriate approximation error for some test problems whose analytical solutions are fortunately known. The numerical experiments discussed include two and three-dimensional spatial domains, emphasizing the effectiveness of this approach in predicting accurate solutions.


Introduction
In this paper, we aim to exploit a novel Artificial Intelligence (AI) methodology, known as Physics-Informed Neural Networks (PINNs), to predict solutions to Cauchy-Dirichlet problems of the type where Ω is a bounded connected open subset of R n (2 ≤ n ≤ 3) with Lipschitz boundary, f and w are given real-valued functions defined over Ω × [0, T ] and ∂ par Ω T respectively, ∇u denotes the spatial gradient of an unknown solution u : Ω × [0, T ) → R, while ( • ) + stands for the positive part.
A motivation for studying problem (1.1) can be found in gas filtration problems (see [1] and [3]).In order to make the paper self-contained, we provide a brief explanation in Section 1.1 below.
As for the parabolic equation (1.1) 1 , the regularity properties of its weak solutions have been recently studied in [2,3] and [8].The main novelty of this PDE is that it exhibits a strong degeneracy, coming from the fact that its modulus of ellipticity vanishes in the region {|∇u| ≤ 1}, and hence its principal part behaves like a Laplace operator only at infinity.
The regularity of solutions to parabolic problems with asymptotic structure of Laplacian type had already been investigated in [11], where a BMO 1 regularity was proved for solutions to asymptotically parabolic systems in the case f = 0 (see also [13], where the local Lipschitz continuity of weak solutions with respect to the spatial variable is established).In addition, we want to mention the results contained in [4], where nonhomogeneous parabolic problems with an asymptotic regularity in divergence form of p-Laplacian type are considered.There, Byun, Oh and Wang establish a global Calderón-Zygmund estimate by converting a given asymptotically regular problem to a suitable regular problem.
Concerning the approach used here, the PINNs are a Scientific Machine Learning (SciML) technique based on Artificial Neural Networks (ANNs) with the feature of adding constraints to make the predicted results more in line with the physical laws of the addressed problem.The concept of PINNs was introduced in [12,15,16] and [17] to solve PDE-based problems.The PINNs predict the solution to a PDE under prescribed initial-boundary conditions by training a neural network to minimize a cost function, called loss function, which penalizes some suitable terms on a set of admissible functions u (for more information, we refer the interested reader to [6]).
The kind of approach we want to propose here can offer effective solutions to real problems such as (1.1) and can be applied in many other different fields: for example, in production and advanced engineering [19], for transportation problems [7], and for virtual thermal sensors using real-time simulations [10].Additionally, it is employed to solve groundwater flow equations [5] and address petroleum and gas contamination [18].
As far as we know, this is one of the first papers demonstrating the effectiveness of the PINN framework for solving strongly degenerate parabolic problems of the type (1.1).

Motivation
Before describing the structure of this paper, we wish to motivate our study by pointing out that, in the physical cases n = 2 and n = 3, degenerate equations of the form (1.1) 1 may arise in gas filtration problems taking into account the initial pressure gradient.
The existence of remarkable deviations from the linear Darcy filtration law has been observed in several systems consisting of a fluid and a porous medium (e.g., the filtration of a gas in argillous rocks).One of the manifestations of this nonlinearity is the existence of a limiting pressure gradient, i.e. the minimum value of the pressure gradient for which fluid motion takes place.In general, fluid motion still occurs for subcritical values of the pressure gradient, but very slowly; when achieving the limiting value of the pressure gradient, there is a marked acceleration of the filtration.Therefore, the limiting-gradient concept provides a good approximation for velocities that are not too low.
In accordance with some experimental results (see [1]), under certain physical conditions one can take the gas filtration law in the very simple form where v = v(x, t) is the filtration velocity, k is the rock permeability, µ is the gas viscosity, p = p(x, t) is the pressure and β is a positive constant.Under this assumption we obtain a particularly simple expression for the gas mass velocity (flux) j, which contains only the gradient of the pressure squared, exactly as in the usual gas filtration problems: where ϱ is the gas density and C is a positive constant.Plugging expression (1.2) into the gas mass-conservation equation, we obtain the basic equation for the pressure: where m is a positive constant.Equation (1.3) implies, first of all, that the steady gas motion is described by the same relations as in the steady motion of an incompressible fluid if we replace the pressure of the incompressible fluid with the square of the gas pressure.In addition, if the gas pressure differs very little from some constant pressure p 0 , or if the gas pressure differs considerably from a constant value only in regions where the gas motion is nearly steady, then the equation for the gas filtration in the region of motion can be "linearized" following L. S. Leibenson, and thus obtaining (see [1] again) Setting u = p 2 and performing a suitable scaling, equation (1.4) turns into which is nothing but equation (1.1) 1 with f ≡ 0. This is why (1.1) 1 is sometimes called the Leibenson equation in the literature.
The paper is organized as follows.Section 2 is devoted to the preliminaries: after a list of some classic notations, we provide details on the strongly degenerate parabolic problem (1.1).In Section 3, we describe the PINN methodology that was employed.Section 4 presents the results that were obtained.Finally, Section 5 provides the conclusions.

Notation and preliminaries
In what follows, the norm we use on R n will be the standard Euclidean one and it will be denoted by | • |.In particular, for the vectors ξ, η ∈ R n , we write ⟨ξ, η⟩ for the usual inner product and |ξ| := ⟨ξ, ξ⟩ for the backward parabolic cylinder with vertex (x 0 , t 0 ) and width ρ.Finally, for a general cylinder Q = A × (t 1 , t 2 ), where A ⊂ R n and t 1 < t 2 , we denote by the usual parabolic boundary of Q.
To give the definition of a weak solution to problem (1.1), we now introduce the function H : R n → R n defined by ) is a weak solution of equation (1.1) 1 if and only if for any test function φ ∈ C ∞ 0 (Ω T ) the following integral identity holds: (2.1) We identify a function as a weak solution of the Cauchy-Dirichlet problem (1.1) if and only if (2.1) holds and, moreover, Therefore, the initial condition u = w on Ω × {0} has to be understood in the usual L 2 -sense (2.2), while the condition u = w on the lateral boundary ∂Ω × (0, T ) has to be meant in the sense of traces, i.e. (u − w) (•, t) ∈ W 1,2 0 (Ω) for almost every t ∈ (0, T ).
Taking p = 2 and ν = 1 in [3, Theorem 1.1], we immediately obtain the following spatial Sobolev regularity result: is a weak solution of equation (1.1) 1 .Then the solution satisfies Furthermore, the following estimate ⋐ Ω T and a positive constant c depending on n, q and R 0 .
From the above result one can easily deduce that u admits a weak time derivative u t , which belongs to the local Lebesgue space L min {2, q} loc (Ω T ).The idea is roughly as follows.Consider equation (1.1) 1 ; since the previous theorem tells us that in a certain pointwise sense the second spatial derivatives of u exist, we may develop the expression under the divergence symbol; this will give us an expression that equals u t , from which we get the desired summability of the time derivative.Such an argument has been made rigorous in [3,Theorem 1.2], from which we can derive the next result.
Theorem 2.4.Under the assumptions of Theorem 2.3, the time derivative of the solution exists in the weak sense and satisfies Furthermore, the following estimate holds true for any parabolic cylinder ⋐ Ω T and a positive constant c depending on n, q and R 0 .
Now, let the assumptions of Theorem 2.3 be in force.For ε ∈ [0, 1] and a couple of standard, non-negative, radially symmetric mollifiers where f is meant to be extended by zero outside Ω T .Observe that f 0 = f and f ε ∈ C ∞ (Ω T ) for every ε ∈ (0, 1]. Next, we consider a domain in space-time denoted by Ω ′ 1,2 := Ω ′ × (t 1 , t 2 ), where Ω ′ ⊆ Ω is a bounded domain with smooth boundary and (t 1 , t 2 ) ⊆ (0, T ).In the following, we will need the definitions below.
) the following integral identity holds: as a weak solution of the Cauchy-Dirichlet problem if and only if (2.4) holds and, moreover, in the usual L 2 -sense and the condition u ε = u on the lateral boundary Due to the strong degeneracy of equation (1.1) 1 , in order to prove Theorems 2.3 and 2.4 above, the authors of [3] resort to the family of approximating parabolic problems (2.5).These problems exhibit a milder degeneracy than (1.1) and the advantage of considering them stems from the fact that the existence of a unique energy solution u ε satisfying the requirements of Definition 2.6 can be ensured by the classic existence theory for parabolic equations (see [14, Chapter 2, Theorem 1.2 and Remark 1.2]).

Physics-informed methodology
PINNs are a type of SciML approach used in neural networks to solve PDEs.Unlike traditional neural networks, PINNs incorporate physics constraints into the model, resulting in predicted outcomes that adhere more closely to the natural laws governing the specific problem being addressed.The general form of the problem involves a PDE along with initial and/or boundary conditions.
In particular, we consider a (well-posed) problem of the type where Ω is a bounded domain in R n , F denotes a nonlinear differential operator, γ is a parameter associated with the physics of the problem, B is an operator defining arbitrary initial-boundary conditions, the functions f and w represent the problem data, while u(x, t) denotes the unknown solution.The objective of PINNs is to predict the solution to (3.1) by training the neural network to minimize a cost function.The neural network's architecture used for PINNs is typically a FeedForward fully-connected Neural Network (FF-DNN), also known as Multi-Layer Perceptron (MLP).In an FF-DNN, information flows only forward direction, in the sense that the neural network does not form a loop.Furthermore, all neurons are interconnected.Once the number N of hidden layers has been chosen, for any i ∈ {1, . . ., N } and set z = (x, t) we define where W i is the weights matrix of the links between the layers i − 1 and i, while b i corresponds to the biases vector.Then, a generic layer of the neural network is defined by for some nonlinear activation function φ i .The output of the FF-DNN, denoted by ûθ (z), can be expressed as a composition of these layers by where θ represents the set of hyperparameters of the neural network and the activation function φ is assumed to be the same for all layers.To solve the differential problem (3.1) using PINNs, the PDE is approximated by finding an optimal set θ * of neural network hyperparameters that minimizes a loss function L. This function consists of two components: the former, denoted by L F , is related to the differential equation, while the latter, here denoted by L B , is connected to the initial-boundary conditions (see Fig. 3.1).In particular, the loss function can be defined as follows where ω F and ω B represent the weights that are usually applied to balance the importance of each component.Hence, we can write The aim of this approach is to approximate the solution of the PDE satisfying the initialboundary conditions.This is known in the literature as the direct problem, which is the only one we will address here.

Numerical results
In this section we evaluate the accuracy and effectiveness of our predictive method, by testing it with five problems of the type (1.1) whose exact solutions are known.For each problem, we will denote the exact solution by u, and the predicted (or approximate) solution by û.Sometimes, by abuse of language, for a given time t ≥ 0 we will refer to the partial maps u(•, t) and û(•, t) as the exact solution and the predicted (or approximate) solution respectively.The meaning will be clear from the context every time.We will deal with each test problem separately, so that no confusion can arise.In the first three problems, Ω will be a bounded domain of R 2 , while, in the last two problems, Ω will denote the open unit sphere of R 3 centered at the origin.
In addition, for each of the test problems, we employed the same neural network architecture.This consists of four layers, each with 20 neurons.We utilized the hyperbolic tangent function as the activation function for both the input layer and the hidden layers, while a linear function served as the activation function for the output layer.Lastly, to train the neural network, we conducted 80000 epochs with a learning rate (lr) of 3×10 −3 and employed the Adaptive Moment Estimation (ADAM) optimizer.The decision to set the lr to the constant value 3 × 10 −3 was based on the observation that this specific hyperparameter led to the optimal convergence of our method.Experimentation with lr set to 1 × 10 −1 highlighted the network's inability to achieve convergence, while using an lr of 1 × 10 −5 allowed the method to converge, albeit requiring a significantly higher number of epochs.The latter scenario, while ensuring convergence, proved to be less computationally efficient.The experiments were performed on a NVIDIA GeForce RTX 3080 GPU with AMD Ryzen 9 5950X 16-Core Processor and 128 GB of RAM.

First test problem
The first test problem that we consider is where The exact solution of this problem is given by Therefore, for any fixed time t ≥ 0 the graph of the function u(•, t) is an elliptic paraboloid.
As time goes on, this paraboloid slides along an oriented vertical axis at a constant velocity,  What has been verified is that the plot of the predicted solution û(•, t) has precisely the same shape and geometric properties as the graph of the exact solution u(•, t), for both short and long times t.Moreover, the time evolution of the approximate solution û exactly mirrors the behavior described for the known solution u.A further interesting aspect that can be noticed is that the level curves of the approximate solution û(•, t) overlap almost perfectly those of u(•, t), provided that t is not very large (see Fig.  We have also noted that, at time t = 0, the approximate solution is basically equal to zero in a very tiny region around the origin (0, 0) of the xy-plane.This means that the said region is composed of "numerical zeros" of the solution predicted at time t = 0, while we know that u(x, y, 0) = 0 if and only if (x, y) = (0, 0).However, this discrepancy is actually negligible, since the order of magnitude of u(x, y, 0) is not greater than 10 −6 within the above region.
To assess the accuracy of our predictive method and the numerical convergence of the solution û toward u in a more quantitative way, we now look at the time behavior of the L 2 -error ∥û(•, t) − u(•, t)∥ L 2 (Ω) by considering the natural quantities and Passing from Cartesian to polar coordinates, one can easily find that and therefore .
The results that we have obtained are shown in Table 1.The estimates of E(1) and E(10) are equal to 2.24 × 10 −5 and 4.30 × 10 −4 respectively, which is very satisfactory, especially considering that the order of magnitude of E rel (1) and E rel (10) is equal to 10 −6 .
In order to get more accurate estimates for larger values of T , for every fixed T ≥ 100 we have used 2.1 × T equispaced points (instead of the initial 21) to discretize the time interval [0, T ].By doing so, we have observed that the variation of the estimate of E(T ) displays a monotonically increasing behavior, in accordance with the definition (4.1).However, even for 100 ≤ T ≤ 500, the order of magnitude of E rel (T ) remains not greater than 10 −6 .Therefore, for this first test problem, we can conclude that our predictive method is indeed very accurate and efficient, on both a short and long time scale.
which is nothing but the approximating problem (2.5) associated with (P1).In what follows, we will denote the exact solution of (4.3) by u ε , while the predicted solution will be denoted by ûε .Throughout our tests, for 10 −9 ≤ ε ≤ 10 −3 , for Ω ′ = Ω and for (t 1 , t 2 ) = (0, T ), we have observed that the plots of the predicted solution ûε (•, t) and the exact solution u(•, t) share the same configurations and geometric peculiarities, on both a short and long time scale (see, e.g. Figure 4.3).Furthermore, we have seen that the evolution over time of ûε reflects the behavior depicted for the solution u quite faithfully.In addition, the contour lines of ûε (•, t) perfectly overlap those of u(•, t), at least for not very long times t (see Fig.   Let us now assume that where z 0 = (0, 2) = (0, 0, 2).Then, the limit in (2.6) suggests that ûε should numerically converge to u as ε ↘ 0. To obtain a numerical evidence of such convergence, we have chosen Ω ′ = B 1/2 (0) and (t 1 , t 2 ) = ( 7 4 , 2) into (4.3) and examined the time behavior of the L 2 -error ∥û ε (•, t) − u(•, t)∥ L 2 (Ω ′ ) , by evaluating the quantities and .
Table 2 shows the results obtained and reveals that the predicted solution ûε converges to u as ε tends to zero, although not very quickly.In fact, the estimates of E ε and E rel (ε) approach zero with a convergence rate much lower than that of ε.Furthermore, they seem to start decreasing monotonically, i.e. without oscillations, for ε ≤ 10

Second test problem
Let α > 0. As a second test problem we consider where The exact solution of problem (P2) is given by u(x, y, t) ≡ u α (x, y, t) := t (x 2 + y 2 ) α .
At any fixed time t > 0, the shape and geometric properties of the graph of u(•, t) strongly depend on the value of the parameter α.
If α = 1 2 , then the graph of u(•, t) is a cone whose vertex coincides with the origin (0, 0, 0) at any given positive time t.As time goes on, the cone in question gets narrower and narrower around the vertical axis.In this case, the plot of the approximate solution û(•, t) has the same form as the graph of the exact solution u(•, t) for both short and long times t > 0, except near the origin, where the tip of the cone appears to have been smoothed out (see Figure 4.5, center).However, this is not a surprise at all, since we already know that for t > 0 the function is not differentiable at the center (0, 0) of Ω.When 0 < α < 1  2 , the graph of u(•, t) is cusp-shaped for any fixed time t > 0, the origin now being a cusp for all positive times.In this case, a loss on convexity occurs, which is also observed in the plot of the predicted solution û(•, t) for all times t > 0 (see, e.g. Figure 4.5, left).
Lastly, when α > 1 2 the graph of u(•, t) is no longer cusp-shaped and becomes increasingly narrow around the vertical axis as t increases.Furthermore, for any fixed t > 0 the exact solution u(•, t) is convex again and its graph gets flatter and flatter near the origin when α >> 1 (see Figure 4.6).
In all three of the above cases, we have noticed that the plot of û(•, t) is basically identical in its shape and geometry to the graph of the exact solution u(•, t), for both short and long periods t.Moreover, also for problem (P2) we have verified that the time evolution of the predicted solution faithfully reflects the trend described for the exact solution in all three previous cases.Therefore, we may conclude that α = 1 2 represents a critical value for the global behavior of both the exact and the predicted solution.Later, we have examined the contour lines of ûα (•, t) for α ∈ {0.3, 0.5, 1.3, 5} and for not very large times t > 0. For every fixed α ∈ {0.3, 0.5, 1.3}, the level curves of ûα (•, t) overlap quite well those of u α (•, t), with some small differences between one case and the other.More precisely, for each α ∈ {0.5, 1.3} the contour lines corresponding to the same level are almost indistinguishable, at least for not very long times t (see, for example, Figure 4.7, where t = 0.5).For α = 5 and t > 0, we have also noted that the approximate solution is essentially equal to zero in a fairly large region Σ t around the origin (0, 0) of the xy-plane (see Fig. 4.8).As already said for problem (P1), this means that such region consists of numerical zeros of û5 (•, t), while for t > 0 we know that u 5 (x, y, t) = 0 if and only if (x, y) = (0, 0).However, this discrepancy is reasonably small for short times, since the order of magnitude of u 5 (x, y, t) does not exceed 10 −2 within Σ t for 0 < t ≤ 10.To evaluate in a more quantitative manner the accuracy of our method in solving problem (P2) and the numerical convergence of the solution û toward u, we may now consider again the quantities (4.1) and (4.2).Passing from Cartesian to polar coordinates, we find so that we now have During the experimental phase, we have estimated E(T ) and E rel (T ) for α ∈ {0.3, 0.5, 1.3, 5} and T ∈ {1, 10, 20, 40, 100}.The results that we have obtained are reported in Tables 3−6 and show that, for any fixed value of α, the estimate of E(T ) follows an increasing trend, as prescribed by (4.1).Furthermore, by analyzing the orders of magnitude of E(T ) and E rel (T ), we may affirm that our approach provides very accurate predictions, on both a short and long-term scale.

Third test problem
We shall now consider the problem where Ω = (−1, 1) × (−1, 1), The exact solution of this problem is given by Therefore, for any fixed time t ≥ 0, the graph of the function u(•, t) is given by the union of the horizontal region and the sliding plane Let us denote by G t := H t ∪ I t the graph of u(•, t).Then, as time goes by, the set G t slides along a vertical axis with a constant velocity and no deformation, since ∂ t u ≡ 1 over Ω T (see Fig. 4.9, above).The plot of the approximate solution û(•, t) roughly resembles that of u(•, t) for both short and long times t ≥ 0, except near the joining line H t ∩ I t , where the graph of the solution appears to have been slightly smoothed (see Figure 4.9, below).However, this is not surprising at all, since we already know that, for any fixed t ≥ 0, the function u(•, t) : Ω → R defined by (4.4) is not differentiable at any point of the open segment S 0 := {(x, y) ∈ Ω : x = 0}.This fact also has repercussions in the comparison between the level curves of u(•, t) and û(•, t), whose superposition is far from being perfect on approaching the segment S 0 from the right, i.e. for x > 0 (see Fig. 4.10).
Furthermore, we have also observed that the evolution of û over time accurately reflects the evolution of the set G t described above.In order to assess in a more quantitative way the accuracy of our method in solving (P3) and the distance between the solutions u and û, we resort again to the quantities defined in (4.1) and (4.2).Through an easy calculation, we get so that we now have Proceeding as for the previous problems, we have estimated E(T ) and E rel (T ) for T ∈ {1, 10, 100, 200, 300}.
Table 7 contains the results obtained and reveals that the estimate of E(T ) exhibits again an increasing behavior, as expected from (4.1).Furthermore, from this table, it seems that the asymptotic trend of the estimate of E rel (T ) may encounter a sort of plateau at T = 100, after which convergence sensibly slows down.We do not know whether this is a typical behavior, since we cannot draw information from (4.5) in this sense.In fact, from the definition of E rel (T ) it is not possible to predict what the combined effect of E(T ) and sup is the ratio of two functions which are both increasing with respect to T and we cannot determine a priori the growth rate of E(T ).Nevertheless, by carefully examining the orders of magnitude of both E(T ) and E rel (T ), we can conclude that our method produces accurate results also in this case, in both short and long-term predictions.

Fourth test problem
We now consider the problem where Ω = {(x, y, z) ∈ R 3 : x 2 + y 2 + z 2 < 1}.This problem is the 3-dimensional version of problem (P1) and its exact solution is given by To evaluate the accuracy of our method in solving problem (P4) and the distance between the predicted solution û and the exact solution u, we confined ourselves to considering the quantities (4.1) and (4.2).Passing from Cartesian to spherical coordinates, one can easily find that and therefore Proceeding as for problem (P1), we have estimated both E(T ) and E rel (T ) for T ∈ {1, 10, 20, 30, 40, 50, 100}.
The data that we have obtained are reported in Table 8 and show that the estimate of E(T ) is monotonically increasing, in agreement with the definition (4.1).From Table 8 it also emerges that the trend of the estimate of E rel (T ) has a sort of plateau between T = 30 and T = 40, after which there is a slight rise.In this regard, the same considerations made for Table 7 apply.However, for every T ≤ 100 the order of magnitude of E rel (T ) is not greater than 10 −5 .Therefore, we may affirm that our predictive method is very accurate and efficient in this case, on both a short and long time scale.

Fifth test problem
Let α > 0 and ω = 2α (2α + 1).The last problem that we consider is where This problem is nothing but the 3-dimensional version of (P2) and its exact solution is given by u(x, y, z, t) ≡ u α (x, y, z, t) := t (x 2 + y 2 + z 2 ) α .
In order to assess the accuracy of our method in solving (P5) and the distance between the approximate solution û and the exact solution u, we again limited ourselves to estimating the quantities defined in (4.1) and (4.2).Switching from Cartesian to spherical coordinates, we can easily obtain This yields Proceeding as for problem (P2), we have estimated E(T ) and E rel (T ) for α ∈ {0.

Conclusions
In this paper, we have explored the ability of PINNs to accurately predict the solutions of some strongly degenerate parabolic problems arising in gas filtration through porous media.
Since there are no general methods for finding analytical solutions to such problems, it is essential to use efficient and accurate numerical methods.blueOne of the most prevalent methods for addressing these problems is the Finite Difference Method (FDM), wherein PDEs are discretized into a system of algebraic equations to be solved numerically.However, the FDM necessitates the discretization of the domain into a grid of cells or nodes, which can become computationally expensive for large and intricate systems.Although the primary objective of this article is not to prove the effectiveness of a PINN compared to a classical numerical method, we engaged in a comparison with the FDM.As established in the literature, for problems characterized by a less complex domain, the FDM typically exhibits a higher level of accuracy compared to PINNs.Nevertheless, in our study, the advantage of using a PINN lies in the ability to test the model on various presented problems (varying the initial/boundary functions and the α parameter), once it has been trained.Additionally, the FDM can be utilized as a benchmark in cases where the solution to the problem is unknown, ensuring a fair comparison under equivalent accuracy conditions.
For the test problems discussed here, whose exact solutions are fortunately known, we have compared the plots of the predicted solutions with those of the analytical solutions.Moreover, to evaluate the accuracy of our predictive method in a purely quantitative way, we have also analyzed the error trends over time.The proposed approach provides accurate results in line with expectations, at least in short-term predictions.However, some issues remain open, such as how to obtain fully reliable plots for the predicted solution when the exact (unknown) one is not differentiable somewhere, and how to reduce or eliminate some slight discrepancies between the contour lines of the predicted solution and those of the analytical solution in the case n = 2.
To the best of our knowledge, this is one of the first papers demonstrating the effectiveness of the PINN framework for solving strongly degenerate parabolic problems with asymptotic structure of Laplacian type.

Figure 3 . 1 :
Figure 3.1: Overall structure of the proposed methodology.An FF-DNN serves as the neural network's architecture.Automatic differentiation is employed to calculate the loss terms associated with the model's dynamics.The loss function comprises two components: the physics loss, represented by L F , and the boundary loss denoted by L B .During the optimization phase, the objective is to minimize the loss function with respect to the set of hyperparameters θ.
without deformation, since ∂ t u ≡ 1 over Ω T (see Fig. 4.1, above).To train the neural network, in each experiment we have initially used 441 points to suitably discretize the domain Ω and its boundary ∂Ω, and 21 equispaced points in the time interval [0, T ].Once the network has been trained, we have made a prediction of the solution to problem (P1) at different times t (Fig. 4.1, below).
Final time T Estimate of E(T ) Estimate of E rel (T )
Final time T Estimate of E(T ) Estimate of E rel (T )

Table 10 :
Estimates of E(T ) and E rel