Enhanced Numerical Solutions for Fractional PDEs Using Monte Carlo PINNs Coupled with Cuckoo Search Optimization

1. Introduction

Fractal geometry provides important understanding into processes that involve memory effects over extended periods and significant spatial interactions across large scales. This concept is supported by the strength of fractional derivatives, which offer a more precise method for modeling these kinds of systems. As a result, fractional partial differential equations (fPDEs) have become widely used to describe phenomena that exhibit memory effects, non-local interactions in space, and power-law behavior [1,2,3,4]. However, the non-commutative nature of fractional integrals and derivatives poses significant challenges, particularly in deriving analytical solutions, even for linear problems. To tackle these issues effectively, various computational Strategies have been created for solving fPDEs [5,6,7,8]. Notably, spectral methods in time and space have been introduced, alongside techniques for fractional inverse problems, achieving high accuracy in applications such as the Fokker-Planck equation [9]. Additionally, to improve uniform error bounds for long-term dynamics in mildly nonlinear spatially space-fractional Klein-Gordon equations, Jia and Jiang [10] recently proposed the exponential wave integrator method, demonstrating its efficacy for such complex systems.

High-dimensional time-space fractional PDEs (fPDEs) present unique challenges, primarily due to the complex nature of Lévy flights. Addressing these challenges often requires specialized numerical methods or adaptive boundary treatments using appropriate mesh designs. Techniques initially developed for lower-dimensional, less complex problems have been extended to higher dimensions [11,12,13,14,15]. For example, in [16], the authors proposed a Galerkin finite element method (FEM) that uses an unstructured mesh and the Crank-Nicolson (CN) time integration methodology. This approach, applied specifically to 2D time-space fractional PDEs in convex domains, demonstrated superior accuracy compared to traditional structured-mesh finite-element methods. Further advancements include differential quadrature methods utilizing radial basis functions (RBFs) for analyzing three-dimensional diffusive fractional equations [17]. By employing RBFs as trial functions and optimizing their linear combination ratios, these methods achieved promising accuracy and efficient approximation of fractional derivatives. Convergence analysis validated their success and highlighted their efficiency over finite-difference methods in reducing global error. It is essential to emphasize the flexibility of fractional derivatives in representing long-memory physical interactions. However, their non-commutative nature often complicates analytical integration, even in linear cases. While various fast algorithms have been developed for small-scale problems to enhance computational efficiency [18,19,20,21,22], solving high-dimensional fPDEs continues to demand more advanced and sophisticated numerical approaches.

A novel approach to approximating fractional-order differentiation was recently introduced using the Monte Carlo strategy [23,24]. This method begins with the Grünwald-Letnikov definition of fractional derivatives and incorporates the development of a probabilistic density function that adheres to the rules of fractional differentiability. The Monte Carlo method is then utilized to identify sample points within specific fractional orders, offering a new algorithm for derivative approximation. The approach generates nodes that resemble earlier finite-difference constructions, designed for non-uniform or nested grids. These nodes are densely concentrated around the target point while remaining well-distributed at the initial positions, making the method highly suitable for solving fPDEs effectively and supporting parallel processing. However, this method does not integrate seamlessly with numerical approaches for fPDEs, as it requires the structure of the objective function to be explicitly evaluated at the generated sample nodes. Physics-Informed Neural Networks (PINNs) [25] have gained popularity as a powerful method for addressing partial differential equations (PDEs) by combining principles of physics with deep learning techniques [26,27,28,29,30]. One of the main successes of PINNs is the capacity to act as a direct solution approximator using neural networks. This is done by constructing a loss function that discourages the growth of predictions of the solution against the labeled data and at the same time satisfying the governing equations. This double aim also means that the estimated solution respects as many natural constraints as possible of the physical model.

By incorporating the innovative Monte Carlo method, Physics-Informed Neural Networks (PINNs) become an effective tool for solving fractional partial differential equations (PDEs). Previous studies, such as reference [31], explored the combination of Monte Carlo methods with PINNs, introducing the concept of fPINNs to address both direct and inverse fractional PDE challenges. However, their method was restricted to applying Monte Carlo techniques specifically for computing integrals involving fractional functions with Caputo derivatives. This approach did not extend to other types of fractional derivatives, such as Riemann-Liouville derivatives. Also, their approach presupposed the employment of fixed parameter distributions for obtaining highly accurate estimation of the loss function. Though this strategy was effective in minimizing computational overhead, it was not very effective when it came to involving specific information concerning the system. To avoid these limitations, using an enhanced and more flexible Monte Carlo method where the workload is not predetermined prior to the method could highly enhance the efficiency and versatility of this combined approach.

In the present study, we present a new method by implementing the Physics-Informed Neural Networks (PINNs) that are based on the Monte Carlo method through integrating it with the Cuckoo Search (CS) algorithm referred as the PINN-CS to solve fractional partial differential equations (fPDEs). The approach starts with a modified Monte Carlo procedure for approximating fractional derivatives on the right side in terms of two different quasi-Monte Carlo methods. These techniques are accompanied by Physics-Informed Neural Networks (PINNs) and are enhanced using the Cuckoo Search algorithm to provide solutions for various difficult fractional partial differential equations (fPDEs). Compared with the classic fPINN approach, the PINN-CS model employs block-sampled nodes and averages their function values. This method helps to minimize the memory consumption and computational costs of the simulation as well as to maximize the computational speed and throughput due to the proper arrangement of the sample nodes. Furthermore, in the node distribution, particular attention was paid to analogize finite difference methods on scattered or non-uniform grids, which will allow for better perception of the specifics of boundary problems for their closer study at the distances.

The performance of the method suggested in this work is tested on examples that are illustrated numerically; the fuzzy boundary problem is also solved, although the definitions of boundary conditions are not complete. Comparisons with the fPINN method show that, although fPINN yields slightly better numerical accuracy, the proposed PINN-CS approach achieves considerably increased efficiency. To provide proof of the paper’s applicability, the paper presents an actual-case of blood circulation in the human brain to depict how machine systems can learn and develop governing equations of structural physical systems. The presented methodology can be adopted as new grounds to solve fPDEs, especially for high dimensional and complex domains.

1.1. Problem Statement

The general structure of fractional partial differential equations (fPDEs) can be expressed as:

$$\begin{array}{c}{}_0{}^{C}D_t^{\alpha }u(x,t)+\mathcal{L}\left(u,u_{xx},(-\mathsf{\Delta }{)}^{{\displaystyle \frac{\beta }{2}}}u,\dots ;x,t\right)=f(x,t),x\in \mathsf{\Omega },t\in [0,T],\\ u\left(x,0\right)=u_0\left(x\right),x\in \mathsf{\Omega },u\left(x,t\right)=u\left(x,t\right),x\in \partial \mathsf{\Omega },t\in \left(0,T\right],\end{array}$$

(1)

In these equations $u(x,t)$ indicates the analytical solution, where $f(x,t)$ represent the forcing term, and ${}_0{}^{C}D_t^{\alpha }$ signifies the Caputo time-fractional derivative operator, which is defined as follows:

$${}_0{}^{C}D_t^{\alpha }u(x,t):={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\int }_0^t (t-s{)}^{-\alpha }{u}^{(n)}(x,s)ds,n-1<\alpha <n,$$

(2)

and $(-\mathsf{\Delta }{)}^{{\displaystyle \frac{\beta }{2}}}$ is the Riesz fractional derivative operator given by:

$$(-\mathsf{\Delta }{)}^{{\displaystyle \frac{\beta }{2}}}u(x,t):=c_{\beta }{( }^{RL}{D}_{x_{\mathrm{l}\mathrm{b}}}^{\beta }u{+ }^{RL}{D}_{x_{\mathrm{u}\mathrm{b}}}^{\beta }u),$$

(3)

with $m-1<\beta <m,c_{\beta }={\displaystyle \frac{-1}{2\mathrm{c}\mathrm{o}\mathrm{s}(\pi \beta /2)}}$, and:

$$\begin{array}{l}{}_{RL}D_{x_{\mathrm{l}\mathrm{b}}}^{\beta }u(x,t):={\displaystyle \frac{1}{\mathsf{\Gamma }(m-\beta )}}{\displaystyle \frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}}{\int }_{x_{\mathrm{l}\mathrm{b}}}^x (x-\tau {)}^{m-\beta -1}u(\tau ,t)d\tau ,\\ {}_{RL}D_{x_{\mathrm{u}\mathrm{b}}}^{\beta }u(x,t):={\displaystyle \frac{(-1{)}^m}{\mathsf{\Gamma }(m-\beta )}}{\displaystyle \frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}}{\int }_x^{x_{\mathrm{u}\mathrm{b}}} (\tau -x{)}^{m-\beta -1}u(\tau ,t)d\tau ,\end{array}$$

