Quintuple and n-Dimensional Laplace Transforms for High-Order Partial Differential Equations: A Unified Analytical and Data-Driven Framework

The higher-order Laplace transform, including quintuple and general n-dimensional variants, extends the classical Laplace transform to functions of multiple independent variables. These transforms are particularly useful for converting complex PDEs into algebraic equations, facilitating analytical and numerical solutions in higher dimensions.

2.1 General definition and properties of the n-dimensional Laplace transform

Let f be a continuous function of n variables, then the Laplace transform of $f\left(t_1, t_2, t_3, \ldots, t_n\right)$ is defined by

$\begin{gathered}L\left\{f\left(t_1, t_2, t_3, \ldots, t_n\right)\right\}=F\left(s_1, s_2, s_3, \ldots, s_n\right)=\int_0^{\infty} \ldots \ldots \ldots \int_0^{\infty} e^{-\sum_{i=1}^n s_i t_i} f\left(t_1, t_2, t_3, \ldots, t_n\right)d t_1 d t_2 d t_3 \ldots d t_n\end{gathered}$          (1)

where, $t_i \geq 0, \quad(i=1,2,3 \ldots n)$ represent independent variables corresponding to time, space, or other physical parameters, and $s_i$ denote the associated Laplace variables and the inverse Laplace transform

$\begin{gathered}f\left(t_1, t_2, t_3, \ldots, t_n\right)=\frac{1}{2 \pi i} \int_{\alpha_1-i \infty}^{\alpha_1+i \infty} e^{s_1 t_1} \ldots \int_{\alpha_n-i \infty}^{\alpha_n+i \infty} e^{s_n t_n}F\left(s_1, s_2, s_3, \ldots, s_n\right) d s_n \ldots d s_3 d s_2 d s_1\end{gathered}$         (2)

which can be written as

$\begin{gathered}f\left(t_1, t_2, t_3, \ldots, t_n\right)=\left(\frac{1}{2 \pi i}\right)^n \int_{\alpha_1-i \infty}^{\alpha_1+i \infty} \ldots \int_{\alpha_n-i \infty}^{\alpha_n+i \infty} e^{\sum_{i=1}^n s_i t_i}F\left(s_1, s_2, s_3, \ldots, s_n\right) d s_1 \ldots d s_n\end{gathered}$         (3)

Each α_i must be chosen such that it lies to the right of all singularities of F with respect to s_i.

The n-th Laplace transforms of certain partial derivatives of a function with n variables are provided below.

1). The first-order partial derivative of the function of n variables can be written as

$\begin{gathered}L\left[\frac{\partial}{\partial t_1} f\left(t_1, t_2, t_3, \ldots, t_n\right)\right]=s_1 F\left(s_1, s_2, s_3, \ldots, s_n\right)-L_{n-1}\left\{f\left(0, t_2, t_3, t_4, \ldots, t_n\right)\right\}\end{gathered}$          (4)

2). The second-order partial derivative of the function of n variables is given by

$\begin{gathered}L\left[\frac{\partial^2}{\partial t_i^2} f\left(t_1, t_2, t_3, \ldots, t_n\right)\right]=s_i^2 F\left(s_1, s_2, s_3, \ldots, s_n\right)-s_i L_{n-1}\left\{f\left(\ldots, t_i=0, \ldots\right)\right\}-L_{n-1}\left\{\left|\frac{\partial f}{\partial t_i}\right|_{t_i=0}\right\}\end{gathered}$         (5)

3). The second-order mixed partial derivative of the function of n variables is defined: Let $f\left(t_1, t_2, t_3, \ldots, t_n\right)$ be a continuous function. For distinct indices $i \neq j$, the mixed second-order partial derivative is $\frac{\partial^2 f}{\partial t_i \partial t_j}$.

The Laplace transform of this mixed derivative is given by

$\begin{gathered}L\left\{\frac{\partial^2 f}{\partial t_i \partial t_j}\right\}=s_i s_j F\left(s_1, s_2, s_3, \ldots, s_n\right)-s_i L_{n-1}\left\{\left|\frac{\partial f}{\partial t_j}\right|_{t_i=0}\right\}-s_i L_{n-1}\left\{\left|\frac{\partial f}{\partial t_i}\right|_{t_j=0}\right\}-L_{n-2}\left|\frac{\partial^2 f}{\partial t_i \partial t_j}\right|_{t_i=0, t_j=0}\end{gathered}$        (6)

4). The n-th order partial derivative of a function of n variables is given by

$\begin{array}{r}L\left[\frac{\partial^n}{\partial t_i^n} f\left(t_1, t_2, t_3, \ldots, t_n\right)\right]=\left(\prod_{i=1}^n s_i\right) F\left(s_1, s_2, s_3, \ldots, s_n\right)-\sum_{\substack{s \subseteq\{1,2,3, \ldots n\} \\ s \neq \emptyset}}\left(\prod_{i \in s} s_i L_{n-|s|}\left\{\left|\frac{\partial^{|s|} f}{\prod_{i \in s} \partial t_i}\right|_{t_i=0 \text { for } i \in s}\right\}\right)\end{array}$          (7)

2.2 General definitions and properties for quintuple (Special case of n-th order)

For $n=5$, the quintuple Laplace transform of $f\left(t_1, t_2, t_3, t_4, t_5\right)$ is given by

$\begin{gathered}L_5\left\{f\left(t_1, t_2, t_3, t_4, t_5\right)\right\}=F\left(s_1, s_2, s_3, s_4, s_5\right) \\ =\iiint \iint_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+s_3 t_3+s_4 t_4+s_5 t_5\right)} \\ f\left(t_1, t_2, t_3, t_4, t_5\right) d t_1 d t_2 d t_3 d t_4 d t_5 \\ =\iiint \iint_0^{\infty} e^{-\Sigma_{i=1}^5 s_i t_i} \\ f\left(t_1, t_2, t_3, t_4, t_5\right) d t_1 d t_2 d t_3 d t_4 d t_5\end{gathered}$          (8)

where, each $t_i \geq 0$, represents independent dimensions of time, space, or other physical parameters and $S_i$ are the Laplace variables and the inverse quintuple Laplace transform