(4)

where $x_{\mathrm{l}\mathrm{b}}$ and $x_{\mathrm{u}\mathrm{b}}$ represent the lower and upper bounds of $\mathsf{\Omega }$ with respect to $x$, defined as:

$$x_{\mathrm{l}\mathrm{b}}:=\mathrm{i}\mathrm{n}\mathrm{f}\left\{x:(x,t)\in \mathsf{\Omega }\subseteq {\mathbb{R}}^d\right\},x_{\mathrm{u}\mathrm{b}}:=\mathrm{s}\mathrm{u}\mathrm{p}\left\{x:(x,t)\in \mathsf{\Omega }\subseteq {\mathbb{R}}^d\right\}$$

(5)

In high-dimensional spaces, ${}_{RL}D_{x_{\mathrm{l}\mathrm{b}}}^{\beta }u$ and ${}_{RL}D_{x_{\mathrm{u}\mathrm{b}}}^{\beta }u$ generalize to expressions such as ${}^{RL}D_{x_{\mathrm{l}\mathrm{b}}}^{\beta }u,{}^{RL}D_{x_{\mathrm{u}\mathrm{b}}}^{\beta }u,{}^{RL}D_{y_{\mathrm{b}}}^{\beta }u,{}^{RL}D_{y_{\mathrm{u}\mathrm{b}}}^{\beta }u$, and so on, allowing for a broad class of operators in high dimensional systems. Here, $\mathcal{L}(\cdot )$ is a spatial differential operator either linear or non-linear, that can include various different orders of derivatives. For clarity, this work focuses on the operator $(-\mathsf{\Delta }{)}^{\beta /2}$ as an example. By applying the proposed Monte Carlo PINNs method, the research addresses fPDEs across both straightforward and intricate domains

2. The Generalized Mante Carlo PINNs Method

In this section, we introduce a new form of the Monte Carlo integrated Physics-Informed Neural Network (PINN) to solve equation (1). Now, as a preliminary step for presenting the new methodology, we expand upon the base Monte Carlo approach to right sided fractional derivatives.

2.1. Mante Carlo Method for Computing Right Sided Fractional Differentiation

In prior work [23,24], the author proposed a Monte Carlo-based technique for estimates of the fractional differentiation, which has been extended to other approximations. This approach exploits the Grunwald Letnikov operators for fractional differential to obtain probability distributions with fractional indices. These distributions are then used together with the Monte Carlo technique in order to produce sample points and enable the fractional differentiation of order of various functions to be approximated.

For instance, Fractional-order derivative defined on the left side of $f(x)$ can be expressed using the formula:

$$\begin{array}{l}{x}_{lb}{D}_x^{\alpha }f(x)\sim A_h^{\alpha }f(x)\\ ={\displaystyle \frac{1}{K}}\left({\displaystyle \frac{1}{h^{\alpha -n}}}\left({\displaystyle \frac{1}{h^n}}\sum _{j=0}^n (-1{)}^j{C}_n^{j}f(x-jh)-{\displaystyle \frac{1}{N}}\sum _{m=1}^N {\displaystyle \frac{1}{h^n}}\sum _{j=0}^n (-1{)}^j{C}_n^{j}f\left(x-\left(Y_m^{(\alpha )}+j\right)h\right)\right)\right)\end{array},$$

(6)

where $x_{Lb}$ represents the domain's left boundary with respect to$x$ and the order of the fractional derivative is $\alpha$ , where $\alpha \in (n,n+1)$. The step size $h$ is determined as $h=(x-$ $\left.x_{Lb}\right)/N$, with $N$ being the total number of points used for the approximation. The parameter $K$ denotes the number of repetitions for sampling. The nodes $\left\{y_m^{(\alpha )},m=1,\dots ,N\right\}$ represent the sampled points in fractional order. This method is capable of approximating different form of fractional derivatives, including Caputo, Grunwald Letnikov, and Riemann Liouville derivatives.

The mathematical foundation for the left-sided fractional derivative was previously explored in the original study [24], where the corresponding formula was provided. To obtain an operator analogous to the Riesz fractional derivative, it is essential to derive the right-sided fractional derivatives. The following derivation illustrates the process for $\alpha \in (\mathrm{1,2})$, though it is adaptable for other values of $\alpha$. The Grunwald Letnikov fractional derivative on $L_1(\mathbb{R})$ is defined as [8]:

$${}_x{}^{GL}D_{x_{ub}}^{\alpha }f(x)=\underset{h⟶0}{\mathrm{l}\mathrm{i}\mathrm{m}} {\mathcal{A}}_{-h}^{\alpha }f(x),\alpha >0,h>0,$$

(7)

where,

$${\mathcal{A}}_{-h}^{\alpha }f(x)={\displaystyle \frac{1}{h^{\alpha }}}{\sum }_{k=0}^{\mathrm{\infty }} (-1{)}^k{\displaystyle \frac{\alpha \left(\alpha -1\right)\cdots \left(\alpha -k+1\right)}{k!}}f\left(x+kh\right),$$

(8)

here, $x_{ub}$ represents the domain's right boundary with respect to x, and $h$ represent the step size, determined as $h=\left(x_{ub}-x\right)/N$, where $N$ represents the total number of approximation points. To proceed, we define:

$$m_k=v(\alpha ,k)=(-1{)}^k{\displaystyle \frac{\mathsf{\Gamma }(\alpha +1)}{k!\mathsf{\Gamma }(\alpha -k+1)}}$$

(9)

From this, we have:

$$m_0=1,m_1=v(\alpha ,1)=-\alpha <0, \mathrm{and} m_k>0, \mathrm{for} k=\mathrm{2,3},\dots$$

(10)

Moreover, the coefficients of variation $m_k$ can be interpreted as the expansion coefficients for the sequence $(1-z{)}^{\alpha }$, satisfying ${\sum }_{k=0}^{\mathrm{\infty }} m_k=0$. Thus, ${\sum }_{k=1}^{\mathrm{\infty }} (-v(\alpha ,k))=1$.

Since:

$$\begin{array}{c}\sum _{k=0}^{\mathrm{\infty }} m_k=0=1-\alpha +\sum _{k=2}^{\mathrm{\infty }} m_k\\ \end{array},$$

(11)

it follows that:

$$1=2-\alpha +\sum _{k=2}^{\mathrm{\infty }} m_k.$$

(12)

Consider $Y$ as an independent random variable satisfying the following conditions:

$$\begin{array}{c}P(Y=1)=p_1=p_1(\alpha )=2-\alpha , \mathrm{where} \alpha \in (\mathrm{0,1}),\hfill \\ P(Y=k)=p_k=p_k(\alpha )=m_k={\displaystyle \frac{(-1{)}^k\mathsf{\Gamma }(\alpha +1)}{k!\mathsf{\Gamma }(\alpha -k+1)}}\mathrm{for} k=\mathrm{2,3},\dots , \mathrm{where} m_k>0.\hfill \end{array}$$

(13)

Thus

$${\sum }_{k=1}^{\mathrm{\infty }} p_k=1$$

(14)

It is noted that:

$E(Y)<\mathrm{\infty }$ (the expected value is finite),

$\mathrm{V}\mathrm{a}\mathrm{r}(Y)=\mathrm{\infty }$ (the variance is infinite).

Next, the summation for $m_k$ is applied to a function $f$:

$${\sum }_{k=0}^{\mathrm{\infty }} m_{k}f(x+kh)=f(x)-\alpha f(x+h)+{\sum }_{k=2}^{\mathrm{\infty }} m_{k}f(x+kh),$$

(15)

which simplifies further to:

$$f(x)-\left(2-p_1(\alpha )\right)f(x+h)+{\sum }_{k=2}^{\mathrm{\infty }} m_{k}f(x+kh),$$

(16)

and finally:

$$f(x)-2f(x+h)+E[f(x+Yh)],$$

(17)

where $E[f(x+Yh)]$ is the expected value over $Y$.

If stochastic process is defined as $\xi (x)=f(x+Yh)$, then for fixed $x,t,h$, the expectation satisfies:

$$E[f(x+Yh)]<\mathrm{\infty }.$$

(18)

Now, let $Y_1,Y_2,\dots ,Y_n$ are to be independent copies of the random variable Y. By the strong law of large numbers:

$${\displaystyle \frac{1}{N}}{\sum }_{n=1}^N f\left(x+Y_{n}h\right)\to E[f(x+Yh)] \mathrm{as} N\to \mathrm{\infty },$$

(19)

with probability 1 for any fixed $t$ and $h$.

Thus, the approximation for the expectation is:

$$A_N^{\alpha ,h}={\displaystyle \frac{1}{h\alpha }}\left[f(x)-2f(x+h)+{\displaystyle \frac{1}{N}}{\sum }_{n=1}^N f\left(x+Y_{n}h\right)\right].$$

(20)

The probability of convergence to $A_h^{\alpha }f(x)$, where $\alpha \in (\mathrm{1,2})$, is also noted. Furthermore, as in the referenced source, we have:

$$A_{N,h}^{\alpha }={\displaystyle \frac{1}{K}}\left({\displaystyle \frac{1}{h\alpha }}\left[f(x)-2f(x+h)+{\displaystyle \frac{1}{N}}{\sum }_{n=1}^N f\left(x+Y_{n}h\right)\right]\right).$$

(21)

Next, sample values $Y_i$ are replaced with Monte Carlo simulations. Define:

$$E_j=\sum _{i=1}^j p_i,$$

(22)

where $p_i=p_i(\alpha )$. The sequence of probabilities is:

$$0=E_0<E_1<E_2<\dots <E_j<\dots , \mathrm{with} p_i=E_i-E_{i-1}.$$

(23)

Let $F$ be a uniformly distributed random variable on the interval $[\mathrm{0,1}]$, then:

$$P\left(E_{j-1}<F\le E_j\right)=p_j,P\left(E_{k-1}<F\le E_k\right)=p_k.$$

(24)

To generate $Y\in \{\mathrm{1,2},\dots \}$, we set $Y=k$ if $E_{k-1}\le F<E_k$.

2.2. Quasi Monte Carlo Approach for Fractional Derivatives

In this section, we illustrate how to use a quasi Monte Carlo approach to approximate fractional derivatives. Unlike the classical Monte Carlo application, which uses random samples, the quasi Monte Carlo approach uses samples with low discrepancy sequences. This method is carried out in two steps:

Step 1: Define ${\tilde{E}}_j={\sum }_{i=1}^j p_i$, where $p_i=p_i(\alpha )$ (as in Equation (13)).

Then:

$$0={\tilde{E}}_0<{\tilde{E}}_1<{\tilde{E}}_2<\cdots ,$$

with:

$$p_i={\tilde{E}}_i-{\tilde{E}}_{i-1}$$

Step 2: Take $\hat{F}$ is a variable chosen at random from the Sobol sequence or the Halton sequence.

Then:

$$P\left({\tilde{E}}_{j-1}<\hat{F}\le {\tilde{E}}_j\right)=p_j.$$

To generate $Y\in \{\mathrm{1,2},\dots \}$, set $Y=k$ if ${\tilde{E}}_{k-1}\le \hat{F}<{\tilde{E}}_k$.

Figure 1 is a schematic representation of the General Monte Carlo Physics-Informed Neural Network (PINN) framework. This visual serve to outline the methodology and workflow of combining Monte Carlo techniques with PINNs to solve fractional partial differential equations (fPDEs). Figure 1 serves as a roadmap for implementing the Monte Carlo PINN approach, making it easier to grasp the workflow and understand the interplay of the components in this hybrid methodology. If you'd like, I can expand further on any specific aspect of this figure.

2.3. The PINN Solver

Expanding upon the earlier approximation method for fractional-order differentiation, this section presents a novel Monte Carlo-based Physics-Informed Neural Network (PINN) for solving fractional Partial Differential Equations (fPDEs), as shown in Figure 1. For example, consider Equation (1) with $\alpha \in (\mathrm{0,1})$ and $\beta =(\mathrm{1,2})$. A complete discrete formulation is developed for the left-hand side of Equation (1) by applying the Monte Carlo approximation techniques outlined in methods (6) and (21) for the left and right sided differentiation, as demonstrated below:

$$\begin{array}{ll}& L^{\mathrm{*}}={\displaystyle \frac{1}{K_t}}\left[{\displaystyle \frac{1}{h_t^{\alpha }}}\left[u(\mathbf{x},t)-{\displaystyle \frac{1}{N_t}}\sum _{m=1}^{N_t} f\left(\mathbf{x},t-Y_m^{(t)}{h}_t\right)\right]\right]\\ & +c_{\beta }\left\{{\displaystyle \frac{1}{{}_{lb}K_{\mathbf{x}}}}\left({\displaystyle \frac{1}{{}_{lb}h_{\mathbf{x}}^{\beta }}}\left[u(\mathbf{x},t)-2u\left(\mathbf{x}{- }^{lb}{h}_{\mathbf{x}},t\right)+{\displaystyle \frac{1}{{}_{lb}N_{\mathbf{x}}}}\sum _{n=1}^{{l }^l{N}_{\mathbf{x}}} u\left(\mathbf{x}-Y_n^{\left({\mathbf{x}}_{lb}\right)lb}{h}_{\mathbf{x}},t\right)\right]\right)\right.\\ & \left.+{\displaystyle \frac{1}{ub}}{K}_{\mathbf{x}}\left({\displaystyle \frac{1}{{}^{bb}h_{\mathbf{x}}^{\beta }}}\left[u(,t)-2u\left(\mathbf{x}{+ }^{ub}{h}_{\mathbf{x}},t\right)+{\displaystyle \frac{1}{ub{N}_{\mathbf{x}}}}\sum _{j=1}^{ub{N}_{\mathbf{x}}} u\left(\mathbf{x}+Y_j^{\left({\mathbf{x}}_{ub}\right)ub}{h}_{\mathbf{x}},t\right)\right]\right)\right\},\end{array}$$

(25)

where $k_t,h_{t}l{b}_{h_x}$, and $u{b}_{h_x}$ are the step sizes with respect to $t$, and:

$$\mathbf{x};h_t=t/N_t{, }^{lb}{h}_{\mathbf{x}}=\left(\mathbf{x}-{\mathbf{x}}_{lb}\right){/ }^{lb}{N}_{\mathbf{x}}, {\mathrm{and} }^{ub}{h}_{\mathbf{x}}=\left({\mathbf{x}}_{ub}-\mathbf{x}\right){/ }^{ub}{N}_{\mathbf{x}};{\left\{Y_m^{(t)}\right\}}_{m=1}^{N_t},{\left\{Y_n^{\left({\mathbf{x}}_{lb}\right)}\right\}}_{n=1}^{l{N}_{\mathbf{x}}}.$$

(26)

Additionally, $N_x$ represents the Monte Carlo method generates samples for ${}^{0}D_x^{\alpha },x_{lb}{D}_x^{\beta }$ and $x_{ub}{D}_x^{\beta }$, respectively. To reduce computational costs, a basic identification is conducted on the data points in the sample set $\left\{Y_i\right\}$, acquired through the Monte Carlo method using a block distribution. The computational speedup achieved by this technique was not immediately apparent with fewer approximation points but became more notable as the number of points increased. The training dataset is defined as:.

$$\mathcal{D}=\left({\mathcal{D}}_{\mathrm{equ} },{\mathcal{D}}_{\mathrm{ini} },{\mathcal{D}}_{\mathrm{bound} }\right),$$

(27)

Here, the network's output acts as an approximate solution.

The corresponding loss function is expressed as:

$$\mathrm{L}\mathrm{o}\mathrm{s}\mathrm{s}(\mathsf{\Theta }=\{ \mathrm{weights}, \mathrm{biases} \})={\omega }_E{\mathrm{M}\mathrm{S}\mathrm{E}}_E+{\omega }_I{\mathrm{M}\mathrm{S}\mathrm{E}}_I+{\omega }_B{\mathrm{M}\mathrm{S}\mathrm{E}}_B,$$

(28)

Where

$$\begin{array}{ll}{\mathrm{M}\mathrm{S}\mathrm{E}}_E& ={\displaystyle \frac{1}{N_E}}\sum _{L=1}^{N_E} {\left[L^{\mathrm{*}}\left(u_{NN}\left(x_i,t_i\right)-f\left(x_i,t_i\right)\right)\right]}^2.\\ {\mathrm{M}\mathrm{S}\mathrm{E}}_B& ={\displaystyle \frac{1}{N_B}}\sum _{L=1}^{N_B} \left[{\left|u_{NN}\left(x_i,t_i\right)-u\left(x_i,t_i\right)\right|}^2\right].\\ {\mathrm{M}\mathrm{S}\mathrm{E}}_I& ={\displaystyle \frac{1}{N_I}}\sum _{L=1}^{N_I} \left[{\left|u_{NN}\left(x_i,0\right)-u_0\left(x_i\right)\right|}^2\right].\end{array}$$