$\begin{gathered}f\left(t_1, t_2, t_3, t_4, t_5\right) \\ =\frac{1}{2 \pi i} \int_{\alpha-i \infty}^{\alpha+i \infty} e^{s_1 t_1} \frac{1}{2 \pi i} \int_{\beta-i \infty}^{\beta+i \infty} e^{s_2 t_2} \\ \frac{1}{2 \pi i} \int_{\gamma-i \infty}^{\gamma+i \infty} e^{s_3 t_3} \frac{1}{2 \pi i} \int_{\delta-i \infty}^{\delta+i \infty} e^{s_4 t_4} \\ \frac{1}{2 \pi i} \int_{\mu-i \infty}^{\mu+i \infty} e^{s_5 t_5} F\left(s_1, s_2, s_3, s_4, s_5\right) d s_5 d s_4 d s_3 d s_2 d s_1\end{gathered}$        (9)

The five-dimensional Laplace transforms corresponding to specific partial derivatives of a function of five variables are presented

1). The first-order partial derivative of a function of five variables is given by

$\begin{gathered}L_5\left[\frac{\partial}{\partial t_1} f\left(t_1, t_2, t_3, t_4, t_5\right)\right]=s_1 F\left(s_1, s_2, s_3, s_4, s_5\right)-L_4\left\{f\left(0, t_2, t_3, t_4, t_5\right)\right\}\end{gathered}$          (10)

2). The second-order partial derivative of a function of five variables is given by

$\begin{gathered}L_5\left[\frac{\partial^2}{\partial t_i^2} f\left(t_1, t_2, t_3, t_4, t_5\right)\right]=s_i^2 F\left(s_1, s_2, s_3, s_4, s_5\right)-s_i L_4\left\{f\left(\ldots, t_i=0, \ldots\right)\right\}-L_4\left[\left|\frac{\partial f}{\partial t_i}\right|_{t_i=0}\left(\ldots, t_i=0, \ldots\right)\right]\end{gathered}$         (11)

3). The fifth-order partial derivative of a function of five variables is given by

$\begin{gathered}L_5\left[\frac{\partial^5}{\partial t_i^5} f\left(t_1, t_2, t_3, t_4, t_5\right)\right]=s_1 s_2 s_3 s_4 s_5 F\left(s_1, s_2, s_3, s_4, s_5\right)-\sum_{\substack{s \subseteq\{1,2,3,4,5\} \\ s \neq \emptyset}}\left(\prod_{i \in s} s_i L_{5-|s|}\left\{\begin{array}{c}\left.\frac{\partial^{|s|} f}{\prod_{i \in s} \partial t_i} \right\rvert\, t_i=0 \\ \text { for } i \in s\end{array}\right\}\right)\end{gathered}$         (12)

The convolution theorem is a powerful tool in the analysis of linear systems and PDEs, extending naturally to multidimensional Laplace transforms, including quintuple (5D) and general n-dimensional cases.

4.1 Theorem: (n-th Laplace transform)

If at the point $\left(s_1, s_2, s_3, \ldots, s_n\right)$, the integrals

$F_1\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\cdots+s_n t_n\right)} f_1\left(t_1, t_2, \ldots, t_n\right) d t_1 d t_2 \ldots d t_n$        (21)

$F_2\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\cdots+s_n t_n\right)} f_2\left(t_1, t_2, \ldots, t_n\right) d t_1 d t_2 \ldots d t_n$          (22)

$F_3\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\cdots+s_n t_n\right)} f_3\left(t_1, t_2, \ldots, t_n\right) d t_1 d t_2 \ldots d t_n$          (23)

$F_{m-1}\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\cdots+s_n t_n\right)} f_{m-1}\left(t_1, t_2, \ldots, t_n\right) d t_1 d t_2 \ldots d t_n$         (24)

are convergent, and in addition, if

$F_m\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\cdots+s_n t_n\right)} f_m\left(t_1, t_2, \ldots, t_n\right) d t_1 d t_2 \ldots d t_n$          (25)

is absolutely convergent, then, the following expression

$\begin{aligned} & F\left(s_1, s_2, \ldots, s_n\right)=F_1\left(s_1, s_2, \ldots, s_n\right) \\ & F_2\left(s_1, s_2, \ldots, s_n\right) \ldots F_m\left(s_1, s_2, \ldots, s_n\right)\end{aligned}$         (26)

is the Laplace transform of the function

$f\left(t_1, t_2, \ldots, t_n\right)=\left(f_1 f_2 \ldots f_m\right)\left(t_1, t_2, \ldots, t_n\right)$         (27)

where, n denotes the dimension of Laplace transform, and m represents the number of functions involved in the convolution and the integral

$\begin{gathered}F\left(s_1, s_2, \ldots, s_n\right)=\int_0^{\infty} \ldots \int_0^{\infty} e^{-\left(s_1 t_1+s_2 t_2+\ldots+s_n t_n\right)} f\left(t_1, t_2, \ldots, t_n\right)d t_1 d t_2 \ldots d t_n\end{gathered}$         (28)

is convergent at $\left(s_1, s_2, \ldots, s_n\right)$, see the detailed proofs in studies [1-4].

4.2 Example: (Quintuple Laplace transform)

Setting $n=5$, let $t=\left(t_1, t_2, t_3, t_4, t_5\right) \quad$ and $\quad s= \left(s_1, s_2, s_3, s_4, s_5\right)$. For function $f_1, f_2, \ldots, f_5$ with convergent Laplace transforms

$F(s)=F_1(s) F_2(s) F_3(s) F_4(s) F_5(s)$         (29)

$f\left(t_1, t_2, t_3, t_4, t_5\right)=\left(f_1 f_2 f_3 f_4 f_5\right)\left(t_1, t_2, t_3, t_4, t_5\right)$         (30)

Explicitly, the convolution integral is

$\begin{aligned} & f\left(t_1, t_2, t_3, t_4, t_5\right) \\ & =\int_0^{t_1} \int_0^{t_2} \int_0^{t_3} \int_0^{t_4} \int_0^{t_5} f_1\left(t_1\right. \\ & -\left(t_{11}+t_{12}+t_{13}+\rho\right), t_2 \\ & -\left(t_{21}+t_{22}+t_{23}+\sigma\right), t_3 \\ & -\left(t_{31}+t_{32}+t_{33}+\mu\right), t_4 \\ & -\left(t_{41}+t_{42}+t_{43}+\phi\right), t_5 \\ & \left.-\left(t_{51}+t_{52}+t_{53}+\psi\right)\right) f_2\left(t_{11}\right. \\ & -\left(t_{12}+t_{13}+\rho\right), t_{21}-\left(t_{22}+t_{23}+\sigma\right), t_{31} \\ & -\left(t_{32}+t_{33}+\mu\right), t_{41}-\left(t_{42}+t_{43}+\phi\right), t_{51} \\ & \left.-\left(t_{52}+t_{53}+\psi\right)\right) f_3\left(t_{12}-\left(t_{13}+\rho\right), t_{22}\right. \\ & -\left(t_{23}+\sigma\right), t_{32}-\left(t_{33}+\mu\right), t_{42}-\left(t_{43}+\phi\right), t_{52} \\ & \left.-\left(t_{53}+\psi\right)\right) f_4\left(t_{13}-\rho, t_{23}-\sigma, t_{33}-\mu, t_{43}\right. \\ & \left.-\phi, t_{53}-\psi\right) f_1(\rho, \sigma, \mu, \phi, \psi) d \rho d \sigma d \mu d \phi d \psi\end{aligned}$          (31)