(29)

The penalty weights ${\omega }_E,{\omega }_I$, and ${\omega }_B$ are manually selected here, but a more efficient method is available in the referenced literature. $N_E,N_I$, and $N_B$ denote the number of collocations, initial conditions, and boundary conditions, respectively. The Cuckoo Search technique is used to find the ideal parameter $\mathsf{\Theta }$ that minimizes the loss function, resulting in an estimated solution that meets the target equation.

2.4. Cuckoo Search (CS) Algorithm

The breeding parasitism behavior of cuckoo birds forms the basis of the Cuckoo Search (CS) algorithm. It applies three main hypotheses:

1. Each bird lays one egg at a time and places it in a randomly selected nest.

2. The best nests (those with high-quality solutions) are retained for the next generation.

3. The number of available nests is limited, and the host bird has a probability of discovering a foreign egg, denoted as $S_p\in [\mathrm{0,1}]$.

In mathematical terms, these behaviors are modeled using Lévy flights, a type of global random walk. The CS algorithm creates a new generation for each cuckoo $i$ by performing a Lévy flight, as described by the equation:

$$Z_i^{t+1}=Z_i^t+\alpha \mathrm{L}\mathrm{e}\mathrm{v}\mathrm{y}(s,\lambda ),$$

(30)

where:

$$\mathrm{Levy} \sim u=t^{-\lambda }.$$

(31)

In most situations, the step size $\alpha >0$ is adjusted to 1. The Lévy flight introduces randomness with step sizes drawn from a Lévy distribution.

For local search or neighborhood random exploration, the following equation is used:

$$Z_i^{t+1}=Z_i^t+\alpha s\oplus H\left(S_p-ϵ\right)\oplus \left(Z_j-Z_k\right).$$

(32)

here

$\oplus$ represents element-wise multiplication.

$Z_j$ and $Z_k$ are two randomly selected solutions (nests) obtained using permutation.

$H$ is the Heaviside step function.

$ϵ$ is a small random value drawn from a normal distribution.

$s$ is the step size.

$S_p$ is a switching probability that balances local and global random walks.

3. Numerical Experiments

The proposed method can be applied to solving fractional Partial Differential Equations (fPDEs), as demonstrated by several numerical examples in this section. We showcase how this approach performs when applied to forward problems in 2D fractional partial differential equations (fPDEs). Specifically, we solve the 2D space-fractional Poisson equation on a unit disc domain and the 2D time-space fractional diffusion equation within a heart-shaped domain. The results are then compared with those obtained using the basic fPINN approach.

Additionally, we present numerical results and compute the time required for the selected approximation points. Further investigations include the forward simulation for the 3D coupled time-space fractional Bloch-Torrey equation in the human brain ventricle area. These results are compared with existing numerical studies. Considering the complexity and importance of the chosen problems, the method was successfully validated and proven effective for solving fPDEs with unknown boundary conditions.

Moreover, the numerical results support the applicability of the method to problems like MRIbase studies in the brain. To assess the discrete data accuracy, we apply $L^2$ relative error and point-wise error measurements:

$${ L}^2 \mathrm{relative} \mathrm{error} =\sqrt{{\displaystyle \frac{\sum _i {\left[u\left({\mathbf{x}}_i,t_i\right)-\widehat{u}\left({\mathbf{x}}_i,t_i\right)\right]}^2}{\sum _i {\left[u\left({\mathbf{x}}_i,t_i\right)\right]}^2}}}, \mathrm{Point}\text{-}\mathrm{wise} \mathrm{error} ={\displaystyle \frac{u\left({\mathbf{x}}_i,t_i\right)-\widehat{u}\left({\mathbf{x}}_i,t_i\right)}{u\left({\mathbf{x}}_i,t_i\right)}},$$

(33)

here, $\mathrm{u}$ represents the exact solution. Numerical experiments were carried out on every case within the two domains outlined below:

Domain ${\mathsf{\Omega }}_1$:

$$\mathrm{Unit} \mathrm{disk} {\mathsf{\Omega }}_1=\left\{x,y:{∥x-x_{\mathrm{center} }∥}^2+{∥y-y_{\mathrm{center} }∥}^2\le 1\right\}$$

Domain ${\mathsf{\Omega }}_2$:

$${\mathsf{\Omega }}_2=\left\{\left(x,y\right):\parallel x{\parallel }^2+\parallel y{\parallel }^2\le x^2+y^2\right\},$$

$\begin{array}{c} \mathrm{where} x=1.6{\sin }^{3}t,y=1.3\mathrm{c}\mathrm{o}\mathrm{s}t-5\mathrm{c}\mathrm{o}\mathrm{s}2t-2\mathrm{c}\mathrm{o}\mathrm{s}3t-\mathrm{c}\mathrm{o}\mathrm{s}4t\\ t\in [\mathrm{0,2}\pi ]\end{array}$

Example 3.1

We examine the following two-dimensional space-fractional PDEs, along with their corresponding initial and boundary conditions:

$$(-△{)}_{xx}^{\alpha /2}u(x,y)+(-△{)}_{yy}^{\beta /2}u(x,y)=f(x,y),(x,y)\in \mathsf{\Omega },$$

(34)

with the boundary conditions:

$$u(x,y)=u_b(x,y),(x,y)\in \partial \mathsf{\Omega }.$$

(35)

The fractional operators ${{(-\Delta )}^{{\displaystyle \frac{\alpha }{2}}}}_{xx}$and ${{(-\Delta )}^{{\displaystyle \frac{\beta }{2}}}}_{yy}$are defined as:

$$(-△{)}_{xx}^{\alpha /2}:={\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}(\pi \alpha /2)}}\left(x_{lb}{D}_x^{\alpha }+x{D}_{x_{ub}}^{\alpha }\right),(-△{)}_{yy}^{\beta /2}:={\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}(\pi \beta /2)}}\left(y_{lb}{D}_y^{\beta }+y{D}_{y_{ub}}^{\beta }\right),$$

(36)

where ${}_{x_{lb}}D_x^{\alpha }$, ${}_{X}D_{x_{ub}}^{\alpha }$$,$${}_{y_{lb}}D_y^{\beta }$, and ${}_{X}D_{y_{ub}}^{\beta }$are the right and left sided Riemann-Liouille fractional derivatives defined over $h^{x_{lb},x_{ub}}$and $y_{lb},y_{ub}$, respectively.

The analytical solution of the above problem is:

$${(x-x_{lb})}^3{(x_{ub}-x)}^3{(y-y_{lb})}^3{(y_{ub}-y)}^3,$$

(37)

where $x_{lb},x_{ub},y_{lb}$ _, and $y_{ub}$ are the upper and lower bounds of $x$ and $y$, respectively. The standard structure of the forcing terms is given as follows

$$f(x,y)=f_1(x,y)+f_2(x,y),$$

(38)

$$\begin{array}{c}{f}_1(x,y)={\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi \alpha }{2}\right)}}\left\{{\displaystyle \frac{{\left|x_{lb}+x_{ub}\right|}^3\mathsf{\Gamma }(4)}{\mathsf{\Gamma }(4-\alpha )}}\left({\left(x-x_{lb}\right)}^{3-\alpha }+{\left(x_{ub}-x\right)}^{3-\alpha }\right)\right.\\ -{\displaystyle \frac{3{\left|x_{lb}+x_{ub}\right|}^2\mathsf{\Gamma }(5)}{\mathsf{\Gamma }(5-\alpha )}}\left({\left(x-x_{lb}\right)}^{4-\alpha }+{\left(x_{ub}-x\right)}^{4-\alpha }\right)\\ +{\displaystyle \frac{3\left|x_{lb}+x_{ub}\right|\mathsf{\Gamma }(6)}{\mathsf{\Gamma }(6-\alpha )}}\left({\left(x-x_{lb}\right)}^{5-\alpha }+{\left(x_{ub}-x\right)}^{5-\alpha }\right)\\ \left.-{\displaystyle \frac{3\mathsf{\Gamma }(7)}{\mathsf{\Gamma }(7-\alpha )}}\left({\left(x-x_{lb}\right)}^{6-\alpha }+{\left(x_{ub}-x\right)}^{6-\alpha }\right)\right\}{\left(y-y_{lb}\right)}^3{\left(y_{ub}-y\right)}^3\end{array}$$

(39)

$$\begin{array}{c}{f}_2(x,y)={\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi \beta }{2}\right)}}\left\{{\displaystyle \frac{{\left|y_{lb}+y_{ub}\right|}^3\mathsf{\Gamma }(4)}{\mathsf{\Gamma }(4-\beta )}}\left({\left(y-y_{lb}\right)}^{3-\beta }+{\left(y_{ub}-y\right)}^{3-\beta }\right)\right.\\ -{\displaystyle \frac{3{\left|y_{lb}+y_{ub}\right|}^2\mathsf{\Gamma }(5)}{\mathsf{\Gamma }(5-\beta )}}\left({\left(y-y_{lb}\right)}^{4-\beta }+{\left(y_{ub}-y\right)}^{4-\beta }\right)\\ +{\displaystyle \frac{3\left|y_{lb}+y_{ub}\right|\mathsf{\Gamma }(6)}{\mathsf{\Gamma }(6-\beta )}}\left({\left(y-y_{lb}\right)}^{5-\beta }+{\left(y_{ub}-y\right)}^{5-\beta }\right)\\ \left.-{\displaystyle \frac{3\mathsf{\Gamma }\left(7\right)}{\mathsf{\Gamma }\left(7-\beta \right)}}\left({\left(y-y_{lb}\right)}^{6-\beta }+{\left(y_{ub}-y\right)}^{6-\beta }\right)\right\}{\left(x-x_{lb}\right)}^3{\left(x_{ub}-x\right)}^3.\end{array}$$

(40)

Example 3.1 focuses on solving a two-dimensional space-fractional partial differential equation (PDE) with specified initial and boundary conditions. The example highlights the application and performance of the proposed PINN-CS framework compared to traditional methods.

Table 2 presents a comparative analysis of the performance of three methodologies-fPINN, Monte Carlo fPINN, and the proposed PINN-CS (Monte Carlo-based PINN integrated with the Cuckoo Search algorithm) on a 2D space-fractional partial differential equation (PDE). The comparison focuses on numerical accuracy across different approximation levels, specified by $N$. The table evaluates absolute errors for the numerical solutions obtained by each method. As $N$ (number of approximation points) increases, the errors for all methods decrease, indicating improved solution accuracy with higher resolution. The improvement rate, however, varies significantly across methods. fPINN, shows moderate accuracy but lags behind the other methods. Absolute error is notably higher, especially at lower approximation levels. Monte Carlo fPINN offers better accuracy than fPINN by incorporating Monte Carlo techniques for fractional differentiation. PINN-CS achieves superior accuracy at all levels of $N$. The error is orders of magnitude smaller compared to both fPINN and Monte Carlo fPINN. For instance, at $N=8$, the error is $3.88\times {10}^{-11}$ for PINN-CS, compared to $3.73\times$ ${10}^{-1}$ for fPINN and $2.91\times {10}^{-1}$ for Monte Carlo fPINN. Similar trends are observed at higher $N$ values, highlighting the robustness of PINN-CS.

Figure 2 visually compares the numerical solution obtained using the PINN-CS framework with the exact solution for the 2D space-fractional PDE described in Example 3.1. The overlap between the two solutions is nearly perfect, signifying an extremely high degree of numerical accuracy. The maximum absolute error, as seen in Table 2, is as low as $3.88\times {10}^{-11}$ when $N=8$. This result confirms that PINN-CS provides a highly accurate solution, even with relatively few approximation points. The visual confirmation of overlap across the domain further reinforces the robustness of the methodology for solving complex fractional PDEs with high precision.

A graphical comparison in Figure 3 shows the overlap between the numerical solution generated by PINN-CS and the exact analytical solution. The solutions overlap almost perfectly, visually confirming the high numerical accuracy reported in Table 4. Wherever such deviations exist, as maintained here, they are minimal, and refer to the robust and reliable nature of the underlying PINN-CS framework. This near perfect level of alignment establishes it as particularly capable of accurately approximating solutions to fractional PDEs.

Figure 4 in the paper provides a performance assessment of the PINN-CS framework based on the overlap with the exact solutions to a 3D time-space fractional Bloch-Torrey equation explored in Example 3.3. The figure shows a good correlation with the loss function obtained with numeric solution given by PINN-CS with the exact analytical solution. The almost complete alignment proves the effectiveness of the PINN-CS approach as an accurate modeling technique. The considered analysis example 3.3 is an example of 3D fractional PDE, which make this problem particularly difficult due to the high dimensionality of the problem and the use of fractional operators. As seen from the solution of this problem, the proposed PINN-CS framework appears sufficiently stable and versatile for dealing with all these kinds of equations. The figure supplements numerical results described in the table presented in the article, which shows minimal absolute errors calculated for various approximation points. The reported error, being of the order of 10^-4, establishes the high accuracy of the proposed PINN-CS framework further. In Figure 4, the performance of the PINN-CS framework is demonstrated when solving the 3D time space fractional Bloch Torrey equation considered in Example 3.3. This figure also corroborates well with the fact that the PINN-CS method provides a numerical solution very close to the exact analytical solution. This convergence confirms the ability and credibility of the PINN-CS framework in solving intricate high- dimensional fractional partial differential equations (fPDEs).

Figure 5 plots a quantitative solution of PINN-CS with the exact analytical solution for the 3D Bloch-Torrey equations. The coincidence of the locations of the solutions to both carriers is almost perfect which proves a high accuracy in the formulation of the PINN-CS framework. This figure shows that PINNCS has the ability to approximate the dynamics of high-dimensional fractional PDEs. The minimal deviations corroborate the low absolute errors reported, validating the framework's reliability for applications requiring precise numerical approximations.

Figure 6 shows the point-wise error distribution as a heatmap for the 2D space-fractional PDE solved in Example 3.1. The errors are minimal across the domain, with most values concentrated in the range of ${10}^{-11}$ to ${10}^{-9}$, as corroborated by numerical results in Table 2. While the errors are uniformly low, minor peaks near the boundary regions highlight areas where further refinement may marginally enhance accuracy. The figure underscores the efficacy of the PINN-CS framework, achieving a point-wise error that is orders of magnitude smaller than the competing methods, such as Monte Carlo fPINN $\left(2.91\times {10}^{-1}\right.$ at $N=8$) and traditional fPINN $(3.73\times$ ${10}^{-1}$ at $N=8$).

Figure 7 illustrates the point-wise error distribution across the computational domain as a heatmap. The errors are minimal, mostly concentrated in the range of ${10}^{-9} to {10}^{-8}$, with slight peaks near boundary regions. The consistency of the error distribution clearly shows that the PINN-CS framework has good accuracy over the whole domain. Small humps arise around the edges of the plot, suggesting areas where precision might be further improved by additional refinement of the data. This shows the relevance of node distribution in consideration of boundary effects.

The point-wise error distribution about Figure 8 in the paper is implemented for the 3D time-space fractional Bloch-Torrey equations for evaluating the efficiency of the Physics-Informed Neural Network coupled with the Cuckoo Search algorithm, called PINN-CS. The heatmap in figure 8 specifies that there exists minimal standard error spread wide and far across the computational domain, with most evaluations clustered closely to the range of ${10}^{-9}$ to ${10}^{-8}$. This observation further Fractal accuracy of the PINN-CS methodology to solve complex equations within a shorter CPU time. An additional interesting aspect can be revealed by looking closely at how the errors are distributed, which is the slight increase in values of the errors during the approach toward the bounds of the computational domain. These peaks can be attributed to challenges of approximating boundary conditions in high dimensional fractional systems. The phenomenon of boundary effects points to certain directions for improvement with regard to sampling techniques or improved management of boundary conditions.

Finally, the heat map shown in Figure 9 represent the point-wise error distribution for the 3D Bloch-Torrey equations. As shown here all the errors are bounded and spread out evenly throughout the computational domain with dominant values between 10^-8 and 10^-9. Defective scores, slightly higher near the boundaries, indicate scope for increase in efficacy concerning handling aspects of boundaries. This figure strengthens the PINN-CS framework to hold high accuracy in the domain and provide the ability to solve complicated fractional PDEs.

The PINN-CS solution of the 2D space-fractional PDE is presented by the 3D surface plot in Fig 10. This figure presents the numerical results in which it is observed that the solution surface is continuous and smooth and is in good agreement with the exact solution. The low point-wise errors ($\sim {10}^{-10}$ for $N=128$) ensure the solution's accuracy across the domain. Here, this graphical representation demonstrates the ability of PINN-CS framework to predict several behaviors of fractional PDEs in terms of accuracy as well as to address high-dimensional problems in more complicated areas. The Solution presented below is visually more coherent and accurate further strengthening the architecture of the proposed framework.