4.3 Theorem

A function $f\left(t_1, t_2, \ldots, t_n\right)$ that is continuous on $[0, \infty)$ and fulfills the growth condition specified in condition (13) can be uniquely reconstructed solely from its Laplace transform $F\left(s_1, s_2, \ldots, s_n\right)$ through the following expression

$\begin{gathered}f\left(t_1, t_2, \ldots, t_n\right)=\lim _{m_1, m_2, \ldots, m_{n \rightarrow \infty}} \frac{(-1)^{m_1+m_2+\cdots+m_n}}{m_{1}!m_{2}!\ldots m_{n}!}\left(\frac{m_1}{t_1}\right)^{m_1+1}\left(\frac{m_2}{t_2}\right)^{m_2+1} \ldots\left(\frac{m_n}{t_n}\right)^{m_n+1}\frac{\partial^{m_1+m_2+\cdots+m_n}}{\partial m_1 \partial m_2 \ldots \partial m_n}\left[\frac{m_1}{t_1}, \frac{m_2}{t_2}, \ldots, \frac{m_n}{t_n}\right]\end{gathered}$          (32)

Evidently, the main challenge in applying in expression (32) to compute the inverse Laplace transform lies in the repeated symbolic differentiation of $F\left(s_1, s_2, \ldots, s_n\right)$. To examine the applicability and efficiency of the theorem, we now consider the following example.

Let the following functions

$f\left(t_1, t_2, \ldots, t_n\right)=e^{-s_1 t_1-s_2 t_2-\cdots-s_n t_n}$         (33)

Then the Laplace transform

$F\left(s_1, s_2, \ldots, s_n\right)=\frac{1}{\left(s_1+a\right)\left(s_2+b\right) \ldots\left(s_n+c\right)}$         (34)

Now, by using the higher-order mixed derivative of the previous expression (34), we derive the result

$\begin{aligned} & \frac{\partial^{m_1+m_2+\cdots+m_n}}{\partial m_1 \partial m_2 \ldots \partial m_n} F\left(s_1, s_2, \ldots, s_n\right)= & \frac{m_{1}!m_{2}!\ldots m_{n}!(-1)^{m_1+m_2+\cdots+m_n}}{\left(s_1+a\right)^{m_1+1}\left(s_2+b\right)^{m_2+1} \ldots\left(s_n+c\right)^{m_n+1}}\end{aligned}$         (35)

Applying theorem (32), the expression (35) becomes

$\begin{gathered}f\left(t_1, t_2, \ldots, t_n\right)=\lim _{m_1, m_2, \ldots, m_{n \rightarrow \infty}} \frac{(-1)^{m_1+m_2+\cdots+m_n}}{m_{1}!m_{2}!\ldots m_{n}!}\left(\frac{m_1}{t_1}\right)^{m_1+1}\left(\frac{m_2}{t_2}\right)^{m_2+1} \ldots\left(\frac{m_n}{t_n}\right)^{m_n+1}\left(a+\frac{m_1}{t_1}\right)^{-m_1-1}\left(b+\frac{m_2}{t_2}\right)^{-m_2-1} \ldots\left(c+\frac{m_n}{t_n}\right)^{-m_n-1}\end{gathered}$        (36)

$\begin{gathered}f\left(t_1, t_2, \ldots, t_n\right)=\lim _{m_1, m_2, \ldots, m_{n \rightarrow \infty}}\left(1+\frac{a t_1}{m_1}\right)^{-m_1-1}\left(1+\frac{b t_2}{m_2}\right)^{-m_2-1} \ldots\left(1+\frac{c t_n}{m_n}\right)^{-m_n-1}\end{gathered}$         (37)

Applying logarithmic properties and L’Hôpital’s Rule to the preceding expression (37) yields the following result

$\begin{gathered}\ln \left(f\left(t_1, t_2, \ldots, t_n\right)\right)=-a t_1-b t_2-\cdots-c t_n\Rightarrow f\left(t_1, t_2, \ldots, t_n\right)=e^{-a t_1-b t_2-\cdots-c t_n}\end{gathered}$         (38)

In this section, we demonstrate the applicability of the quintuple and higher-order Laplace transform operators to the solution of higher-dimensional PDEs. These examples illustrate how the transforms reduce complex partial differential operators into algebraic forms in the transform domain, thereby simplifying the process of obtaining exact or semi-analytical solutions. MATLAB is then employed to visualize representative solutions in lower-dimensional slices.

6.1 Application: Fifth-order partial differential equations

Consider the following fifth-order linear PDEs

$\begin{array}{r}\frac{\partial^5}{\partial x_1 \partial x_2 \partial x_3 \partial x_4 \partial x_5} u\left(x_1, x_2, x_3, x_4, x_5\right)-u\left(x_1, x_2, x_3, x_4, x_5\right)=0\end{array}$        (56)

subjected to the following conditions:

$\begin{aligned} & u\left(0, x_2, x_3, x_4, x_5\right)=e^{x_2+x_3+x_4+x_5} \\ & u\left(x_1, 0, x_3, x_4, x_5\right)=e^{x_1+x_3+x_4+x_5} \\ & u\left(x_1, x_2, 0, x_4, x_5\right)=e^{x_1+x_2+x_4+x_5} \\ & u\left(x_1, x_2, x_3, 0, x_5\right)=e^{x_1+x_2+x_3+x_5} \\ & u\left(x_1, x_2, x_3, x_4, 0\right)=e^{x_1+x_2+x_3+x_4}\end{aligned}$       (57)

Applying the quintuple Laplace transform on both sides of Eq. (56), and using the derivative property (Eq. (12)), the transformed PDE becomes

$\begin{gathered}s_1 s_2 s_3 s_4 s_5 U\left(s_1, s_2, s_3, s_4, s_5\right)-G\left(s_1, s_2, s_3, s_4, s_5\right)-U\left(s_1, s_2, s_3, s_4, s_5\right)=0\end{gathered}$      (58)