Table 3 provides parameter settings for Example 3.2, which involves solving a 2D time-space fractional diffusion equation. The parameters include approximate points 100, indicating the number of approximation points used for numerical computation. Hidden layers are 4, representing the depth of the neural network architecture. Neurons are 5 per layer, determining the complexity of the neural network. Activation function is $tanh$, chosen for its smooth and differentiable nature, suitable for physics-informed neural networks. Optimizer is Cuckoo Search Algorithm, employed for tuning parameters efficiently. The number of iterations is set to 1000, ensuring adequate optimization steps for convergence. The fractional orders $\alpha =0.3$ and $\beta =1.7$ highlight the characteristics of the diffusion process in the time and spatial domains, respectively.

Example 3.2

We move forward by examining the following 2D time-space fractional diffusion equation, along with its initial and boundary conditions, with in three distinct regions ${\Omega }_1$.

$$\begin{array}{c}{}_c{}^{0}D_t^{\alpha }u(x,y,t)=-K_1(-\Delta {)}_{xx}^{{\beta }_1}u(x,y,t)-K_2(-△{)}_{yy}^{{\beta }_2}u(x,y,t)\\ +f(x,y,t),(x,y,t)\in {\mathsf{\Omega }}_i\times (0,T]\end{array},$$

(41)

with the initial conditions:

$$u\left(x,y,0\right)=u_0\left(x,y\right),\left(x,y\right)\in {\mathsf{\Omega }}_i,$$

(42)

with the boundary condition:

$$u(x,y,t)=0,(x,y)\in \partial {\mathsf{\Omega }}_i\times (0,T].$$

(43)

The exact solution considered for the above problem is given by:

$$u\left(x, y, t\right)= {(x-x_{lb})}^3{(x_{ub}-x)}^3{(y-y_{lb})}^3{(y_{ub}-y)}^3\cos \left(t\right).$$

(44)

Given this exact solution, the forcing term f (x, y, t) can be expressed as:

$$f(x,y,t)=f_1(x,y,t)+f_2(x,y,t)+f_3(x,y,t).$$

(45)

Where $f$₁$(x,y,t),f$₂$(x,y,t),$ and$f$₃$(x,y,t)$ are expressed as:

$$f_1(x, y, t) = {(x-x_{lb})}^3{(x_{ub}-x)}^3{(y-y_{lb})}^3{(y_{ub}-y)}^3\cos (t+{\displaystyle \frac{п\alpha }{2}}),$$

(46)

$$\begin{array}{ll}{f}_2(x,y,t)=& {\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}(\pi \alpha /2)}}\left\{{\displaystyle \frac{{\left|x_{lb}+x_{ub}\right|}^3\mathsf{\Gamma }(4)}{\mathsf{\Gamma }(4-\alpha )}}\left({\left(x-x_{lb}\right)}^{3-\alpha }+{\left(x_{ub}-x\right)}^{3-\alpha }\right)\right.\\ & -{\displaystyle \frac{3{\left|x_{lb}+x_{ub}\right|}^2\mathsf{\Gamma }(5)}{\mathsf{\Gamma }(5-\alpha )}}\left({\left(x-x_{lb}\right)}^{4-\alpha }+{\left(x_{ub}-x\right)}^{4-\alpha }\right)\\ & +{\displaystyle \frac{3\left|x_{lb}+x_{ub}\right|\mathsf{\Gamma }(6)}{\mathsf{\Gamma }(6-\alpha )}}\left({\left(x-x_{lb}\right)}^{5-\alpha }+{\left(x_{ub}-x\right)}^{5-\alpha }\right)\\ & \left.-{\displaystyle \frac{3\mathsf{\Gamma }\left(7\right)}{\mathsf{\Gamma }\left(7-\alpha \right)}}\left({\left(x-x_{lb}\right)}^{6-\alpha }+{\left(x_{ub}-x\right)}^{6-\alpha }\right)\right\}{\left(y-y_{lb}\right)}^3{\left(y_{ub}-y\right)}^{3}cos(t),\end{array}$$

(47)

$$\begin{array}{ll}{f}_3(x,y,t)=& {\displaystyle \frac{1}{2\mathrm{c}\mathrm{o}\mathrm{s}(\pi \beta /2)}}\left\{{\displaystyle \frac{{\left|y_{lb}+y_{ub}\right|}^3\mathsf{\Gamma }(4)}{\mathsf{\Gamma }(4-\beta )}}\left({\left(y-y_{lb}\right)}^{3-\beta }+{\left(y_{ub}-y\right)}^{3-\beta }\right)\right.\\ & -{\displaystyle \frac{3{\left|y_{lb}+y_{ub}\right|}^2\mathsf{\Gamma }(5)}{\mathsf{\Gamma }(5-\beta )}}\left({\left(y-y_{lb}\right)}^{4-\beta }+{\left(y_{ub}-y\right)}^{4-\beta }\right)\\ & +{\displaystyle \frac{3\left|y_{lb}+y_{ub}\right|\mathsf{\Gamma }(6)}{\mathsf{\Gamma }(6-\beta )}}\left({\left(y-y_{lb}\right)}^{5-\beta }+{\left(y_{ub}-y\right)}^{5-\beta }\right)\\ & \left.-{\displaystyle \frac{3\mathsf{\Gamma }\left(7\right)}{\mathsf{\Gamma }\left(7-\beta \right)}}\left({\left(y-y_{lb}\right)}^{6-\beta }+{\left(y_{ub}-y\right)}^{6-\beta }\right)\right\}{\left(x-x_{lb}\right)}^3{\left(x_{ub}-x\right)}^{3}cos(t).\end{array}$$

(48)

For the sake of simplicity, we will assume $K_1$$=1, =1,$ and $T=1$. The Table 3 shows the parameters setting for the 2D time-space fractional diffusion equation. Example 3.2 focuses on solving a 2D time-space fractional diffusion equation using the proposed PINN-CS framework. The governing equation is characterized by fractional orders $\alpha =0.3$ (time) and $\beta =1.7$ (space), reflecting the non-local temporal and spatial dynamics inherent in fractional systems. The problem includes specified initial and boundary conditions over a given domain, with the exact solution provided for validation purposes.

Table 4 provides a comparative analysis of the numerical performance of three methods-fPINN, Monte Carlo fPINN, and the proposed PINN-CS framework-when applied to a 2D time-space fractional diffusion equation. The comparison is based on different values of $N$, the number of approximation points, for fractional orders $\alpha =0.3$ and $\beta =1.7$. It is observed that fPINN shows relatively lower accuracy, with the error remaining in the range of ${10}^{-1}$ even at higher approximation points. Monte Carlo fPINN improves upon fPINN by using Monte Carlo techniques, achieving errors in the order of ${10}^{-2}$ for larger $N$. On the other hand, PINN-CS outperforms both, with errors consistently in the range of ${10}^{-8}$ to ${10}^{-9}$ across all $N$, demonstrating a dramatic improvement in accuracy. Across all methods, increasing $N$ leads to lower errors, indicating improved numerical resolution. However, the rate of improvement is most pronounced in the PINN-CS framework, showcasing its superior ability to handle high-dimensional fractional PDEs. Moreover, at $N=128$, the error for fPINN remains significantly higher compared to PINN-CS, with a difference of several orders of magnitude. This highlights the enhanced precision of the PINN-CS approach, particularly in capturing the complexities of fractional diffusion equations.

Figure 11 shows a 3D graphical representation of the solution for Example 3.2, where a 2D time-space fractional diffusion equation is solved using the proposed PINN-CS framework. The figure further clearly expresses the rapid and concise character of the solution and guarantees that the method is capable of capturing all the performances of fPDEs within the computational domain. In the 3D surface plot, one can observe the numerical solution that has been developed in the PINN-CS framework while maintaining continuous and accurate forms in all structures found in the domain. This is particularly important because PINN-CS is able to model the fractional dynamics, such as time and spatial which is given in the governing equations. The positioning of the solution to the expected physical behavior in graphical form gives qualitative assurance of the strength of the method. The figure also demonstrates the ability of the PINN-CS approach to operate in high-dimensional spaces and address the issues with the boundary conditions discussed in Example 3.2.

Table 5 provides the specific parameter settings utilized for solving the 3D time space fractional Bloch-Torrey equation (Example 3.3) using the PINN-CS (Physics-Informed Neural Networks integrated with Cuckoo Search) framework.

Example 3.3

Examine the 3D time-space fractional Bloch-Torrey equation defined with in a spherical domain Ω= (x, y, z):$\sqrt{x^2+y^2+z^2}$ < r = 1/2 considering the homogeneous initial and boundary conditions [34].

$${K_{\mathsf{\alpha }}}^{\mathrm{C}}{D^{\mathsf{\alpha }}}_{\mathrm{t}}u(x,y,z,t)=K_{\mathsf{\beta }}{R}^{2\mathsf{\beta }}u(x,y,z,t)+f(x,y,z,t),$$

$$u(x,y,z,0)=0,(x,y,z)\in \Omega ,$$

(49)

$$u(x,y,z,t)=0,(x,y,z,t)\in R^3\Omega \times (0,T],$$

Where

$$R^{2\mathsf{\beta }}:{=}^{\mathrm{Rz}} D^{{2\mathsf{\beta }}_{\mathrm{x}}}{+}^{\mathrm{Rz}} D^{{2\mathsf{\beta }}_{\mathrm{y}}}{+}^{\mathrm{Rz}} D^{{2\mathsf{\beta }}_{\mathrm{z}}},$$

(50)

where^Rz$D$^2β_xis the Riesz fractional derivative operator,

$${}^{\mathrm{Rz}}D^{{2\mathsf{\beta }}_{\mathrm{x}}}f(x):={\displaystyle \frac{-1}{2cos(\pi \beta )}}(\mathrm{xlb}{D}^{{2\mathsf{\beta }}_{\mathrm{x}}}f(x)+xD{2}^{\mathsf{\beta }}{x}_{\mathrm{ub}}f(x)).$$

(51)

where:

$x_{ib}{D}_x^{2\beta }f(x)$ is the left sided Riemann-Liouville fractional derivative of order $2\beta$,

${}_{x}D_{x_m}^{2\beta }f(x)$ is the right sided Riemann-Liouville fractional derivative of order $2\beta$,

$x_{lb}$ and $x_{ub}$ are the lower and upper bounds for the domain of $x$,

$\mathrm{c}\mathrm{o}\mathrm{s}(\pi \beta )$ is a normalization factor to ensure symmetry.

The exact solution of the given problem is:

$$u(x,y,z,t)=2^4{t}^2(x^2+y^2+z^2-r^2{)}^2.$$

(52)

The forcing term of the above problem can be expressed as:

$$\begin{gathered} f(x,y,z,t)=\mathrm{h} (\mathrm{x}, \sqrt{r^2-y^2-z^2}) + \mathrm{h} (\mathrm{y}, \sqrt{r^2-x^2-z^2}) + \mathrm{h} (\mathrm{z}, \sqrt{r^2-x^2-y^2} \\ +{\displaystyle \frac{32k_{\alpha }}{\mathrm{┌}(3-\alpha )}}{t}^{2\text{-}\mathsf{\alpha }}(x^2+y^2+z^2-r^2{)}^2,) \end{gathered}$$

(53)

With

$$\begin{gathered} h(x,s)={\displaystyle \frac{8K_{\beta t^2}}{cos(\frac{\beta \pi }{2})}}({\displaystyle \frac{8s^2}{\mathrm{┌}(3-2\beta )}}((x+s{)}^{2\text{-}2\mathsf{\beta }}+(s-x{)}^{2\text{-}2\mathsf{\beta }})-{\displaystyle \frac{24s}{\mathrm{┌}(4-2\beta )}}((x+s{)}^{3\text{-}2\mathsf{\beta }}+(s-x{)}^{3\text{-}2\mathsf{\beta }}) \\ +{\displaystyle \frac{24s}{\mathrm{┌}(5-2\beta )}}((x+s{)}^{4\text{-}2\mathsf{\beta }}+(s-x{)}^{4\text{-}2\mathsf{\beta }}, \end{gathered}$$

(54)

for simplicity we take K_α=1, K_β=1 and T=1.

Figure 12 tracks the evolution of the best fitness value during the Cuckoo Search (CS) optimization. The initial fitness value starts relatively high, reflecting the significant discrepancy between the initial model parameters and the exact solution. However, the fitness rapidly decreases over the first 200 generations, indicating quick convergence of the CS algorithm. Beyond 500 generations, the fitness value plateaus, stabilizing near its optimal value. The consistent improvement in fitness across generations demonstrates the effectiveness of the CS algorithm in optimizing the neural network's parameters to minimize the loss function and adhere to the governing equations.

Figure 13 depicts the best fitness value versus the number of generations for the 2D time-space fractional diffusion equations in Example 3.2. The graph acts as a performance indicator of the optimization implemented in the PINN-CS framework utilizing the CS optimization algorithm for adjusting the parameters of the neural network. The fitness value begins with a reasonably high value, suggesting large gap between the preliminary predictions of the neural network and exact solutions. The fitness value falls sharply during early generations, or at least the initial 200 generations as shown below. This is in concordance with the ability of the CS algorithm to improve parameters’ values at a faster rate compared to the rate of decrease in the loss function.

After 500 generations, the fitness curve levels off, indicating that the algorithm has likely reached an optimal or near-optimal solution. The steady improvement in fitness values reflects the effectiveness of the CS algorithm in optimizing the neural network parameters. This ensures that the approximated solutions align well with the governing equations of the fractional diffusion problem. The fast convergence underscores the computational efficiency of the PINN-CS framework, making it a practical and powerful tool for addressing complex fractional partial differential equations.

Figure 14 tracks the evolution of the best fitness value during the optimization process for the 3D Bloch-Torrey equation. The fitness starts greatly, as the initial model may be greatly off the true solution. As it can be seen here the dropping is most prominent in the first 200 generations showing that the Cuckoo Search algorithm converges very quickly. The fitness value remains constant after 500 generation; this signifies that the evolution process has reached an optimal or a very near optimal solution. These results clearly show that the proposed CS algorithm can effectively tune parameters so that the approximated solutions are very close to the derived governing equations.

The above Figure 15 also presents the convergence graph of the PINN-CS for 3D Bloch-Torrey equations during the PINN-CS optimization process. Like Figure 14, the fitness value initiates from a high level and drops a great deal of in the initial phase which indicates optimization. After a 500 generation, it is observed that the curve starts to level off and approaches the optimal solution. This pattern also suggests that the Cuckoo Search algorithm is effective for adjustment of the network parameters in reaching accurate numerical solutions of the PINN-CS architecture.

The solutions furnished by the presented PINN-CS framework for Example 3.3 are given in Table 6. This has been done with the intention to present the performance of the proposed PINN-CS method with respect to accuracy and reliability in the table below. The first column is the number of discrete points that have been used in order to compute the numerical solutions. These points refer to the size of your computational domain where larger values of ‘N’ nets a greater level of detail. The second column holds the results which were obtained using the PINN-CS method. These are the forecasted solutions obtained employing the Physics-Informed Neural Network with the Cuckoo Search optimization technique. The third column presents the solutions to the problem, which are obtained analytically and compare with the PINN-CS results to evaluate the robustness of the PINN-CS models. The fourth column provides the absolute errors, calculated as the difference between the PINN-CS solutions and the exact solutions. These errors quantify the accuracy of the numerical method. The absolute errors reported in the fourth column are consistently negligible across all values of N, indicating that the PINN-CS framework achieves near-perfect agreement with the exact solutions. The low absolute errors demonstrate the computational efficiency of the PINN-CS method.

Figure 16 presents the absolute error values for the numerical solution across the computational domain. For $N=128$, the absolute error drops to $5.25\times {10}^{-10}$, as per Table 2. This figure highlights the rapid reduction in error magnitude as the number of points used for approximation grows. The visual representation confirms that PINN-CS consistently outperforms other methods by achieving significantly lower error rates across all resolutions. For example, at $N=$ 32, the PINN-CS error $\left(1.12\times {10}^{-9}\right)$ is orders of magnitude smaller than both Monte Carlo fPINN ($5.97\times {10}^{-2}$) and fPINN $\left(6.71\times {10}^{-2}\right)$.

Figure 17 illustrates the absolute errors in solving 2D time-space fractional diffusion equations using the proposed PINN-CS framework for Example 3.2. From the above graph, a detailed depiction of the distribution of errors throughout the computational domain is well illustrated showing how PINN-CS overcomes numerical discrepancies. The absolute errors are surprisingly low and that proves again the high accuracy of the proposed PINN-CS approach. For most computational points, the error is in the range of ${10}^{-9}$ to ${10}^{-10}$, significantly outperforming other methods like fPINN and Monte Carlo fPINN. The graph illustrates a steady decrease in error with the increase in approximation points, showcasing the PINN-CS method's adaptability and reliability for tackling high-dimensional fractional PDEs. Although the overall error remains minimal, minor spikes near the boundaries indicate potential improvements in node placement or boundary condition strategies to boost accuracy. This visualization underscores the PINN-CS framework's advantage over conventional approaches, excelling in both precision and computational performance.

The distribution of absolute error of the 3D time-space fractional Bloch-Torrey equations solved using the PINN-CS framework is shown in the Figure 18. Their absolute error values are very small in all cases at 10^-8 to 10^-9 and are relatively small across the computational domain. Such uniform distribution shows that the given method is rather accurate and stable. Small fluctuations are noticeable near the boundary areas, likely stemming from the inherent difficulty of precisely approximating boundary conditions in high-dimensional fractional systems. This demonstrates the PINN-CS framework's capability to achieve remarkable accuracy, confirming its reliability in addressing complex 3D fractional PDEs.

Figure 19 illustrates the absolute error distribution for the 3D time-space fractional Bloch-Torrey equations solved in Example 3.4. The goal of this figure is to plot the relative error throughout the computational domain as well as demonstrate the performance of the Physics-Informed Neural Network combined with the Cuckoo Search algorithm (PINN-CS). The absolute errors are small and are largely within the order of 10^-9 to 10^-8. Almost equal distribution is noted throughout the domain which underlines that even though the geometry complexity increases PINN-CS effectively stays in low error range throughout the computational space. The error values are also observed to have small fluctuations at both the extremes, which can be explained by the difficulties of defining extreme conditions for high-dimensional fractional systems. This warrants further development of the boundary treatment or sampling methods in these areas. Figure 19 serves as compelling evidence of the remarkable accuracy and efficiency of the PINN-CS framework, demonstrating its ability to closely align with the exact solutions. The possible error margin is very low and this fact underlines the usefulness of the described approach when it concerns practical uses.

Example 3.4

The three-term time-space fractional Bloch-Torrey equation in three dimensions on a finite domain is considered:

$$\begin{gathered} K_{{\alpha }_0}{\displaystyle \frac{\partial u\left(x,y,z,t\right)}{\partial t}}+K_{{\alpha }_1}{\displaystyle \frac{{\partial }^{{\alpha }_1}u\left(x,y,z,t\right)}{\partial t^{{\alpha }_1}}}+K_{{\alpha }_2}{\displaystyle \frac{{\partial }^{{\alpha }_2}u\left(x,y,z,t\right)}{\partial t^{{\alpha }_2}}}=K_{\beta }({\displaystyle \frac{{\partial }^{\beta }u\left(x,y,z,t\right)}{\partial |x{|}^{\beta }}}+{\displaystyle \frac{{\partial }^{\beta }u\left(x,y,z,t\right)}{\partial |y{|}^{\beta }}}+ \\ {\displaystyle \frac{{\partial }^{\beta }u\left(x,y,z,t\right)}{\partial |z{|}^{\beta }}})+ \\ f(x,y,z,t)\left(0\le x\le \mathrm{1,0}\le y\le \mathrm{1,0}\le z\le 1,t\ge 0\right), \end{gathered}$$

(55)

$$f\left(x,y,z,t\right)={\sum }_{i=1}^6 f_i\left(x,y,z,t\right),$$

(56)

$$f_1(x,y,z,t)={\displaystyle \frac{K_{{\alpha }_0}\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(2\right)}}t{x}^2(1-x{)}^2{y}^2(1-y{)}^2{z}^2(1-z{)}^2,$$

(57)

$$f_2(x,y,z,t)={\displaystyle \frac{K_{{\alpha }_1}\mathsf{\Gamma }(3)}{\mathsf{\Gamma }\left(3-{\alpha }_1\right)}}{t}^{2-{\alpha }_1}{x}^2(1-x{)}^2{y}^2(1-y{)}^2{z}^2(1-z{)}^2,$$

(58)

$$f_3(x,y,z,t)={\displaystyle \frac{K_{{\alpha }_2}\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(3-{\alpha }_2\right)}}{t}^{2-{\alpha }_2}{x}^2(1-x{)}^2{y}^2(1-y{)}^2{z}^2(1-z{)}^2 ,$$

(59)

$$\begin{gathered} f_4(x,y,z,t)={\displaystyle \frac{K_{\beta }}{2\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\beta \pi }{2}\right)}}{t}^2{y}^2(1-y{)}^2{z}^2(1-z{)}^2\left[{\displaystyle \frac{\mathsf{\Gamma }(3)}{\mathsf{\Gamma }(3-\beta )}}\left(x^{2-\beta }+(1-x{)}^{2-\beta }\right)\right. \\ \left.+{\displaystyle \frac{\mathsf{\Gamma }(5)}{\mathsf{\Gamma }(5-\beta )}}\left(x^{4-\beta }+(1-x{)}^{4-\beta }\right)-{\displaystyle \frac{2\mathsf{\Gamma }(4)}{\mathsf{\Gamma }(4-\beta )}}\left(x^{3-\beta }+(1-x{)}^{3-\beta }\right)\right], \end{gathered}$$