$\begin{gathered}\left(s_1 s_2 s_3 s_4 s_5-1\right) U\left(s_1, s_2, s_3, s_4, s_5\right) =G\left(s_1, s_2, s_3, s_4, s_5\right)\end{gathered}$       (59)

By applying the Laplace transformation to the boundary conditions (57) and applying the exponential shifting property (Eq. (52)) gives

$\begin{aligned} & G\left(s_1, s_2, s_3, s_4, s_5\right)=\frac{1}{\left(s_2-1\right)\left(s_3-1\right)\left(s_4-1\right)\left(s_5-1\right)}+ \\ & \frac{1}{\left(s_1-1\right)\left(s_3-1\right)\left(s_4-1\right)\left(s_5-1\right)}+\frac{1}{\left(s_1-1\right)\left(s_2-1\right)\left(s_4-1\right)\left(s_5-1\right)}+ \\ & \frac{1}{\left(s_1-1\right)\left(s_2-1\right)\left(s_3-1\right)\left(s_5-1\right)}+\frac{1}{\left(s_1-1\right)\left(s_2-1\right)\left(s_3-1\right)\left(s_4-1\right)}\end{aligned}$       (60)

where the five terms correspond to the five boundary conditions.

Then the Laplace domain solution is 

$U\left(s_1, s_2, s_3, s_4, s_5\right)=\frac{G\left(s_1, s_2, s_3, s_4, s_5\right)}{\left(s_1 s_2 s_3 s_4 s_5-1\right)}$       (61)

Now taking the inverse quintuple Laplace transform on both sides of Eq. (61) yields

$L^{-1}\left[U\left(s_1, s_2, s_3, s_4, s_5\right)\right]=L^{-1}\left[\frac{G\left(s_1, s_2, s_3, s_4, s_5\right)}{\left(s_1 s_2 s_3 s_4 s_5-1\right)}\right]$      (62)

$\begin{aligned} u\left(x_1, x_2, x_3, x_4, x_5\right) & =L^{-1}\left\lceil\frac{G\left(s_1, s_2, s_3, s_4, s_5\right)}{\left(s_1 s_2 s_3 s_4 s_5-1\right)}\right\rceil =e^{x_1+x_2+x_3+x_4+x_5}\end{aligned}$       (63)

Thus, the explicit solution satisfies the governing PDE (56) together with the boundary conditions (57).

Table 6 presents the approximated numerical values of the function $u\left(x_1, x_2, 0.5,0.5,0.5\right) \times 10^5$ on grid $[0,2] \times [0,2](\Delta x=0.2)$, fixed $x_3=x_4=x_5=0.5$. The values are computed directly from the exact analytical solution $u=e^{x_1+x_2+x_3+x_4+x_5}=e^{x_1+x_2+1.5}$ and then scaled by $10^5$ and rounded to the nearest integer. The numerical results demonstrate a consistent monotonic increase in both coordinate directions. For any fixed $x_2$, the values increase steadily as $x_1$ increases from 0 to 2 . Similarly, for any fixed $x_1$, the values increase as $x_2$ increases. This behavior reflects the exponential dependence of the solution on the sum $x_1+ x_2$. The smallest value occurs at $\left(x_1, x_2\right)=(0,0)$, while the largest value appears at $(2,2)$, confirming that the function grows continuously across the domain without any local maxima or symmetric peak. The surface implied by the table is therefore smooth and strictly increasing toward the upperright corner of the domain.

Table 6. Approximated values of $u\left(x_1, x_2, 0.5,0.5,0.5\right) \times 10^5$ on grid $[0,2] \times[0,2](\Delta x=0.2)$, fixed $x_3=x_4=x_5=0.5$

Overall, Table 6 is fully consistent with the analytical expression and correctly captures the exponential growth pattern of the solution, ensuring agreement between theoretical formulation and computed data.

Figure 5 illustrates the three-dimensional surface of $u\left(x_1, x_2, 0.5,0.5,0.5\right) \times 10^5$ over the domain $[0,2] \times [0,2](\Delta x=0.2)$, fixed $x_3=x_4=x_5=0.5$. The surface is generated from the analytical solution $u=e^{x_1+x_2+1.5}$. It is evident that the surface increases smoothly and monotonically along both $x_1$ and $x_2$ directions, with the minimum value at (0,0) and the maximum at (2,2). The absence of any central peak or symmetry confirms that the solution exhibits pure exponential growth. The color gradient, transitioning from blue (lower values) to red (higher values), further highlights the continuous increase of the function across the domain. This visualization is fully consistent with the analytical formulation and supports the correctness of the computed numerical results.

5.png

Figure 5. 3D surface plot of the function $u\left(x_1, x_2, x_3, x_4, x_5\right)$ with fixed values $x_3=x_4=x_5=0.5$

6.1.1 Deep Operator Network fits Laplace inversion

The classical numerical inversion of a multi-dimensional Laplace transforms $F\left(s_1, \ldots, s_5\right)$, F relies on iterated Bromwich integrals. This becomes computationally expensive in high dimensions due to nested complex contour integrations, sensitivity to branch cuts, and the need for many evaluations. DeepONet learns the operator G mapping input functions (here, F evaluated on a grid of s-points or parametrized in s-domain) to output functions (f evaluated at t-points). Trained on pairs (F_i, f_i) from known examples (e.g., generated via forward Laplace on test functions or simulations), it generalizes to unseen F, enabling fast, data-driven inversion without repeated contour integrals.

Branch net encodes the input function F (e.g., samples of F at fixed "sensor" points in s-space); trunk net encodes query locations $\left(t_1, \ldots, t_5\right)$; output is their (possibly nonlinear) combination. 

It leverages the universal operator approximation theorem, making it suitable for learning transforms like Laplace (explicitly mentioned in early DeepONet works as a demonstrable case, i.e., learning Laplace on Gaussian processes or bounded supports).

6.1.2 Relevant variants for Laplace-related problems

The most directly relevant extension is the Laplace Neural Operator (LNO) [9], explicitly designed to improve on FNO for time-dependent problems by working in the Laplace domain:

• LNO decomposes kernels via pole residue representations in the s-domain: $K(s)=\sum \beta_n /\left(s-\mu_n\right)$, with trainable poles μₙ and residues βₙ.
• It naturally captures transient responses, exponential decays, and non-periodic signals, exactly the behavior of inverse Laplace transforms (damped oscillations, decaying exponentials, etc.).
• Forward: Input → lifting → Laplace-layer (pole-residue multiplication + convolution) → projection → output.
• For inversion: Train LNO to map F(s)→f(t) by treating the inverse transform as the learned operator (physics-informed or supervised on known pairs).
• Advantages over pure DeepONet/FNO: Better handling of transients (your Gaussian-like ridge in the s-slice suggests multi-exponential time behavior); exponential convergence in some cases; avoids periodicity assumptions of Fourier methods.

6.1.3 Comparisons from recent studies (physics-informed neural operators for control/ODEs)

• DeepONet: Strong at enforcing initial conditions, good generalization, but slower training in high dimensions.
• FNO: Fast on grids, excels in periodic cases, but is weaker on transients.
• LNO: Competitive on transient-rich problems (e.g., damped systems), but can be unstable or require tuning (more poles/layers) for highly nonlinear cases.

6.1.4 Practical application to 5D quintuple Laplace setting

Analytical uniqueness: DeepONet/LNO provides practical numerical uniqueness and fast evaluation:

1. Data generation to use known f (sums of multi-exponentials, products separable in tᵢ) to compute F via numerical quadrature of the 5D integral (3.5), or symbolic if possible.
2. Training supervised pairs (F samples on s-grid, f on t-grid). Or physics-informed (PI-DeepONet): Enforce Bromwich-like residuals or exponential-order bounds.
3. Inversion workflow input new $F\left(s_1, \ldots, s_5\right) \rightarrow$ trained DeepONet/LNO → query f at desired $\left(t_1, \ldots, t_5\right)$.
4. For slice example: $\left(F\left(s_1, s_2 \mid s_3=s_4=s_5=1\right)\right.$ Gaussian ridge) A 2D+ trunk variant could learn the inverse on that slice quickly; extend to full 5D with separable tricks.
5. Benefits in high dimensions avoid the curse of dimensionality in classical nested Bromwich (5 nested contours → huge cost). DeepONet generalizes across parametric.

DeepONet (especially the LNO variant) is one of the strongest modern machine learning approaches for numerical multi-dimensional Laplace inversion, and it turns a uniqueness theorem into a fast, generalizable computational tool for recovering f from F in practice, particularly suited to transient behaviors visible in various perspectives in Figure 6. 

6.png

Figure 6. Two-layer feed-forward deep neural network model (Deep Operator Networks (DeepONets))