(60)

$$f_5(x,y,z,t)=f_4(y,x,z,t),f_6(x,y,z,t)=f_4(z,y,x,t).$$

(61)

With the initial boundary conditions:

$$u(x,y,z,0)=├ u(x,y,z,t)┤|\_\Gamma =0.$$

(62)

The analytical solution of the given problem is:

$$u(x,y,z,t)=t^2{x}^2(1-\mathrm{x}{)}^2{y}^2(1-y{)}^2{z}^2(1-z{)}^2.$$

(63)

Table 7 shows the parameter setting for the three-term time-space fractional Bloch–Torrey equation in 3-D on a finite domain.

Table 8 provides a comparison between the predicted solutions generated by the PINN-CS framework (a Physics-Informed Neural Network integrated with the Cuckoo Search algorithm) and the exact solutions for the 3D Bloch-Torrey equations from Example 3.4. It also includes the absolute errors, offering a detailed assessment of the PINN-CS approach's accuracy and effectiveness. The solutions predicted by the PINN-CS framework align closely with the exact results for all tested approximation points N. The absolute errors are remarkably low, ranging between 10⁻¹⁰ and 10⁻⁸. This level of precision highlights the reliability of the PINN-CS method in delivering accurate results across the entire computational domain, even as the problem complexity varies.

4. Conclusion

In this study, we introduced a novel framework combining Monte Carlo-based Physics-Informed Neural Networks (PINNs) with the Cuckoo Search (CS) algorithm, referred to as PINN-CS, for efficiently solving fractional partial differential equations (fPDEs). The proposed framework capitalizes on the advantages of applying QMC approaches for approximations of the fractional derivatives and the optimization ability of the CS algorithm to reach high total precision and fast convergence rate. Several numerical simulations were performed and the algorithm was tested to confirm that PINN-CS method works for the 2D space-fractional Poisson equations, the 2D time-space fractional diffusion equations, and the 3D time-space fractional Bloch-Torrey equations. The results consist that the proposed PINN-CS framework is more accurate and efficient than the fPINN and Monte Carlo fPINN techniques. The usage of the structured sampling nodes and block-sampled approximations saves computer space and reduces computational time considerably however, the precision of the solution estimate is still high. Further, the performance of proposed PINN-CS in approximating the solutions of irregular domains and the corresponding fuzzy boundary conditions clearly shows that the deep learning based PINNs are suited well for applied problems. The convergence analysis indicates fast convergence of the network parameters with the help of the CS algorithm thereby making the approximated solution to be in correct physical approximation with governing equation. The best absolute errors and near perfect overlaps with exact solutions testify towards the effectiveness and efficiency of the proposed framework. Further, the results show the significance of proper distribution of nodes and further advanced optimization for achieving increased accuracy while solving fractional PDEs of high dimensions.

The PINN-CS framework defines new base lines for solving fractional PDEs especially for issues that arise with irregular geometries, fuzzy boundaries, and high dimensional systems. The approach described in this research is easily transferable to computational science and engineering including biomechanics, materials design, and heated fluids. The PINN-CS framework defines new base lines for solving fractional PDEs especially for issues that arise with irregular geometries, fuzzy boundaries, and high dimensional systems. The approach described in this research is easily transferable to computational science and engineering including biomechanics, materials design, and heated fluids.