Figure 6 shows that DeepONet is a powerful neural architecture for learning nonlinear operators, particularly effective for solving parametric PDEs like the Laplace equation in high dimensions. In practical high-dimensional applications, such as a 5D quintuple Laplace setting (Laplace's equation in 5 spatial dimensions, possibly with quintuple referring to five-fold parametric or coordinate complexity), DeepONet approximates the solution operator mapping input functions to field solutions. It excels at breaking the curse of dimensionality compared to traditional solvers or standard PINNs, enabling efficient inference for multi-query scenarios (varying boundary conditions or parameters in 5D domains). The architecture uses a branch net (encoding inputs) and trunk net (query locations), combined via dot products, often with physics-informed training for data scarcity. This approach supports fast surrogate modeling in complex physics, quantum, or continuum systems involving high-D Laplace/Poisson problems.

DeepONet is a deep learning framework designed to approximate nonlinear operators, i.e., mappings from functions to functions. Table 7 summarizes typical PDE examples where DeepONet excels:

1. Linear parametric ODEs serve as classical benchmarks to validate operator learning.
2. Advection-diffusion and Burgers’ equations demonstrate DeepONet’s ability to capture nonlinear and shock-dominated dynamics efficiently.
3. Elliptic PDEs such as Poisson or Helmholtz show DeepONet’s capability for solving spatially complex problems, especially on irregular meshes.
4. Navier-Stokes equations illustrate surrogate modeling for fluid flows, where combining DeepONet with CNN branches enhances spatial feature extraction.
5. Diffusion-reaction and Eikonal equations highlight DeepONet’s capacity to handle parametric inputs, coefficient variations, and discontinuities, making it useful in chemical kinetics, geophysics, and optics.

DeepONet’s strength lies in learning mappings across function spaces, making it ideal for parameterized PDEs, multi-dimensional domains, and physics-informed modeling.

Table 7. Representative PDE problems demonstrating Deep Operator Networks (DeepONets) performance

Table 8. Deep Operator Network (DeepONet) vs. Fourier Neural Operator (FNO) for partial differential equations (PDEs) learning

Table 8 summarizes key differences between DeepONet and FNO for PDE surrogate modeling. DeepONet excels at arbitrary queries, irregular domains, and small data or physics-informed scenarios, while FNO often trains faster due to Fast Fourier Transform (FFT) based layers and performs well on periodic or smooth problems on regular grids. The choice depends on problem type, domain geometry, and data availability.

6.2 Application: General n-th order partial differential equation

Consider the following general n-dimensional PDEs

$\begin{gathered}\frac{\partial^n}{\partial x_1 \partial x_2 \partial x_3 \ldots \partial x_n} u\left(x_1, x_2, x_3, \ldots, x_n\right)-u\left(x_1, x_2, x_3, \ldots, x_n\right)=0\end{gathered}$       (64)

where, $u=u\left(x_1, x_2, x_3, \ldots, x_n\right)$, and each variable $x_i \in[0, \infty]$ for $i=1,2,3 \ldots, n$.

With exponential boundary conditions

$\begin{gathered}u\left(x_1, x_2, \ldots, x_{i-1}, 0, x_{i+1}, \ldots, x_n\right)=e^{\sum_{j \neq i} x_j}, \\ \text { for all } i=1,2,3 \ldots, n\end{gathered}$      (65)

Applying the n-dimensional Laplace transform and the n-th-order derivative property (Eq. (7)) on both sides of Eq. (64), the transformed PDE becomes

$\begin{gathered}\left(\prod_{i=1}^n s_i\right) U\left(s_1, s_2, \ldots, s_n\right)-G\left(s_1, s_2, \ldots, s_n\right)-U\left(s_1, s_2, \ldots, s_n\right)=0\end{gathered}$      (66)

This implies

$\left(\left(\prod_{i=1}^n s_i\right)-1\right) U\left(s_1, s_2, \ldots, s_n\right)=G\left(s_1, s_2, \ldots, s_n\right)$       (67)

By taking the Laplace transform of the boundary conditions (65) gives

$L\left[e^{\sum_{j \neq i} x_j}\right]=\prod_{j \neq i} \frac{1}{s_j-1}$       (68)

So, summing over all n such expressions

$G\left(s_1, s_2, \ldots, s_n\right)=\sum_{i=1}^n \prod_{\substack{j=1 \\ j \neq i}} \frac{1}{s_j-1}$      (69)

Then the Laplace domain solution is

$U\left(s_1, s_2, \ldots, s_n\right)=\frac{\sum_{i=1}^n \prod_{\substack{j=1 \\ j \neq i}} \frac{1}{s_j-1}}{\left(\prod_{i=1}^n s_i\right)-1}$       (70)

Now, taking the inverse n-th Laplace transform on both sides of Eq. (70)

$L^{-1}\left[U\left(s_1, s_2, \ldots, s_n\right)\right]=L^{-1}\left[\frac{\sum_{i=1}^n \prod_{\substack{j=1 \\ j \neq i}} \frac{1}{s_j-1}}{\left(\prod_{i=1}^n s_i\right)-1}\right]$      (71)

$\begin{gathered}u\left(x_1, x_2, x_3, \ldots, x_n\right)=L^{-1}\left[\frac{G\left(s_1, s_2, \ldots, s_n\right)}{\left(\prod_{i=1}^n s_i\right)-1}\right]=e^{x_1+x_2+\cdots+x_n}\end{gathered}$       (72)

Hence the result.

6.3 Application: n-dimensional heat equation

Consider the heat equation in n-dimensional is given by

$\frac{\partial u\left(x_1, x_2, \ldots, x_n, t\right)}{\partial t}=\sum_{i=1}^n \frac{\partial^2 u}{\partial x_i{ }^2}, \quad t>0$        (73)

With the boundary and initial conditions specified below

$\begin{gathered}u\left(x_1, \ldots, x_{i-1}, 0, x_{i+1}, \ldots, x_n, t\right)=e^{\sum_{j \neq i} x_j+n t} \\ \text { for all } i=1,2, \ldots, n \\ u\left(x_1, x_2, \ldots, x_n, 0\right)=e^{x_1+x_2+\cdots+x_n}\end{gathered}$       (74)

Now, by applying the (n + 1)-fold Laplace transform and using the second-order derivative property (Eq. (5)) for spatial terms and the first-order derivative property (Eq. (4)) for time, the result is

$\begin{aligned} & s U\left(p_1, \ldots, p_n, s\right)-F\left(p_1, \ldots, p_n\right)=\sum_{i=1}^n\left(p_i{ }^2 U\left(p_1, \ldots, p_n, s\right)-B_i\right)\end{aligned}$        (75)

Simplifying Eq. (75), we get

$\left(s-\sum_{i=1}^n p_i{ }^2\right) U=F-\sum_{i=1}^n B_i$      (76)

where, $B_i$ denote the boundary terms and let, $G\left(p_1, \ldots, p_n, s\right)=F-\sum_{i=1}^n B_i$.

Then the Eq. (76) becomes

$U\left(p_1, \ldots, p_n, s\right)=\frac{G\left(p_1, \ldots, p_n, s\right)}{s-\sum_{i=1}^n p_i^2}$       (77)

From initial and boundary conditions (74)

$F\left(p_1, \ldots, p_n\right)=\frac{1}{\left(p_1-1\right)\left(p_2-1\right) \ldots\left(p_n-1\right)}$       (78)

$B_i=\frac{1}{(s-n)} \prod_{j \neq i} \frac{1}{\left(p_j-1\right)}$      (79)

$G\left(p_1, \ldots, p_n, s\right)=\frac{1}{\prod_{i=1}^n\left(p_i-1\right)} \frac{s-\sum_{i=1}^n p_i}{s-n}$      (80)

Then Eq. (77) becomes

$\begin{gathered}U\left(p_1, \ldots, p_n, s\right)=\frac{1}{\prod_{i=1}^n\left(p_i-1\right)} \frac{s-\sum_{i=1}^n p_i}{(s-n)\left(s-\sum_{i=1}^n p_i^2\right)}\end{gathered}$       (81)

Taking the inverse Laplace transform in all (n + 1) variables on both sides of Eq. (81), we obtain 

$u\left(x_1, \ldots, x_n, t\right)=e^{x_1+x_2+\cdots+x_n+n t}=e^{\sum_{i=1}^n e^{x_i+n t}}$      (82)

To demonstrate the potential of data-driven approximation, a DeepONet model can be trained using solution data generated from Eq. (82) to learn the operator mapping between spatial–temporal coordinates and the temperature field.

Table 9. Exact analytical values (×10⁶) of the n-dimensional heat equation solution $u\left(x_1, x_3, t\right)$ on a 21 × 21 grid $\left(x_1, x_3 \in[0,2], \Delta x=0.1\right)$, with all other coordinates fixed at 0.5

Table 9 reports the exact analytical values of $u\left(x_1, x_3, t\right)$ obtained from the closed-form solution $u\left(x_1, \ldots, x_n, t\right)=e^{x_1+x_2+\cdots+x_n+n t}$. A two-dimensional slice is considered by varying $x_1$ and $x_3$ on a 21 × 21 grid over [0,2] with step $\Delta x=0.1$, while the remaining spatial variables are fixed at 0.5. The values are presented as (u × 10⁶) for scaling convenience. The results show a monotonic increase along both spatial directions, which is consistent with the exponential analytical solution derived using the higher-order Laplace transform method.

Figure 7 illustrates the surface visualization of the analytical solution $u\left(x_1, x_3, t\right)$ obtained from the n-dimensional heat equation for a two-dimensional spatial slice where the remaining spatial coordinates are fixed. The surfaces are plotted for five time instances: t = 0.00,   0.25,   0.50,   0.75. The results show that the magnitude of the solution increases progressively with time, which reflects the exponential structure of the analytical expression derived using the higher-order Laplace transform method. The color scale indicates the magnitude of $u\left(x_1, x_3, t\right)$, where warmer colors correspond to larger values. The smooth surface and consistent growth pattern confirm the analytical behavior predicted by the exact solution.

7.png

Figure 7. Time evolution of the analytical solution $u\left(x_1, x_3, t\right)$ for a two-dimensional slice of the n-dimensional heat equation at different time levels

6.3.1 Machine learning approximation of the n-dimensional heat equation solution using deep neural networks or Deep Operator Network

In addition to the analytical solution obtained through the higher-order Laplace transform method, modern machine learning techniques can be employed to approximate the solution of PDEs. Among these approaches, DNNs and operator learning architectures such as DeepONet have shown significant potential in learning complex mappings between input variables and solution fields.

For the n-dimensional heat equation, the solution can be interpreted as a mapping between spatial–temporal variables and the corresponding temperature distribution. In particular, the variables $\left(x_1, \ldots, x_n, t\right)$ serve as inputs, while the function $u\left(x_1, \ldots, x_n, t\right)$ represents the output solution. A machine learning model can therefore be trained to approximate this relationship using data generated from the analytical solution derived in Eq. (82). DeepONet is particularly suitable for such problems because it learns operators rather than simple numerical functions. The architecture typically consists of two subnetworks: a branch network, which encodes input functions or parameters, and a trunk network, which represents the spatial–temporal coordinates at which the solution is evaluated. The outputs of these subnetworks are combined through an inner product to obtain the final approximation of the solution. In the present framework, training data can be generated by sampling points within the computational domain $x_i \in[0,2]$ and $t \in[0,1]$. The corresponding solution values are computed using the analytical expression obtained from the Laplace transform formulation. These input–output pairs are then used to train the neural network so that it learns the underlying relationship governing the heat diffusion process.

After training, the machine learning model can rapidly predict the temperature distribution for new spatial–temporal inputs without repeatedly solving the differential equation. This significantly reduces computational cost when evaluating solutions over large multidimensional domains or when performing repeated simulations. Consequently, combining the analytical Laplace transform framework with machine learning models such as DNNs or DeepONet provides a hybrid methodology for solving high-dimensional PDEs. The transform-based analysis ensures mathematical accuracy, while the data-driven model enables efficient approximation and fast evaluation of the solution.

8.png

Figure 8. Two-layer feed-forward deep neural network model (deep neural networks (DNN) or Deep Operator Network (DeepONet))

Figure 8 represents a two-layer feed-forward DNN architecture used for learning the relationship between input variables and the predicted output. First, the input layer receives several input variables (in this figure, 4 inputs). These inputs represent the features or data points that are provided to the neural network. Next, the inputs are passed to the hidden layer, where each input is multiplied by a weight (W) and combined with a bias (b). After this linear transformation, an activation function (Sigmoid) is applied. The hidden layer contains multiple neurons (about 10 neurons in the diagram) that help the model capture nonlinear patterns in the data.

Finally, the processed information from the hidden layer is passed to the output layer, where another weighted combination is performed to generate the final output (1 output value).

Overall, this feed-forward network processes data sequentially from input → hidden layer → output, enabling the model to approximate complex mathematical relationships and make predictions.

9.png

Figure 9. Neural network approximation of the analytical solution of the n-dimensional heat equation

Figure 9 illustrates the performance of a neural network model trained to approximate the analytical solution of the n-dimensional heat equation derived using the higher-order Laplace transform method. The horizontal axis represents the exact solution values obtained from the analytical expression in Eq. (82), while the vertical axis shows the corresponding predictions produced by the trained neural network model. Each point in the scatter plot corresponds to a sampled spatial-temporal input $\left(x_1, \ldots, x_n, t\right)$ within the computational domain. The clustering of points close to the diagonal line indicates strong agreement between the analytical solution and the neural network predictions. This demonstrates that the machine learning model successfully captures the underlying functional relationship governing the heat diffusion process.

The results confirm that neural network–based surrogate models can effectively approximate the analytical solution of high-dimensional PDEs. Such data-driven models provide a computationally efficient alternative for evaluating PDE solutions over large domains or for performing repeated simulations.

6.4 Application: n-dimensional diffusion equation

Finally, consider the diffusion equation, which has the same mathematical form as the heat equation but is interpreted in terms of concentration rather than temperature.

$\begin{gathered}\frac{\partial u\left(x_1, x_2, \ldots, x_n, t\right)}{\partial t}=\sum_{i=1}^n \frac{\partial^2 u}{\partial x_i{ }^2}, \\ \text { on }[0, \infty]^n, t>0\end{gathered}$       (83)

With the boundary and initial conditions specified

$\begin{gathered}u\left(x_1, x_2, \ldots, x_n, 0\right)=\prod_{i=1}^n \sin \left(m_i \pi x_i\right) \\ u\left(x_1, \ldots, x_i=0, \ldots, x_n, t\right)=0\end{gathered}$      (84)

Now, by applying the (n + 1)-fold Laplace transform and using the second-order derivative property (Eq. (5)) for spatial terms and the first-order derivative property (Eq. (4)) for time, the result is 

$\begin{aligned} & s U\left(p_1, \ldots, p_n, s\right)-G\left(p_1, \ldots, p_n\right)=-\sum_{i=1}^n p_i{ }^2 U\left(p_1, \ldots, p_n, s\right)\end{aligned}$      (85)

Simplifying the previous equation, we get,

$U\left(p_1, \ldots, p_n, s\right)=\frac{G\left(p_1, \ldots, p_n\right)}{\left(s+\sum_{i=1}^n p_i^2\right)}$       (86)

Applying the Laplace transform to the initial condition (84) leads to

$L\left[u\left(x_1, x_2, \ldots, x_n, 0\right)\right]=\prod_{i=1}^n L\left[\sin \left(m_i \pi x_i\right)\right]$       (87)

$G\left(p_1, \ldots, p_n\right)=\prod_{i=1}^n \frac{m_i \pi}{p_i{ }^2+\left(m_i \pi\right)^2}$       (88)

Substituting this value into Eq. (86) gives the result as 

$U\left(p_1, \ldots, p_n, s\right)=\frac{\prod_{i=1}^n \frac{m_i \pi}{p_i{ }^2+\left(m_i \pi\right)^2}}{\left(s+\sum_{i=1}^n p_i{ }^2\right)}$      (89)

Applying the inverse Laplace transform with the linearity property to the previous equation yields

$u\left(x_1, \ldots, x_n, t\right)=\prod_{i=1}^n \sin \left(m_i \pi x_i\right) e^{-\pi^2\left(\sum_{i=1}^n m_i^2\right) t}$      (90)

Table 10 presents the computed values of the diffusion solution $u\left(x_1, x_2, t ; x_3=0.5\right) \times 10$ at selected spatial points for five different time levels t = 0, 0.25,0.50, 0.75, and 1.00. The analytical solution is based on the product of sine functions multiplied by an exponential decay term $e^{-3 \pi^2 t}$. Because the sine function vanishes at the boundaries of the domain $[0,1]$, the solution values remain zero at points where $x_1=0$ or $x_2=0$ or $x_1=1$ or $x_2=1$.

At the initial time t = 0, the solution reaches its maximum value at the interior point (0.5,0.5), corresponding to the peak of the sine distribution. As time increases, the exponential factor rapidly reduces the magnitude of the solution at all spatial points. Consequently, the values decrease steadily with time while preserving the spatial structure determined by the sine terms. This behavior reflects the characteristic property of diffusion processes, where an initial concentration gradually dissipates over time.

Table 10. Temporal evolution of the diffusion solution $u\left(x_1, x_2, t ; x_3=0.5\right) \times 10$ at selected spatial points for t = 0, 0.25, 0.50, 0.75, and 1.0, based on the analytical solution $u=\prod_{i=1}^3 \sin \left(\pi x_i\right) e^{-3 \pi^2 t}$ on the domain [0,1]

Figure 10 illustrates the time evolution of the diffusion solution $u\left(x_1, x_2, t ; x_3=0.5\right) \times 10$ using three-dimensional surface plots for several time levels. At t = 0, the solution forms a smooth sine-shaped surface with a maximum at the center of the domain (0.5,0.5) and zero values along the boundaries. This represents the initial distribution of the diffusing quantity.

As time progresses, the amplitude of the surface decreases significantly due to the exponential decay term in the analytical solution. The peak height gradually diminishes while the overall shape remains centered in the domain. By larger time values, the surface becomes nearly flat, indicating that the magnitude of the solution approaches zero throughout the region. The color variation in the surface plots highlights this rapid reduction, demonstrating the diffusion-driven decay of the initial pulse.

10.png

6.4.1 Machine learning application for the n-dimensional diffusion model

Beyond solving the n-dimensional diffusion equation analytically using higher-order Laplace transforms, modern machine learning techniques offer a complementary approach for approximating its solution. Data generated from the analytical solution (Eq. (90)) can serve as a training set, where the input consists of the spatial temporal coordinates $\left(x_1, \ldots, x_n, t\right)$ and the output is the corresponding solution $u\left(x_1, \ldots, x_n, t\right)$.

Among various options, DNNs, support vector machines, and operator-learning architectures like DeepONet are particularly effective for capturing complex functional dependencies inherent in PDEs. The training data is typically obtained by sampling points from the domain and computing their exact solution values using the Laplace transform-based formula.

Once trained, these models can quickly predict the diffusion behavior at new points without the need to repeatedly solve the differential equation, offering a substantial reduction in computational cost for large multidimensional domains or multiple simulation runs. This hybrid methodology, combining analytical and data-driven techniques, provides an efficient and scalable framework for modeling high-dimensional diffusion phenomena and other intricate PDE systems.

Figure 11 illustrates the architecture of a DNN used to approximate a mapping between input variables and the predicted output. The network begins with an input layer containing three variables, which represent the features or parameters provided to the model. These inputs are passed to the first hidden layer (Hidden Layer 1), which consists of 50 neurons. In this layer, each neuron performs a weighted linear combination of the inputs using weight parameters (W) and bias terms (b). The result is then passed through a nonlinear activation function, allowing the network to capture complex relationships in the data. The transformed outputs are forwarded to Hidden Layer 2, which also contains 50 neurons. Similar to the first hidden layer, each neuron computes a weighted sum of the incoming signals and adds a bias term, followed by a nonlinear activation function. This second layer enables the model to learn deeper and more abstract representations of the input features. Finally, the processed information reaches the output layer, which contains a single neuron. This neuron performs another weighted combination of the previous layer’s outputs and produces the final predicted value. The output layer typically uses a linear activation when the network is applied to regression problems, such as approximating solutions of PDEs.

11.png

Figure 11. Three-layer feed-forward deep neural network model (deep neural network (DNN) or Deep Operator Network (DeepONet))

Overall, the architecture represents a fully connected feedforward neural network, where information flows sequentially from the input layer through two hidden layers to generate the final prediction. Such networks are commonly used for function approximation and scientific computing tasks, including the learning of PDE solution mappings.

Figure 12 shows a 3D scatter comparison between the analytical solution and the prediction produced by a DNN for a two-dimensional diffusion equation evaluated at a fixed time within a square domain. The horizontal axes represent the spatial variables $x_1$ and $x_2$, while the vertical axis represents the solution value uuu, which may correspond to temperature, concentration, or another diffusive quantity.

12.png

Figure 12. Deep neural network (DNN) approximation of the 2D diffusion equation compared with the analytical solution

Blue points indicate the exact analytical solution, whereas red points represent the DNN-predicted values at the same spatial locations. Both sets of points form a smooth, hill-shaped surface with the highest values near the center of the domain and gradually decreasing toward the boundaries, reflecting the typical behavior of diffusion from a localized source. Visually, the predicted values closely follow the analytical solution, particularly near the central region. Slight differences appear near the boundaries and at lower solution values, which is common when neural networks approximate PDE solutions. Overall, the figure demonstrates that the DNN can effectively reproduce the main features of the diffusion process and provide a reliable approximation of the analytical solution.

This study generalizes the classical Laplace transform to address quintuple (five-dimensional) and general n-dimensional PDEs. A formal definition of the n-th order Laplace transform and its operational properties for higher-order and mixed partial derivatives are developed, allowing complex multidimensional PDEs to be converted into algebraic forms in the transform domain. This transformation simplifies analytical solution procedures and enables reconstruction of the original functions. Theoretical results establish existence and uniqueness conditions, confirming that the transforms are well-defined under exponential growth assumptions. In addition, the generalized Bromwich inversion integral guarantees accurate and unique recovery of solutions in the original domain.

To illustrate the practical implementation of the method, several computational examples and visualizations were generated using MATLAB. One representative example considered a scaled Laplace slice F(s₁, s₂), which exhibited a Gaussian-like peak of approximately 5.8 near the point (1,7). Surface plots revealed symmetric patterns and resonance-type structures in the transform domain, highlighting how multidimensional variable interactions influence the transformed solution space. Additional tabulated results for selected bivariate functions and their partial derivatives further validated the consistency of the transformation and inversion processes.

The results show that the generalized Laplace transform framework is an effective analytical approach for simplifying and solving high-dimensional PDEs. Unlike traditional numerical methods such as the finite difference method and finite element method, the transform approach provides clearer insight into the structural behavior of solutions in multidimensional systems. However, computational difficulty increases as the dimension n exceeds five, particularly during the numerical inversion stage. Overall, higher-order Laplace transforms form a strong analytical basis for studying complex systems, and future work will aim to improve inversion techniques, investigate fractional extensions, and integrate machine learning to broaden practical applications.